# # (gH)   -_-  metauto.pl  ;  TimeStamp (unix) : 27 Octobre 2014 vers 18:30

##################################################################
##################################################################
##                                                              ##
##  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 metauto.pl (gH) version 1.59\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() ;
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 (long)",
  "run abundance computations (long)",
  "run alpha-diversity analysis, determine factors, make R beanplots and taxonomic clustering",
  "run beta-diversity analysis",
  "run amova analysis") ;

# path to all script files including this metauto.pl

my $pathf = substr($0,0,length($0)-length("metauto.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 : metauto.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) ;  # 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) ; }
        &writeOkStep(6) ;
      } ; # 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

# 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") ;

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] ;

# 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) ;

# run R script beanplot.r

&rscript("beanplots",$factor) ;

# run R script taxonomy.r

&rscript("taxonomy",$factor) ;

} ; # 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

} ; # fin sub rscript

####################################################################
####################################################################

# end of file metauto.pl

####################################################################
####################################################################
