@@ -8,6 +8,7 @@ use lib "$Bin";
88use lib " $RealBin /../lib" ;
99use lib " $RealBin /../ext/lib/perl5" ;
1010use File::Basename;
11+ use Time::Piece;
1112use PhaME;
1213use Cwd;
1314
@@ -49,13 +50,20 @@ AUTHORS
4950 Los Alamos National Laboratory
5051=cut
5152
53+ # ...............................................................................
54+ # global variables
55+ my $PACKAGENAME = " PhaME" ;
56+ my $PACKAGEVERSION = " 1.0.0" ;
57+ my $AUTHORS = " Migun Shakya, Chienchi Lo, Sanaa Ahmed" ;
58+ my $URL = ' https://github.com/LANL-Bioinformatics/PhaME' ;
59+ my $vcheck = ' ' ;
5260my $refdir ; # variable storing directory with reference sequence
5361my $workdir ; # variable storing directory where all analysis will be done
5462my $outdir ;
5563my $project ;
5664my $rsignal = 0;
5765my $reference ;
58- my $aligner = " bowtie" ;
66+ my $aligner = " bowtie" ;
5967my $name ;
6068my $path ;
6169my $suffix ;
@@ -91,17 +99,80 @@ my $mappingGaps;
9199my $ptree ;
92100my $buildSNPdb = 0;
93101my $message ;
94- my $PACKAGENAME =" PhaME" ;
95- my $PACKAGEVERSION =" 1.0.0" ;
96102my $bindir = getBinDirectory();
97103
104+ # ...............................................................................
105+ # command line options
106+
98107GetOptions(
108+ ' vcheck' => sub {check_version()},
99109 ' version|v' => sub { version() },
100- ' help|h' => sub { printusage() },
110+ ' help|h' => sub { printusage() }
111+
101112);
102113
114+ if ( !@ARGV ) {
115+ print (@ARGV );
116+ print " $0 : Argument required.\n " ;
117+ print " For help message use:\n " ;
118+ print " runPhaME -h.\n " ;
119+ exit 1;
120+ }
121+
103122my $control = $ARGV [0] || " phame.ctl" ;
104- # #########################################################################################
123+
124+ # ...............................................................................
125+ # check for depenencies
126+
127+ sub parse_version {
128+
129+ # from snippy
130+ my ( $cmd , $min , $re ) = @_ ;
131+ open my $ver , ' -|' , $cmd . ' 2>&1'
132+ or err (" Could not open pipe to '$cmd '" );
133+ my $blob = join ' ' , <$ver >;
134+ err (" Could not read output of '$cmd '" ) unless $blob ;
135+
136+ # msg("parse=[$blob]");
137+ $blob =~ $re ;
138+ $ver = defined $1 ? $1 : 0;
139+ err (" Need $cmd >= $min but you have $ver - please upgrade it." )
140+ if $ver < $min ;
141+ msg(" Checking version: $cmd is >= $min - ok, have $ver " );
142+ return $ver ;
143+ }
144+
145+ sub check_version {
146+ # quit now if --vcheck was provided
147+ parse_version( ' samtools --version' , 1.3, qr / samtools\s (\d +\.\d +)/ ms );
148+ parse_version( ' bcftools --version' , 1.3, qr / bcftools\s (\d +\.\d +)/ ms );
149+ parse_version( ' nucmer --version' , 3.1, qr / version\s (\d +\.\d +)/ ms );
150+ parse_version( ' bowtie2 --version' , 2.3, qr / version\s (\d +\.\d +)/ ms );
151+ parse_version( ' bwa' , 0.7, qr / Version:\s (\d +\.\d +)/ ms );
152+ parse_version( ' FastTree' , 2.1, qr / version\s (\d +\.\d +)/ ms );
153+ parse_version( ' bbmap.sh -v' , 38.06, qr / Version\s (\d +\.\d +)/ ms );
154+ parse_version( ' raxmlHPC-PTHREADS -version' , 8.2, qr / version\s (\d +\.\d +)/ ms );
155+ msg(" Dependences look good!" );
156+ exit (0);
157+ }
158+
159+ # ----------------------------------------------------------------------
160+
161+ sub msg {
162+ my $t = localtime ;
163+ my $line = " [" . $t -> hms . " ] @_ \n " ;
164+ print STDERR $line ;
165+ }
166+
167+ # ----------------------------------------------------------------------
168+
169+ sub err {
170+ msg(@_ );
171+ exit (2);
172+ }
173+
174+ # ...............................................................................
175+
105176# # Read in control file
106177open ( CTL, " $control " ) || die " Please provide a control file" ;
107178
@@ -124,16 +195,18 @@ while (<CTL>) {
124195 if (/ project\s *=\s *(\S +)\s *#{0,1}.*$ / ) { $project = $1 ; }
125196
126197 if (/ reference\s *=\s *(0|1|2)\s *#{0,1}.*$ / ) { $rsignal = $1 ; }
127- if (/ aligner\s *=\s *(\S +)\s *#{0,1}.*$ / ){ $aligner = $1 ;}
198+ if (/ aligner\s *=\s *(\S +)\s *#{0,1}.*$ / ) { $aligner = $1 ; }
128199 if ( $rsignal == 1 && / reffile\s *=\s *(\S +)\s *#{0,1}.*$ / ) {
129200 $reference = " $refdir /$1 " ;
130201 }
131- elsif ($rsignal ==0){
132- opendir (DIR, $refdir );
133- my @reffiles = grep { / .[fna|fa|fasta|fsa]$ / && -f " $refdir /$_ " } readdir (DIR);
134- $reference = " $refdir /" . $reffiles [int (rand (scalar (@reffiles )))];
135- closedir DIR;
136- }
202+ elsif ( $rsignal == 0 ) {
203+ opendir ( DIR, $refdir );
204+ my @reffiles
205+ = grep { / .[fna|fa|fasta|fsa]$ / && -f " $refdir /$_ " } readdir (DIR);
206+ $reference
207+ = " $refdir /" . $reffiles [ int ( rand ( scalar (@reffiles ) ) ) ];
208+ closedir DIR;
209+ }
137210
138211 if (/ cdsSNPS\s *=\s *(0|1)\s *#{0,1}.*$ / ) { $gsignal = $1 ; }
139212
@@ -214,7 +287,6 @@ elsif ( $rsignal == 1 ) {
214287
215288( $name , $path , $suffix ) = fileparse( " $reference " , qr /\. [^.]*/ );
216289
217-
218290print (" \n " );
219291&print_timeInterval( " \n $runtime \n " , " Loading information\n " );
220292print " \t Refd:\t $refdir \n " ;
@@ -226,18 +298,17 @@ if ( !-d $refdir || !-d $workdir || !-d $outdir ) {
226298}
227299
228300print " \t Reference:\t $reference \n " ;
229- if ( !-e " $reference " )
230- { print " File $reference does not exist.\n " ;
301+ if ( !-e " $reference " ) {
302+ print " File $reference does not exist.\n " ;
231303 exit ;
232- }
233-
304+ }
234305
235306if ( $gsignal == 1 ) {
236307 if ( -e " $refdir /$name .gff" ) {
237308 $annotation = " $refdir /$name .gff" ;
238309 print " \t Annotation:\t $annotation \n " ;
239310 }
240- elsif (-e " $refdir /$name .gff3" ) {
311+ elsif ( -e " $refdir /$name .gff3" ) {
241312 $annotation = " $refdir /$name .gff3" ;
242313 print " \t Annotation:\t $annotation \n " ;
243314 }
@@ -257,11 +328,11 @@ END_MESSAGE
257328}
258329
259330if ( $pselection > 0 ) {
260- if ( -e " $refdir /$name .gff" ){
331+ if ( -e " $refdir /$name .gff" ) {
261332 $genefile = " $refdir /$name .gff" ;
262333 print " \t Genes:\t $genefile \n " ;
263334 }
264- elsif ( -e " $refdir /$name .gff3" ){
335+ elsif ( -e " $refdir /$name .gff3" ) {
265336 $genefile = " $refdir /$name .gff3" ;
266337 print " \t Genes:\t $genefile \n " ;
267338 }
@@ -519,7 +590,7 @@ if ( $check == 0 ) {
519590 $contig_list {$list } = $genome_size ;
520591 }
521592 }
522- if ( $read_mapping == 1 ) {
593+ if ( $read_mapping == 1 ) {
523594 if ( $files =~ / .+\. f{1}a?s?t?q$ /
524595 && $files !~ / .+\. f{1}n|s?a?s?t?a$ /
525596 && $files !~ / .+\. contig$| \. contigs$| \. ctg$] $ / )
@@ -620,16 +691,20 @@ if ( $read_mapping == 1 ) {
620691 " map" , $project );
621692 PhaME::removeGaps( $bindir , $reference , $mappingGaps );
622693 &print_timeInterval( $runtime , " Mapping reads to reference\n " );
623- my $end = PhaME::readsMapping( $workdir , $bindir ,
624- " $workdir /reads_list.txt" , $threads , $name , $error , $aligner , $logfile );
694+ my $end
695+ = PhaME::readsMapping( $workdir , $bindir ,
696+ " $workdir /reads_list.txt" , $threads , $name , $error , $aligner ,
697+ $logfile );
625698 &print_timeInterval( $runtime , " $end \n " );
626699 }
627700 elsif ( $time == 1 ) {
628701 my $tempdir = $outdir . ' /temp/' ;
629702 ` mkdir -p $tempdir ; cp $reference $tempdir ` ;
630703 &print_timeInterval( $runtime , " Mapping reads to reference\n " );
631- my $end = PhaME::readsMapping( $workdir , $bindir ,
632- " $workdir /reads_list.txt" , $threads , $name , $error , $aligner , $logfile );
704+ my $end
705+ = PhaME::readsMapping( $workdir , $bindir ,
706+ " $workdir /reads_list.txt" , $threads , $name , $error , $aligner ,
707+ $logfile );
633708 &print_timeInterval( $runtime , " $end \n " );
634709 }
635710}
@@ -638,7 +713,7 @@ if ( $buildSNP == 1 ) {
638713 &print_timeInterval( $runtime , " Identifying SNPs\n " );
639714 &print_timeInterval( $runtime , " Starting with gaps\n " );
640715 PhaME::identifyGaps( $outdir , " $workdir /working_list.txt" , $name , " snp" ,
641- $project );
716+ $project );
642717 &print_timeInterval( $runtime , " Gap identification complete\n " );
643718 if ( $gsignal == 1 ) {
644719 &print_timeInterval( $runtime , " Preparing to identify SNPs\n " );
@@ -654,31 +729,30 @@ if ( $buildSNP == 1 ) {
654729 " $workdir /working_list.txt" , $project , $gsignal , $error , $logfile ,
655730 $cutoff );
656731 if ( $gsignal == 1 ) {
657- print
658- " \t SNPs from contigs and/or raw reads will be parsed\n " ;
659- print
660- " \t to synonymous and non-synonymous if given\n\n " ;
732+ print " \t SNPs from contigs and/or raw reads will be parsed\n " ;
733+ print " \t to synonymous and non-synonymous if given\n\n " ;
661734 my ( $genname , $genpath , $gensuffix )
662735 = fileparse( " $annotation " , qr /\. [^.]*/ );
663- my $cds_gff = $outdir . " /$genname \_ cds.gff" ;
736+ my $cds_gff = $outdir . " /$genname \_ cds.gff" ;
664737 my $snps_folder = $outdir . " /snps" ;
665- my $ref_fasta = " $workdir /files/$genname .fna" ;
666- PhaME::SNPsAnalysis($cds_gff , $snps_folder , $ref_fasta , $logfile , $error );
738+ my $ref_fasta = " $workdir /files/$genname .fna" ;
739+ PhaME::SNPsAnalysis( $cds_gff , $snps_folder , $ref_fasta , $logfile ,
740+ $error );
667741 }
668-
742+
669743 &print_timeInterval( $runtime , " $end \n " );
670744}
671745
672746if ( $buildtree == 1 || $bs == 1 ) {
673- &print_timeInterval($runtime , " Preparing all SNP phylogeny\n " );
747+ &print_timeInterval( $runtime , " Preparing all SNP phylogeny\n " );
674748 my $end = PhaME::buildTree( $bindir , $outdir , $threads , $tree ,
675749 " $project \_ all" , $error , $logfile );
676- &print_timeInterval( $runtime , " DONE\n " );
750+ &print_timeInterval( $runtime , " DONE\n " );
677751 if ( $gsignal == 1 ) {
678- &print_timeInterval( $runtime , " Preparing CDS phylogeny\n " );
752+ &print_timeInterval( $runtime , " Preparing CDS phylogeny\n " );
679753 PhaME::buildTree( $bindir , $outdir , $threads , $tree , " $project \_ cds" ,
680754 $error , $logfile );
681- &print_timeInterval( $runtime , " DONE\n " );
755+ &print_timeInterval( $runtime , " DONE\n " );
682756 }
683757 &print_timeInterval( $runtime , " $end \n " );
684758 if ( $bsignal == 1 ) {
@@ -854,13 +928,13 @@ sub version {
854928}
855929
856930sub printusage {
857- print " \n $PACKAGENAME (version $PACKAGEVERSION )\n " ;
858- print " \n USAGE: runPhaME <control file>\n\n " ;
859- print " -h show this help message and exit\n " ;
860- print " -v show version number and exit\n " ;
861- exit ;
862- print " \n " ;
863- print " \n " ;
931+ print " \n $PACKAGENAME (version $PACKAGEVERSION )\n " ;
932+ print " \n USAGE: runPhaME <control file>\n\n " ;
933+ print " -h show this help message and exit\n " ;
934+ print " -v show version number and exit\n " ;
935+ print " --vcheck checks if all depenencies are correct and in path.\n " ;
936+ exit ;
937+ print " \n " ;
938+ print " \n " ;
864939}
865940
866-
0 commit comments