Mercurial > repos > nml > abacas
diff abacas.1.1.pl @ 0:a1fdc6925620 draft default tip
"planemo upload for repository https://github.com/phac-nml/abacas commit f6856961094e89e4cad0ee7df6c2a49bf005e4bf"
| author | nml |
|---|---|
| date | Thu, 21 Nov 2019 12:53:20 -0500 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/abacas.1.1.pl Thu Nov 21 12:53:20 2019 -0500 @@ -0,0 +1,1793 @@ +#!/usr/bin/env perl +# Copyright (C) 2008 Genome Research Limited. All Rights Reserved. +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 2 +# of the License, or (at your option) any later version. +#This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + + +use strict; +use warnings; +use POSIX qw(ceil floor); +use Getopt::Std; + +#------------------------------------------------------------------------------- + +if (@ARGV < 1) { usage();} + +my ($help, $reference, $query_file, $choice, $sen, $seed, $mlfas, $fasta_bin, $avoid_Ns, + $tbx, $min_id, $min_cov, $diff_cov, $min_len,$add_bin_2ps, $pick_primer, + $flank, $chk_uniq,$redoMummer, $is_circular,$escapeToPrimers, $debug, $gaps_2file, $prefix) + =checkUserInput( @ARGV ); + +my $ref_inline; +if ($escapeToPrimers ==1) +{ + pickPrimers ($reference, $query_file, $flank, $chk_uniq); + exit; +} +#BEGIN +#------------------------------------------------------------------------------- +print_header(); + +$ref_inline = Ref_Inline($reference); +#Get length of the reference sequence +my $ref_len = length($ref_inline); +print "PREPARING DATA FOR $choice \n"; + +################### +# Running MUMmer # +################### +my ($path_dir, $run_mum, $path_toPass); +if ($debug) +{ + print "the seed is $seed \n"; + print "RedoMummer= ",$redoMummer."\n"; +} +my @do_mum_return; +if ($redoMummer) +{ + @do_mum_return = doMummer($reference, $query_file, $choice,$sen,$seed,$min_id, $min_cov, $diff_cov,$min_len, $debug, $is_circular) or die "Couldn't run MUMmer\n"; +} + +#################################### +# Processing tiling output # +#################################### +if ($debug) { + print "Do tiling...\n"; +} + +#-------------------------------- +my $mummer_tiling = $do_mum_return[0]; +$path_dir = $do_mum_return[2]; +$path_toPass = $do_mum_return[1]; +################## +#Do Tiling +#------------------------------------------- +doTiling ($mummer_tiling, $path_toPass, $path_dir,$reference, $query_file, $choice, $prefix,$mlfas, $fasta_bin, $avoid_Ns, $ref_len, $gaps_2file, $ref_inline, $add_bin_2ps, $pick_primer, $flank, $chk_uniq, $tbx); + + + + + +########################################################################################################################################################## +#################################Contig ordering ########################################################## +########################### + +sub help +{ + +die <<EOF + +*********************************************************************************** +* ABACAS: Algorithm Based Automatic Contiguation of Assembled Sequences * +* * +* * +* Copyright (C) 2009 The Wellcome Trust Sanger Institute, Cambridge, UK. * +* All Rights Reserved. * +* * +*********************************************************************************** + +USAGE +abacas.pl -r <reference file: single fasta> -q <query sequence file: fasta> -p <nucmer/promer> [OPTIONS] + + -r reference sequence in a single fasta file + -q contigs in multi-fasta format + -p MUMmer program to use: 'nucmer' or 'promer' +OR +abacas.pl -r <reference file: single fasta> -q <pseudomolecule/ordered sequence file: fasta> -e +OPTIONS + -h print usage + -d use default nucmer/promer parameters + -s int minimum length of exact matching word (nucmer default = 12, promer default = 4) + -m print ordered contigs to file in multifasta format + -b print contigs in bin to file + -N print a pseudomolecule without "N"s + -i int mimimum percent identity [default 40] + -v int mimimum contig coverage [default 40] + -V int minimum contig coverage difference [default 1] + -l int minimum contig length [default 1] + -t run tblastx on contigs that are not mapped + -g string (file name) print gaps on reference to file name + -a append contigs in bin to the pseudomolecule + -o prefix output files will have this prefix + -P pick primer sets to close gaps + -f int number of flanking bases on either side of a gap for primer design (default 350) + -R Redo mummer [default 1, use -R 0 to avoid running mummer] + -e Escape contig ordering i.e. go to primer design + -c Reference sequence is circular + +EOF +} +######## +sub usage +{ + +die <<EOF + +USAGE +abacas.pl -r <reference file: single fasta> -q <query sequence file: fasta> -p <nucmer/promer> [OPTIONS] + -r reference sequence in a single fasta file + -q contigs in multi-fasta format + -p MUMmer program to use: 'nucmer' or 'promer' +for contig ordering and primer design + +OR +abacas.pl -r <reference file: single fasta> -q <pseudomolecule/ordered sequence file: fasta> -e 1 +to escape contig ordering and go directly to primer design + +OR +abacas.pl -h for help + +EOF +} +######## +################## +sub print_header +{ + print " +*********************************************************************************** +* ABACAS: Algorithm Based Automatic Contiguation of Assembled Sequences * +* * +* * +* Copyright (C) 2009 The Wellcome Trust Sanger Institute, Cambridge, UK. * +* All Rights Reserved. * +* * +*********************************************************************************** +\n"; +} +######################################### +sub checkUserInput{ + my %options; + getopts('hr:q:p:ds:mbNi:v:V:l:tg:ao:Pf:R:u:ecD', \%options); + + my ($help, $reference, $query_file, $choice, $sen, $seed, $mlfas, $fasta_bin, $avoid_Ns, + $tbx, $min_id, $min_cov, $diff_cov, $min_len,$add_bin_2ps, $pick_primer, + $flank, $chk_uniq,$redoMummer, $is_circular,$escapeToPrimers, $debug, $gaps_2file, $prefix); + + if($options{h}) { + $help = $options{h}; + help(); + } + else{ + $help =0; + } + if ($options{r} && $options{q} ){ + ($reference, $query_file) = ($options{r},$options{q}); + + } + else + { + usage() unless $options{e}; + } + + if ($options{p}) + { + $choice = $options{p}; + unless ($choice eq "nucmer" || $choice eq "promer"){ + print "Unknown MuMmer function\n Please use nucmer or promer\n"; + exit; + } + } + else + { + usage() unless $options{e}; + } + if ($options{e}) #$escapeToPrimers) + { + print_header(); + print "Primer design selected,... escaping contig ordering\n"; + $escapeToPrimers = 1; + $chk_uniq = "nucmer"; + $choice = ""; + } + else + { + $escapeToPrimers = 0; + } + if ($options{d}) {$sen =1;} else {$sen =0;} #print $sen , " ---sen\n"; + #print $options{t}, "\n"; exit; + if($options{t}) {$tbx = 1;} else {$tbx = 0;} #print $tbx, " ---tbx\n"; # + unless ($options{s}){ + if ($choice eq "nucmer"){ + $seed = 12; + } + else{ + $seed =4; + } + } + if ($options{m}){$mlfas =1;} else { $mlfas =0; } #print $mlfas , " ---mlfasta\n"; + if ($options{b}) { $fasta_bin =1;} else {$fasta_bin =0;} + if ($options {N}) {$avoid_Ns =1;} else {$avoid_Ns=0;} + + unless($options{i}) {$min_id =40;} + unless ($options{v}) {$min_cov =40;} + unless ($options{V}) {$diff_cov =1;} + unless ($options {l}) {$min_len = 1;} + if ($options{a}) {$add_bin_2ps = 1; }else {$add_bin_2ps =0;} + if ($options{P}) {$pick_primer=1;} else {$pick_primer =0;} + if ($options{f}) {$flank = $options{f};} else {$flank = 1000;} + if ($options{u}) {$chk_uniq = $options{u}; } else {$chk_uniq = "nucmer";} + + if($options{R}) {$redoMummer = $options{R}; }else {$redoMummer = 1;} + unless ($options{c}) {$is_circular = 0;} + if ($options{g}) {$gaps_2file = $options{g};} else {$gaps_2file ="";} + if($options{o}) {$prefix = $options{o};} else {$prefix = "";} + unless ($options{D}) {$debug =0}; + if ($tbx ==1 && $fasta_bin !=1) + { + print "ERROR: Please use -t -b if you want to run tblastx on contigs in bin\n"; + exit; + } + #print $tbx, "\n"; exit; + return ($help, $reference, $query_file, $choice, $sen, $seed, $mlfas, $fasta_bin, $avoid_Ns, + $tbx, $min_id, $min_cov, $diff_cov, $min_len,$add_bin_2ps, $pick_primer, + $flank, $chk_uniq,$redoMummer, $is_circular,$escapeToPrimers, $debug, $gaps_2file, $prefix); + +} +############# +## get the reference sequence in one line +#-------------------------------------------------- +sub Ref_Inline +{ + my $ref = shift; + open (refFH, $ref) or die "Could not open file\n"; + my $seq =""; + my @r = <refFH>; + my $num_chr =0; + foreach(@r){ + if ($_ =~ /\>/){ + $num_chr +=1; + } + } + if ($num_chr > 1){ + print "\nERROR: Please use a single fasta reference file. You can simply merge chromosomes in to a union fasta file.\n\n"; + exit; + } + shift @r; + foreach(@r){ + chomp; + $seq = $seq.$_; + } + return $seq; +} +################ +# Run mummer +#-------------------------------------------- +sub doMummer +{ + my ($reference, $query_file, $choice, $sen,$seed,$min_id, $min_cov, $diff_cov,$min_len, $debug, $is_circular ) = @_; + + my $df = 'delta-filter'; + my $st = 'show-tiling'; + my $ask = 'which'; + my ($path_toPass, $run_mum); # params to return... + my ($command, $Path, $dir) = checkProg($choice); + my ($run_df, $df_path, $df_dir) = checkProg($df); + my ($run_st, $st_path, $st_dir) = checkProg($st); + my (@running, @deltaRes, @coordsRes); + if ($choice eq "nucmer") + { + if ($sen ==0) + { + @running = `$command --maxmatch -l $seed -p $choice $reference $query_file &> /dev/null`; + @deltaRes = `$run_df -q $choice.delta >$choice.filtered.delta`; + if ($is_circular == 1) + { + @coordsRes = `$run_st -c -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`; + } + else + { + @coordsRes = `$run_st -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`; + + } + } + else + { + @running = `$command -p $choice $reference $query_file &> /dev/null`; + @deltaRes = `$run_df -q $choice.delta >$choice.filtered.delta`; + if ($is_circular ==1) {@coordsRes = `$run_st -c -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`;} + else { @coordsRes = `$run_st -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`;} + } + + } + else + { + if ($sen ==0) + { + @running = `$command --maxmatch -l $seed -x 1 -p $choice $reference $query_file &> /dev/null`; + @deltaRes = `$run_df -q $choice.delta >$choice.filtered.delta`; + if ($is_circular == 1) + { + @coordsRes= `$run_st -c -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`; + } + else + { + @coordsRes= `$run_st -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`; + } + } + else + { + @running = `$command -l $seed -p $choice $reference $query_file &> /dev/null`; + @deltaRes = `$run_df -q $choice.delta >$choice.filtered.delta`; + if ($is_circular == 1) { + @coordsRes= `$run_st -c -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`; + } + else + { + @coordsRes= `$run_st -i $min_id -v $min_cov -V $diff_cov -l $min_len -R -u unused_contigs.out $choice.filtered.delta > $choice.tiling`; + } + } + } + my $Coordsfull= "$choice.tiling"; + return ($Coordsfull,$Path, $dir); +} + +## ############################################################################# +sub checkProg{ #checks if a given excutable is in the path... + my $prog = shift; + my $ask = 'which'; + my @check_prog = `$ask $prog`; + my $path_toPass; + my $path_dir; + my $command; + if (defined $check_prog[0] && $check_prog[0] =~ /$prog$/) + { + $path_toPass = $prog; + $command = $prog; + } + else + { + print "\nENTER the directory for your ", $prog, " executables [default ./]: "; + my $path=<STDIN>; + chomp $path; + $path_dir = $path; + if ($path_dir =~/\/$/) + { + $path_dir = $path_dir; + } + else + { + $path_dir = $path_dir."/"; + } + my @final_check = `$ask $command`; + if (exists $final_check[0] && $final_check[0] =~ /$prog$/) + { + $command = $path_dir.$prog; + $path_toPass = $command; + } + else + { + print "ERROR: Could not run ", $prog, ", please check if it is installed in your path\n or provide the directory \n"; + exit; + } + } + + return ($command, $path_toPass, $path_dir); +} + +############################## +# converts a fasta file to an ordered single line +#-------------------------------------------------- +sub Fasta2ordered +{ + if( @_ != 1 ) + { + print "Usage: Fasta2ordered <fasta_file>\n"; + exit; + } + my $fast = shift; #print $fast; exit; + my @fasta = split (/\n/, $fast); + if ($fasta[0] =~ /\>/ ) #remove chromosome name if exists in the onput sequence. + { + my $ch_name = $fasta[0]; + shift @fasta; + } + #print $fasta[0]; exit; + foreach(@fasta){chomp;} + my $num_lines = scalar(@fasta); + my $dna = ''; + for(my $i=0; $i< $num_lines; $i+=1) + { + $dna = $dna.$fasta[$i]; + } + my $ordered_dna = $dna; + return $ordered_dna; +} +############################################################## +# Hash input contigs +#------------------------------------------------ + +sub hash_contigs { + if( @_ != 1 ) + { + print "Usage: hash_contigs contigs_file"; + exit; + } + my $contigs_file = shift; + if( $contigs_file =~ /\.gz$/ ) + { + open(CONTIGS, $contigs_file, "gunzip -c $contigs_file |" ) or die "Cant open contigs file: $!\n"; + } + else + { + open( CONTIGS, $contigs_file) or die "Cant open contigs file: $!\n"; + } + + my %contigs_hash; # hash to store contig names and sequences + my $contigName; + + + while (<CONTIGS>) ##add error checking... + { + if (/^>(\S+)/) { + $contigName=$1; + } + else { + chomp; + $contigs_hash{$contigName} .= $_; + } + } + close(CONTIGS); + #tdo + ## check if qual exists + my %contigs_qual_hash; + + + if (-r "$contigs_file.qual" or -r "$contigs_file.qual.gz") { + if( -r "$contigs_file.qual.gz" ) + { + open(CONTIGS, "$contigs_file.qual.gz", "gunzip -c $contigs_file.qual.gz |" ) or die "Cant open contigs file: $!\n"; + } + else + { + open( CONTIGS, "$contigs_file.qual" ) or die "Cant open contigs file: $!\n"; + } + + + while (<CONTIGS>) { + if (/^>(\S+)/) { + $contigName=$1; + } + else { + chomp; + $contigs_qual_hash{$contigName} .= $_; + } + } + + } # end tdo # end if exist + + + return (\%contigs_hash,\%contigs_qual_hash); +} +####################### +############################## +### it gets a delta name +sub getMummerComparison{ + my $deltaName = shift; + + + ### transform the delta file to coords + my $call ="show-coords -H -T -q $deltaName > $deltaName.coords "; + !system(" $call") or die "Problems doing the show-coords comparison: $call $!\n"; + + + ### willh old results + my %h; + + ### has as index the postion with the max hits + my %h_max; + my $tmp=0; + my $tmp_index; + my $key=''; + my $is_promer =0; + if ($deltaName =~/^promer/) + { + $is_promer =1; + } + + open (F,"$deltaName.coords") or die "Problem in getComparisonFile to open file $deltaName.coords\n"; + my @File=<F>; + my @a; + @a=split(/\s+/,$File[0]); + $tmp=$a[5]; + $tmp_index=$a[0]; + if ($is_promer ==1) + { + $key=$a[12]; ## nucmer: $key = $a[8] + } + else + { + $key = $a[8]; + } + + foreach (@File) { + @a=split(/\s+/); + + if ($is_promer ==1) + { + push @{ $h{$a[12]}}, "$a[12]\t$a[11]\t$a[7]\t$a[5]\t0\t0\t$a[2]\t$a[3]\t$a[0]\t$a[1]\t1\t$a[5]\n"; + if ($key eq $a[12] and $a[5]>$tmp) + { + $tmp=$a[5]; # length + $tmp_index=$a[0]; # position reference + } + elsif ($key ne $a[12]) { + ### here possible bugg... + $h_max{$tmp_index}=$key; + $key=$a[12]; + $tmp=$a[5]; # length + $tmp_index=$a[0]; # position reference + + } + } #end if + #else i.e. if nucmer + else + { + push @{ $h{$a[8]}}, "$a[8]\t$a[6]\t$a[6]\t$a[5]\t0\t0\t$a[2]\t$a[3]\t$a[0]\t$a[1]\t1\t$a[5]\n"; + if ($key eq $a[8] and $a[5]>$tmp) + { + $tmp=$a[5]; # length + $tmp_index=$a[0]; # position reference + } + elsif ($key ne $a[8]) { + ### here possible bugg... + $h_max{$tmp_index}=$key; + $key=$a[8]; + $tmp=$a[5]; # length + $tmp_index=$a[0]; # position reference + + } + } + + } + $h_max{$tmp_index}=$key; +# print Dumper %h_max; + + return (\%h,\%h_max); +} +################################## +########################### +sub writeBinContigs2Ref{ + my $nameBin = shift; + my $name = shift; + + open (F, "$nameBin") or die "Couldn't find file $nameBin: $!\n"; + + my @ar; + + my $count=0; + + while (<F>) { + push @ar, $_; + $count++; + } + #### sa4: added error checking:- if file is empty + if (scalar(@ar) < 1) + { + print "No contigs in unusedcontigs file\n"; + $count = 0; + + } + else + { + open (F, "> $name.notMapped.contigs.tab") or die "Couldn't write file $name.tab: $!\n"; + print F doArt(\@ar); + close(F); + } + return $count; + +} + + +############################## +sub doArt{ + my ($ref) = @_; + + + ## hash of array with all positions of the contig + my %Pos; + + ## Hash with note of result line of nucmer + my %lines; + + foreach (@$ref) { + chomp; + my @ar=split(/\t/); + push @{ $Pos{$ar[12]}}, "$ar[0]..$ar[1]"; + $lines{$ar[12]} .= "FT $_\n"; + } + + my $res; + + foreach my $contig (keys %lines) { + + if (scalar(@{ $Pos{$contig} } >1)) { + my $tmp; + + foreach (@{ $Pos{$contig} }) { + $tmp.="$_,"; + } + $tmp =~ s/,$//g; # get away last comma + $res .= "FT contig join($tmp)\n"; + } + else { + $res .= "FT contig $Pos{$contig}[0]\n"; + } + $res .= "FT /systematic_id=\"$contig\"\n"; + $res .= "FT /note=\"Contig $contig couldn't map perfectly.\n"; + $res .= $lines{$contig}; + $res .= "FT \"\n"; + + } + + return $res; +} + +########################################################################## +#--------------------------- +sub makeN #creates a string of Ns +{ + my $n = shift; + my $Ns= "N"; + for (my $i =1; $i < $n; $i+=1) + { + $Ns = $Ns."N"; + } + return $Ns; +} + +########################################################################### +## reverse complement a sequence +#--------------------------------------------- +sub revComp { + my $dna = shift; + my $revcomp = reverse($dna); + + $revcomp =~ tr/ACGTacgt/TGCAtgca/; + + return $revcomp; +} + +################################################################################ +### function to visualize rep. regions in Reference genome. +sub findRepeats +{ + my $reference = shift; + my $name = shift; + my $path_prog = shift; + + # get path + my ($path_coords) = $path_prog; + + $path_coords =~ s/nucmer/show-coords/; + + my $call = "$path_prog --maxmatch -c 100 -b 40 -p $name.repeats -l 25 $reference $reference &> /dev/null "; + !system("$call") or die "Problems doing the nucmer comparison: $call $!\n"; + $call ="$path_coords -r -c -l $name.repeats.delta > $name.repeats.coords "; + !system(" $call") or die "Problems doing the show-coords comparison: $call $!\n"; + + my @Res; + open (F, "$name.repeats.coords" ) or die "Problems to open file $reference.repeats.coords: Is MUMmer installed correctly and inserted in the PATH environment variable? ($!)\n"; + $_=<F>; $_=<F>; $_=<F>; $_=<F>;$_=<F>; + while (<F>) { + my @ar = split(/\s+/); + if (!($ar[1] == $ar[4] or $ar[2] == $ar[5] or $ar[7] > 100000)) { # to exclude self alignment + + foreach ($ar[1]..$ar[2]) { + $Res[($_-1)]++; + } + } + } + + ### write the result to the plot file + + my $res; + foreach (@Res){ + if (defined($_)) { + $res .= "$_\n"; + } + else { + $res .= "0\n"; + } + } + + open (F, "> $name.Repeats.plot") or die "Couldn't open file $name.plot to write: $! \n"; + print F $res; + close(F); + + ### delete files + unlink("$name.repeats.delta"); + unlink("$name.repeats.coords"); + unlink("$name.repeats.cluster"); +} + +################################################################################ +# reverse a list of qualities +#------------------------------ +sub reverseQual{ + my $str = shift; + + $str =~ s/\s+$//g; + + my @ar = split(/\s/,$str); + my $res; + + for (my $i=(scalar(@ar)-1);$i>=0;$i--) { + $res.="$ar[$i] "; + } + return $res; +} + +############################################################################## +# ---------------- +sub getPosCoords{ + my $ref_ar = shift; + my $contig = shift; + + my $offset = shift; + my $res; + # print Dumper $$ref_ar{$contig}; + foreach (@{$$ref_ar{$contig}}) { + print "in getPos Coords: $_\n";; + + my @ar=split(/\t/); + $ar[6]+=$offset; + $ar[7]+=$offset; + $res .= join("\t",@ar);; + } + + return $res; + +} +############################################################################# +# ----------------------------- +sub getPosCoordsTurn{ + my $ref_ar = shift; + my $contig = shift; + + my $offset = shift; + + + my $res; + + # print Dumper $$ref_ar{$contig}; + foreach (@{$$ref_ar{$contig}}) { + my @ar=split(/\t/); + my $tmp_8=$ar[8]; + my $tmp_9=$ar[9]; + + $ar[8]=$ar[6]+$offset; + $ar[9]=$ar[7]+$offset; + $ar[6]=$tmp_8; + $ar[7]=$tmp_9; + + ## change query subject + + $res .= join("\t",@ar);; + } + + return $res; + +} + +############################################################################ +#------------------------ +sub printStats{ + + my ($num_fortillingcontigs,$num_notsetTilling,$num_mapped, $num_contigs, $num_inComparisoncontigs, $ref_len, $total_bases_mpd) = @_; + $num_fortillingcontigs=$num_notsetTilling+$num_mapped; + my $res; + $res.= "Short statistics of run.\n"; + $res.= "$num_contigs\t\tcontigs entered to be mapped against the reference.\n"; + $res.= sprintf("$num_inComparisoncontigs\t%0.2f \%\tmade a hit against the reference to the given parameter (-s -d etc)\n",($num_inComparisoncontigs*100/$num_contigs)); + $res.= sprintf("$num_fortillingcontigs\t%0.2f \%\twere considered for the tilling graph (-s -d etc)\n",($num_fortillingcontigs*100/$num_contigs)); + $res.= sprintf("$num_mapped\t%0.2f \%\tare mapped in the tilling graph (-s -d etc)\n",($num_mapped*100/$num_contigs)); + $res.= sprintf("\nCoverage: The reference is $ref_len long. Up to $total_bases_mpd bp (%0.2f \%) are covered by the contigs (This doesn't mean that these regions are really similar...)\n",($total_bases_mpd*100/$ref_len)); + + print $res; + +} +################################################## +#### Do Tiling +#---------------------------------------------------- +sub doTiling { + + my ($mummer_tiling, $path_toPass, $path_dir,$reference, $query_file, $choice, $prefix,$mlfas, $fasta_bin, $avoid_Ns, $ref_len, $gaps_2file, $ref_inline, $add_bin_2ps, $pick_primer, $flank, $chk_uniq, $run_blast) = @_; + + ### these are also defined in the main script.... to be changed!! + my ($num_contigs , $num_inbincontigs, $avg_cont_size,$num_overlaps , $num_gaps, + $num_mapped,$total_bases_mpd,$p_ref_covered,$num_ambigus,$num_inComparisoncontigs, + $num_fortillingcontigs,$num_notsetTilling)= (0,0,0,0,0,0,0,0,0,0,0,0); + + my ($href, $ref_contigs_qual) = hash_contigs($query_file); + my $qualAvailable=0; + my %contigs_hash = %{$href}; + my @c_keys = keys %contigs_hash; + $num_contigs = scalar(@c_keys); + my @cont_lens; + my (@ids ,$id, $id_len); + my (@Rs, @Re, @G, @Len, @cov, @pid, @orient, @Cid); + my (@Ps, @Pe); + my ($total); + my $g; #define gap size between contigs + my $tiling_gap; #gap size from tiling graph + open (TIL, $mummer_tiling) or die "Could not open $mummer_tiling: $!"; + while (<TIL>) + { + chomp; + if ($_ =~/^>/) + { + my $line = substr $_, 1; + my @splits = split /\s+/, $line; + $id = $splits[0]; + push @ids, $id; + $id_len= $splits[1]; + } + else + { + my @splits = split /\s+/, $_; + push @Rs, $splits[0]; + push @Re, $splits[1]; + push @G, $splits[2]; + push @Len, $splits[3]; + push @cov, $splits[4]; + push @pid, $splits[5]; + push @orient, $splits[6]; + push @Cid, $splits[7]; + } + + } + close (TIL); + if (scalar(@Rs) != scalar(@Re)) + { + print "ERROR: unequal array size\n"; + exit; + } + else + { + $total = scalar(@Rs); + $num_mapped = scalar(@Rs); + } + my $ref_loc = $reference; # get locations of reference and query files + my $qry_loc = $query_file; + my $dif_dir =0; #assume query and reference are in the working directory + my @splits_reference = split (/\//, $reference); + my $new_reference_file = $splits_reference[(scalar(@splits_reference)-1)]; + my @splits_query = split (/\//, $query_file); + my $new_query_file = $splits_query[(scalar(@splits_query)-1)]; + if ($prefix eq "") + { + $prefix = $new_query_file."_".$new_reference_file; + } + #------------------------------------------------------------------- + #define file handles for output files and open files to write output + #------------------------------------------------------------------- + my ($seqFH,$tabFH,$binFH,$crunchFH, $gapFH, $mlFH, $dbinFH, $avoidNFH); + open ($seqFH, '>', $prefix . '.fasta') or die "Could not open file $prefix.fasta for write: $!\n"; + open ($tabFH, '>', $prefix . '.tab') or die "Could not open file $prefix.tab for write: $!\n"; + open ($binFH, '>', $prefix . '.bin') or die "Could not open file $prefix.bin for write: $!\n"; + open ($crunchFH, '>', $prefix . '.crunch') or die "Could not open file $prefix.crunch for write: $!\n"; + open ($gapFH, '>', $prefix . '.gaps') or die "Could not open file $prefix.gaps for write: $!\n"; + if ($mlfas ==1) + { + open ($mlFH, '>', $prefix . '.contigs.fas') or die "Could not open file $prefix.contigs.fas for write: $!\n"; + } + if ($fasta_bin ==1) + { + open ($dbinFH, '>', $prefix . '.contigsInbin.fas') or die "Could not open file $prefix.contigsInbin.fas for write: $!\n"; + } + if ($avoid_Ns ==1) + { + open ($avoidNFH, '>', $prefix .'.NoNs.fasta') or die "Could not open file $prefix.NoNs.fasta for write: $!\n"; + } + #------------------------------------------------------------------------- + # Writing tiling graph and generating a pseudomolecule + # Note use use ps for pseudomolecule + #@Ps = start of ps, and @Pe = end of ps + my $ps_start =1; + $Ps[0] = 1; + $Pe[0] = $Ps[0] + $Len[0] -1; + my $tmp_qual; + my $tmp_nqual; + my $tmp_seq =""; + my $tmp_nseq =""; + print "Total contigs = $total \n"; + + #------------------------------------------------------------ + # The 'for loop' loops over each contig in the Tiling output + #Writing to file is done for each contig to speed up the process + #This part could potentially be a separate subroutine + + print $tabFH "ID ",$id, "\n"; + print $seqFH ">", "ordered_", $id, "\n"; + for (my $i=1; $i <= $total; $i+=1) + { + my $covv =sprintf("%.0f",$cov[$i -1]); #ROUNDING + my $pidd = sprintf("%.0f", $pid[$i -1]); + my ($contig_coord, $color, $contig_seq); + my $contig_qual=''; + $tiling_gap = $G[$i -1]; + if ($tiling_gap <= 5){ #insert 100Ns for overlaps and gaps of size less than or equal to 5bp + $g = 99; # default gap size to + } + else{ + $g = $tiling_gap; + } + if (defined($Len[$i])) + { + $Ps[$i] = $Pe[$i-1] +$g +1; + $Pe[$i] = $Ps[$i] + $Len[$i] -1; + $total_bases_mpd+=$Len[$i]; + } + + if ($Rs[$i -1] <0) #check if a reference starting position is less than 0 + { + $Rs[$i -1] =1; + } + + if($orient[$i-1] eq "+") + { + $contig_coord = $Ps[$i -1]."..".$Pe[$i-1]; + $color = 4; + $contig_seq = $contigs_hash{$Cid[$i-1]}; + } + else + { + $contig_coord = "complement(".$Ps[$i -1]."..".$Pe[$i-1].")"; + $color =3; + $contig_seq = revComp($contigs_hash{$Cid[$i-1]}); #REVERSE COMPLEMENT A SEQUENCE + + } + push (@cont_lens, length($contig_seq)); + + # tdo + if (defined($$ref_contigs_qual{$Cid[$i-1]})) { + ## flag to know, that the qual exists + $qualAvailable=1; + $contig_qual = $$ref_contigs_qual{$Cid[$i-1]}; + } + $tmp_qual .= $contig_qual; + $tmp_seq .= $contig_seq; + if ($avoid_Ns ==1) + { + $tmp_nseq.= $contig_seq; + #tdo + $tmp_nqual .= $contig_qual; + } + if ($mlfas ==1) + { + my $head = $Cid[$i-1]; + my $multifasta_seq = write_Fasta_headers ($contig_seq,$head); + print $mlFH $multifasta_seq; + } + if ($Re[$i -1] > $ref_len) + { + $Re[$i -1] = $ref_len -1; + } + if ($Pe[$i -1] > length($tmp_seq)) + { + $Pe[$i -1] = length($tmp_seq); + } + + #----------------------------------------- + print $crunchFH $covv, " ", $pidd, " ", $Ps[$i -1], " ", $Pe[$i -1], " ", $Cid[$i -1], " ", $Rs[$i -1], " ", $Re[$i-1], " ", "unknown NONE\n"; + + #WRITE FEATURE FILE + print $tabFH "FT contig ",$contig_coord, "\n"; + print $tabFH "FT ", "/systematic_id=\"", $Cid[$i-1],"\"","\n"; + print $tabFH "FT ", "/method=\"", "mummer\"", "\n"; + print $tabFH "FT ", "/Contig_coverage=\"",$cov[$i -1], "\"", "\n"; + print $tabFH "FT ", "/Percent_identity=\"",$pid[$i -1], "\"", "\n"; + if ($tiling_gap > 1) #WRITE GAP LOCATIONS AND SIZE TO FILE + { + my $gap_start = $Pe[$i -1] +1; + my $gap_end = ""; + if (defined $Ps[$i]) + { + $gap_end = $Ps[$i] -1; + } + my $ref_start = $Re[$i -1] +1; + my $ref_end; + if (defined $Rs[$i]) + { + $ref_end =$Rs[$i]-1; + } + else + { + $ref_end = "END"; + } + print $gapFH "Gap\t",$tiling_gap, "\t", $gap_start, "\t", $gap_end, "\t", $ref_start, "\t", $ref_end,"\n"; + if ($gaps_2file ne "" && $ref_start < $ref_len) + { + my $ref_gapsFH; + open ($ref_gapsFH, '>', $gaps_2file.'.Gaps_onRef') or die "Could not open file $gaps_2file.Gaps_onRef for write: $!\n"; + my $gapOnref = substr ($ref_inline, $ref_start, $g); + print $ref_gapsFH ">",$g,"_",$ref_start, "\n"; + my $file_toPrint = write_Fasta ($gapOnref); + print $ref_gapsFH $file_toPrint; + } + } + else + { + $color = 5; + print $tabFH "FT ", "/Overlapping=\"", "YES\"", "\n"; + } + my $ns = makeN($g); + $tmp_seq = $tmp_seq.$ns; + #tdo + for (1..length($ns)) { + $tmp_qual .= "0 "; + } + print $tabFH "FT ", "/colour=\"",$color, "\"", "\n"; + } + #------------------------------------------------------------------ + #tdo + my @Quality_Array; + if ($qualAvailable) { + @Quality_Array = split(/\s/,$tmp_qual); + my $res; + foreach (@Quality_Array) { + $res .= "$_\n"; + } + ## get name + my @splits_query = split (/\//, $query_file); + $new_query_file = $splits_query[(scalar(@splits_query)-1)]; + open (F,"> $new_query_file.qual.plot") or die "problems\n"; + print F $res; + close(F); + } + ##WRITE PSEUDOMOLECULE WITHOUT 'N's + #-------------------------------------- + if ($avoid_Ns ==1) + { + print $avoidNFH ">", "ordered_", $id, "without 'N's","\n"; + my $toWrite = write_Fasta ($tmp_nseq); + print $avoidNFH $toWrite; + } + #################################### + #WRITE CONTIGS WITH NO HIT TO FILE # + ################################# + my %Cids; + + foreach(@Cid) + { + chomp; + $Cids{$_} = 1; + } + my @contigs_2bin = (); + my %h_contigs_2bin; + + foreach (@c_keys) + { + push(@contigs_2bin, $_) unless exists $Cids{$_}; + + } + foreach(@contigs_2bin) + { + $h_contigs_2bin{$_}=1; + + print $binFH "$_ \n"; + + } + $num_inbincontigs= scalar(@contigs_2bin); + + ######## + # WRITE PSEUDOMOLECULE TO FILE + #---------------------------------- + my $new_seq = $tmp_seq; + my $prev_len = length($tmp_seq); + my $total_len = $prev_len; + foreach (@contigs_2bin) + { + chomp; + #my $binseq = $contigs_hash{$contigs_2bin[$i]}; + my $l = length ($contigs_hash{$_}); + $total_len +=$l; + } + if ($add_bin_2ps ==1) #appending unmapped contigs to pseudomolecule + { + + for (my $i =0; $i < scalar(@contigs_2bin); $i+=1) + { + my $binseq = $contigs_hash{$contigs_2bin[$i]}; + $new_seq .=$contigs_hash{$contigs_2bin[$i]}; + my $len_current_contig = length($binseq); + my $start = $prev_len +1; + my $end = $start + $len_current_contig -1; + my $col = 7; + if ($start > $total_len) + { + $start = $total_len; + } + if ($end >$total_len) + { + $end = $total_len; + } + my $co_cord = $start."..".$end; + my $note = "NO_HIT"; + print $tabFH "FT contig ",$co_cord, "\n"; + print $tabFH "FT ", "/systematic_id=\"", $contigs_2bin[$i],"\"","\n"; + print $tabFH "FT ", "/method=\"", "mummer\"", "\n"; + print $tabFH "FT ", "/colour=\"",$col, "\"", "\n"; + print $tabFH "FT ", "/", $note, "\n"; + $prev_len= $end; + } + } + my $to_write = write_Fasta ($new_seq); + print $seqFH $to_write; + ######## + #WRITE CONTIGS IN BIN TO FILE # + #------------------------------------------------------ + if ($fasta_bin ==1) + { + foreach(@contigs_2bin) + { + print $dbinFH ">", $_, "\n"; + my $to_write = write_Fasta($contigs_hash{$_}); + print $dbinFH $to_write; + } + } + #unlink ("$choice.delta"); + #unlink ("$choice.filtered.delta"); + #unlink ("$choice.cluster"); + #unlink ("$choice.tiling"); + #PRINT FINAL MESSAGE + print " FINISHED CONTIG ORDERING\n"; + print "\nTo view your results in ACT\n\t\t Sequence file 1: $new_reference_file\n\t\t Comparison file 1: $prefix.crunch\n\t\t Sequence file 2: $prefix.fasta\n + \t\tACT feature file is: $prefix.tab\n + \t\tContigs bin file is: $prefix.bin\n + \t\tGaps in pseudomolecule are in: $prefix.gaps\n\n"; + + #Run tblastx.... + if ($run_blast ==1) + { + print "Running tblastx on contigs in bin...\nThis may take several minutes ...\n"; + my $formatdb = 'formatdb -p F -i' ; +# my @formating = ` + !system("$formatdb $new_reference_file") or die "ERROR: Could not find 'formatdb' for blast\n"; + my $blast_opt = 'blastall -m 9 -p tblastx -d '; + my $contigs_inBin = $prefix.'.contigsInbin.fas'; +# my @bigger_b = ` + !system("$blast_opt $new_reference_file -i $contigs_inBin -o blast.out") or die "ERROR: Could not find 'blastall' , please install blast in your working path (other dir==0)\n$blast_opt $new_reference_file -i $contigs_inBin -o blast.out\n \n"; + } + + + if ($pick_primer == 1) + { + print " DESIGNING PRIMERS FOR GAP CLOSURE...\n"; + my $qq = "$prefix.fasta"; + pickPrimers($qq, $reference, $flank, $path_toPass, $chk_uniq,$qualAvailable,@Quality_Array); + } +} + + +#------------------------------ +sub write_Fasta { + my $sequence = shift; + my $fasta_seq =""; + my $length = length($sequence); + if ($length <= 60) + { + $fasta_seq = $sequence ; + } + elsif ($length> 60 ) + { + for (my $i =0; $i < $length; $i+=60) + { + my $tmp_s = substr $sequence, $i, 60; + $fasta_seq .= $tmp_s."\n"; + } + } + + return $fasta_seq; +} + +sub write_Fasta_headers { + my $sequence = shift; + my $header = shift; + my $fasta_seq =">$header\n"; + my $length = length($sequence); + if ($length <= 60) + { + $fasta_seq.= $sequence ; + } + elsif ($length> 60 ) + { + for (my $i =0; $i < $length; $i+=60) + { + my $tmp_s = substr $sequence, $i, 60; + $fasta_seq .= $tmp_s."\n"; + } + } + + return $fasta_seq; +} + +#---------------------------------------------- END OF CONTIG ORDERING SUBROUTINES---------------------------------------------------- + +#----------------------------------------------- PRIMER DESIGN --------------------------------------------------- +sub pickPrimers +{ + #$ps = pseudo molecule,$rf = reference, $flan = flanking region size + my ($ps,$rf, $flan, $passed_path, $chk_uniq,$qualAvailable, @Quality_Array); + if (@_==4){ + ($rf,$ps, $flan, $chk_uniq) = @_; + print "Primers without ordering..\n"; + print $rf; + $passed_path = "nucmer"; + $qualAvailable =0; + @Quality_Array = []; + } + else #(@_ == 7) + { + ($ps,$rf, $flan, $passed_path, $chk_uniq, $qualAvailable, @Quality_Array) = @_; + } + + my $dna=''; + my @gappedSeq; + my $records=''; + my @sequence; + my $input=''; + #tdo + my @gappedQual; + #my $quality=''; + my $path_toPass = $passed_path; + my @fasta; + my $query=''; + my @exc_regions; + my $ch_name; + #my $flank = $flan; + open (FH, $rf) or die "Could not open reference file\n"; + open (FH2, $ps) or die "Could not open query/pseudomolecule file\n"; + my $ref; #print ".... ", $rf; exit; + my @r = <FH>; + my @qry = <FH2>; + my $dn = join ("", @qry); + $ref = join ("", @r); + $dna = Fasta2ordered ($dn); + #check if primer3 is installed + my $pr3 = "primer3_core"; + my ($pr3_path, $pr3_Path, $pr3_path_dir) = checkProg ($pr3); + #my @check_prog = `which primer3_core`; + open (PRI, '>primer3.summary.out') or die "Could not open file for write\n"; + + # + + #print $ref; exit; + #PARSING FOR PRIMER3 INPUT + my ($opt,$min,$max,$optTemp,$minTemp,$maxTemp,$flank,$lowRange,$maxRange,$gcMin,$gcOpt,$gcMax,$gclamp,$exclude,$quality) = getPoptions($qualAvailable); + my ($gap_position,@positions, %seq_hash); + + my $exc1 = $flank -$exclude; #start of left exclude + print "Please wait... extracting target regions ...\n"; + #regular expression extracts dna sequence before and after gaps in sequence (defined by N) + while($dna=~ /([atgc]{$flank,$flank}N+[atgc]{$flank,$flank})/gi) + { + $records= $1; + push (@gappedSeq, $records); + $gap_position = index($dna, $records); + push @positions, $gap_position; + $seq_hash{$gap_position}=$records; + #dna + if ($qualAvailable) { + my $res; + for (my $nn=($gap_position-1); $nn <= ($gap_position-1+length($records)-1); $nn++) { + $res.="$Quality_Array[$nn] "; + } + push @gappedQual, $res; + } + } + #loop prints out primer targets into a file format accepted by primer3 + my $count=1; + my $identify=''; + my $seq_num = scalar @gappedSeq; + my $name= " "; + + my ($totalp, @left_name, @right_name, @left_seq, @right_seq); + + my ($leftP_names, $rightP_names, $leftP_seqs, $rightP_seqs, $leftP_start, $leftP_lens, $rightP_ends, $rightP_lens,$left_Exclude,$right_Exclude, $primers_toExc, $prod_size)= + ("","","","","","","","","","", "", ""); + + print $seq_num, " gaps found in target sequence\n"; + print "Please wait...\nLooking for primers...\n"; + print "Running Primer3 and checking uniquness of primers...\nCounting left and right primer hits from a nucmer mapping (-c 15 -l 15)\n"; + + for (my $i=0; $i<$seq_num; $i+=1) + { + $identify = $count++; + if (defined $ch_name) + { + $name = $ch_name; + } + my $len = length($gappedSeq[$i]); + my $exc2 = $len - $flank; + open(FILE, '>data') or die "Could not open file\n"; + #tdo + my $qual=''; + if ($qualAvailable) { + $qual="PRIMER_SEQUENCE_QUALITY=$gappedQual[$i]\nPRIMER_MIN_QUALITY=$quality\n"; + } + +#WARNING: indenting the following lines may cause problems in primer3 +print FILE "PRIMER_SEQUENCE_ID=Starting_Pos $positions[$i] +SEQUENCE=$gappedSeq[$i] +PRIMER_OPT_SIZE=$opt +PRIMER_MIN_SIZE=$min +PRIMER_MAX_SIZE=$max +PRIMER_OPT_TM=$optTemp +PRIMER_MIN_TM=$minTemp +PRIMER_MAX_TM=$maxTemp +PRIMER_NUM_NS_ACCEPTED=1 +PRIMER_PRODUCT_SIZE_RANGE=$lowRange-$maxRange +PRIMER_MIN_GC=$gcMin +PRIMER_GC_CLAMP =$gclamp +PRIMER_OPT_GC_PERCENT=$gcOpt +PRIMER_MAX_GC=$gcMax +PRIMER_INTERNAL_OLIGO_EXCLUDED_REGION=$exc1,$exclude $exc2,$exclude +".$qual."Number To Return=1 +=\n"; +close FILE; + + #runs primer3 from commandline + + ################# NOTE: PRIMER3 SHOULD BE IN YOUR WORKING PATH ######### + + my @Pr3_output = `$pr3_path -format_output <data`; + #print $positions[$i], "\t", $i, " ", $path_toPass, " ", $rf, $exc1, " ",$exc2, "\n"; + my $fil = join (":%:", @Pr3_output); + my ($uniq_primer, $string,$left_nm,$right_nm,$left_sq, $right_sq,$left_strt,$left_ln, $right_End,$right_ln,$primers_toExclude, $product_size) + = check_Primers ($fil, $positions[$i], $i,$path_toPass, $rf, $exc1, $exc2); + + print PRI $string; + if ($uniq_primer ==1) + { + $leftP_names.=$left_nm."\n"; + $rightP_names.=$right_nm."\n"; + $leftP_seqs.=$left_sq."\n"; + $rightP_seqs.=$right_sq."\n"; + $leftP_start.=$left_strt."\n"; + $leftP_lens.=$left_ln."\n"; + $rightP_ends.=$right_End."\n"; + $rightP_lens.=$right_ln."\n"; + $left_Exclude.=$exc1."\n"; + $right_Exclude.=$exc2."\n"; + $prod_size.=$product_size."\n"; + } + if ($primers_toExclude ne "") + { + $primers_toExc.= $primers_toExclude; #."\n"; + } + + } + write_Primers ($leftP_names, $rightP_names, $leftP_seqs, $rightP_seqs, $leftP_start, $leftP_lens, $rightP_ends, $rightP_lens,$primers_toExc,$left_Exclude,$right_Exclude, $prod_size); + #write_Primers (@left_name, @right_name, @left_seq, @right_seq,@left_start, @left_len, @right_end, @right_len, @left_exclude, @right_exclude, $primers_toExclude); + +} + +#checks the uniqueness of primers +#input an array with promer3 output for each gap +sub check_Primers +{ + + my ($fil, $position, $index,$path_toPass, $rf, $exc1, $exc2) = @_; + my @Pr3_output = split /:%:/, $fil; + my ($left_name, $right_name, $left_seq, $right_seq, $left_start,$left_len,$right_end,$right_len,$left_exclude,$right_exclude) = ("", "", "", "", "", "", "", "", "", ""); + my $primers_toExclude =""; + my $product_size =""; + my $string =""; + my $uniq_primer = 0; + $string.="=========================================================================================\n"; + $string.="Primer set for region starting at ".$position."\n"; + + if (defined $Pr3_output[5] && defined $Pr3_output[6]) + { + if ($Pr3_output[5]=~ /LEFT PRIMER/) + { + # print $Pr3_output[5]; + #check uniquness of primer against the genome + my @splits_1 = split (/\s+/, $Pr3_output[5]); + my $left_primer = $splits_1[8]; + my $left_st = $splits_1[2]; + my $left_length = $splits_1[3]; + + my @splits_2 = split (/\s+/, $Pr3_output[6]); + my $right_primer = $splits_2[8]; + my $right_st = $splits_2[2]; + my $right_length = $splits_2[3]; + + open (QRY_1, '>./left_query'); # open a file for left primers + print QRY_1 ">left_01\n"; # + print QRY_1 $left_primer,"\n"; + open (QRY_2, '>./right_query'); + print QRY_2 ">right_01\n"; + print QRY_2 $right_primer,"\n"; + + my ($left_count, $right_count); + #if ($chk_uniq eq "nucmer") + #{ + my $options = "-c 15 --coords -l 15 "; + my $rq = "right_query"; + my $lq = "left_query"; + my (@right_ps, @left_ps); + # print $path_toPass, "\t", $options, "\n"; + + + my @Rrun = `$path_toPass $options -p R $rf $rq &> /dev/null`; + print "."; + my $f1 = "R.coords"; + open (RP, $f1) or die "Could not open file $f1 while checking uniqueness of right primer\n"; + while (<RP>) + { + chomp; + if ($_ =~ /right_01$/) + { + push @right_ps, $_; + } + } + close (RP); + my @Lrun = `$path_toPass $options -p L $rf $lq &> /dev/null`; + print "."; + my $f2 = "L.coords"; + open (LQ, $f2) or die "Could not open file $f2\n"; + while (<LQ>) + { + chomp; + if ($_ =~ /left_01$/) + { + push @left_ps, $_; + } + } + close (LQ); + $right_count = scalar (@right_ps); + $left_count = scalar(@left_ps); + #check if a primer is not in the excluded region:: + my $primer_NearEnd =0; + if ($left_st > $exc1 || $right_st < $exc2) + { + $primer_NearEnd = 1; + } + + if ($left_count < 2 && $right_count<2 && $primer_NearEnd ==0) + { + $string.=$left_count."\t".$Pr3_output[5]."\n"; + $string.=$right_count."\t".$Pr3_output[6]."\n"; + $string.="***************************** PRIMER3 OUTPUT **************************\n"; + foreach (@Pr3_output) {$string.=$_;} + + my @prod_size_split = split /\s+/, $Pr3_output[10]; + + $product_size = substr($prod_size_split[2], 0, -1); + $left_name = $position; + $right_name = $position; + my $lp_uc = uc ($left_primer); + my $rp_uc = uc($right_primer); + #print $left_count, "..", $right_count, "\t"; + $left_seq = $lp_uc; + $right_seq= $rp_uc; + + $left_start= $left_st; + $left_len = $left_length; + + $right_end = $right_st; + $right_len = $right_length; + + $left_exclude = $exc1; + $right_exclude =$exc2; + $uniq_primer =1; + } + else + { + if ($primer_NearEnd ==1) + { + $string.="One of the oligos is near the end of a contig\n"; + } + else + { + $string.="Primer set not unique\n"; + } + $primers_toExclude.=">L.".$position."\n".$left_primer."\n"; + $primers_toExclude.=">R.".$position."\n".$right_primer."\n"; + } + + } + else + { + $string.="No Primers found\n"; + } + } + + return ($uniq_primer, $string,$left_name,$right_name,$left_seq, $right_seq,$left_start,$left_len, $right_end,$right_len,$primers_toExclude, $product_size); + + +} + +###------------------------------------ +# Writes primers and their regions to file +sub write_Primers { + my ($leftP_names, $rightP_names, $leftP_seqs, $rightP_seqs, $leftP_start, $leftP_lens, $rightP_ends, $rightP_lens,$primers_toExclude,$left_Exclude,$right_Exclude, $product_sizes) = @_; + my (@left_name, @right_name, @left_seq, @right_seq, @left_start, @left_len, @right_end, @right_len, @left_exclude, @right_exclude, @product_size); + + #open files to read + @left_name = split /\n/, $leftP_names; + @right_name= split /\n/, $rightP_names; + @left_seq = split /\n/, $leftP_seqs; + @right_seq = split /\n/, $rightP_seqs; + + @left_start = split /\n/, $leftP_start; + @left_len = split /\n/, $leftP_lens; + @right_end = split/\n/, $rightP_ends; + @right_len = split /\n/, $rightP_lens; + @left_exclude = split /\n/, $left_Exclude; + @right_exclude = split /\n/,$right_Exclude; + @product_size = split /\n/, $product_sizes; + + my $primers_withSize =""; + open (SEN, '>sense_primers.out') or die "Could not open file for write\n"; + open (ASEN, '>antiSense_primers.out') or die "Could not open file for write\n"; + open (REG_1, '>sense_regions.out') or die "Could not open file for write\n"; + open (REG_2, '>antiSense_regions.out') or die "Could not open file for write\n"; + + if ($primers_toExclude ne "") + { + open (PEX, '>primers_toExclude.out') or die "Could not open file for write\n"; + print PEX $primers_toExclude; + } + + + my $totalp = scalar (@left_name); + + #print $totalp, "\n"; exit; + + my $well_pos; + my $max_plates = ceil($totalp/96); + #print "MAX Ps ", $max_plates, "\n"; + my $plate=1; + my $sen =""; + my $asen =""; + my $plate_counter =0; + my $wells = 96; + for (my $index =0; $index < $totalp; $index += $wells) + { + my $do = $index; + my $upper_bound= $index + $wells; + if ($upper_bound > $totalp) + { + $upper_bound = $totalp; + } + + for (my $j=$index; $j <= ($upper_bound-1); $j+=1) + { + my $i = $j; + if ($j < 96) + { + $well_pos = get_WellPosition ($j); + } + else + { + $well_pos = get_WellPosition ($j - $wells) + } + + #$primers_withSize.=$product_size[$i]."\t"."Plate_".$plate. "\t\tS.".$i."\tS.".$left_name[$i]."\t". + print SEN "Plate_".$plate, "\t\t","S.", $i, "\tS.", $left_name[$i], "\t", $left_seq[$i], "\t\t+", "\t", $well_pos, "\n"; + print ASEN "Plate_".$plate, "\t\t","AS.", $i, "\tAS.", $right_name[$i], "\t", $right_seq[$i], "\t\t-","\t", $well_pos,"\n"; + print REG_1 "Plate_".$plate, "\t\t","S.", $i, "\t", $left_start[$i], "\t", $left_len[$i], "\n"; + print REG_2 "Plate_".$plate, "\t\t","AS.", $i, "\t", $right_end[$i], "\t",$right_len[$i], "\n"; + + } + $plate +=1; + } + + #delete tmp. files + #my $rm = "rm -f"; + system ("rm -f data left_query right_query R.delta R.cluster R.coords L.delta L.cluster L.coords"); + print "\nPRIMER DESIGN DONE\n\n"; + # end of primer design program +}#// +##### +# returns a well position for oligos +sub get_WellPosition{ + + my $j = shift; + my $well_pos; + if ($j < 12) + { + $well_pos = "a".($j+1); + } + elsif ($j>11 && $j<24) { + $well_pos = "b". (($j+1) -12); + } + elsif ($j>23 && $j<36) { + $well_pos = "c". (($j+1) -24); + } + elsif ($j>35 && $j<48) { + $well_pos = "d". (($j+1) - 36); + } + elsif($j>47 && $j<60) { + $well_pos = "e". (($j+1) -48); + } + elsif ($j>59 && $j<72) + { + $well_pos = "f". (($j+1) - 60); + } + elsif ($j>71 && $j< 84) + { + $well_pos = "g". (($j+1) - 72); + } + elsif ($j>83 && $j<96) + { + $well_pos = "h". (($j+1) - 84); + } + return $well_pos; +} + + +#################################################################### +#get options for primer design +#---------------------- +sub getPoptions{ + + my $qualAvailable = shift; + #### USER INPUTS ########## + #ask for optimum primer size + print "\nEnter Optimum Primer size (default 20 bases):"; + my $opt=<STDIN>; + chomp $opt; + if($opt eq '') + { + $opt = 20; + } + #ask for minimum primer size + print "\nEnter Minimum Primer size (default 18 bases):"; + my $min=<STDIN>; + chomp $min; + if($min eq '') + { + $min = 18; + } + #ask for maximum primer size + print "\nEnter Maximum Primer size (default 27 bases):"; + my $max= <STDIN>; + chomp $max; + if($max eq '') + { + $max= 27; + } + #ask for optimum primer temperature + print "\nEnter Optimum melting temperature (Celcius) for a primer oligo (default 60.0C):"; + my $optTemp=<STDIN>; + chomp $optTemp; + if($optTemp eq '') + { + $optTemp = 60.0; + } + #ask for minimum primer temperature + print "\nEnter Minimum melting temperature (Celcius) for a primer oligo (default 57.0C):"; + my $minTemp=<STDIN>; + chomp $minTemp; + if($minTemp eq '') + { + $minTemp = 57.0; + } + #ask for maximum primer temperature + print "\nEnter Maximum melting temperature (Celcius) for a primer oligo (default 63.0C):"; + my $maxTemp=<STDIN>; + chomp $maxTemp; + if($maxTemp eq '') + { + $maxTemp = 63.0; + } + print "\nEnter flanking region size (default 1000 bases): "; + my $flank=<STDIN>; + chomp $flank; + if ($flank eq '') + { + $flank = 1000; + } + #ask for primer product range + print "\nEnter minimum product size produced by primers (default =flanking size):"; + my $lowRange=<STDIN>; + chomp $lowRange; + if($lowRange eq '') + { + $lowRange = $flank; + } + print "\nEnter maxmimum product size produced by primers (default 7000):"; + my $maxRange=<STDIN>; + chomp $maxRange; + if($maxRange eq '') + { + $maxRange = 7000; + } + #ask for minimum GC content in primers + print "\nEnter minimum GC content in primers (default 20%):"; + my $gcMin=<STDIN>; + chomp $gcMin; + if($gcMin eq '') + { + $gcMin = 20.0; + } + #ask for optimum GC content in primers + print "\nEnter optimum GC content in primers (default 50%):"; + my $gcOpt=<STDIN>; + chomp $gcOpt; + if($gcOpt eq '') + { + $gcOpt = 50.0; + } + #ask for maximum GC content in primers + print "\nEnter maximum GC content in primers (default 80%):"; + my $gcMax=<STDIN>; + chomp $gcMax; + if($gcMax eq '') + { + $gcMax = 80.0; + } + print "\nEnter GC clamp (default 1):"; + my $gclamp=<STDIN>; + chomp $gclamp; + if($gclamp eq '') + { + $gclamp = 1; + } + print "\nEnter size of region to exclude at the end of contigs (default 100 bases):"; + my $exclude=<STDIN>; + chomp $exclude; + if ($exclude eq '') + { + $exclude = 100; + } + + + #tdo + my $quality=''; + if ($qualAvailable) + { + + print "\nEnter minimum quality for primer pick (default 40):"; + $quality=<STDIN>; + chomp $quality; + if($quality eq '') + { + $quality = 40; + } + } + + +return ($opt,$min,$max,$optTemp,$minTemp,$maxTemp,$flank,$lowRange,$maxRange,$gcMin,$gcOpt,$gcMax,$gclamp,$exclude, $quality); + +} +############### +#-----------------------------------------------------END of PRIMER DESIGN ---------------------------------------------------------------- +#-----------------------------------------------------END OF ABACAS ----------------------------------------------------------------------- + +