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 : 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 à la page principale de   (gH)