Valid XHTML     Valid CSS2    

meta 

Structure

Evolution

Et

Diversité des microbiotes des semences

1. Un script pour automatiser les traitements en métagénomique avec mothur et R

Le script metauto16S.pl et ses fichiers annexes dont ceux pour mothur et ses programmes perl sont disponibles dans l'archive metauto16S.zip.

Le but du script est d'aider à la réalisation d'une analyse métagénomique sur des séquences fastq associées à divers "facteurs" (jours, pays, échantillons...).

1.1 Préparation des données et recopie du script

Avant d'utiliser le script, il est nécessaire de créer un répertoire d'analyse structuré de la façon suivante :

  • il y a un sous-répertoire Input/ qui contient les séquences fastq 
  • il y a un sous-répertoire Output/ (vide au départ) qui contiendra les fichiers temporaires de calcul 
  • il y a un sous-répertoire Results/ (vide au départ) qui contiendra les fichiers PDF et CSV générés à la fin de l'analyse 
  • les fichiers de l'archive ont été décompressés dans le répertoire de l'analyse au-dessus de Input, Output et Results.

Par exemple si votre analyse se nomme Metaseed, le répertoire Metaseed doit ressembler à :


     Metaseed/            # répertoire principal de l'analyse avec le script et ses annexes
     Metaseed/Input/      # données fastq
     Metaseed/Output/     # fichiers temporaires
     Metaseed/Results/    # fichiers résultats
     
     

Pour que le script puisse détecter des facteurs (échantillons, dates, pays...) il faut que les fichiers fastq utilisent le séparateur _ (soit : underscore ou encore : espace souligné, "tiret du 8"). Par exemple un fichier identifié par S01_H24_GTAAGT_L001_R1.fastq fait partie potentiellement de deux facteurs, détectés par S01 et par H24. Si par contre les informations sont collées, comme dans S01H24_GTAAGT_L001_R1.fastq le script ne sera pas capable de reconnaitre les bons facteurs.

1.2 Contenu de l'archive et utilisation du script

Pour utiliser le script, il suffit d'aller dans le répertoire d'analyse et d'y décompresser tous les fichiers de l'archive. Ensuite, il faut taper perl metauto16S.pl puis la ou les options désirées. L'analyse comporte 7 étapes qu'il faut réaliser dans l'ordre.

Si vous disposez du programme wget il vous suffit de copier/coller les instructions suivantes pour que vous n'ayez plus qu'à recopier vos données dans le répertoire Input/ avant d'exécuter le script :


     # on est dans le répertoire d'analyse
     # création des sous-répertoires
     
     mkdir Input
     mkdir Output
     mkdir Results
     
     # rapatriement de l'archive du script
     
     wget http://forge.info.univ-angers.fr/~gh/Metaseed/metauto16S.zip
     
     # décompression de l'archive
     
     unzip metauto16S.zip
     
     # il n'y a plus qu'à recopier les données dans Input/
     # avant d'exécuter le script...
     

Le script est écrit en perl mais vous n'avez pas besoin de connaitre perl pour vous en servir. Par contre, perl doit être installé. De même, le script utilise des programmes R mais, là encore, vous n'avez pas besoin de connaitre R pour vous servir du script. Par contre, R doit être installé et les packages suivants doivent avoir été installés (faites-vous aider si vous ne savez pas comment installer un package) :


     beanplot
     car
     ggplot2
     plyr
     reshape2
     scales
     

Voici le contenu de l'archive du script :

Archive:  metauto16S.zip
 Length   Method    Size  Cmpr    Date    Time   CRC-32   Name
--------  ------  ------- ---- ---------- ----- --------  ----
   36530  Defl:N     8320  77% 2015-01-14 09:52 71571fbb  metauto16S.pl
    6944  Defl:N     1981  72% 2014-10-13 14:13 b21d1b71  elim_LowFreqOTU_from_shared_make_database.pl
    5489  Defl:N     1578  71% 2014-10-20 17:30 5f1c74b3  calc_abundance_count.pl
    4441  Defl:N      926  79% 2014-07-27 11:40 50f1311c  quality.mothur
     244  Defl:N      141  42% 2015-01-13 16:27 02c37986  amova.mothur
     596  Defl:N      250  58% 2015-01-13 16:30 a9d9173b  beta-div.mothur
      54  Defl:N       52   4% 2013-07-18 09:51 bdee38e1  test16S.oligos
200975572  Defl:N  1505886  99% 2013-11-19 09:55 0af82ba4  silva.v4.fasta
15907100  Defl:N  2536620  84% 2012-07-30 16:19 7b35f1d5  trainset9_032012.pds.fasta
 1123514  Defl:N   145259  87% 2012-07-30 16:19 207c5b46  trainset9_032012.pds.tax
    1444  Defl:N      736  49% 2015-01-13 17:18 21c4e424  beanplots.r
    7113  Defl:N     2849  60% 2015-01-13 17:19 784cf0be  taxonomy.r
    4597  Defl:N     1027  78% 2015-01-13 17:52 776ecd4e  nmdspcoa.r
     243  Defl:N      185  24% 2014-10-20 18:10 816f081c  alpha-divMOCK.mothur
     233  Defl:N      178  24% 2014-10-20 17:56 44113438  alpha-divNOMOCK.mothur
    1355  Defl:N      607  55% 2014-10-20 17:24 4bc85b3d  groupsMOCK.mothur
    1266  Defl:N      573  55% 2014-10-18 19:31 9ea012d9  groupsNOMOCK.mothur
--------          -------  ---                            -------
218076735          4207168  98%                            17 files

Le script comporte des options, visibles via l'option --help. Voici quelques exemples d'appel du script :

1.3 Rappel de l'aide

$> perl metauto16S.pl --help

 metauto16S.pl (gH) version 1.63

    syntax : metauto16S.pl OPTION where OPTION can be

     --help     # displays this help
     --steps    # details the steps
     --status   # shows completed steps
     --clean    # removes temporary and log files
     --step 1 | 2 | 3 | 4 | 5 | 6 | 7 --to 2 | 3 | 4 | 5 | 6 | 7
     --continue # starts step 1 or proceeds to next step
     --factors  # to see what the possible factors are (step 4)
     --step 5 --factor 1 | 2 | 3... # use the appropriate 16S.design*

1.4 Description des étapes

$> perl metauto16S.pl --steps

 metauto16S.pl (gH) version 1.63

 --steps of metauto script:

  1. build stability.files
  2. check quality of fragments
  3. build groups (may be long)
  4. run abundance computations (may be long)
  5. run alpha-diversity analysis, determine factors
     + use of R for beanplots and taxonomic clustering
  6. run beta-diversity analysis
     + use of R for NMDS and PCOA
  7. run amova analysis

1.5 Statut en cours d'analyse


     $>  perl metauto16S.pl --step 1
     
     [...]
     
     $>  perl metauto16S.pl --continue
     
     [...]
     
     $>  perl metauto16S.pl --status
     
      metauto16S.pl (gH) version 1.63
     
      --status of metauto script:
     
       step 1 (build stability.files and determine factors):
          completed
       step 2 (check quality of fragments):
          completed
       step 3 (build groups (long)):
          UNCOMPLETED
       step 4 (abundance computations):
          UNCOMPLETED
       step 5 (alpha-diversity step with R beanplots and taxonomic clustering):
          UNCOMPLETED
       step 6 (beta-diversity step):
          UNCOMPLETED
       step 7 (amova analysis):
          UNCOMPLETED
     
     
     

2. Un exemple naif d'exécution pour les étapes 1 à 7

La version longue, non complètement expurgée des sorties mothur, peut être lue dans le fichier metauto_demo_long.txt.


     $> metauto --step 1 --to 7
     
     metauto.pl (gH) version 1.47
     
     Step 1:  build stability.files
     
      required file alpha-divMOCK.mothur is present
      required file alpha-divNOMOCK.mothur is present
      required file amova.mothur is present
      required file beanplots.r is present
      required file beta-div.mothur is present
      required file calc_abundance.pl is present
      required file elim_LowFreqOTU_from_shared_make_database.pl is present
      required file groupsMOCK.mothur is present
      required file groupsNOMOCK.mothur is present
      required file quality.mothur is present
      required file silva.v4.fasta is present
      required file test16S.oligos is present
      required file trainset9_032012.pds.fasta is present
      required file trainset9_032012.pds.tax is present
     
      no MOCK files found.
     
     
      2 files written in stability.files for *R1.fastq files.
     
      -- 2 stability.files *R1* and *R2* used.
     
      -- step 1 completed.
     
     Step 2:  check quality of fragments (long)
     
       -- mothur quality.mothur
     
       mothur v.1.33.3
       [...]
     
       mothur > quit()
     
      -- step 2 completed.
     
     Step 3:  build groups
     
       -- mothur groupsNOMOCK.mothur
       [...]
     
       -- step 3 completed.
     
     Step 4:  abundance computations
     
       -- elim_LowFreqOTU
     
       perl elim_LowFreqOTU_from_shared.pl -s 16S.an.unique_list.shared -f 0.1 -o 16S.an.unique_list.abund.shared
     
       -- calc_abundance
     
       perl calc_abundance.pl -s 16S.an.unique_list.abund.shared -f 1000 -o 16S.an.unique_list.abund.proportion.shared
     
       -- step 4 completed.
     
     Step 5:  alpha-diversity step with R beanplots and taxonomic clustering
     
       -- mothur alpha-divNOMOCK.mothur
       [...]
     
       let's use factor 1 to build Output/16S.design ;
       values: S01 S03
     
       -- R beanplots.r
       [...]
     
       -- R taxonomy.r
       [...]
     
       -- step 5 completed.
     
     Step 6:  beta-diversity step
     
       -- mothur beta-div.mothur
       [...]
     
       -- step 6 completed.
     
     Step 7:  amova analysis
     
        -- mothur amova.mothur
        [...]
     
        -- step 7 completed.
     
     Start: 20/10/2014 18:11 ; end: 20/10/2014 18:31 ; elapsed time: 20 min 03 sec.
     

3. Choix de facteur et étape 5

Le script essaie de détecter les facteurs avec l'option --factors. Les facteurs possibles sont alors affichés, il suffit de choisir celui que l'on veut avant d'exécuter ou de ré-éxécuter l'étape 5. Par défaut, l'étape 5 utilise --factor 1.

Attention : il faut avoir exécuté l'étape 4 (bien sûr !) pour que les facteurs puissent être déterminés.


     $> metauto --factors
     
      metauto.pl (gH) version 1.35
     
      Looking for factors in 16S.an.unique_list.abund.shared...
     
      for factors, you may use the option(s)
      --factor 1 : S01 S03
      --factor 2 : H0 H24
     
     $> metauto --step 5 # équivalent à --step 5 --factor 1
     
        Step 5:  run alpha-diversity analysis, determine factors, make R beanplots and taxonomic clustering
     
        let's use factor 1 to build Output/16S.design ;
        values: S01 S03
     
        -- mothur alpha-div.mothur
        [...]
     
        -- R beanplots.r
        [...]
     
        -- R taxonomy.r
        [...]
     
     -- step 5 completed.
     
     
     $> metauto --step 5 --factor 2
     
        metauto.pl (gH) version 1.35
     
        Step 5:  run alpha-diversity analysis, determine factors, make R beanplots and taxonomic clustering
     
        let's use factor 2 to build Output/16S.design ;
        values: H0 H24
     
        -- mothur alpha-div.mothur
        [...]
     
        -- R beanplots.r
        [...]
     
        -- R taxonomy.r
        [...]
     
     -- step 5 completed.
     

4. Exemple de résultats

Ce script a notamment permis de traiter, dans le cadre d'une analyse amplicon 16S 172 fichiers de reads, de taille moyenne 47 Mo soit une taille totale de 8,2 Go (1,5 Go une fois compressés). en à peu près 16 h sur un petit serveur dédié (2 processeurs, 4 coeurs, 32 Go RAM) pour 4 facteurs d'intérêt. Vous pouvez consulter l'archive des résultats.

5. Contenu du script metauto16S.pl


     # # (gH)   -_-  metauto16S.pl  ;  TimeStamp (unix) : 13 Janvier 2015 vers 18:26
     
     ##################################################################
     ##################################################################
     ##                                                              ##
     ##  metauto.pl: a script for an automated treatement of         ##
     ##              fastq sequences                                 ##
     ##                                                              ##
     ##  authors: gilles.hunault@univ-angers.fr                      ##
     ##           matthieu.barret@angers.inra.fr                     ##
     ##                                                              ##
     ##  with the help of martial.briand.@angers.inra.fr             ##
     ##           for the external perl scripts                      ##
     ##                                                              ##
     ##################################################################
     ##                                                              ##
     ## see:                                                         ##
     ##                                                              ##
     ##    http://forge.info.univ-angers.fr/~gh/Metaseed/metauto.php ##
     ##                                                              ##
     ## for a Web documentation                                      ##
     ##                                                              ##
     ##################################################################
     ##################################################################
     
     print "\n metauto16S.pl (gH) version 1.63\n" ;
     
     use strict   ;
     use warnings ;
     
     my $in  = "Input/" ;
     my $out = "Output/" ;
     my $res = "Results/" ;
     
     # requires mothur 1.33
     
     my $mothur = "mothur" ; # you may write here the path to mothur
     
     my $date1 = &dateEtHeure() ;
     my $time1 = time() ;
     if (-f $out) { system("echo $date1 > ".$out."time.start") ; } ;
     
     ################################################################
     #
     # show help if no arguments
     #
     ################################################################
     
     if ($#ARGV==-1) { # no argument, display the options
     
      &help() ;
      exit(-1) ;
     
     } ; # fin de test sur les arguments
     
     # constant global variables
     
     my @stepsArray = (
       "build stability.files",
       "check quality of fragments",
       "build groups (may be long)",
       "run abundance computations (may be long)",
       "run alpha-diversity analysis, determine factors\n     + use of R for beanplots and taxonomic clustering",
       "run beta-diversity analysis\n     + use of R for NMDS and PCOA",
       "run amova analysis") ;
     
     # path to all script files including this metauto.pl
     
     my $pathf = substr($0,0,length($0)-length("metauto16S.pl")) ;
     
     ################################################################
     #
     # check arguments and decide what to do
     #
     ################################################################
     
     my $firstArg = $ARGV[0] ;
     
     if ($firstArg eq '--help') {
      &help()
     } elsif ($firstArg eq '--steps') {
      &steps() ;
     } elsif ($firstArg eq '--status') {
      &status() ;
     } elsif ($firstArg eq '--statut') { # for the french people
      &status() ;
     } elsif ($firstArg eq '--clean') {
      &clean() ;
     } elsif ($firstArg eq '--continue') {
      &continue() ;
     } elsif ($firstArg eq '--factors') {
      &possibleFactors() ;
     #} elsif ($firstArg eq '--factor') {
     # &buildForFactor(@ARGV) ;
     } elsif ($firstArg eq '--step') {
     
        if ($#ARGV>0) {
          &step(@ARGV) ;
        } else {
          print " nbargs $#ARGV \n\n" ;
          print "\n after --step you must supply an integer from 1 to 7\n\n " ;
          exit(-2) ;
        } ; # fin si
     
     } else {
       print "\n unknown option \"$firstArg\". use metauto --help to list the options.\n\n" ;
     } ; # fin si
     
     exit(0) ;
     
     # end of main program
     
     ####################################################################
     ####################################################################
     
     
     ####################################################################
     #
     # subprograms
     #
     ####################################################################
     
     sub help {
     
     ####################################################################
     
     # display all options
     
     print << 'FIN_HELP' ;
     
         syntax : metauto16S.pl OPTION where OPTION can be
     
          --help     # displays this help
          --steps    # details the steps
          --status   # shows completed steps
          --clean    # removes temporary and log files
          --step 1 | 2 | 3 | 4 | 5 | 6 | 7 --to 2 | 3 | 4 | 5 | 6 | 7
          --continue # starts step 1 or proceeds to next step
          --factors  # to see what the possible factors are (step 4)
          --step 5 --factor 1 | 2 | 3... # use the appropriate 16S.design*
     
     FIN_HELP
     exit(-1) ;
     
     } ; # end of help
     
     ####################################################################
     
     sub steps {
     
     ####################################################################
     
     # details what the steps do
     
     print "\n".' --'."steps of metauto script:\n\n" ;
     
     foreach my $istep (0..$#stepsArray) {
       my $jstep = $istep + 1 ;
       print "  $jstep. ".$stepsArray[$istep]."\n" ;
     } ; # fin de foreach
     
     print "\n" ;
     exit(-1) ;
     
     } ; # end of steps
     
     ####################################################################
     
     sub status {
     
     ####################################################################
     
     # tells if steps have been completed or not
     
     print "\n".' --'."status of metauto script:\n\n" ;
     
     foreach my $istep (0..$#stepsArray) {
       my $jstep = $istep + 1 ;
       print "  step $jstep (".$stepsArray[$istep]."):\n" ;
       my $fn = $out."step".$jstep.".ok" ;
       if (-e $fn) {
           print "     completed\n" ;
       } else {
           print "     UNCOMPLETED \n" ;
       } ; # fin de si
     } ; # fin de foreach
     
     print "\n" ;
     exit(-1) ;
     
     } ; # end of status
     
     ####################################################################
     
     sub clean {
     
     ####################################################################
     
     # remove some files to free disk space
     # and the *ok files when starting step 1
     
     print "\n".' --'."clean the temporary files:\n\n" ;
     print " 1. remove *.logfile            files\n" ;
     system("rm -rf ".$out."*.logfile") ;
     print " 2. remove ".$out."*.temp               files\n" ;
     system("rm -rf ".$out."*.temp") ;
     print " 3. remove 16S.an.unique_list*  files\n" ;
     system("rm -rf ".$out."16S.an.unique_list.abund.shared") ;
     system("rm -rf ".$out."16S.an.unique_list.abund.proportion.shared") ;
     print " 4. remove step*.ok             files\n" ;
     system("rm -rf ".$out."step*.ok") ;
     print " 5. remove mothur*.ok           files\n" ;
     system("rm -rf ".$out."mothur*.ok") ;
     
     print "\n" ;
     exit(-1) ;
     
     } ; # end of clean
     
     ####################################################################
     
     sub continue {
     
     ####################################################################
     
     # starts step 1 or determines the last completed step
     # and tries to execute the next one
     
     my $istep = 0 ;
     my $jstep = $istep + 1 ;
     my $fn = $out."step".$jstep.".ok" ;
     while (-e $fn) {
       $istep++ ;
       if ($istep>$#stepsArray) {
         print "\n Greetings! All steps completed. nothing to do.\n\n" ;
         exit(0) ,
       } ; # fin de si
       $jstep = $istep + 1 ;
       $fn = $out."step".$jstep.".ok" ;
     } ; # fin de tant que
     
     &step("--step $jstep") ;
     
     } ; # end of continue
     
     ####################################################################
     
     sub step {
     
     ####################################################################
     
     # executes one or more steps depending on the --continue option
     # or the --step option and and eventually the --to option
     
     my $args   = join(" ",@_) ;
     my $from   = 0 ;
     my $to     = 0 ;
     my $factor = 0 ;
     
     # first determine what step(s) to execute
     
     if ($args =~ /--to/) { # use --step and --to
     
        $args =~ /--step *([1-7])/ ;
        if ($1) { $from = $1 ; } ;
     
        $args =~ /--to *([2-7])/ ;
        if ($1) { $to   = $1 ; } ;
     
        if ($to==0) {
          print "\n".' if you use the syntax --step N1 --to N2,' ;
          print "\n".' after --step you must supply an integer from 1 to 7 and' ;
          print "\n".' after --to   you must supply an integer from 2 to 7.'."\n\n" ;
          exit(-2) ;
        } # fin si
     
     } else {
     
     # use only --step
     
        $args =~ /--step *([1-7])/ ;
        if ($1) {
           $from = $1 ;
           $to = $from ;
        } else {
          print "\n".' after --step you must supply an integer from 1 to 7.'."\n\n" ;
          exit(-2) ;
        } # fin si
     
     } # fin si
     
     # if step 5 determine the factor number
     
     if (($from>=5) | ($to<=5)) {
     
      $factor = 1 ;
      if ($args =~ /--factor/) { # which factor to use
     
        $args =~ /--factor *([1-7]?).*$/ ;
        if ($1) { $factor= $1 ; } ;
        if ($factor =~ /[^1-7]/) {
          print "\n".' if you use the syntax --factor #F,' ;
          print "\n".' after --factor you must supply an integer from 1 to 4.\n\n' ;
          exit(-2) ;
        } # fin si
     
      } # fin si
     
     } # fin si
     
     # and now, let's do it: call the right subprograms
     
     my $istep = $from ;
     while ($istep<=$to) {
       my $jstep = $istep -1 ;
       print "\n Step $istep:  ".$stepsArray[$jstep]."\n" ;
     
       if ($istep==1) {
     
     ################################################################
     #
     # step 1 : verification of the needed files
     #          construction of stability.files
     #
     ################################################################
     
         my @log0 = &check_mothur_version() ;
         my @log1 = &check_Files() ;
         my @log2 = &build_stability_files() ;
         &writeOkStep(1) ;
     
       } elsif ($istep==2) {
     
         if (-e $out."step1.ok") {
     
     ################################################################
     #
     # step 2 : script mothur 1 for quality of fragments
     #
     ################################################################
     
           my $cmd = "$mothur ../quality.mothur" ;
           print "\n-- $cmd \n\n" ;
           #print 'system("(cd Output ; $cmd | tee quality_sor.txt )")'."\n" ;
           my $mo = "quality_sor.txt" ;
           my $rc2 = system("(cd Output ; $cmd | tee $mo )") ;
           if (-e $out."mothur1.ok") {
             if (-e $out.$mo) { system("cp ".$out.$mo." ".$res.$mo) ; }
             &writeOkStep(2) ;
           } ; # fin si
     
         } else { # tell step 1 was not completed
     
           &noStep(1) ;
     
         } ; # fin si
     
       } elsif ($istep==3) {
     
     ################################################################
     #
     # step 3 : script mothur 2 for groups -- may lead to "kernel panic"
     #
     ################################################################
     
         if (-e $out."step2.ok") {
     
           # depending on MOCK files via the contents of mf.option -- Mock Files option
           # we use groupsMOCK.mothur or groupsNOMOCK.mothur file
     
           my $mf = `cat mf.option` ;
           chomp($mf) ;
           my $gm  = "?" ;
           if ($mf eq "NOMOCK") {
              $gm  = "groupsNOMOCK.mothur" ;
           } else {
              $gm  = "groupsMOCK.mothur" ;
           } ; # fin si
           my $cmd = "$mothur ../$gm" ;
           # print " mf $mf DONC gm $gm \n\n" ;
           print "\n-- $cmd \n\n" ;
           #print 'system("(cd Output ; $cmd | tee groups_sor.txt )")'."\n" ;
           my $mo = "groups_sor.txt" ;
           my $rc3 = system("(cd Output ; $cmd | tee $mo )") ;
           if (-e $out."mothur2.ok") {
             if (-e $out.$mo) { system("cp ".$out.$mo." ".$res.$mo) ; }
             &writeOkStep(3) ;
           } ; # fin si
     
         } else { # tell step 2 was not completed
     
           &noStep(2) ;
     
         } ; # fin si
     
       } elsif ($istep==4) {
     
     ################################################################
     #
     # step 4 : script perls for abundance
     #
     ################################################################
     
         if (-e $out."step3.ok") {
     
           &elim_LowFreqOTU() ;
           &calc_abundance()  ;
           &possibleFactors() ;
           &writeOkStep(4)    ;
     
         } else {
     
           &noStep(3) ; # tell step 3 was not completed
     
         } ; # fin si
     
       } elsif ($istep==5) {
     
     ################################################################
     #
     # step 5 : script mothur 3 for alpha-diversity and R beanplots and taxonomy
     #
     ################################################################
     
         if (-e $out."step4.ok") {
     
           # depending on MOCK files via the contents of mf.option -- Mock Files option
           # we use alpha-divMOCK.mothur or alpah-divNOMOCK.mothur file
     
           my $mf = `cat mf.option` ;
           chomp($mf) ;
           my $am  = "?" ;
           if ($mf eq "NOMOCK") {
              $am  = "alpha-divNOMOCK.mothur" ;
           } else {
              $am  = "alpha-divMOCK.mothur" ;
           } ; # fin si
           my $cmd = "$mothur ../$am" ;
           my $mo  = "alpha-div_sor.txt" ;
           my $rc5 = system("(cd Output ; $cmd | tee $mo) ") ;
     
           &buildForFactor($factor) ;
     
           my $designs = $out."16S.design" ;         # standard name
           my $designf = $out."16S.design$factor" ;  # name with factor number
     
           system("cp $designs $designf") ;
           unless (-e $designs) {
             print " no $designs file. Impossible to complete step 4.\n" ;
             exit(-1) ;
           } ; # fin unless
     
           rscripts($factor,"step5") ;  # run R scripts beanplots.r and taxonomy.r
     
           if (-e $out."mothur3.ok") {
                if (-e $out.$mo) { system("cp ".$out.$mo." ".$res.$mo) ; }
                if (-e $out."rscript1.ok") {
                   if (-e $out."rscript2.ok") {
                     &writeOkStep(5) ;
                   } ; # fin si
                } ; # fin si
           } ; # fin si
     
         } else {
     
           &noStep(4) ; # tell step 4 was not completed
     
         } ; # fin si
     
       } elsif ($istep==6) {
     
     ################################################################
     #
     # step 6 : script mothur 4 for beta-diversity
     #
     ################################################################
     
         if (-e $out."step5.ok") {
     
           my $cmd = "$mothur ../beta-div.mothur" ;
           my $mo  = "beta-div_sor.txt" ;
           my $rc6 = system("(cd Output ; $cmd | tee $mo )") ;
           if (-e $out."mothur4.ok") {
             if (-e $out.$mo) { system("cp ".$out.$mo." ".$res.$mo) ; }
             # ajout janvier 2015
             rscripts($factor,"step6") ;  # run R script nmdspcoa.r
             if (-e $out."rscript3.ok") {
               &writeOkStep(6) ;
             } ; # fin si
           } ; # fin si
     
         } else {
     
           &noStep(5) ; # tell step 5 was not completed
     
         } ; # fin si
     
       } elsif ($istep==7) {
     
     ################################################################
     #
     # step 7 : script mothur 5 for amova
     #
     ################################################################
     
         if (-e $out."step6.ok") {
     
           my $cmd = "$mothur ../amova.mothur" ;
           my $mo  = "amova_sor.txt" ;
           my $rc7 = system("(cd Output ; $cmd | tee $mo )") ;
           print "\n";
     
           if (-e $out."mothur3.ok") {
             if (-e $out.$mo) { system("cp ".$out.$mo." ".$res.$mo) ; }
             &writeOkStep(7) ;
           } ; # fin si
     
         } else {
     
           &noStep(6) ; # tell step 6 was not completed
     
         } ; # fin si
     
       } else {
     
         print "\n unknown step number $istep. stop.\n\n" ;
         exit(-3) ;
     
       } ; # finsi
     
       $istep++ ;
     
     } ; # fin de tant que
     
     my $date2 = &dateEtHeure() ;
     system("echo $date2 >  ".$out."time.stop") ;
     my $time2 = time() ;
     my $duree = $time2 - $time1 ;
     print " Start: $date1 ; end: $date2 ; elapsed time: $duree sec" ;
     if ($duree>60) {
       print " (" ;
       my $hrs = int($duree/(24*60.0)) ;
       if ($hrs>0) {
         print "$hrs h " ;
       } ; # fin si
       my $rst = $duree - (24*60.0)*$hrs ;
       my $min = int($rst / 60.0) ;
       my $sec = $rst - 60* $min ;
       print "$min min $sec sec)" ;
     } ; # fin si
     print ".\n" ;
     exit(-1) ;
     
     } ; # end of step
     
     ####################################################################
     
     sub logfiles {
     
     ####################################################################
     
     return( glob('*.logfile') ) ;
     
     } ; # end of logfles
     
     ####################################################################
     
     sub noStep {
     
     ####################################################################
     
     my $stepn1 = $_[0] ;
     my $stepn2 = $stepn1 + 1 ;
     
     print "\n  !! step $stepn1 was uncomplete." ;
     print "\n  impossible to execute step $stepn2. STOP.\n\n" ;
     exit(-3) ;
     
     } ; # end of noStep
     
     ####################################################################
     
     sub writeOkStep {
     
     ####################################################################
     
     my $step = $_[0] ;
     
     my $fn   = $out."step".$step.".ok" ;
     
     open(FILE,">$fn") or die ("\n\n! Unable to write to $fn\n\n") ;
     print FILE "ok\n" ;
     close(FILE) ;
     
     print "\n -- step $step completed.\n\n" ;
     
     } ; # end of writeOkStep
     
     ####################################################################
     
     sub check_mothur_version {
     
     ####################################################################
     
     my $required = "1.33" ;
     my $moth     = 'mothur v' ;
     my $cmd      = '(cd $out ; mothur version.mothur | grep "mothur v" )  ' ;
     
     print $cmd."\n" ;
     
     my $mothurv = `$cmd` ;
     print $mothurv."\n" ;
     my $lngthmv = 1+length($moth) ;
     $mothurv = substr($mothurv,$lngthmv,4) ;
     if ($mothurv<$required) {
         print "\n" ;
         print " mothur version >= $required required. \n" ;
         print " your version is $mothurv. please update mothur. STOP.\n\n" ;
         exit(-4) ;
     } # fin si
     
     } ; # end of check_mothur_version
     
     ####################################################################
     
     sub check_Files {
     
     ####################################################################
     
     # create the directories if not present
     
     system("mkdir Input")   unless (-e "Input/") ;
     system("mkdir Output")  unless (-e "Output/") ;
     system("mkdir Results") unless (-e "Results/") ;
     
     # first delete the 8 step files if present
     
     for (my $file=0;$file<=7;$file++) {
       my $fn = $out."step".$file.".ok" ;
       if (-e $fn) {
         unlink($fn) ;
       } ; # end if
     } ; # end for file
     
     # next delete the 4 mothur step files if present
     
     for (my $file=1;$file<=4;$file++) {
       my $fn = $out."mothur".$file.".ok" ;
       if (-e $fn) {
         unlink($fn) ;
       } ; # end if
     } ; # end for file
     
     # next delete the 3 R  step files if present
     
     for (my $file=1;$file<=3;$file++) {
       my $fn = $out."rscript".$file.".ok" ;
       if (-e $fn) {
         unlink($fn) ;
       } ; # end if
     } ; # end for file
     
     # check if the files needed by mothur are present
     
     my $nbe   = 0 ; # number of errors
     
     my @listf = ("test16S.oligos","silva.v4.fasta","trainset9_032012.pds.fasta","trainset9_032012.pds.tax",
                  "quality.mothur","groupsMOCK.mothur","groupsNOMOCK.mothur",
                  "alpha-divMOCK.mothur","alpha-divNOMOCK.mothur",
                  "beta-div.mothur","amova.mothur",
                  "elim_LowFreqOTU_from_shared_make_database.pl","calc_abundance_count.pl",
                  "beanplots.r","taxonomy.r","nmdspcoa.r") ;
     
     foreach my $filename (sort @listf) {
       my $fn = "$pathf"."$filename" ;
       print " required file $filename is " ;
       if (-e $fn) {
         print "present" ;
         unless (-e $filename) {
           my $cmd = "cp ".$pathf."$filename ." ;
           system($cmd) ;
         } ; # end if
       } else {
         print "absent. ERROR." ;
         print "(absolute file pathname : $fn)\n" ;
         $nbe++ ;
       } ; # end if
       print "\n" ;
     } ; # end of foreach
     
     # if files are missing, report it and stop
     
     if ($nbe>0) {
       print "\nsome files are missing, stop.\n" ;
       exit(-1) ;
     } ; # end if
     
     # we have to copy some files to the Output directory for motur
     
     my @listf2 = ("test16S.oligos","silva.v4.fasta","trainset9_032012.pds.fasta","trainset9_032012.pds.tax") ;
     
     foreach my $filename (sort @listf2) {
       my $fn = "$pathf"."$filename" ;
       if (-e $fn) {
         unless (-e $out.$filename) {
           my $cmd = "cp ".$pathf."$filename $out" ;
           system($cmd) ;
         } ; # end unless
       } ; # end if
     } ; # end of foreach
     
     # let's detect MOCK Mock and mock files
     
     my $mock = 0 ;
     my @fastq = glob($in."*fastq") ;
     foreach my $nameFq (sort @fastq) {
       if ($nameFq =~ /mock/i) {
         $mock++ ;
         print " MOCK file #mock found: $nameFq\n" ;
       } ; # fin si
     } ; # end of foreach
     
     if ($mock==0) {
         print " no MOCK files found.\n" ;
         system("echo NOMOCK > mf.option" )
     } else {
         system("echo MOCKFILES > mf.option" )
     } ; # fin si
     
     } ; # end of check_Files
     
     ####################################################################
     
     sub build_stability_files {
     
     ####################################################################
     
     print "\n" ;
     
     # then build a list for "*R1.fastq" files
     # if a "*R2.fastq" correspond, make an entry in stability.files
     
     my $fs = $out."stability.files" ;
     open(FILE,">$fs") or die ("\n\n! Unable to write to $fs\n\n") ;
     print " writing to $fs \n" ;
     
     # *.R1.fastq and *.R2.fastq OR *.R1_001.fastq and *.R2_001.fastq ?
     
     my $spec1 = "R1_001.fastq" ;
     my $spec2 = "R2_001.fastq" ;
     my $lensp = length($spec1) ;
     my @fastq = glob($in."*$spec1") ;
     
     my $nbs = 0 ; # number of files
     foreach my $nameR1 (sort @fastq) {
     #print " A - nameR1 : $nameR1 \n" ;
       my $common  = substr($nameR1,0,length($nameR1)-$lensp) ;
       my $commonS = substr($common,length($in)) ;
       if ($nameR1 =~ /mock/i) { $commonS = "MOCK" ; } ;
       my $nameR2  = $common.$spec2 ;
     #print " A - nameR2 $nameR2  \n\n" ;
     #print " A - common $commonS \n\n" ;
       if (-e $nameR2) {
         $nbs++ ;
         print FILE "$commonS "."../".$nameR1." ../".$nameR2."\n" ;
       } ; # end if
     } ; # end of foreach
     
     close(FILE) ;
     if ($nbs>0) { print "\n $nbs files written in $fs for *$spec1 files.\n" ; } ;
     
     if ($nbs==0)  { # sometimes, files are named R1.fastq instead of R1_001.fastq
     
     open(FILE,">$fs") or die ("\n\n! Unable to write to $fs\n\n") ;
     
     $spec1 = "R1.fastq" ;
     $spec2 = "R2.fastq" ;
     $lensp = length($spec1) ;
     @fastq = glob($in."*$spec1") ;
     
     $nbs = 0 ; # number of files
     foreach my $nameR1 (sort @fastq) {
     #print " B - nameR1 : $nameR1 \n" ;
       my $common  = substr($nameR1,0,length($nameR1)-$lensp) ;
       my $commonS = substr($common,length($in)) ;
       if ($nameR1 =~ /mock/i) { $commonS = "MOCK" ; } ;
       my $nameR2 = $common.$spec2 ;
     #print " B - nameR2 $nameR2  \n\n" ;
     #print " B - common $commonS \n\n" ;
       if (-e $nameR2) {
         $nbs++ ;
         print FILE "$commonS "."../".$nameR1." ../".$nameR2."\n" ;
       } ; # end if
     } ; # end of foreach
     
     close(FILE) ;
     if ($nbs>0) { print " $nbs files written in $fs for *$spec1 files.\n" ; } ;
     
     } # fin si
     
     # if stability.files is empty, report it and stop
     
     if ($nbs==0) {
       print " stability.files not found. Check if *.fastq files are present. STOP.\n" ;
       exit(-1) ;
     } else {
       print "\n -- $nbs stability.files *R1* and *R2* used.\n" ;
     
       # now, build the file 16S.design which corresponds to the factors
       # ignore the part before the first minus or _  sign and form ATCAG etc.
       # use stability.files as the input source
     
     ##  open(FILE,"<$fs") or die ("\n\n! Unable to read from $fs\n\n") ;
     ##  my %series   ;
     ##  my %facteurs ;
     ##  while (my $lig=<FILE>) {
     ##    # don't use MOCK* files
     ##    unless ($lig =~ /MOCK/i) {
     ##    print " LIGNE $lig " ;
     ##       my @mots = split(/\s+/,$lig) ;
     ##       my $prefactor = $mots[0] ;
     ##       $prefactor =~ /^(.*?)[_-](.*?)_/ ;
     ##    print " PREFACTEUR  $prefactor \n " ;
     ##       my $serie   = $1 ;
     ##       my $facteur = $2 ;
     ##       $facteurs{$facteur}++ ;
     ##       $series{$serie}++ ;
     ##       #print " FACTEUR  $facteur \n " ;
     ##    } ; # fin si
     ##  } ; # fin de tant que
     ##  close(FILE) ;
     ##
     
     } ; # end if
     
     # useful ? at last, make and return a list of all files named *.logfiles
     
     return( &logfiles() ) ;
     
     } ; # end of build_stability_files
     
     ####################################################################
     
     sub elim_LowFreqOTU {
     
     ####################################################################
     
     #perl elim_LowFreqOTU_from_shared_make_database.pl -s 16S.an.unique_list.0.03.pick.shared -f 0.1 -id 16S.an.unique_list.database -o 16S.an.unique_list.abund.shared -od 16S.an.unique_list.abund.database
     #
     print "\n-- elim_LowFreqOTU \n\n" ;
     
     my $inputf  ;
     my $outfile ;
     my $odfile  ;
     
     my $idfile   = "16S.an.unique_list.database" ;
     
     # input
     my $inputf1  = "16S.an.unique_list.shared" ;
     my $inputf2  = "16S.an.unique_list.0.03.pick.shared" ;
     
     # output
     my $outfile1 = "16S.an.unique_list.abund.shared" ;
     my $odfile1  = "16S.an.unique_list.abund.database " ;
     
     my $outfile2 = "16S.an.unique_list.0.03.pick.abund.shared" ;
     #my $odfile2  = "16S.an.unique_list.0.03.pick.abund.database " ;
     my $odfile2  = "16S.an.unique_list.abund.database " ;
     
     my $mf = `cat mf.option` ;
     chomp($mf) ;
     my $gm  = "?" ;
     if ($mf eq "NOMOCK") {
        $inputf  = $inputf1 ;
        $outfile = $outfile1 ;
        $odfile  = $odfile1  ;
     } else {
        $inputf  = $inputf2 ;
        $outfile = $outfile2 ;
        $odfile  = $odfile2  ;
     } ; # fin si
     
     unless (-e  $out.$idfile) {
       print "\n file $idfile is missing, end of step.\n" ;
       exit(-1) ;
     } ; # end if
     
     #perl elim_LowFreqOTU_from_shared_make_database.pl -s 16S.an.unique_list.0.03.pick.shared -f 0.1 -id 16S.an.unique_list.database -o 16S.an.unique_list.abund.shared -od 16S.an.unique_list.abund.database
     #R
     
     if (-e  $out.$inputf) {
       my $cmd  = "perl ../elim_LowFreqOTU_from_shared_make_database.pl -s $inputf -f 0.1 -o $outfile -id $idfile -od $odfile " ;
       print "$cmd\n";
       system("( cd Output ; $cmd ) ") ;
     } else {
       print "\n file $inputf is missing, end of step.\n" ;
       exit(-1) ;
     } ; # end if
     
     # at last, make and return a list of all files named *.logfiles
     
     return( &logfiles() ) ;
     
     } ; # end of elim_LowFreqOTU
     
     ####################################################################
     
     sub calc_abundance {
     
     ####################################################################
     
     print "\n-- calc_abundance  \n\n" ;
     
     my $inputf  = "" ;
     my $outfil  = "" ;
     
     my $inputf1 = "16S.an.unique_list.abund.shared" ;
     my $outfil1 = "16S.an.unique_list.abund.proportion.shared" ;
     my $inputf2 = "16S.an.unique_list.0.03.pick.abund.shared" ;
     my $outfil2 = "16S.an.unique_list.0.03.pick.abund.proportion.shared" ;
     
     my $mf = `cat mf.option` ;
     chomp($mf) ;
     my $gm  = "?" ;
     if ($mf eq "NOMOCK") {
        $inputf = $inputf1 ;
        $outfil = $outfil1 ;
     } else {
        $inputf = $inputf2 ;
        $outfil = $outfil2 ;
     } ; # fin si
     
     #perl calc_abundance_count.pl -s 16S.an.unique_list.abund.shared -o 16S.an.unique_list.abund.proportion.shared -c 16S.an.unique_list.abund.proportion.count_table -f 10000
     
     if (-e  $out.$inputf) {
       my $cmd ;
       #$cmd = "perl ../calc_abundance.pl -s $inputf -f 1000 -o $outfil" ;
       $cmd = "perl ../calc_abundance_count.pl -s $inputf -o 16S.an.unique_list.abund.proportion.shared -c 16S.an.unique_list.abund.proportion.count_table -f 10000 " ;
     
       print "$cmd\n";
       system("( cd Output ; $cmd ) ") ;
     } else {
       print "\n file $inputf is missing, end of step.\n" ;
       exit(-1) ;
     } ; # end if
     
     } ; # end of calc_abundance
     
     ####################################################################
     
     sub report {
     
     ####################################################################
     
     print "\n-- report \n" ;
     
     my @log2 = &logfiles() ;
     
     #system("ls -alt step*ok") ;
     #system("ls -alt mothur*ok") ;
     system("ls -alt *ok") ;
     system("rm -rf stability.trim*") ;
     
     } ; # end of report
     
     ####################################################################
     
     sub findFactors {
     
     ####################################################################
     
     # use Group values of file 16S.an.unique_list.abund.shared
     # to determine which factors could be used
     
     ##my $fn = $out."16S.an.unique_list.abund.shared" ;
     
     my $fn ;
     my $fn1 = $out."16S.an.unique_list.abund.shared" ;
     my $fn2 = $out."16S.an.unique_list.0.03.pick.abund.shared" ;
     
     my $mf = `cat mf.option` ;
     chomp($mf) ;
     my $gm  = "?" ;
     if ($mf eq "NOMOCK") {
        $fn = $fn1 ;
     } else {
        $fn = $fn2  ;
     } ; # fin si
     
     open( FIC ,"<$fn") || die "\n impossible to read  $fn \n\n" ;
     
     # get all the file names
     
     my %pf0  ; # contiendra les associations lettres/facteurs numéro 1
     my %pf1  ; # contiendra les associations lettres/facteurs        2
     my %pf2  ; # contiendra les associations lettres/facteurs        3
     my %pf3  ; # contiendra les associations lettres/facteurs        4
     my %pf4  ; # contiendra les associations lettres/facteurs        5
     
     my $ldf = "" ; # liste des fichiers
     my $lefic ;    # fichier courant
     while (my $lig=<FIC>) {
       my @tdv = split(/\s+/,$lig) ;
     #print  " TDV : ".$tdv[1]."\n" ;
       $lefic  = $tdv[1] ;
     #print "lefic $lefic \n" ;
       if ($tdv[1] ne "Group") {
         $tdv[1] =~ /^(.*?)[_-][ACGT]{4,10}/ ;
         if ($1) {
             my $alldeb = $1 ;
             if (length($alldeb)==1) { $alldeb = $lefic  ; } ;
     #print "donc $alldeb \n" ;
             if ($1 ne "MOCK") { $ldf .= " $alldeb " ; } ;
     
         } else {
             #print "rien ? \n" ;
             $ldf .= " ".$lefic." " ;
         } ; # finsi
       } ; # fin si
     #print " au final : $ldf\n\n" ;
     } ; # fin de tant que
     
     close(FIC) ;
     
     # now try to split each file name and find the differents words used
     
     #print "ldf : $ldf\n" ;
     my @ldf = split(/\s+/,$ldf) ;
     foreach my $fn (@ldf) {
       my @parts = split(/[_-]/,$fn) ;
       if ($parts[0]) {
         $pf0{$parts[0]}++ ;
       } ; # fin si
       if ($parts[1]) {
         $pf1{$parts[1]}++ ;
       } ; # fin si
       if ($parts[2]) {
         $pf2{$parts[2]}++ ;
       } ; # fin si
       if ($parts[3]) {
         $pf3{$parts[3]}++ ;
       } ; # fin si
       if ($parts[4]) {
         $pf4{$parts[4]}++ ;
       } ; # fin si
     } ; # fin pour chaque
     
     my $nbf = -1 ;
     my @rl  = () ; # returned list
     my $lf ;
     
     if ((keys %pf0)>1) {
       $nbf++ ;
       $lf = join(" ",sort keys %pf0) ;
       $rl[$nbf] = $lf ;
     } ; # fin si
     
     if ((keys %pf1)>1) {
       $nbf++ ;
       $lf = join(" ",sort keys %pf1) ;
       $rl[$nbf] = $lf ;
     } ; # fin si
     
     if ((keys %pf2)>1) {
       $nbf++ ;
       $lf = join(" ",sort keys %pf2) ;
       $rl[$nbf] = $lf ;
     } ; # fin si
     
     if ((keys %pf3)>1) {
       $nbf++ ;
       $lf = join(" ",sort keys %pf3) ;
       $rl[$nbf] = $lf ;
     } ; # fin si
     
     if ((keys %pf4)>1) {
       $nbf++ ;
       $lf = join(" ",sort keys %pf4) ;
       $rl[$nbf] = $lf ;
     } ; # fin si
     
     return(@rl)
     
     } ; # end of findFactors
     
     ####################################################################
     
     sub showFactors {
     
     ####################################################################
     
     print "\n" ;
     
     my $nbf = 0 ;
     foreach my $lf (@_) {
       $nbf++ ;
       if ($nbf==1) {
         print " for factors in --step 5, you may use the option(s) \n" ;
       } ; # fin si
       print " --factor $nbf : $lf \n" ;
     } ; # fin pour chaque
     
     if ($nbf==0) {
     print " No factor found. Stop.\n\n" ;
     exit(-2) ;
     } ; # fin si
     
     } ; # end of showFactors
     
     ####################################################################
     
     sub possibleFactors {
     
     ####################################################################
     
     # use Group values of file 16S.an.unique_list.abund.shared
     # to determine which factors could be used
     
     print "\n" ;
     my $fn ;
     my $fn1 = $out."16S.an.unique_list.abund.shared" ;
     my $fn2 = $out."16S.an.unique_list.0.03.pick.abund.shared" ;
     
     my $mf = `cat mf.option` ;
     chomp($mf) ;
     my $gm  = "?" ;
     if ($mf eq "NOMOCK") {
        $fn = $fn1 ;
     } else {
        $fn = $fn2  ;
     } ; # fin si
     
     print " Looking for factors in $fn ...\n" ;
     
     unless(-e $fn) {
       print "\n file $fn is missing, impossible to determine factors.\n" ;
       print " Maybe you have not yet executed step 4.\n" ;
       exit(-1) ;
     } else {
       my @lfp = &findFactors() ;
     #print " findFactors renvoie : ".join(" ! ",@lfp)."\n" ;
       &showFactors(@lfp) ;
     } ; # fin si
     
     } ; # end of possibleFactors
     
     ####################################################################
     
     sub buildForFactor {
     
     ####################################################################
     
     my $factor ;
     
     # which factor to use ? default is 1
     
       $_[0]  =~ /([1-7])/ ;
       if ($1) {
         $factor = $1 ;
       } else {
         $factor = 1 ;
       } ; # fin si
     
     #  if unless ($1) {
     #    print "\n no factor number specified. Stop.\n\n" ;
     #    exit(-5) ;
     #
     #  } else {
     
         # look for the possible factors
     
         my @lfp = &findFactors() ;
     #print " findFactors renvoie : ".join(" * ",@lfp)."\n" ;
     
         # check the factor number is OK
     
         my $ndf = $factor ;
         my $nbfa = 1+ $#lfp ;
         if ($ndf>$nbfa) {
            print "\n" ;
            print " there are at most $nbfa factor(s), so you can not use factor #$ndf.\n\n" ;
            exit(-5) ;
         } ; # fin si
     
         # then proceed
     
         my $fnd = $out."16S.design" ;
         print "\n " ;
         print "let's use factor $ndf to build $fnd ; \n" ;
         my $fd = $out."16S.design" ;
         open(FILE,">$fd") or die ("\n\n! Unable to write to $fd\n\n") ;
         print FILE "File                   Factor    \n" ;
         my @mdf = split(" ",$lfp[$ndf-1]) ;
         print " values: ".$lfp[$ndf-1]."\n" ;
     
         # my $fn = $out."16S.an.unique_list.abund.shared" ;
         my $fn ;
         my $fn1 = $out."16S.an.unique_list.abund.shared" ;
         my $fn2 = $out."16S.an.unique_list.0.03.pick.abund.shared" ;
     
         my $mf = `cat mf.option` ;
         chomp($mf) ;
         my $gm  = "?" ;
         if ($mf eq "NOMOCK") {
            $fn = $fn1 ;
         } else {
            $fn = $fn2  ;
         } ; # fin si
     
     
         open( FIC ,"<$fn") || die "\n impossible to read  $fn \n\n" ;
         my $ldf = "" ; # liste des fichiers
         my $nmoda = "" ;
         my $alldeb ;
         my $vu ;
         while (my $lig=<FIC>) {
     #print "\nline $lig \n" ;
           $vu = 0 ;
           my @tdv = split(/\s+/,$lig) ;
           my $fichier = $tdv[1]." " ;
           if ($fichier eq "Group ") {
             $vu = 1 ;
           } else {
              $vu = 0 ;
              $fichier =~ /^(.*?)[_-][ACGT]{4,10}/ ;
              if ($1) {
               if ($1 =~ /MOCK/i) {
                 $vu = 1 ;
               } ; # finsi
               $alldeb = $1 ;
               if (length($alldeb)==1) { $alldeb = $fichier ;  } ;
              } else {
               $alldeb = $fichier ;
              } ; # finsi
     
     #print "\nprocessing $alldeb \n" ;
                 # boucle à réécrire en tantr que
              my $nmoda = "" ;
              foreach my $moda (@mdf) {
     #print "processing $alldeb $moda \n" ;
                   if (index($alldeb,$moda)>-1) {
                     $vu++ ;
                     $nmoda = $moda ;
     #print  "VU !! $fichier       $nmoda\n" ;
                     print FILE "$fichier       $nmoda\n" ;
                     last ;
                   } ; # fin si
              } ; # fin pour chaque modalité
     
     
              if ($vu==0) {
                print " the file *$fichier* has no referenced value for your factor.\n" ;
     } else {
     #print " for $fichier referenced value = $nmoda \n" ;
              } ; # fin de si
           } ; # finsi
     
         } ; # fin de tant que
     
         close(FIC) ;
         close(FILE) ;
         #print "\n you may now use $fnd within --step 4 \n" ;
     
         # write factor number to file Output/current.factor for R to use it
         system("echo $ndf > ".$out."current.factor")
     
     #  } ; # fin si
     
     } ; # end of buildForFactor
     
     ####################################################################
     
     sub dateEtHeure {
     
     ####################################################################
     
         my ($sec,$min,$hour,$mday,$mmon,$year)=localtime();
         $mmon = $mmon + 1 ;
         $year = $year + 1900 ;
         if (length($sec)<2)  { $sec = "0$sec"   } ;
         if (length($mday)<2) { $mday = "0$mday" } ;
         if (length($mmon)<2) { $mmon = "0$mmon" } ;
         my $now  = $mday."/".$mmon."/".$year;
         $now .= " ".$hour.":".$min ;
     
         return $now ;
     
     } ; # fin sub dateEtHeure
     
     ####################################################################
     
     sub rscripts {
     
     ####################################################################
     
     my $factor = $_[0] ;
     my $rstepr = $_[1] ;
     
     # process 16S.an.unique_list.0.03.cons.taxonomy
     # in order to read it with R
     
     my $fin = $out."16S.an.unique_list.0.03.cons.taxonomy" ;
     open( FIN ,"<$fin") || die "\n impossible to read file $fin \n\n" ;
     
     my $fout = $out."16Staxons.txt" ;
     open( FOUT ,">$fout") || die "\n impossible to write file $fout \n\n" ;
     
     my $nbl = 0 ;
     while (my $lig=<FIN>) {
       $nbl++ ;
       my @mots = split(/;/,$lig) ;
       my $nbmots = $#mots ;
       #print " ligne $nbl avec $nbmots mots \n" ;
       if ($#mots<1) { next ; } ;
       #print "AVANT $lig" ;
       my $ligApres = $lig ;
       $ligApres =~ s/\([0-9]+\)//g ;
       $ligApres =~ s/"//g ;
       $ligApres =~ s/;/\t/g ;
       #print "APRES $ligApres" ;
       #print "\n" ;
       # if ($nbl>5) { exit(-1) ; } ;
       print FOUT $ligApres ;
     } ; # fin tant que
     
     close(FIN) ;
     close(FOUT) ;
     
     
     if ($rstepr eq "step5") {
     
        # run R script beanplot.r
     
        &rscript("beanplots",$factor) ;
     
        # run R script taxonomy.r
     
        &rscript("taxonomy",$factor) ;
     } ; # fin si
     
     if ($rstepr eq "step6") {
     
        # run R script nmdspcoa.r
     
        &rscript("nmdspcoa",$factor) ;
     } ; # fin si
     
     } ; # fin sub rscripts
     
     ####################################################################
     
     sub rscript {
     
     ####################################################################
     
     my $nomscr = $_[0] ;
     my $factor = $_[1] ;
     
     my $cmd  ;
     my $resr ;
     my $flog ;
     
     $flog = "../".$res.$nomscr."_F".$factor."_log.txt" ;
     $cmd  = " R --vanilla < ../$nomscr.r > $flog 2> $flog " ;
     print   " R --vanilla < ../$nomscr.r  > $flog 2> $flog \n" ;
     $resr = `(cd Output ; $cmd)` ;
     
     ##
     ##if ($nomscr eq "beanplots") {
     ##    $cmd = " echo ok > ".$out."rscript1.ok" ;
     ##    print($cmd."\n") ;
     ##    system($cmd) ;
     ##} ; # fin si
     ##
     ##if ($nomscr eq "taxonomy") {
     ##    $cmd = " echo ok > ".$out."rscript2.ok" ;
     ##    print($cmd."\n") ;
     ##    system($cmd) ;
     ##} ; # fin si
     ##
     ##if ($nomscr eq "nmdspcoa") {
     ##    $cmd = " echo ok > ".$out."rscript3.ok" ;
     ##    print($cmd."\n") ;
     ##    system($cmd) ;
     ##} ; # fin si
     
     
     } ; # fin sub rscript
     
     ####################################################################
     ####################################################################
     
     # end of file metauto.pl
     
     ####################################################################
     ####################################################################
     

 

 

retour gH    Retour à la page principale de   (gH)