diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample
index 085e75a158b..d30119c9c16 100644
--- a/tool_conf.xml.sample
+++ b/tool_conf.xml.sample
@@ -176,6 +176,7 @@
s
+
diff --git a/tools/regVariation/microsatellite_birthdeath.pl b/tools/regVariation/microsatellite_birthdeath.pl
new file mode 100644
index 00000000000..3bf909a4428
--- /dev/null
+++ b/tools/regVariation/microsatellite_birthdeath.pl
@@ -0,0 +1,3979 @@
+#!/usr/bin/perl -w
+use strict;
+use warnings;
+use Term::ANSIColor;
+use Pod::Checker;
+use File::Basename;
+use IO::Handle;
+use Cwd;
+use File::Path qw(make_path remove_tree);
+use File::Temp qw/ tempfile tempdir /;
+my $tdir = tempdir( CLEANUP => 0 );
+chdir $tdir;
+my $dir = getcwd;
+#print "current dit=$dir\n";
+
+use vars qw (%treesToReject %template $printer $interr_poscord $interrcord $no_of_interruptionscord $stringfile @tags
+$infocord $typecord $startcord $strandcord $endcord $microsatcord $motifcord $sequencepos $no_of_species
+$gapcord %thresholdhash $tree_decipherer @sp_ident %revHash %sameHash %treesToIgnore %alternate @exactspecies @exacttags);
+use FileHandle;
+use IO::Handle; # 5.004 or higher
+
+#my @ar = ("/Users/ydk/work/rhesus_microsat/results/galay/chr22_5sp.maf.txt", "/Users/ydk/work/rhesus_microsat/results/galay/dataset_11.dat",
+#"/Users/ydk/work/rhesus_microsat/results/galay/chr22_5spec.maf.summ","hg18,panTro2,ponAbe2,rheMac2,calJac1","((((hg18, panTro2), ponAbe2), rheMac2), calJac1)","9,10,12,12",
+#"10","0.8");
+my @ar = @ARGV;
+my ($maf, $orth, $summout, $species_set, $tree_definition, $thresholds, $FLANK_SUPPORT, $SIMILARITY_THRESH) = @ar;
+$SIMILARITY_THRESH=$SIMILARITY_THRESH/100;
+#########################
+$SIMILARITY_THRESH = $SIMILARITY_THRESH/100;
+my $EDGE_DISTANCE = 10;
+my $COMPLEXITY_SUPPORT = 20;
+load_thresholds("9_10_12_12");
+#########################
+
+my $complexity=int($COMPLEXITY_SUPPORT * (1/40));
+
+#print "complexity=$complexity\n";
+#;
+
+#$printer = 1;
+
+my $rando = int(rand(1000));
+my $localdate = `date`;
+$localdate =~ /([0-9]+):([0-9]+):([0-9]+)/;
+my $info = $rando.$1.$2.$3;
+
+#---------------------------------------------------------------------------
+# GETTING INPUT INFORMATION AND OPENING INPUT AND OUTPUT FILES
+
+
+my @thresharr = (0, split(/,/,$thresholds));
+
+my $megamatch = $orth.".megamatch.net.axt"; #"/gpfs/home/ydk104/work/rhesus_microsat/axtNet/hg18.panTro2.ponAbe2.rheMac2.calJac1/chr1.hg18.panTro2.ponAbe2.rheMac2.calJac1.net.axt";
+my $megamatchlck = $megamatch.".lck";
+unlink $megamatchlck;
+
+#my $selected= $orth;
+#my $eventfile = $orth;
+#$selected = $selected."_SELECTED";
+#$selected = $selected."_".$SIMILARITY_THRESH;
+#my $runtime = $selected.".runtime";
+
+my $inputtags = "H:C:O:R:M";
+$inputtags = $ARGV[3] if exists $ARGV[3] && $ARGV[3] =~ /[A-Z]:[A-Z]/;
+
+my @all_tags = split(/:/, $inputtags);
+my $inputsp = "hg18:panTro2:ponAbe2:rheMac2:calJac1";
+$inputsp = $ARGV[4] if exists $ARGV[4] && $ARGV[3] =~ /[0-9]+:/;
+@sp_ident = split(/:/,$inputsp);
+my $junkfile = $orth."_junk";
+
+my $sh = load_sameHash(1);
+my $rh = load_revHash(1);
+
+#print "inputs are : \n"; foreach(@ARGV){print $_,"\n";}
+#open (SELECT, ">$selected") or die "Cannot open selected file: $selected: $!";
+open (SUMMARY, ">$summout") or die "Cannot open summout file: $summout: $!";
+#open (RUN, ">$runtime") or die "Cannot open orth file: $runtime: $!";
+#my $ctlfile = "baseml\.ctl"; #$ARGV[4];
+#my $treefile = "/gpfs/home/ydk104/work/rhesus_microsat/codes/lib/"; #1 THIS IS THE THE TREE UNDER CONSIDERATION, IN NEWICK
+my %registeredTrees = ();
+my @removalReasons =
+("microsatellite is compound",
+"complex structure",
+"if no. if micros is more than no. of species",
+"if more than one micro per species ",
+"if microsat contains N",
+"different motif than required ",
+"more than zero interruptions",
+"microsat could not form key ",
+"orthologous microsats of different motif size ",
+"orthologous microsats of different motifs ",
+"microsats belong to different alignment blocks altogether",
+"microsat near edge",
+"microsat in low complexity region",
+"microsat flanks dont align well",
+"phylogeny not informative");
+my %allowedhash=();
+#---------------------------------------------------------------------------
+# WORKING ON MAKING THE MEGAMATCH FILE
+my $chromt=int(rand(10000));
+my $p_chr=$chromt;
+
+$tree_definition=~s/,/, /g;
+$tree_definition =~ s/, +/, /g;
+my @exactspeciesset_unarranged = split(/,/,$species_set);
+my $largesttree = "$tree_definition;";
+$tree_definition=~s/[\)\(, ]/\t/g;
+
+my @treespecies=split(/\t+/,$tree_definition);
+
+foreach my $spec (@treespecies){
+ foreach my $espec (@exactspeciesset_unarranged){
+ push @exactspecies, $spec if $spec eq $espec;
+ }
+}
+#print "exactspecies=@exactspecies\n";
+my $focalspec = $exactspecies[0];
+my $arranged_species_set=join(".",@exactspecies);
+@exacttags=@exactspecies;
+foreach my $extag (@exacttags){
+ $extag =~ s/hg18/H/g;
+ $extag =~ s/panTro2/C/g;
+ $extag =~ s/ponAbe2/O/g;
+ $extag =~ s/rheMac2/R/g;
+ $extag =~ s/calJac1/M/g;
+}
+my $chr_name = join(".",("chr".$p_chr),$arranged_species_set, "net", "axt");
+#print "sending to maftoAxt_multispecies: $maf, $tree_definition, $chr_name, $species_set .. focalspec=$focalspec \n";
+maftoAxt_multispecies($maf, $tree_definition, $chr_name, $species_set);
+my @filterseqfiles= ($chr_name);
+ $largesttree =~ s/hg18/H/g;
+ $largesttree =~ s/panTro2/C/g;
+ $largesttree =~ s/ponAbe2/O/g;
+ $largesttree =~ s/rheMac2/R/g;
+ $largesttree =~ s/calJac1/M/g;
+#---------------------------------------------------------------------------
+
+my ($lagestnodes, $largestbranches) = get_nodes($largesttree);
+shift (@$lagestnodes);
+my @extendedtitle=();
+
+my $title = ();
+my $parttitle = ();
+my @titlearr = ();
+my @firsttitle=($focalspec."chrom", $focalspec."start", $focalspec."end", $focalspec."motif", $focalspec."motifsize", $focalspec."threshold");
+
+my @finames= qw(chr start end motif motifsize microsat mutation mutation.position mutation.from mutation.to insertion.details deletion.details);
+
+my @fititle=();
+
+foreach my $spec (split(",",$species_set)){
+ push @fititle, $spec;
+ foreach my $name (@finames){
+ push @fititle, $spec.".".$name;
+ }
+}
+
+
+my @othertitle=qw(somechr somestart somened event source);
+
+my @fnames = ();
+push @fnames, qw(insertions_num deletions_num motinsertions_num motinsertionsf_num motdeletions_num motdeletionsf_num noninsertions_num nondeletions_num) ;
+push @fnames, qw(binsertions_num bdeletions_num bmotinsertions_num bmotinsertionsf_num bmotdeletions_num bmotdeletionsf_num bnoninsertions_num bnondeletions_num) ;
+push @fnames, qw(dinsertions_num ddeletions_num dmotinsertions_num dmotinsertionsf_num dmotdeletions_num dmotdeletionsf_num dnoninsertions_num dnondeletions_num) ;
+push @fnames, qw(ninsertions_num ndeletions_num nmotinsertions_num nmotinsertionsf_num nmotdeletions_num nmotdeletionsf_num nnoninsertions_num nnondeletions_num) ;
+push @fnames, qw(substitutions_num bsubstitutions_num dsubstitutions_num nsubstitutions_num indels_num subs_num);
+
+my @fullnames = ();
+
+foreach my $lnode (@$lagestnodes){
+ my @pair = @$lnode;
+ my @nodemutarr = ();
+ for my $p (@pair){
+ # print "p = $p\n";
+ $p =~ s/[\(\), ]+//g;
+ $p =~ s/H/hg18/g;
+ $p =~ s/C/panTro2/g;
+ $p =~ s/O/ponAbe2/g;
+ $p =~ s/R/rheMac2/g;
+ $p =~ s/M/calJac1/g;
+ foreach my $n (@fnames) { push @fullnames, $p.".".$n;}
+ }
+}
+print SUMMARY "#",join("\t", @firsttitle, @fititle, @othertitle);
+
+print SUMMARY "\t",join("\t", @fullnames);
+#$title = $title."\t".join("\t", @fullnames);
+
+print SUMMARY "\t",join("\t", @fnames);
+#$title= $title."\t".join("\t", @fnames);
+
+print SUMMARY "\t","tree","\t", "cleancase", "\n";
+#$title= $title."\t"."tree"."\t"."cleancase". "\n";
+
+#print $title; #;
+
+#print "all_tags = @all_tags\n";
+
+for my $no (3 ... $#all_tags+1){
+# print "no=$no\n"; #;
+ @tags = @all_tags[0 ... $no-1];
+ #print "tags = = @tags\n" if $printer == 1;
+ %template=();
+ my @nextcounter = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ #next if scalar(@tags) < 4;
+
+ #print "now doing tags = @tags, no = $no\n";
+ open (ORTH, "<$orth") or die "Cannot open orth file: $orth: $!";
+
+# print SUMMARY join "\t", qw (species chr start end branch motif microsat mutation position from to insertion deletion);
+
+
+ ##################### T E M P O R A R Y #####################
+ my @finaltitle=();
+ my @singletitle = qw (species chr start end motif motifsize microsat strand microsatsize col10 col11 col12 col13);
+ my $endtitle = ();
+ foreach my $tag (@tags){
+ my @tempsingle = ();
+
+ foreach my $single (@singletitle){
+ push @tempsingle, $tag.$single;
+ }
+ @finaltitle = (@finaltitle, @tempsingle);
+ }
+
+# print SUMMARY join("\t",@finaltitle),"\n";
+
+ #############################################################
+
+ #---------------------------------------------------------------------------
+ # GET THE TREE FROM TREE FILE
+ my $tree = ();
+ $tree = "((H, C), O)" if $no == 3;
+ $tree = "(((H, C), O), R)" if $no == 4;
+ $tree = "((((H, C), O), R), M)" if $no == 5;
+# $tree=~s/;$//g;
+# print "our tree = $tree\n";
+ #---------------------------------------------------------------------------
+ # LOADING HASH CONTAINING ALL POSSIBLE TREES:
+ $tree_decipherer = "/gpfs/home/ydk104/work/rhesus_microsat/codes/lib/tree_analysis_".join("",@tags).".txt";
+ load_allPossibleTrees($tree_decipherer, \%template, \%alternate);
+
+ #---------------------------------------------------------------------------
+ # LOADING THE TREES TO REJECT FOR BIRTH ANALYSIS
+ %treesToReject=();
+ %treesToIgnore=();
+ load_treesToReject(@tags);
+ load_treesToIgnore(@tags);
+ #---------------------------------------------------------------------------
+ # LOADING INPUT DATA INTO HASHES AND ARRAYS
+
+
+ #1 THIS IS THE POINT WHERE WE CAN FILTER OUT LARGE MICROSAT CLUSTERS
+ #2 AS WELL AS MULTIPLE-ALIGNMENT-BLOCKS-SPANNING MICROSATS (KIND OF
+ #3 IMPLICIT IN THE FIRST PART OF THE SENTENCE ITSELF IN MOST CASES).
+
+ my %orths=();
+ my $counterm = 0;
+ my $loaded = 0;
+ my %seen = ();
+ my @allowedchrs = ();
+# print "no = $no\n"; #;
+
+ while (my $line = ){
+ #print "line=$line\n";
+ $line =~ s/>hg18/>H/g;
+ $line =~ s/>panTro2/>C/g;
+ $line =~ s/>ponAbe2/>O/g;
+ $line =~ s/>rheMac2/>R/g;
+ $line =~ s/>calJac1/>M/g;
+ my @micros = split(/>/,$line); # LOADING ALL THE MICROSAT ENTRIES FROM THE CLUSTER INTO @micros
+ #print "micros=",printarr(@micros),"\n"; #;
+ shift @micros; # EMPTYING THE FIRST, EMTPY ELEMENT OF THE ARRAY
+ $no_of_species = adjustCoordinates($micros[0]);
+ next if $no_of_species != $no;
+ $counterm++;
+ #------------------------------------------------
+ $nextcounter[0]++ if $line =~ /compound/;
+ next if $line =~ /compound/; # GETTING RID OF COMPOUND MICROSATS
+ #------------------------------------------------
+ #next if $line =~ /[A-Za-z]>[a-zA-Z]/;
+ #------------------------------------------------
+ chomp $line;
+ my $match_count = ($line =~ s/>/>/g); # COUNTING THE NUMBER OF MICROSAT ENTRIES IN THE CLUSTER
+ #print "number of species = $match_count\n";
+ my $stopper = 0;
+ foreach my $mic (@micros){
+ my @local = split(/\t/,$mic);
+ if ($local[$typecord] =~ /\./ || exists($local[$no_of_interruptionscord+2])) {$stopper = 1; $nextcounter[1]++;
+ last; }
+ # REMOVING CLUSTERS WITH THE CYRPTIC, (UNRESOLVABLY COMPLEX) MICROSAT ENTRIES IN THEM
+ }
+ next if $stopper ==1;
+ #------------------------------------------------
+ $nextcounter[2]++ if (scalar(@micros) >$no_of_species);
+
+ next if (scalar(@micros) >$no_of_species); #1 REMOVING MICROSAT CLUSTERS WITH MORE NUMBER OF MICROSAT ENTRIES THAN THE NUMBER OF SPECIES IN THE DATASET.
+ #2 THIS IS SO BECAUSE SUCH CLUSTERS IMPLY THAT IN AT LEAST ONE SPECIES, THERE IS MORE THAN ONE MICROSAT ENTRY
+ #3 IN THE CLUSTER. THUS, HERE WE ARE GETTING RID OF MICROSATS CLUSTERS THAT INCLUDE MULTUPLE, NEIGHBORING
+ #4 MICROSATS, AND STICK TO CLEAN MICROSATS THAT DO NOT HAVE ANY MICROSATS IN NEIGHBORHOOD.
+ #5 THIS 'NEIGHBORHOOD-RANGE' HAD BEEN DECIDED PREVIOUSLY IN OUR CODE multiSpecies_orthFinder4.pl
+ my $nexter = 0;
+ foreach my $tag (@tags){
+ my $tagcount = ($line =~ s/>$tag\t/>$tag\t/g);
+ if ($tagcount > 1) { $nexter =1; #print colored ['red'],"multiple entires per species : $tagcount of $tag\n" if $printer == 1;
+ next;
+ }
+ }
+
+ if ($nexter == 1){
+ $nextcounter[3]++;
+ next;
+ }
+ #------------------------------------------------
+ foreach my $mic (@micros){ #1 REMOVING MICROSATELLITES WITH ANY 'N's IN THEM
+ my @local = split(/\t/,$mic);
+ if ($local[$microsatcord] =~ /N/) {$stopper =1; $nextcounter[4]++;
+ last;}
+ }
+ next if $stopper ==1;
+ #print "till here 1\n"; #;
+ #------------------------------------------------
+ my @micros_copy = @micros;
+
+ my $tempmicro = shift(@micros_copy); #1 CURRENTLY OBTAINING INFORMATION FOR THE FIRST
+ #2 MICROSAT IN THE CLUSTER.
+ my @tempfields = split(/\t/,$tempmicro);
+ my $prevtype = $tempfields[$typecord];
+ my $tempmotif = $tempfields[$motifcord];
+
+ my $tempfirstmotif = ();
+ if (scalar(@tempfields) > $microsatcord + 2){
+ if ($tempfields[$no_of_interruptionscord] >= 1) { #1 DISCARDING MICROSATS WITH MORE THAN ZERO INTERRUPTIONS
+ #2 IN THE FIRST MICROSAT OF THE CLUSTER
+ $nexter =1; #print colored ['blue'],"more than one interruptions \n" if $printer == 1;
+ }
+ }
+ if ($nexter == 1){
+ $nextcounter[6]++;
+ next;
+ } #1 DONE OBTAINING INFORMATION REGARDING
+ #2 THE FIRST MICROSAT FROM THE CLUSTER
+
+ if ($tempmotif =~ /^\[/){
+ $tempmotif =~ s/^\[//g;
+ $tempmotif =~ /([a-zA-Z]+)\].*/;
+ $tempfirstmotif = $1; #1 OBTAINING THE FIRTS MOTIF OF MICROSAT
+ }
+ else {$tempfirstmotif = $tempmotif;}
+ my $prevmotif = $tempfirstmotif;
+
+ my $key = ();
+ if ($tempmicro =~ /([0-9]+)\s+($focalspec)\s(chr[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)/ ) {
+ $key = join("\t",$1, $2, $4, $5);
+ }
+ else{
+ # print "counld not form a key \n" if $printer == 1;
+ $nextcounter[7]++;
+ next;
+ }
+ #----------------- #1 NOW, AFTER OBTAINING INFORMATION ABOUT
+ #2 THE FIRST MICROSAT IN THE CLUSTER, THE
+ #3 FOLLOWING LOOP GOES THROUGH THE OTHER MICROSATS
+ #4 TO SEE IF THEY SHARE THE REQUIRED FEATURES (BELOW)
+
+ foreach my $micro (@micros_copy){
+ my @fields = split(/\t/,$micro);
+ #-----------------
+ if (scalar(@fields) > $microsatcord + 2){ #1 DISCARDING MICROSATS WITH MORE THAN ONE INTERRUPTIONS
+ if ($fields[$no_of_interruptionscord] >= 1) {$nexter =1; #print colored ['blue'],"more than one interruptions \n" if $printer == 1;
+ $nextcounter[6]++;
+ last; }
+ }
+ #-----------------
+ if (($prevtype ne "0") && ($prevtype ne $fields[$typecord])) {
+ $nexter =1; #print colored ['yellow'],"microsat of different type \n" if $printer == 1;
+ $nextcounter[8]++;
+ last; } #1 DISCARDING MICROSAT CLUSTERS WHERE MICROSATS BELONG
+ #----------------- #2 TO DIFFERENT TYPES (MONOS, DIS, TRIS ETC.)
+ $prevtype = $fields[$typecord];
+
+ my $motif = $fields[$motifcord];
+ my $firstmotif = ();
+
+ if ($motif =~ /^\[/){
+ $motif =~ s/^\[//g;
+ $motif =~ /([a-zA-Z]+)\].*/;
+ $firstmotif = $1;
+ }
+ else {$firstmotif = $motif;}
+
+ my $motifpattern = $firstmotif.$firstmotif;
+ my $prevmotifpattern = $prevmotif.$prevmotif;
+
+ if (($prevmotif ne "0")&&(($motifpattern !~ /$prevmotif/i)||($prevmotifpattern !~ /$firstmotif/i)) ) {
+ $nexter =1; #print colored ['green'],"different motifs used \n$line\n" if $printer == 1;
+ $nextcounter[9]++;
+ last;
+ } #1 DISCARDING MICROSAT CLUSTERS WHERE MICROSATS BELONG
+ #2 TO DIFFERENT MOTIFS
+ my $prevmotif = $firstmotif;
+ #-----------------
+
+ for my $t (0 ... $#tags){ #1 DISCARDING MICROSAT CLUSTERS WHERE MICROSAT ENTRIES BELONG
+ #2 DIFFERENT ALIGNMENT BLOCKS
+ if ($micro =~ /([0-9]+)\s+($focalspec)\s(chr[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)/ ) {
+ my $key2 = join("\t",$1, $2, $4, $5);
+ if ($key2 ne $key){
+# print "microsats belong to diffferent alignment blocks altogether\n" if $printer == 1;
+ $nextcounter[10]++;
+ $nexter = 1; last;
+ }
+ }
+ else{
+ # print "counld not form a key \n" if $printer == 1;
+ $nexter = 1; last;
+ }
+ }
+
+ }
+ #####################
+ if ($nexter == 1){
+ # print "nexting\n" if $printer == 1;
+ next;
+ }
+ else{
+# print "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n$key:\n$line\nvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv\n" if $printer == 1;
+ push (@{$orths{$key}},$line);
+ $loaded++;
+ if ($line =~ /($focalspec)\s([a-zA-Z0-9]+)\s([0-9]+)\s([0-9]+)/ ) {
+# print "$line\n" if $printer == 1; #if $line =~ /Contig/;
+# print "################ ################\n" if $printer == 1;
+ push @allowedchrs, $2 if !exists $allowedhash{$2};
+ $allowedhash{$2} = 1;
+ my $key = join("\t",$1, $2, $3, $4);
+ #print "print the shit: $key\n" if $printer == 1;
+ $seen{$key} = 1;
+ }
+ else { #print "Key could not be formed in SPUT for ($org) ($title) ([0-9]+) ([0-9]+)\n";
+ }
+ }
+ }
+ close ORTH;
+
+# print "now studying where we lost microsatellites: @nextcounter\n";
+ for my $reason (0 ... $#nextcounter){
+# print $removalReasons[$reason]."\t".$nextcounter[$reason],"\n";
+ }
+# print "\ntotal number of keys formed = ", scalar(keys %orths), " = \n";
+# print "done filtering .. counterm = $counterm and loaded = $loaded\n";
+
+ #----------------------------------------------------------------------------------------------------------------
+ # NOW GENERATING THE ALIGNMENT FILE WITH RELELEVENT ALIGNMENTS STORED ONLY.
+
+ while (1){
+ if (-e $megamatchlck){
+# print "waiting to write into $megamatchlck\n";
+ sleep 10;
+ }
+ else{
+ open (MEGAMLCK, ">$megamatchlck") or die "Cannot open megamatchlck file $megamatchlck: $!";
+ open (MEGAM, ">$megamatch") or die "Cannot open megamatch file $megamatch: $!";
+ last;
+ }
+ }
+
+ foreach my $seqfile (@filterseqfiles){
+ my $fullpath = $seqfile;
+
+# print "opening file: $fullpath\n";
+ open (MATCH, "<$fullpath") or die "Cannot open MATCH file $fullpath: $!";
+ my $matchlines = 0;
+
+ while (my $line = ) {
+ if ($line =~ /($focalspec)\s([a-zA-Z0-9]+)\s([0-9]+)\s([0-9]+)/ ) {
+ my $key = join("\t",$1, $2, $3, $4);
+ if (exists $seen{$key}){
+ while (1){
+ $matchlines++;
+ print MEGAM $line;
+ $line = ;
+ print MEGAM "\n" if $line !~/[0-9a-zA-Z]/;
+ last if $line !~/[0-9a-zA-Z]/;
+ }
+ }
+ }
+ }
+# print "matchlines = $matchlines\n";
+ close MATCH;
+ }
+ close MEGAMLCK;
+
+ unlink $megamatchlck;
+ close MEGAM;
+ undef %seen;
+ #----------------------------------------------------------------------------------------------------------------
+
+ #---------------------------------------------------------------------------
+ # NOW, AFTER FILTERING MANY MICROSATS, AND LOADING THE FILTERED ONES INTO
+ # THE HASH %orths , WE GO THROUGH THE ALIGNMENT FILE, AND STUDY THE
+ # FLANKING SEQUENCES OF ALL THESE MICROSATS, TO FILTER THEM FURTHER
+ #$printer = 1;
+
+ my $microreadcounter=0;
+ my $contigsentered=0;
+ my $contignotrightcounter=0;
+ my $keynotformedcounter=0;
+ my $keynotfoundcounter= 0;
+ my $dotcounter = 0;
+
+ open (BO, "<$megamatch") or die "Cannot open alignment file: $megamatch: $!";
+
+ while (my $line = ){
+# print "." if $dotcounter % 100 ==0;
+# print "\n" if $dotcounter % 5000 ==0;
+# print "dotcounter = $dotcounter\n " if $printer == 1;
+ next if $line !~ /^[0-9]+/;
+ $dotcounter++;
+# print colored ['green'], "~" x 60, "\n" if $printer == 1;
+# print colored ['green'], $line;# if $printer == 1;
+ chomp $line;
+ my @fields2 = split(/\t/,$line);
+ my $key2 = ();
+ my $alignment_no = (); #1 TEMPORARY
+ if ($line =~ /([0-9]+)\s+($focalspec)\s(chr[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)/ ) {
+ $key2 = join("\t",$1, $2, $4, $5);
+ $alignment_no=$1;
+ }
+ else {print "seq line $line incompatible\n"; $keynotformedcounter++; next;}
+
+ $no_of_species = adjustCoordinates($line);
+ $contignotrightcounter++ if $no_of_species != $no;
+ # print "contignotrightcounter=$contignotrightcounter\n";
+ # print "no_of_species=$no_of_species\n";
+ # print "no=$no\n";
+
+ next if $no_of_species != $no;
+
+ # print "key = $key2\n" if $printer == 1;
+ my @clusters = (); #1 EXTRACTING MICROSATS CORRESPONDING TO THIS
+ #2 ALIGNMENT BLOCK
+ if (exists($orths{$key2})){
+ @clusters = @{$orths{$key2}};
+ $contigsentered++;
+ delete $orths{$key2};
+ }
+ else{
+ #print "orth does not exist\n";
+ $keynotfoundcounter++;
+ next;
+ }
+
+ my %sequences=(); #1 WILL STORE SEQUENCES IN THE CURRENT ALIGNMENT BLOCK
+ my $humseq = ();
+ foreach my $tag (@tags){ #1 READING THE ALIGNMENT FILE AND CAPTURING SEQUENCES
+ my $seq = ; #2 OF ALL SPECIES.
+ chomp $seq;
+ $sequences{$tag} = " ".$seq;
+ #print "sequences = $sequences{$tag}\n" if $printer == 1;
+ $humseq = $seq if $tag =~ /H/;
+ }
+
+
+ foreach my $cluster (@clusters){ #1 NOW, GOING THROUGH THE CLUSTER OF MICROSATS
+ #print "x" x 60, "\n" if $printer == 1;
+ #print colored ['red'],"cluster = $cluster\n";
+ $largesttree =~ s/hg18/H/g;
+ $largesttree =~ s/panTro2/C/g;
+ $largesttree =~ s/ponAbe2/O/g;
+ $largesttree =~ s/rheMac2/R/g;
+ $largesttree =~ s/calJac1/M/g;
+
+ $microreadcounter++;
+ my @micros = split(/>/,$cluster);
+ shift @micros;
+
+ my $edge_microsat=0; #1 THIS WILL HAVE VALUE "1" IF MICROSAT IS FOUND
+ #2 TO BE TOO CLOSE TO THE EDGES OF ALIGNMENT BLOCK
+
+ my @starts= (); my %start_hash=(); #1 STORES THE START AND END COORDINATES OF MICROSATELLITES
+ my @ends = (); my %end_hash=(); #2 SO THAT LATER, WE WILL BE ABLE TO FIND THE EXTREME
+ #3 COORDINATE VALUES OF THE ORTHOLOGOUS MIROSATELLITES.
+
+ my %microhash=();
+ my %microsathash=();
+ my %nonmicrosathash=();
+ my $motif=(); #1 BASIC MOTIF OF THE MICROSATELLITE.. THERE'S ONLY 1
+ #print "tags=@tags\n";
+ for my $i (0 ... $#tags){ #1 FINDING THE MICROSAT, AND THE ALIGNMENT SEQUENCE
+ #2 CORRESPONDING TO THE PARTICULAR SPECIES (AS PER
+ #3 THE VARIABLE $TAG;
+ my $tag = $tags[$i];
+ # print $seq;
+ my $locus="NULL"; #1 THIS WILL STORE THE MICROSAT OF THIS SPECIES.
+ #2 IF THERE IS NO MICROSAT, IT WILL REMAIN "NULL"
+
+ foreach my $micro (@micros){
+ # print "micro=$micro, tag=$tag\n";
+ if ($micro =~ /^$tag/){ #1 MICROSAT OF THIS SPECIES FOUND..
+ $locus = $micro;
+ my @fields = split(/\t/,$micro);
+ $motif = $fields[$motifcord];
+ $microsathash{$tag}=$fields[$microsatcord];
+ # print "fields=@fields, and startcord=$startcord = $fields[$startcord]\n";
+ push(@starts, $fields[$startcord]);
+ push(@ends, $fields[$endcord]);
+ $start_hash{$tag}=$fields[$startcord];
+ $end_hash{$tag}=$fields[$endcord];
+ last;
+ }
+ else{$microsathash{$tag}="NULL"}
+ }
+ $microhash{$tag}=$locus;
+
+ }
+
+
+
+ my $extreme_start = smallest_number(@starts); #1 THESE TWO ARE THE EXTREME COORDINATES OF THE
+ my $extreme_end = largest_number(@ends); #2 MICROSAT CLUSTER ACCROSS ALL THE SPECIES IN
+ #3 WHOM IT IS FOUND TO BE ORTHOLOGOUS.
+
+ #print "starts=@starts... ends=@ends\n";
+
+ my %up_flanks = (); #1 CONTAINS UPSTEAM FLANKING REGIONS FOR EACH SPECIES
+ my %down_flanks = (); #1 CONTAINS DOWNDTREAM FLANKING REGIONS FOR EACH SPECIES
+
+ my %up_largeflanks = ();
+ my %down_largeflanks = ();
+
+ my %locusandflanks = ();
+ my %locusandlargeflanks = ();
+
+ my %up_internal_flanks=(); #1 CONTAINS SEQUENCE BETWEEN THE $extreme_start and the
+ #2 ACTUAL START OF MICROSATELLITE IN THE SPECIES
+ my %down_internal_flanks=(); #1 CONTAINS SEQUENCE BETWEEN THE $extreme_end and the
+ #2 ACTUAL end OF MICROSATELLITE IN THE SPECIES
+
+ my %alignment=(); #1 CONTAINS ACTUAL ALIGNMENT SEQUENCE BETWEEN THE TWO
+ #2 EXTEME VALUES.
+
+ my %microsatstarts=(); #1 WITHIN EACH ALIGNMENT, IF THERE EXISTS A MICROSATELLITE
+ #2 THIS HASH CONTAINS THE START SITE OF THE MICROSATELLITE
+ #3 WIHIN THE ALIGNMENT
+
+
+ for my $i (0 ... $#tags){ #1 NOW THAT WE HAVE GATHERED INFORMATION REGARDING
+ #2 SEQUENCE ALIGNMENT AND MICROSATELLITE COORDINATES
+ #3 AS WELL AS THE EXTREME COORDINATES OF THE
+ #4 MICROSAT CLUSTER, WE WILL PROCEED TO EXTRACT THE
+ #5 FLANKING SEQUENCE OF ALL ORGS, AND STUDY IT IN
+ #6 MORE DETAIL.
+ my $tag = $tags[$i];
+ # print "tag=$tag.. seqlength = ",length($sequences{$tag})," extreme_start=$extreme_start and extreme_end=$extreme_end\n";
+ my $upstream_gaps = (substr($sequences{$tag}, 0, $extreme_start) =~ s/\-/-/g); #1 NOW MEASURING THE NUMBER OF GAPS IN THE UPSTEAM
+ #2 AND DOWNSTREAM SEQUENCES OF THE MICROSATs IN THIS
+ #3 CLUSTER.
+
+ my $downstream_gaps = (substr($sequences{$tag}, $extreme_end) =~ s/\-/-/g);
+ if (($extreme_start - $upstream_gaps )< $EDGE_DISTANCE || (length($sequences{$tag}) - $extreme_end - $downstream_gaps) < $EDGE_DISTANCE){
+ $edge_microsat=1;
+
+ last;
+ }
+ else{
+ $up_flanks{$tag} = substr($sequences{$tag}, $extreme_start - $FLANK_SUPPORT, $FLANK_SUPPORT);
+ $down_flanks{$tag} = substr($sequences{$tag}, $extreme_end+1, $FLANK_SUPPORT);
+
+ $up_largeflanks{$tag} = substr($sequences{$tag}, $extreme_start - $COMPLEXITY_SUPPORT, $COMPLEXITY_SUPPORT);
+ $down_largeflanks{$tag} = substr($sequences{$tag}, $extreme_end+1, $COMPLEXITY_SUPPORT);
+
+
+ $alignment{$tag} = substr($sequences{$tag}, $extreme_start, $extreme_end-$extreme_start+1);
+ $locusandflanks{$tag} = $up_flanks{$tag}."[".$alignment{$tag}."]".$down_flanks{$tag};
+ $locusandlargeflanks{$tag} = $up_largeflanks{$tag}."[".$alignment{$tag}."]".$down_largeflanks{$tag};
+
+ if ($microhash{$tag} ne "NULL"){
+ $up_internal_flanks{$tag} = substr($sequences{$tag}, $extreme_start , $start_hash{$tag}-$extreme_start);
+ $down_internal_flanks{$tag} = substr($sequences{$tag}, $end_hash{$tag} , $extreme_end-$end_hash{$tag});
+ $microsatstarts{$tag}=$start_hash{$tag}-$extreme_start;
+# print "tag = $tag, internal flanks = $up_internal_flanks{$tag} and $down_internal_flanks{$tag} and start = $microsatstarts{$tag}\n" if $printer == 1;
+ }
+ else{
+ $nonmicrosathash{$tag}=substr($sequences{$tag}, $extreme_start, $extreme_end-$extreme_start+1);
+
+ }
+ # print "up flank for species $tag = $up_flanks{$tag} \ndown flank for species $tag = $down_flanks{$tag} \n" if $printer == 1;
+
+ }
+
+ }
+ $nextcounter[11]++ if $edge_microsat==1;
+ next if $edge_microsat==1;
+
+
+ my $low_complexity = 0; #1 VALUE WILL BE 1 IF ANY OF THE FLANKING REGIONS
+ #2 IS FOUND TO BE OF LOW COMPLEXITY, BY USING THE
+ #3 FUNCTION sub test_complexity
+
+
+ for my $i (0 ... $#tags){
+# print "i = $tags[$i]\n" if $printer == 1;
+ if (test_complexity($up_largeflanks{$tags[$i]}, $COMPLEXITY_SUPPORT) eq "LOW" || test_complexity($down_largeflanks{$tags[$i]}, $COMPLEXITY_SUPPORT) eq "LOW"){
+# print "i = $i, low complexity regions: $up_largeflanks{$tags[$i]}: ",test_complexity($up_largeflanks{$tags[$i]}, $COMPLEXITY_SUPPORT), " and $down_largeflanks{$tags[$i]} = ",test_complexity($down_largeflanks{$tags[$i]}, $COMPLEXITY_SUPPORT),"\n" if $printer == 1;
+ $low_complexity =1; last;
+ }
+ }
+
+ $nextcounter[12]++ if $low_complexity==1;
+ next if $low_complexity == 1;
+
+
+ my $sequence_dissimilarity = 0; #1 THIS VALYE WILL BE 1 IF THE SEQUENCE SIMILARITY
+ #2 BETWEEN ANY OF THE SPECIES AGAINST THE HUMAN
+ #3 FLANKING SEQUENCES IS BELOW A CERTAIN THRESHOLD
+ #4 AS DESCRIBED IN FUNCTION sub sequence_similarity
+ my %donepair = ();
+ for my $i (0 ... $#tags){
+ # print "i = $tags[$i]\n" if $printer == 1;
+# next if $i == 0;
+ # print colored ['magenta'],"THIS IS UP\n" if $printer == 1;
+
+ for my $b (0 ... $#tags){
+ next if $b == $i;
+ my $pair = ();
+ $pair = $i."_".$b if $i < $b;
+ $pair = $b."_".$i if $b < $i;
+ next if exists $donepair{$pair};
+ my ($up_similarity,$upnucdiffs, $upindeldiffs) = sequence_similarity($up_flanks{$tags[$i]}, $up_flanks{$tags[$b]}, $SIMILARITY_THRESH, $info);
+ my ($down_similarity,$downnucdiffs, $downindeldiffs) = sequence_similarity($down_flanks{$tags[$i]}, $down_flanks{$tags[$b]}, $SIMILARITY_THRESH, $info);
+ $donepair{$pair} = $up_similarity."_".$down_similarity;
+
+# print RUN "$up_similarity $upnucdiffs $upindeldiffs $down_similarity $downnucdiffs $downindeldiffs\n";
+
+ if ( $up_similarity < $SIMILARITY_THRESH || $down_similarity < $SIMILARITY_THRESH){
+ $sequence_dissimilarity =1;
+ last;
+ }
+ }
+ }
+ $nextcounter[13]++ if $sequence_dissimilarity==1;
+
+ next if $sequence_dissimilarity == 1;
+ my ($simplified_microsat, $Hchrom, $Hstart, $Hend, $locusmotif, $locusmotifsize) = summarize_microsat($cluster, $humseq);
+ # print "simplified_microsat=$simplified_microsat\n"; ;
+ my ($tree_analysis, $alternative_trees, $conformation) = treeStudy($simplified_microsat);
+
+ if (exists $treesToReject{$tree_analysis}){
+ $nextcounter[14]++;
+ next;
+ }
+
+# my $adjuster=();
+# if ($no_of_species == 4){
+# my @sields = split(/\t/,$simplified_microsat);
+# my $somend = pop(@sields);
+# my $somestart = pop(@sields);
+# my $somechr = pop(@sields);
+# $adjuster = "NA\t" x 13 ;
+# $simplified_microsat = join ("\t", @sields, $adjuster).$somechr."\t".$somestart."\t".$somend;
+# }
+# if ($no_of_species == 3){
+# my @sields = split(/\t/,$simplified_microsat);
+# my $somend = pop(@sields);
+# my $somestart = pop(@sields);
+# my $somechr = pop(@sields);
+# $adjuster = "NA\t" x 26 ;
+# $simplified_microsat = join ("\t", @sields, $adjuster).$somechr."\t".$somestart."\t".$somend;
+# }
+#
+ $registeredTrees{$tree_analysis} = 1 if !exists $registeredTrees{$tree_analysis};
+ $registeredTrees{$tree_analysis}++ if exists $registeredTrees{$tree_analysis};
+
+ if (exists $treesToIgnore{$tree_analysis}){
+ my @appendarr = ();
+
+ print SUMMARY $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize], "\t", $simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
+ #print "SUMMARY ",$Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize], "\t", $simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
+# print SELECT $Hchrom,"\t",$Hstart,"\t",$Hend,"\t","NOEVENT", "\t\t", $cluster,"\n";
+
+ foreach my $lnode (@$lagestnodes){
+ my @pair = @$lnode;
+ my @nodemutarr = ();
+ for my $p (@pair){
+ my @mutinfoarray1 = ();
+ for (1 ... 38){
+ push (@mutinfoarray1, "NA")
+ }
+ print SUMMARY join ("\t", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
+ }
+
+ }
+ for (1 ... 38){
+ push (@appendarr, "NA")
+ }
+ print SUMMARY join ("\t", @appendarr,"NULL", "NULL"),"\n";
+ # print "SUMMARY ",join ("\t", @appendarr,"NULL", "NULL"),"\n"; #;
+ next;
+ }
+
+ my ($mutations_array, $nodes, $branches_hash, $alivehash, $primaryalignment) = peel_onion($tree, \%sequences, \%alignment, \@tags, \%microsathash, \%nonmicrosathash, $motif, $tree_analysis, $thresholdhash{length($motif)}, \%microsatstarts);
+
+ if ($mutations_array eq "NULL"){
+ my @appendarr = ();
+
+ print SUMMARY $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize],"\t",$simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
+ # print "SUMMARY ", $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize],"\t",$simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
+# print SELECT $Hchrom,"\t",$Hstart,"\t",$Hend,"\t","EVENT", "\t\t", $cluster,"\n";
+
+ foreach my $lnode (@$lagestnodes){
+ my @pair = @$lnode;
+ my @nodemutarr = ();
+ for my $p (@pair){
+ my @mutinfoarray1 = ();
+ for (1 ... 38){
+ push (@mutinfoarray1, "NA")
+ }
+ print SUMMARY join ("\t", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
+ # print join ("\t", "SUMMARY", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
+ }
+
+ }
+ for (1 ... 38){
+ push (@appendarr, "NA")
+ }
+ print SUMMARY join ("\t", @appendarr,"NULL", "NULL"),"\n";
+ # print join ("\t","SUMMARY", @appendarr,"NULL", "NULL"),"\n"; #;
+ next;
+ }
+
+
+# print "sent: \n" if $printer == 1;
+# print "nodes = @$nodes, branches array:\n" if $mutations_array ne "NULL" && $printer == 1;
+
+ my ($newstatus, $newmutations_array, $newnodes, $newbranches_hash, $newalivehash, $finalalignment) = fillAlignmentGaps($tree, \%sequences, \%alignment, \@tags, \%microsathash, \%nonmicrosathash, $motif, $tree_analysis, $thresholdhash{length($motif)}, \%microsatstarts);
+# print "newmutations_array returned = \n",join("\n",@$newmutations_array),"\n" if $newmutations_array ne "NULL" && $printer == 1;
+ my @finalmutations_array= ();
+ @finalmutations_array = selectMutationArray($mutations_array, $newmutations_array, \@tags, $alivehash, \%alignment, $motif) if $newmutations_array ne "NULL";
+ @finalmutations_array = selectMutationArray($mutations_array, $mutations_array, \@tags, $alivehash, \%alignment, $motif) if $newmutations_array eq "NULL";
+
+ my ($besttree, $treescore) = selectBetterTree($tree_analysis, $alternate{$tree_analysis}, \@finalmutations_array);
+ my $cleancase = "UNCLEAN";
+
+ $cleancase = checkCleanCase($besttree, $finalalignment) if $treescore > 0 && $finalalignment ne "NULL" && $finalalignment =~ /\!/;
+ $cleancase = checkCleanCase($besttree, $primaryalignment) if $treescore > 0 && $finalalignment eq "NULL" && $primaryalignment =~ /\!/ && $primaryalignment ne "NULL";
+ $cleancase = "CLEAN" if $finalalignment eq "NULL" && $primaryalignment !~ /\!/ && $primaryalignment ne "NULL";
+ $cleancase = "CLEAN" if $finalalignment ne "NULL" && $finalalignment !~ /\!/ ;
+ $besttree = "NULL" if $treescore <= 0;
+ print SUMMARY $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize],"\t",$simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
+ # print "SUMMARY ", $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize],"\t",$simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
+
+# print SELECT $Hchrom,"\t",$Hstart,"\t",$Hend,"\t","EVENT", "\t\t", $cluster,"\n";
+
+ my @mutinfoarray =();
+
+ foreach my $lnode (@$lagestnodes){
+ my @pair = @$lnode;
+ my $joint = "(".join(", ",@pair).")";
+ my @nodemutarr = ();
+
+ for my $p (@pair){
+ foreach my $mut (@finalmutations_array){
+ $mut =~ /node=([A-Z, \(\)]+)/;
+ push @nodemutarr, $mut if $p eq $1;
+ }
+# print "from pair @pair, p=$p\n";
+ @mutinfoarray = summarizeMutations(\@nodemutarr, $besttree);
+ print SUMMARY join ("\t", @mutinfoarray[0...($#mutinfoarray-1)] ),"\t";
+ # print "SUMMARY ",join ("\t", @mutinfoarray[0...($#mutinfoarray-1)] ),"\t";
+ }
+
+ }
+
+ @mutinfoarray = summarizeMutations(\@finalmutations_array, $besttree);
+ print SUMMARY join ("\t", @mutinfoarray ),"\t";
+ print SUMMARY $cleancase, "\n";
+ # print "SUMMARY ",join ("\t", @mutinfoarray,$cleancase ),"\n"; #;
+ # print "summarized\n"; ;
+
+ my %indelcatch = ();
+ my %substcatch = ();
+ my %typecatch = ();
+ my %nodescatch = ();
+ my $mutconcat = join("\t", @finalmutations_array)."\n";
+ my %indelposcatch = ();
+ my %subsposcatch = ();
+
+ foreach my $fmut ( @finalmutations_array){
+# next if $fmut !~ /indeltype=[a-zA-Z]+/;
+ #print RUN $fmut, "\n";
+ $fmut =~ /node=([a-zA-Z, \(\)]+)/;
+ my $lnode = $1;
+ $nodescatch{$1}=1;
+
+ if ($fmut =~ /type=substitution/){
+ # print "fmut=$fmut\n";
+ $fmut =~ /from=([a-zA-Z\-]+)\tto=([a-zA-Z\-]+)/;
+ my $from=$1;
+ # print "from=$from\n";
+ my $to=$2;
+ # print "to=$to\n";
+ push @{$substcatch{$lnode}} , ("from:".$from." to:".$to);
+ $fmut =~ /position=([0-9]+)/;
+ push @{$subsposcatch{$lnode}}, $1;
+ }
+
+ if ($fmut =~ /insertion=[a-zA-Z\-]+/){
+ $fmut =~ /insertion=([a-zA-Z\-]+)/;
+ push @{$indelcatch{$lnode}} , $1;
+ $fmut =~ /indeltype=([a-zA-Z]+)/;
+ push @{$typecatch{$lnode}}, $1;
+ $fmut =~ /position=([0-9]+)/;
+ push @{$indelposcatch{$lnode}}, $1;
+ }
+ if ($fmut =~ /deletion=[a-zA-Z\-]+/){
+ $fmut =~ /deletion=([a-zA-Z\-]+)/;
+ push @{$indelcatch{$lnode}} , $1;
+ $fmut =~ /indeltype=([a-zA-Z]+)/;
+ push @{$typecatch{$lnode}}, $1;
+ $fmut =~ /position=([0-9]+)/;
+ push @{$indelposcatch{$lnode}}, $1;
+ }
+ }
+
+ # print $simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t" if $printer == 1;
+ # print join ("<\t>", @mutinfoarray),"\n" if $printer == 1;
+ # print "where mutinfoarray = @mutinfoarray\n" if $printer == 1;
+ # #print RUN ".";
+
+ # print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
+ # print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
+
+ # print colored ['red'],"finalmutations_array=\n" if $printer == 1;
+ foreach (@finalmutations_array) {
+# print colored ['red'], "$_\n" if $_ =~ /type=substitution/ && $printer == 1 ;
+# print colored ['yellow'], "$_\n" if $_ !~ /type=substitution/ && $printer == 1 ;
+
+ }# if $line =~ /cal/;# && $line =~ /chr4/;
+
+# print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
+# print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
+# print "tree analysis = $tree_analysis\n" if $printer == 1;
+
+ # my $mutations = "@$mutations_array";
+
+
+ next;
+ for my $keys (@$nodes) {foreach my $key (@$keys){
+ #print "key = $key, => $branches_hash->{$key}\n";
+ }
+ # print "x" x 50, "\n";
+ }
+ my ($birth_steps, $death_steps) = decipher_history($mutations_array,join("",@tags),$nodes,$branches_hash,$tree_analysis,$conformation, $alivehash, $simplified_microsat);
+ }
+ }
+ close BO;
+# print "now studying where we lost microsatellites:";
+# print "x" x 60,"\n";
+ for my $reason (0 ... $#nextcounter){
+# print $removalReasons[$reason]."\t".$nextcounter[$reason],"\n";
+ }
+# print "x" x 60,"\n";
+# print "In total we read $microreadcounter microsatellites after reading through $contigsentered contigs\n";
+# print " we lost $keynotformedcounter contigs as they did not form the key, \n";
+# print "$contignotrightcounter contigs as they were not of the right species configuration\n";
+# print "$keynotfoundcounter contigs as they did not contain the microsats\n";
+# print "... In total we went through a file that had $dotcounter contigs...\n";
+# print join ("\n","remaining orth keys = ", (keys %orths),"");
+
+
+# print "now printing counted trees: \n";
+ if (scalar(keys %registeredTrees) > 0){
+ foreach my $keyb ( sort (keys %registeredTrees) )
+ {
+# print "$keyb : $registeredTrees{$keyb}\n";
+ }
+ }
+
+
+}
+
+my @summarizarr = ("+C=+C +R.+C -HCOR,+C",
+"+H=+H +R.+H -HCOR,+H",
+"-C=-C -R.-C +HCOR,-C",
+"-H=-H -R.-H +HCOR,-H",
+"+HC=+HC",
+"-HC=-HC",
+"+O=+O -HCOR,+O",
+"-O=-O +HCOR,-O",
+"+HCO=+HCO",
+"-HCO=-HCO",
+"+R=+R +R.+C +R.+H",
+"-R=-R -R.-C -R.-H");
+
+foreach my $line (@summarizarr){
+ next if $line !~ /[A-Za-z0-9]/;
+# print $line;
+ chomp $line;
+ my @fields = split(/=/,$line);
+# print "title = $fields[0]\n";
+ my @parts=split(/ +/, $fields[1]);
+ my %partshash = ();
+ foreach my $part (@parts){$partshash{$part}=1;}
+ my $count=0;
+ foreach my $key ( sort keys %registeredTrees ){
+ next if !exists $partshash{$key};
+# print "now adding $registeredTrees{$key} from $key\n";
+ $count+=$registeredTrees{$key};
+ }
+# print "$fields[0] : $count\n";
+}
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+sub largest_number{
+ my $counter = 0;
+ my($max) = shift(@_);
+ foreach my $temp (@_) {
+
+ #print "finding largest array: $maxcounter \n";
+ if($temp > $max){
+ $max = $temp;
+ }
+ }
+ return($max);
+}
+
+sub smallest_number{
+ my $counter = 0;
+ my($min) = shift(@_);
+ foreach my $temp (@_) {
+ #print "finding largest array: $maxcounter \n";
+ if($temp < $min){
+ $min = $temp;
+ }
+ }
+ return($min);
+}
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+sub baseml_parser{
+ my $outputfile = $_[0];
+ open(BOUT,"<$outputfile") or die "Cannot open output of upstream baseml $outputfile: $!";
+ my @info = ();
+ my @branchields = ();
+ my @distanceields = ();
+ my @bout = ;
+ #print colored ['red'], @bout ,"\n";
+ for my $b (0 ... $#bout){
+ my $bine=$bout[$b];
+ #print colored ['yellow'], "sentence = ",$bine;
+ if ($bine =~ /TREE/){
+ $bine=$bout[$b++];
+ $bine=$bout[$b++];
+ $bine=$bout[$b++];
+ #print "FOUND",$bine;
+ chomp $bine;
+ $bine =~ s/^\s+//g;
+ @branchields = split(/\s+/,$bine);
+ $bine=$bout[$b++];
+ chomp $bine;
+ $bine =~ s/^\s+//g;
+ @distanceields = split(/\s+/,$bine);
+ #print "LASTING..............\n";
+ last;
+ }
+ else{
+ }
+ }
+
+ close BOUT;
+# print "branchfields = @branchields and distanceields = @distanceields\n" if $printer == 1;
+ my %distance_hash=();
+ for my $d (0 ... $#branchields){
+ $distance_hash{$branchields[$d]} = $distanceields[$d];
+ }
+
+ $info[0] = $distance_hash{"9..1"} + $distance_hash{"9..2"};
+ $info[1] = $distance_hash{"9..1"} + $distance_hash{"8..9"}+ $distance_hash{"8..3"};
+ $info[2] = $distance_hash{"9..1"} + $distance_hash{"8..9"}+$distance_hash{"7..8"}+$distance_hash{"7..4"};
+ $info[3] = $distance_hash{"9..1"} + $distance_hash{"8..9"}+$distance_hash{"7..8"}+$distance_hash{"6..7"}+$distance_hash{"6..5"};
+
+# print "\nsending back: @info\n" if $printer == 1;
+
+ return join("\t",@info);
+
+}
+
+
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+sub test_complexity{
+ my $printer = 0;
+ my $sequence = $_[0];
+ my $COMPLEXITY_SUPPORT = $_[1];
+ my $complexity=int($COMPLEXITY_SUPPORT * (1/40)); #1 THIS IS AN ARBITRARY THRESHOLD SET FOR LOW COMPLEXITY.
+ #2 THE INSPIRATION WAS WEB MILLER'S MAIL SENT ON
+ #3 19 Apr 2008 WHERE HE CLASSED AS HIGH COMPLEXITY
+ #4 REGION, IF 40 BP OF SEQUENCE HAS AT LEAST 3 OF
+ #5 EACH NUCLEOTIDE. HENCE, I NORMALIZE THIS PARAMETER
+ #6 FOR THE ACTUAL LENGTH OF $FLANK_SUPPORT SET BY
+ #7 THE USER.
+ #8 WEB MILLER SENT THE MAIL TO YDK104@PSU.EDU
+
+
+
+ my $As = ($sequence=~ s/A/A/gi);
+ my $Ts = ($sequence=~ s/T/T/gi);
+ my $Gs = ($sequence=~ s/G/G/gi);
+ my $Cs = ($sequence=~ s/C/C/gi);
+ #print "seq = $sequence, As=$As, Ts=$Ts, Gs=$Gs, Cs=$Cs\n" if $printer == 1;
+
+ my $ans = ();
+
+ return "HIGH" if $As >= $complexity && $Ts >= $complexity && $Cs >= $complexity && $Gs >= $complexity;
+
+ my @nts = ("A","T","G","C","-");
+
+ my $lowcomplex = 0;
+
+ foreach my $nt (@nts){
+ $lowcomplex =1 if $sequence =~ /(($nt\-*){10,})/i;
+# print "caught with a mono of $nt : $1 in $sequence\n" if $sequence =~ /(($nt\-*){10,})/i;
+ $lowcomplex =1 if $sequence =~ /(($nt[A-Za-z]){10,})/i;
+ $lowcomplex =1 if $sequence =~ /(([A-Za-z]$nt){10,})/i;
+# print "caught with a di with $nt : $2 in $sequence\n" if $sequence =~ /(($nt[A-Za-z]){10,})/i || $sequence =~ /(([A-Za-z]$nt){10,})/i;
+ my $nont = ($sequence=~ s/$nt/$nt/gi);
+
+ }
+# print "leaving for now.. $sequence\n" if $printer == 1 && $lowcomplex == 0;
+ #;
+ return "HIGH" if $lowcomplex == 0;
+ return "LOW" ;
+}
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+sub sequence_similarity{
+ my $printer = 0;
+ my @seq1 = split(/\s*/, $_[0]);
+ my @seq2 = split(/\s*/, $_[1]);
+ my $similarity_thresh = $_[2];
+ my $info = $_[3];
+# print "input = @_\n" if $printer == 1;
+ my $seq1str = $_[0];
+ my $seq2str = $_[1];
+ $seq1str=~s/\-//g; $seq2str=~s/\-//g;
+ my $similarity=0;
+
+ my $nucdiffs=0;
+ my $nucsims=0;
+ my $indeldiffs=0;
+
+ for my $i (0...$#seq1){
+ $similarity++ if $seq1[$i] =~ /$seq2[$i]/i ; #|| $seq1[$i] =~ /\-/i || $seq2[$i] =~ /\-/i ;
+ $nucsims++ if $seq1[$i] =~ /$seq2[$i]/i && ($seq1[$i] =~ /[a-zA-Z]/i && $seq2[$i] =~ /[a-zA-Z]/i);
+ $nucdiffs++ if $seq1[$i] !~ /$seq2[$i]/i && ($seq1[$i] =~ /[a-zA-Z]/i && $seq2[$i] =~ /[a-zA-Z]/i);
+ $indeldiffs++ if $seq1[$i] !~ /$seq2[$i]/i && $seq1[$i] =~ /\-/i || $seq2[$i] =~ /\-/i;
+ }
+ my $sim = $similarity/length($_[0]);
+ return ( $sim, $nucdiffs, $indeldiffs ); #<= $similarity_thresh;
+}
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+
+sub load_treesToReject{
+ my @rejectlist = ();
+ my $alltags = join("",@_);
+ @rejectlist = qw (-HCOR +HCOR) if $alltags eq "HCORM";
+ @rejectlist = qw ( -HCO|+R +HCO|-R) if $alltags eq "HCOR";
+ @rejectlist = qw ( -HC|+O +HC|-O) if $alltags eq "HCO";
+
+ %treesToReject=();
+ $treesToReject{$_} = $_ foreach (@rejectlist);
+ #print "loaded to reject for $alltags; ", $treesToReject{$_},"\n" foreach (@rejectlist); #;
+}
+#--------------------------------------------------------------------------------------------------------
+sub load_treesToIgnore{
+ my @rejectlist = ();
+ my $alltags = join("",@_);
+ @rejectlist = qw (-HCOR +HCOR +HCORM -HCORM) if $alltags eq "HCORM";
+ @rejectlist = qw ( -HCO|+R +HCO|-R +HCOR -HCOR) if $alltags eq "HCOR";
+ @rejectlist = qw ( -HC|+O +HC|-O +HCO -HCO) if $alltags eq "HCO";
+
+ %treesToIgnore=();
+ $treesToIgnore{$_} = $_ foreach (@rejectlist);
+ #print "loaded ", $treesToIgnore{$_},"\n" foreach (@rejectlist);
+}
+#--------------------------------------------------------------------------------------------------------
+sub load_thresholds{
+ my @threshold_array=split(/[,_]/,$_[0]);
+ unshift @threshold_array, "0";
+ for my $size (1 ... 4){
+ $thresholdhash{$size}=$threshold_array[$size];
+ }
+}
+#--------------------------------------------------------------------------------------------------------
+sub load_allPossibleTrees{
+ #1 THIS FILE STORES ALL POSSIBLE SCENARIOS OF MICROSATELLITE
+ #2 BIRTH AND DEATH EVENTS ON A 5-PRIMATE TREE OF H,C,O,R,M
+ #3 IN FORM OF A TEXT FILE. THIS WILL BE USED AS A TEMPLET
+ #4 TO COMPARE EACH MICROSATELLITE CLUSTER TO UNDERSTAND THE
+ #5 EVOLUTION OF EACH LOCUS. WE WILL THEN DISCARD SOME
+ #6 MICROSATS ACCRODING TO THEIR EVOLUTIONARY BEHAVIOUR ON
+ #7 THE TREE. MOST PROBABLY WE WILL REMOVE THOSE MICROSATS
+ #8 THAT ARE NOT SUFFICIENTLY INFORMATIVE, LIKE IN CASE OF
+ #9 AN OUTGROUP MICROSATELLITE BEING DIFFERENT FRON ALL OTHER
+ #10 SPECIES IN THE TREE.
+ my $tree_list = $_[0];
+# print "file to be loaded: $tree_list\n";
+
+ my @trarr = ();
+ @trarr = ("#H C O CONCLUSION ALTERNATE",
+"+ + + +HCO NA",
+"+ _ _ +H NA",
+"_ + _ +C NA",
+"_ _ + -HC|+O NA",
+"+ _ + -C +H",
+"_ + + -H +C",
+"+ + _ +HC|-O NA",
+"_ _ _ -HCO NA") if $tree_list =~ /_HCO\.txt/;
+ @trarr = ("#H C O R CONCLUSION ALTERNATE",
+"_ _ _ _ -HCOR NA",
+"+ + + + +HCOR NA",
+"+ + + _ +HCO|-R +H.+C.+O",
+"+ + _ _ +HC +H.+C;-O",
+"+ _ _ _ +H +HC,-C",
+"_ + _ _ +C +HC,-H",
+"_ _ + _ +O -HC|-H.-C",
+"_ _ + + -HC -H.-C",
+"+ _ _ + +H|-C.-O +HC,-C",
+"_ + _ + +C -H.-O",
+"_ + + _ -H +C.+O",
+"_ _ _ + -HCO|+R NA",
+"+ _ + _ +H.+O|-C NA",
+"_ + + + -H -HC,+C",
+"+ _ + + -C -HC,+H",
+"+ + _ + -O +HC") if $tree_list =~ /_HCOR\.txt/;
+
+ @trarr = ("#H C O R M CONCLUSION ALTERNATE",
+"+ + + + _ +HCOR NA",
+"+ + + _ + -R +HCO;+HC.+O;+H.+C.+O",
+"+ + _ + + -O -HCO,+HC|-HCO,+HC;-HCO,(+H.+C)",
+"+ _ + + + -C -HC,+H;+HCO,(+H.+O)",
+"_ + + + + -H -HC,+C;-HCO,(+C.+O)",
+"_ _ _ _ + -HCOR NA",
+"_ _ _ + _ +R -HC.-O;-H.-C.-O",
+"_ _ + _ _ +O +HCO,-HC;+HCO,(-H.-C)",
+"_ + _ _ _ +C +HC,-H;+HCO,(-H.-O)",
+"+ _ _ _ _ +H +HC,-C;+HCO,(-C.-O)",
+"+ + + _ _ +HCO +H.+C.+O",
+"+ + _ + _ -O +R.+HC|-HCO,+HC;+H.+C.+R|-HCO,(+H.+C)",
+"+ _ + + _ -C -HC,+H;+H.+O.+R|-HCO,(+H.+O)",
+"_ + + + _ -H -HC,+C;+C.+O.+R|-HCO,(+C.+O)",
+"_ _ _ + + -HCO -HC.-O;-H.-C.-O",
+"_ _ + _ + +O +HCO,-HC;+HCO,(-H.-C)",
+"_ + _ _ + +C +HC,-H;+HCO,(-H.-O)",
+"+ _ _ _ + +H -HC,+H;+HCO,(-C.-O)",
+"+ + _ _ + +HC -R.-O|+HCO,-O|+H.+C;-HCO,+HC;-HCO,(+H.+C)",
+"+ _ + _ + -R.-C|+HCO,-C|+H.+O NA",
+"_ + + _ + -R.-H|+HCO,-H|+C.+O NA",
+"_ _ + + _ -HC +R.+O|-HCO,+O|+HCO,-HC",
+"_ + _ + _ +R.+C|-HCO,+C|-HC,+C +HCO,(-H.-O)",
+"+ _ _ + _ +R.+H|-C.-O +HCO,(-C.-O)",
+"+ _ _ + + -O.-C|-HCO,+H +R.+H;-HCO,(+R.+H)",
+"_ + _ + + -O.-H|-HCO,+C +R.+C;-HCO,(+R.+C)",
+"_ + + _ _ +HCO,-H|+O.+C NA",
+"+ _ + _ _ +HCO,-C|+O.+H NA",
+"_ _ + + + -HC -H.-C|-HCO,+O",
+"+ + _ _ _ +HC +H.+C|+HCO,-O|-HCO,+HC;-HCO,(+H.+C)",
+"+ + + + + +HCORM NA") if $tree_list =~ /_HCORM\.txt/;
+
+
+ my $template_p = $_[1];
+ my $alternate_p = $_[2];
+ #1 THIS IS THE HASH IN WHICH INFORMATION FROM THE ABOVE FILE
+ #2 GETS STORED, USING THE WHILE LOOP BELOW. HERE, THE KEY
+ #3 OF EACH ROW IS THE EVOLUTIONARY CONFIGURATION OF A LOCUS
+ #4 ON THE PRIMATE TREE, BASED ON PRESENCE/ABSENCE OF A MICROSAT
+ #5 AT THAT LOCUS, LIKE SAY "+ + + _ _" .. EACH COLUMN BELONGS
+ #6 TO ONE SPECIES; HERE THE COLUMN NAMES ARE "H C O R M".
+ #7 THE VALUE FOR EACH ENTRY IS THE MEANING OF THE ABOVE
+ #8 CONFIGURATION (I.E., CONFIGURAION OF THE KEY. HERE, THE
+ #9 VALUE WILL BE +HCO, SIGNIFYING A BIRTH IN HUMAN-CHIMP-ORANG
+ #10 COMMON ANCESTOR. THIS HASH HAS BEEN LOADED HERE TO BE USED
+ #11 LATER BY THE SUBROUTINE sub treeStudy{} THAT STUDIES
+ #12 EVOLUTIONARY CONFIGURAION OF EACH MICROSAT LOCUS, AS
+ #13 MENTIONED ABOVE.
+ my @keys_array=();
+ foreach my $line (@trarr){
+ next if $line =~ /^#/;
+ chomp $line;
+ my @fields = split("\t", $line);
+ push @keys_array, $fields[0];
+# print "loading: $fields[0]\n";
+ $template_p->{$fields[0]}[0] = $fields[1];
+ $template_p->{$fields[0]}[1] = 0;
+ $alternate_p->{$fields[1]} = $fields[2];
+
+ }
+# print "loaded the trees with keys: @keys_array\n";
+ return $template_p, \@keys_array, $alternate_p;
+}
+
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+sub checkCleanCase{
+ my $printer = 0;
+ my $tree = $_[0];
+ my $finalalignment = $_[1];
+
+ #print "IN checkCleanCase: @_\n";
+ #;
+ my @indivspecies = $tree =~ /[A-Z]/g;
+ $finalalignment =~ s/\./_/g;
+ my @captured = $finalalignment =~ /[A-Za-z, \(\):]+\![:A-Za-z, \(\)]/g;
+
+ my $unclean = 0;
+
+ foreach my $sp (@indivspecies){
+ foreach my $cap (@captured){
+ $cap =~ s/:[A-Za-z\-]+//g;
+ my @sps = $cap =~ /[A-Z]+/g;
+ my $spsc = join("", @sps);
+# print "checking whether imp species $sp is present in $cap i.e, in $spsc\n " if $printer == 1;
+ if ($spsc =~ /$sp/){
+# print "foind : $sp\n";
+ $unclean = 1; last;
+ }
+ }
+ last if $unclean == 1;
+ }
+ #;
+ return "CLEAN" if $unclean == 0;
+ return "UNCLEAN";
+}
+
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+#--------------------------------------------------------------------------------------------------------
+
+
+sub adjustCoordinates{
+ my $line = $_[0];
+ my $no_of_species = $line =~ s/(chr[0-9a-zA-Z]+)|(Contig[0-9a-zA-Z\._\-]+)/x/g;
+ $infocord = 2 + (4*$no_of_species) - 1;
+ $typecord = 2 + (4*$no_of_species) + 1 - 1;
+ $motifcord = 2 + (4*$no_of_species) + 2 - 1;
+ $gapcord = $motifcord+1;
+ $startcord = $gapcord+1;
+ $strandcord = $startcord+1;
+ $endcord = $strandcord + 1;
+ $microsatcord = $endcord + 1;
+ $sequencepos = 2 + (5*$no_of_species) + 1 -1 ;
+ $interr_poscord = $microsatcord + 3;
+ $no_of_interruptionscord = $microsatcord + 4;
+ $interrcord = $microsatcord + 2;
+ #print "$line\n startcord = $startcord, and endcord = $endcord and no_of_species = $no_of_species\n" if $printer == 1;
+ return $no_of_species;
+}
+
+
+sub printhash{
+ my $alivehash = $_[0];
+ my @tags = @$_[1];
+# print "print hash\n";
+ foreach my $tag (@tags){
+# print "$tag=",$alivehash->{$tag},"\n" if exists $alivehash->{$tag};
+ }
+
+ return "\n"
+}
+sub peel_onion{
+ my $printer = 0;
+# print "received: @_\n" ; #;
+ $printer = 0;
+ my ($tree, $sequences, $alignment, $tagarray, $microsathash, $nonmicrosathash, $motif, $tree_analysis, $threshold, $microsatstarts) = @_;
+# print "in peel onion.. tree = $tree \n" if $printer == 1;
+ my %sequence_hash=();
+
+
+# for my $i (0 ... $#sequences){ $sequence_hash{$species[$i]}=$sequences->[$i]; }
+
+
+ my %node_sequences=();
+
+ my %node_alignments = (); #NEW, Nov 28 2008
+ my @tags=();
+ my @locus_sequences=();
+ my %alivehash=();
+ foreach my $tag (@$tagarray) {
+ #print "adding: $tag\n";
+ push(@tags, $tag);
+ $node_sequences{$tag}=join ".",split(/\s*/,$microsathash->{$tag}) if $microsathash->{$tag} ne "NULL";
+ $alivehash{$tag}= $tag if $microsathash->{$tag} ne "NULL";
+ $node_sequences{$tag}=join ".",split(/\s*/,$nonmicrosathash->{$tag}) if $microsathash->{$tag} eq "NULL";
+ $node_alignments{$tag}=join ".",split(/\s*/,$alignment->{$tag}) ;
+ push @locus_sequences, $node_sequences{$tag};
+ #print "adding to node_seq: $tag = ",$node_alignments{$tag},"\n";
+ }
+
+
+
+ my ($nodes_arr, $branches_hash) = get_nodes($tree);
+ my @nodes=@$nodes_arr;
+# print "recieved nodes = " if $printer == 1;
+# foreach my $key (@nodes) {print "@$key " if $printer == 1;}
+
+# print "\n" if $printer == 1;
+
+ #POPULATE branches_hash WITH INFORMATION ABOUT LIVESTATUS
+ foreach my $keys (@nodes){
+ my @pair = @$keys;
+ my $joint = "(".join(", ",@pair).")";
+ my $copykey = join "", @pair;
+ $copykey =~ s/[\W ]+//g;
+# print "for node: $keys, copykey = $copykey and joint = $joint\n" if $printer == 1;
+ my $livestatus = 1;
+ foreach my $copy (split(/\s*/,$copykey)){
+ $livestatus = 0 if !exists $alivehash{$copy};
+ }
+ $alivehash{$joint} = $joint if !exists $alivehash{$joint} && $livestatus == 1;
+# print "alivehash = $alivehash{$joint}\n" if exists $alivehash{$joint} && $printer == 1;
+ }
+
+ @nodes = reverse(@nodes); #1 THIS IS IN ORDER TO GO THROUGH THE TREE FROM LEAVES TO ROOT.
+
+ my @mutations_array=();
+
+ my $joint = ();
+ foreach my $node (@nodes){
+ my @pair = @$node;
+# print "now in the nodes for loop, pair = @pair\n and sequences=\n" if $printer == 1;
+ $joint = "(".join(", ",@pair).")";
+ my @pair_sequences=();
+
+ foreach my $tag (@pair){
+# print "$tag: $node_alignments{$tag}\n" if $printer == 1;
+ print $node_alignments{$tag},"\n" if $printer == 1;
+ push @pair_sequences, $node_alignments{$tag};
+ }
+# print "ppeel onion joint = $joint , pair_sequences=>@pair_sequences< , pair=>@pair<\n" if $printer == 1;
+
+ my ($compared, $substitutions_list) = base_by_base_simple($motif,\@pair_sequences, scalar(@pair_sequences), @pair, $joint);
+ $node_alignments{$joint}=$compared;
+ push( @mutations_array,split(/:/,$substitutions_list));
+# print "newly added to node_sequences: $node_alignments{$joint} and list of mutations =\n", join("\n",@mutations_array),"\n" if $printer == 1;
+ }
+
+# print "now sending for analyze_mutations: mutation_array=@mutations_array, nodes=@nodes, branches_hash=$branches_hash, alignment=$alignment, tags=@tags, alivehash=\%alivehash, node_sequences=\%node_sequences, microsatstarts=$microsatstarts, motif=$motif\n" if $printer == 1;
+ ## if $printer == 1;
+
+ my $analayzed_mutations = analyze_mutations(\@mutations_array, \@nodes, $branches_hash, $alignment, \@tags, \%alivehash, \%node_sequences, $microsatstarts, $motif);
+
+# print "returning: ", $analayzed_mutations, \@nodes, $branches_hash,"\n" if scalar @mutations_array > 0 && $printer == 1;
+# print "returning: NULL, NULL, NULL " if scalar @mutations_array == 0 && $printer == 1;
+# print "final node alignment = $node_alignments{$joint}\n" if $printer == 1;
+
+ # if $printer == 1;
+
+
+ return ($analayzed_mutations, \@nodes, $branches_hash, \%alivehash, $node_alignments{$joint}) if scalar @mutations_array > 0;
+ return ("NULL",\@nodes,$branches_hash, \%alivehash, "NULL") if scalar @mutations_array == 0;
+}
+
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+
+sub get_nodes{
+ my $printer = 0;
+
+ my $tree=$_[0];
+ #$tree =~ s/ +//g;
+ $tree =~ s/\t+//g;
+ $tree=~s/;//g;
+ print "tree=$tree\n" if $printer == 1;
+ my @nodes = ();
+ my @onions=($tree);
+ my %branches=();
+ foreach my $bite (@onions){
+ $bite=~ s/^\(|\)$//g;
+ chomp $bite;
+# print "tree = $bite \n";
+# ;
+ $bite=~ /([ ,\(\)A-Z]+)\,\s*([ ,\(\)A-Z]+)/;
+ #$tree =~ /(\(\(\(H, C\), O\), R\))\, (M)/;
+ my @raw_nodes = ($1, $2);
+ print "raw nodes = $1 and $2\n" if $printer == 1;
+ push(@nodes, [@raw_nodes]);
+ foreach my $node (@raw_nodes) {push (@onions, $node) if $node =~ /,/;}
+ foreach my $node (@raw_nodes) {$branches{$node}="(".$bite.")"; print "adding to branches: $node = ($bite)\n" if $printer == 1;}
+ print "onions = @onions\n" if $printer == 1; if $printer == 1;
+ }
+ $printer = 0;
+ return \@nodes, \%branches;
+}
+
+
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+sub analyze_mutations{
+ my ($mutations_array, $nodes, $branches_hash, $alignment, $tags, $alivehash, $node_sequences, $microsatstarts, $motif) = @_;
+ my $locuslength = length($alignment->{$tags->[0]});
+ my $printer = 0;
+
+
+# print " IN analyzed_mutations....\n" if $printer == 1; # \n mutations array = @$mutations_array, \nAND locuslength = $locuslength\n" if $printer == 1;
+ my %mutation_hash=();
+ my %froms_megahash=();
+ my %tos_megahash=();
+ my %position_hash=();
+ my @solutions_array=();
+ foreach my $mutation (@$mutations_array){
+# print "loadin mutation: $mutation\n" if $printer == 1;
+ my %localhash= $mutation =~ /([\S ]+)=([\S ]+)/g;
+ $mutation_hash{$localhash{"position"}} = {%localhash};
+ push @{$position_hash{$localhash{"position"}}},$localhash{"node"};
+# print "feeding position hash with $localhash{position}: $position_hash{$localhash{position}}[0]\n" if $printer == 1;
+ $froms_megahash{$localhash{"position"}}{$localhash{"node"}}=$localhash{"from"};
+ $tos_megahash{$localhash{"position"}}{$localhash{"node"}}=$localhash{"to"};
+# print "just a trial: $mutation_hash{$localhash{position}}{position}\n" if $printer == 1;
+# print "loadin in tos_megahash: $localhash{position} {$localhash{node} = $localhash{to}\n" if $printer == 1;
+# print "loadin in from: $localhash{position} {$localhash{node} = $localhash{from}\n" if $printer == 1;
+ }
+
+# print "now going through each position in loculength:\n" if $printer == 1;
+ ## if $printer == 1;
+
+ for my $pos (0 ... $locuslength-1){
+# print "at position: $pos\n" if $printer == 1;
+
+ if (exists($mutation_hash{$pos})){
+ my @local_nodes=@{$position_hash{$pos}};
+# print "found mutation: @{$position_hash{$pos}} : @local_nodes\n" if $printer == 1;
+
+ foreach my $local_node (@local_nodes){
+# print "at local node: $local_node ... from state = $froms_megahash{$pos}{$local_node}\n" if $printer == 1;
+ my $open_insertion=();
+ my $open_deletion=();
+ my $open_to_substitution=();
+ my $open_from_substitution=();
+ if ($froms_megahash{$pos}{$local_node} eq "-"){
+ # print "here exists a microsatellite from $local_node to $branches_hash->{$local_node}\n" if $printer == 1 && exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};;
+ # print "for localnode $local_node, amd the realated branches_hash:$branches_hash->{$local_node}, nexting as exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}}\n" if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}} && $printer == 1;
+ #next if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
+ $open_insertion=$tos_megahash{$pos}{$local_node};
+ for my $posnext ($pos+1 ... $locuslength-1){
+# print "in first if .... studying posnext: $posnext\n" if $printer == 1;
+ last if !exists ($froms_megahash{$posnext}{$local_node});
+# print "for posnext: $posnext, there exists $froms_megahash{$posnext}{$local_node}.. already, open_insertion = $open_insertion.. checking is $froms_megahash{$posnext}{$local_node} matters\n" if $printer == 1;
+ $open_insertion = $open_insertion.$tos_megahash{$posnext}{$local_node} if $froms_megahash{$posnext}{$local_node} eq "-";
+# print "now open_insertion=$open_insertion\n" if $printer == 1;
+ delete $mutation_hash{$posnext} if $froms_megahash{$posnext}{$local_node} eq "-";
+ }
+ print "1 Feeding in: ", join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_insertion", "deletion="),"\n" if $printer == 1;
+ push (@solutions_array, join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_insertion", "deletion="));
+ }
+ elsif ($tos_megahash{$pos}{$local_node} eq "-"){
+ # print "here exists a microsatellite to $local_node from $branches_hash->{$local_node}\n" if $printer == 1 && exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};;
+ # print "for localnode $local_node, amd the realated branches_hash:$branches_hash->{$local_node}, nexting as exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}}\n" if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
+ #next if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
+ $open_deletion=$froms_megahash{$pos}{$local_node};
+ for my $posnext ($pos+1 ... $locuslength-1){
+ print "in 1st elsif studying posnext: $posnext\n" if $printer == 1;
+ print "nexting as nextpos does not exist\n" if !exists ($tos_megahash{$posnext}{$local_node}) && $printer == 1;
+ last if !exists ($tos_megahash{$posnext}{$local_node});
+ print "for posnext: $posnext, there exists $tos_megahash{$posnext}{$local_node}\n" if $printer == 1;
+ $open_deletion = $open_deletion.$froms_megahash{$posnext}{$local_node} if $tos_megahash{$posnext}{$local_node} eq "-";
+ delete $mutation_hash{$posnext} if $tos_megahash{$posnext}{$local_node} eq "-";
+ }
+ print "2 Feeding in:", join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_deletion"), "\n" if $printer == 1;
+ push (@solutions_array, join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_deletion"));
+ }
+ elsif ($tos_megahash{$pos}{$local_node} ne "-"){
+ # print "here exists a microsatellite from $local_node to $branches_hash->{$local_node}\n" if $printer == 1 && exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};;
+ # print "for localnode $local_node, amd the realated branches_hash:$branches_hash->{$local_node}, nexting as exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}}\n" if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
+ #next if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
+ # print "microsatstart = $microsatstarts->{$local_node} \n" if exists $microsatstarts->{$local_node} && $pos < $microsatstarts->{$local_node} && $printer == 1;
+ next if exists $microsatstarts->{$local_node} && $pos < $microsatstarts->{$local_node};
+ $open_to_substitution=$tos_megahash{$pos}{$local_node};
+ $open_from_substitution=$froms_megahash{$pos}{$local_node};
+ print "open from substitution: $open_from_substitution \n" if $printer == 1;
+ for my $posnext ($pos+1 ... $locuslength-1){
+ #print "in last elsif studying posnext: $posnext\n";
+ last if !exists ($tos_megahash{$posnext}{$local_node});
+ print "for posnext: $posnext, there exists $tos_megahash{$posnext}{$local_node}\n" if $printer == 1;
+ $open_to_substitution = $open_to_substitution.$tos_megahash{$posnext}{$local_node} if $tos_megahash{$posnext}{$local_node} ne "-";
+ $open_from_substitution = $open_from_substitution.$froms_megahash{$posnext}{$local_node} if $tos_megahash{$posnext}{$local_node} ne "-";
+ delete $mutation_hash{$posnext} if $tos_megahash{$posnext}{$local_node} ne "-" && $froms_megahash{$posnext}{$local_node} ;
+ }
+ print "open from substitution: $open_from_substitution \n" if $printer == 1;
+
+ #IS THE STRETCH OF SUBSTITUTION MICROSATELLITE-LIKE?
+ my @motif_parts=split(/\s*/,$motif);
+ #GENERATING THE FLEXIBLE LEFT END
+ my $left_query=();
+ for my $k (1 ... $#motif_parts) {
+ $left_query= $motif_parts[$k]."|)";
+ $left_query="(".$left_query;
+ }
+ $left_query=$left_query."?";
+
+ print "left_quewry = $left_query\n" if $printer == 1;
+ #GENERATING THE FLEXIBLE RIGHT END
+ my $right_query=();
+ for my $k (0 ... ($#motif_parts-1)) {
+ $right_query= "(|".$motif_parts[$k];
+ $right_query=$right_query.")";
+ }
+ $right_query=$right_query."?";
+ print "right_query = $right_query\n" if $printer == 1;
+ print "Hence, searching for: ^$left_query($motif)+$right_query\$\n" if $printer == 1;
+
+ my $motifcomb=$motif x 50;
+ print "motifcomb = $motifcomb\n" if $printer == 1;
+ if ( ($motifcomb =~/$open_to_substitution/i) && (length ($open_to_substitution) >= length($motif)) ){
+ print "sequence microsat-like\n" if $printer == 1;
+ my $all_microsat_like = 0;
+ print "3 feeding in: ", join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_from_substitution"), "\n" if $printer == 1;
+ push (@solutions_array, join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_from_substitution"));
+ print "4 feeding in: ", join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_to_substitution", "deletion="), "\n" if $printer == 1;
+ push (@solutions_array, join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_to_substitution", "deletion="));
+
+ }
+ else{
+ print "5 feeding in: ", join("\t", "node=$local_node","type=substitution" ,"position=$pos", "from=$open_from_substitution", "to=$open_to_substitution", "insertion=", "deletion="), "\n" if $printer == 1;
+ push (@solutions_array, join("\t", "node=$local_node","type=substitution" ,"position=$pos", "from=$open_from_substitution", "to=$open_to_substitution", "insertion=", "deletion="));
+ }
+ #IS THE FROM-SEQUENCE MICROSATELLITE-LIKE?
+
+ }
+ # if $printer ==1;
+ }
+ # if $printer ==1;
+ }
+ }
+
+ print "\n", "#" x 50, "\n" if $printer == 1;
+ foreach my $tag (@$tags){
+ print "$tag: $alignment->{$tag}\n" if $printer == 1;
+ }
+ print "\n", "#" x 50, "\n" if $printer == 1;
+
+ print "returning SOLUTIONS ARRAY : \n",join("\n", @solutions_array),"\n" if $printer == 1;
+ #print "end\n";
+ # if
+ return \@solutions_array;
+}
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
+
+sub base_by_base_simple{
+ my $printer = 0;
+ my ($motif, $locus, $no, $pair0, $pair1, $joint) = @_;
+ my @seq_array=();
+ print "IN SUBROUTUNE base_by_base_simple.. information received = @_\n" if $printer == 1;
+ print "pair0 = $pair0 and pair1 = $pair1\n" if $printer == 1;
+
+ my @example=split(/\./,$locus->[0]);
+ print "example, for length = @example\n" if $printer == 1;
+ for my $i (0...$no-1){push(@seq_array, [split(/\./,$locus->[$i])]); print "for $i, from $locus->[$i], seq_array = >@{$seq_array[$i]}<\n" if $printer == 1;}
+
+ my @compared_sequence=();
+ my @substitutions_list;
+ for my $i (0...scalar(@example)-1){
+
+ #print "i = $i\n" if $printer == 1;
+ #print "comparing $seq_array[0][$i] and $seq_array[1][$i] \n" ;#if $printer == 1;
+ if ($seq_array[0][$i] =~ /!/ && $seq_array[1][$i] !~ /!/){
+
+ my $resolution= resolve_base($seq_array[0][$i],$seq_array[1][$i], $pair1 ,"keep" );
+ # print "ancestral = $resolution\n" if $printer == 1;
+
+ if ($resolution =~ /$seq_array[1][$i]/i && $resolution !~ /!/){
+ push @substitutions_list, add_mutation($i, $pair0, $seq_array[0][$i], $resolution );
+ }
+ elsif ( $resolution !~ /!/){
+ push @substitutions_list, add_mutation($i, $pair1, $seq_array[1][$i], $resolution);
+ }
+ push @compared_sequence,$resolution;
+ }
+ elsif ($seq_array[0][$i] !~ /!/ && $seq_array[1][$i] =~ /!/){
+
+ my $resolution= resolve_base($seq_array[1][$i],$seq_array[0][$i], $pair0, "invert" );
+ # print "ancestral = $resolution\n" if $printer == 1;
+
+ if ($resolution =~ /$seq_array[0][$i]/i && $resolution !~ /!/){
+ push @substitutions_list, add_mutation($i, $pair1, $seq_array[1][$i], $resolution);
+ }
+ elsif ( $resolution !~ /!/){
+ push @substitutions_list, add_mutation($i, $pair0, $seq_array[0][$i], $resolution);
+ }
+ push @compared_sequence,$resolution;
+ }
+ elsif($seq_array[0][$i] =~ /!/ && $seq_array[1][$i] =~ /!/){
+ push @compared_sequence, add_bases($seq_array[0][$i],$seq_array[1][$i], $pair0, $pair1, $joint );
+ }
+ else{
+ if($seq_array[0][$i] !~ /^$seq_array[1][$i]$/i){
+ push @compared_sequence, $pair0.":".$seq_array[0][$i]."!".$pair1.":".$seq_array[1][$i];
+ }
+ else{
+ # print "perfect match\n" if $printer == 1;
+ push @compared_sequence, $seq_array[0][$i];
+ }
+ }
+ }
+ print "returning: comared = @compared_sequence \nand substitutions list =\n", join("\n",@substitutions_list),"\n" if $printer == 1;
+ return join(".",@compared_sequence), join(":", @substitutions_list) if scalar (@substitutions_list) > 0;
+ return join(".",@compared_sequence), "" if scalar (@substitutions_list) == 0;
+}
+
+
+sub resolve_base{
+ my $printer = 0;
+
+ print "IN SUBROUTUNE resolve_base.. information received = @_\n" if $printer == 1;
+ my ($optional, $single, $singlesp, $arg) = @_;
+ my @options=split(/!/,$optional);
+ foreach my $option(@options) {
+ $option=~s/[A-Z\(\) ,]+://g;
+ if ($option =~ /$single/i){
+ print "option = $option , returning single: $single\n" if $printer == 1;
+ return $single;
+ }
+ }
+ print "returning ",$optional."!".$singlesp.":".$single. "\n" if $arg eq "keep" && $printer == 1;
+ print "returning ",$singlesp.":".$single."!".$optional. "\n" if $arg eq "invert" && $printer == 1;
+ return $optional."!".$singlesp.":".$single if $arg eq "keep";
+ return $singlesp.":".$single."!".$optional if $arg eq "invert";
+
+}
+
+sub same_length{
+ my $printer = 0;
+ my @locus = @_;
+ my $temp = shift @locus;
+ $temp=~s/-|,//g;
+ foreach my $l (@locus){
+ $l=~s/-|,//g;
+ return 0 if length($l) != length($temp);
+ $temp = $l;
+ }
+ return 1;
+}
+sub treeStudy{
+ my $printer = 0;
+# print "template DEFINED.. received: @_\n" if defined %template;
+# print "only received = @_" if !defined %template;
+ my $stopper = 0;
+ if (!defined %template){
+ $stopper = 1;
+ %template=();
+ print "tree decipherer = $tree_decipherer\n" if $printer == 1;
+ my ( $template_ref, $keys_array)=load_allPossibleTrees($tree_decipherer, \%template);
+ print "return = $template_ref and @{$keys_array}\n" if $printer == 1;
+ foreach my $key (@$keys_array){
+ print "addding : $template_ref->{$key} for $key\n" if $printer == 1;
+ $template{$key} = $template_ref->{$key};
+ }
+ }
+
+ for my $templet ( keys %template ) {
+ # print "$templet => @{$template{$templet}}\n";
+ }
+ if !defined %template;
+
+ my $strict = 0;
+
+ my $H = 0;
+ my $Hchr = 1;
+ my $Hstart = 2;
+ my $Hend = 3;
+ my $Hmotif = 4;
+ my $Hmotiflen = 5;
+ my $Hmicro = 6;
+ my $Hstrand = 7;
+ my $Hmicrolen = 8;
+ my $Hinterpos = 9;
+ my $Hrelativepos = 10;
+ my $Hinter = 11;
+ my $Hinterlen = 12;
+
+ my $C = 13;
+ my $Cchr = 14;
+ my $Cstart = 15;
+ my $Cend = 16;
+ my $Cmotif = 17;
+ my $Cmotiflen = 18;
+ my $Cmicro = 19;
+ my $Cstrand = 20;
+ my $Cmicrolen = 21;
+ my $Cinterpos = 22;
+ my $Crelativepos = 23;
+ my $Cinter = 24;
+ my $Cinterlen = 25;
+
+ my $O = 26;
+ my $Ochr = 27;
+ my $Ostart = 28;
+ my $Oend = 29;
+ my $Omotif = 30;
+ my $Omotiflen = 31;
+ my $Omicro = 32;
+ my $Ostrand = 33;
+ my $Omicrolen = 34;
+ my $Ointerpos = 35;
+ my $Orelativepos = 36;
+ my $Ointer = 37;
+ my $Ointerlen = 38;
+
+ my $R = 39;
+ my $Rchr = 40;
+ my $Rstart = 41;
+ my $Rend = 42;
+ my $Rmotif = 43;
+ my $Rmotiflen = 44;
+ my $Rmicro = 45;
+ my $Rstrand = 46;
+ my $Rmicrolen = 47;
+ my $Rinterpos = 48;
+ my $Rrelativepos = 49;
+ my $Rinter = 50;
+ my $Rinterlen = 51;
+
+ my $Mchr = 52;
+ my $Mstart = 53;
+ my $Mend = 54;
+ my $M = 55;
+ my $Mmotif = 56;
+ my $Mmotiflen = 57;
+ my $Mmicro = 58;
+ my $Mstrand = 59;
+ my $Mmicrolen = 60;
+ my $Minterpos = 61;
+ my $Mrelativepos = 62;
+ my $Minter = 63;
+ my $Minterlen = 64;
+
+ #-------------------------------------------------------------------------------#
+ my @analysis=();
+
+
+ my %speciesOrder = ();
+ $speciesOrder{"H"} = 0;
+ $speciesOrder{"C"} = 1;
+ $speciesOrder{"O"} = 2;
+ $speciesOrder{"R"} = 3;
+ $speciesOrder{"M"} = 4;
+ #-------------------------------------------------------------------------------#
+
+ my $line = $_[0];
+ chomp $line;
+
+ my @f = split(/\t/,$line);
+ print "received array : @f.. recieved tags = @tags\n" if $printer == 1;
+
+ # collect all motifs
+ my @motifs=();
+ @motifs = ($f[$Hmotif], $f[$Cmotif], $f[$Omotif], $f[$Rmotif], $f[$Mmotif]) if $tags[$#tags] =~ /M/;
+ @motifs = ($f[$Hmotif], $f[$Cmotif], $f[$Omotif], $f[$Rmotif]) if $tags[$#tags] =~ /R/;
+ @motifs = ($f[$Hmotif], $f[$Cmotif], $f[$Omotif]) if $tags[$#tags] =~ /O/;
+# print "motifs in the array = $f[$Hmotif], $f[$Cmotif], $f[$Omotif], $f[$Rmotif]\n" if $tags[$#tags] =~ /R/;;
+ print "motifs = @motifs\n" if $printer == 1;
+ my @translation = ();
+ foreach my $motif (@motifs){
+ push(@translation, "_") if $motif eq "NA";
+ push(@translation, "+") if $motif ne "NA";
+ }
+ my $translate = join(" ", @translation);
+# print "translate = >$translate< and analysis = $template{$translate}[0].. on the other hand, ",$template{"- - +"}[0],"\n";
+ my @analyses = split(/\|/,$template{$translate}[0]);
+
+ print "motifs = @motifs, analyses = @analyses\n" if $printer == 1;
+
+ if (scalar(@analyses) == 1) {
+ #print "analysis = $analyses[0]\n";
+ if ($analyses[0] !~ /,|\./ ){
+ if ($analyses[0] =~ /\+/){
+ my $analysis = $analyses[0];
+ $analysis =~ s/\+|\-//g;
+ my @species = split(/\s*/,$analysis);
+ my @currentMotifs = ();
+ foreach my $specie (@species){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); print "pushing into currentMotifs: $speciesOrder{$specie}: $motifs[$speciesOrder{$specie}]\n" if $printer == 1;}
+ print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
+ $template{$translate}[1]++ if $strict == 1 && consistency(@currentMotifs) ne "NULL";
+ $template{$translate}[1]++ if $strict == 0;
+ print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
+ }
+ else{
+ my $analysis = $analyses[0];
+ $analysis =~ s/\+|\-//g;
+ my @species = split(/\s*/,$analysis);
+ my @currentMotifs = ();
+ my @complementarySpecies = ();
+ my $allSpecies = join("",@tags);
+ foreach my $specie (@species){ $allSpecies =~ s/$specie//g; }
+ foreach my $specie (split(/\s*/,$allSpecies)){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); print "pushing into currentMotifs: $speciesOrder{$specie}: $motifs[$speciesOrder{$specie}]\n" if $printer == 1;;}
+ print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
+ print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
+ }
+ }
+
+ elsif ($analyses[0] =~ /,/) {
+ my @events = split(/,/,$analyses[0]);
+ print "events = @events \n " if $printer == 1;
+ if ($events[0] =~ /\+/){
+ my $analysis1 = $events[0];
+ $analysis1 =~ s/\+|\-//g;
+ my $analysis2 = $events[1];
+ $analysis2 =~ s/\+|\-//g;
+ my @nSpecies = split(/\s*/,$analysis2);
+ print "original anslysis = $analysis1 " if $printer == 1;
+ foreach my $specie (@nSpecies){ $analysis1=~ s/$specie//g;}
+ print "processed anslysis = $analysis1 \n" if $printer == 1;
+ my @currentMotifs = ();
+ foreach my $specie (split(/\s*/,$analysis1)){push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
+ print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
+ print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
+ }
+ else{
+ my $analysis1 = $events[0];
+ $analysis1 =~ s/\+|\-//g;
+ my $analysis2 = $events[1];
+ $analysis2 =~ s/\+|\-//g;
+ my @pSpecies = split(/\s*/,$analysis2);
+ my @currentMotifs = ();
+ foreach my $specie (@pSpecies){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
+ print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
+ print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
+
+ }
+
+ }
+ elsif ($analyses[0] =~ /\./) {
+ my @events = split(/\./,$analyses[0]);
+ foreach my $event (@events){
+ print "event = $event \n" if $printer == 1;
+ if ($event =~ /\+/){
+ my $analysis = $event;
+ $analysis =~ s/\+|\-//g;
+ my @species = split(/\s*/,$analysis);
+ my @currentMotifs = ();
+ foreach my $specie (@species){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
+ #print consistency(@currentMotifs),"<- \n";
+ print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
+ print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
+ }
+ else{
+ my $analysis = $event;
+ $analysis =~ s/\+|\-//g;
+ my @species = split(/\s*/,$analysis);
+ my @currentMotifs = ();
+ my @complementarySpecies = ();
+ my $allSpecies = join("",@tags);
+ foreach my $specie (@species){ $allSpecies =~ s/$specie//g; }
+ foreach my $specie (split(/\s*/,$allSpecies)){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
+ #print consistency(@currentMotifs),"<- \n";
+ print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
+ $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
+ print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
+ }
+ }
+
+ }
+ }
+ else{
+ my $finalanalysis = ();
+ $template{$translate}[1]++;
+ foreach my $analysis (@analyses){ ;}
+ }
+ # test if motifs where microsats are present, as indeed of same the motif composition
+
+
+
+ for my $templet ( keys %template ) {
+ # print "now returning: $templet\n";
+ if (@{ $template{$templet} }[1] > 0){
+ print "returning in the end: $templet and $translate\n" if $printer == 1;
+ $template{$templet}[1] = 0;
+ return (@{$template{$templet}}[0], $translate);
+ }
+ }
+ undef %template;
+ print "sending NULL\n" if $printer == 1;
+ return ("NULL", $translate);
+
+}
+
+
+sub consistency{
+ my @motifs = @_;
+ print "in consistency \n" if $printer == 1;
+ print "motifs sent = >",join("|",@motifs),"< \n" if $printer == 1;
+ return $motifs[0] if scalar(@motifs) == 1;
+ my $prevmotif = shift(@motifs);
+ my $stopper = 0;
+ for my $i (0 ... $#motifs){
+ next if $motifs[$i] eq "NA";
+ my $templet = $motifs[$i].$motifs[$i];
+ if ($templet !~ /$prevmotif/i){
+ $stopper = 1; last;
+ }
+ }
+ return $prevmotif if $stopper == 0;
+ return "NULL" if $stopper == 1;
+}
+sub summarize_microsat{
+ my $printer = 0;
+ my $line = $_[0];
+ my $humseq = $_[1];
+
+ my @gaps = $line =~ /[0-9]+\t[0-9]+\t[\+\-]/g;
+ my @starts = $line =~ /[0-9]+\t[\+\-]/g;
+ my @ends = $line =~ /[\+\-]\t[0-9]+/g;
+ print "starts = @starts\tends = @ends\n" if $printer == 1;
+ for my $i (0 ... $#gaps) {$gaps[$i] =~ s/\t[0-9]+\t[\+\-]//g;}
+ for my $i (0 ... $#starts) {$starts[$i] =~ s/\t[\+\-]//g;}
+ for my $i (0 ... $#ends) {$ends[$i] =~ s/[\+\-]\t//g;}
+
+ my $minstart = array_smallest_number(@starts);
+ my $maxend = array_largest_number(@ends);
+
+ my $humupstream_st = substr($humseq, 0, $minstart);
+ my $humupstream_en = substr($humseq, 0, $maxend);
+ my $no_of_gaps_to_start = 0;
+ my $no_of_gaps_to_end = 0;
+ $no_of_gaps_to_start = ($humupstream_st =~ s/\-/x/g) if $humupstream_st=~/\-/;
+ $no_of_gaps_to_end = ($humupstream_en =~ s/\-/x/g) if $humupstream_en=~/\-/;
+
+ my $locusmotif = ();
+
+
+ print "IN SUB SUMMARIZE_MICROSAT $line\n" if $printer == 1;
+ #return "NULL" if $line =~ /compound/;
+ my $Hstart = "NA";
+ my $Hend = "NA";
+ chomp $line;
+ my $match_count = ($line =~ s/>/>/g);
+ #print "number of species = $match_count\n";
+ my @micros = split(/>/,$line);
+ shift @micros;
+ my $stopper = 0;
+
+
+ foreach my $mic (@micros){
+ my @local = split(/\t/,$mic);
+ if ($local[$microsatcord] =~ /N/) {$stopper =1; last;}
+ }
+ return "NULL" if $stopper ==1;
+
+ #------------------------------------------------------
+
+ my @arranged = ();
+ for my $arr (0 ... $#exacttags) {$arranged[$arr] = '0';}
+
+ foreach my $micro (@micros){
+ for my $i (0 ... $#exacttags){
+ if ($micro =~ /^$exacttags[$i]/){
+ $arranged[$i] = $micro;
+ last;
+ }
+ }
+ }
+# print "arranged = @arranged \n" ; ;;
+
+ my @endstatement = ();
+ my $turn = 0;
+ my $species_counter = 0;
+ # print scalar(@arranged),"\n";
+
+ my $species_no=0;
+
+ my $orthHchr = 0;
+
+ foreach my $micro (@arranged) {
+ $micro =~ s/\t\t/\t \t/g;
+ $micro =~ s/\t,/\t ,/g;
+ $micro =~ s/,\t/, \t/g;
+ print "------------------------------------------------------------------------------------------\n" if $printer == 1;
+ chomp $micro;
+ if ($micro eq '0'){
+ push(@endstatement, join("\t",$exacttags[$species_counter],"NA","NA","NA","NA",0 ,"NA", "NA", 0,"NA","NA","NA", "NA" ));
+ $species_counter++;
+ print join("|","ENDSTATEMENT:",@endstatement),"\n" if $printer == 1;
+ next;
+ }
+ # print $micro,"\n";
+ print "micro = $micro \n" if $printer == 1;
+ my @fields = split(/\t/,$micro);
+ my $microcopy = $fields[$microsatcord];
+ $microcopy =~ s/\[|\]|-//g;
+ my $microsatlength = length($microcopy);
+ print "microsat = $fields[$microsatcord] and microsatlength = $microsatlength\n" if $printer == 1;
+# print "sp_ident = @sp_ident.. species_no=$species_no\n";
+ $micro =~ /$sp_ident[$species_no]\s(\S+)\s([0-9]+)\s([0-9]+)/;
+
+
+ my $sp_chr=$1;
+ my $sp_start=$2 + $fields[$startcord] - $fields[$gapcord];
+ my $sp_end= $sp_start + $microsatlength - 1;
+
+ $species_no++;
+
+ $micro =~ /$focalspec\s(\S+)\s([0-9]+)\s([0-9]+)/;
+ $orthHchr=$1;
+ $Hstart=$2+$minstart-$no_of_gaps_to_start;
+ $Hend=$2+$maxend-$no_of_gaps_to_end;
+
+ print "Hstart = $Hstart = $fields[4] + $fields[$startcord] - $fields[$gapcord]\n" if $printer == 1;
+
+ my $motif = $fields[$motifcord];
+ my $firstmotif = ();
+ my $strand = $fields[$strandcord];
+ # print "strand = $strand\n";
+
+
+ if ($motif =~ /^\[/){
+ $motif =~ s/^\[//g;
+ $motif =~ /([a-zA-Z]+)\].*/;
+ $firstmotif = $1;
+ }
+
+ else {$firstmotif = $motif;}
+ print "firstmotif =$firstmotif : \n" if $printer == 1;
+ $firstmotif = allCaps($firstmotif);
+
+ if (exists $revHash{$firstmotif} && $turn == 0) {
+ $turn=1 if $species_counter==0;
+ $firstmotif = $revHash{$firstmotif};
+ }
+
+ elsif (exists $revHash{$firstmotif} && $turn == 1) {$firstmotif = $revHash{$firstmotif}; $turn = 1;}
+ print "changed firstmotif =$firstmotif\n" if $printer == 1;
+ # ;
+ $locusmotif = $firstmotif;
+
+ if (scalar(@fields) > $microsatcord + 2){
+ print "fields = @fields ... interr_poscord=$interr_poscord=$fields[$interr_poscord] .. interrcord=$interrcord=$fields[$interrcord]\n" if $printer == 1;
+
+ my @interposes = ();
+ @interposes = split(",",$fields[$interr_poscord]) if $fields[$interr_poscord] =~ /,/;
+ $interposes[0] = $fields[$interr_poscord] if $fields[$interr_poscord] !~ /,/ ;
+ print "interposes=@interposes\n" if $printer == 1;
+ my @relativeposes = ();
+ my @interruptions = ();
+ @interruptions = split(",",$fields[$interrcord]) if $fields[$interrcord] =~ /,/;
+ $interruptions[0] = $fields[$interrcord] if $fields[$interrcord] !~ /,/;
+ my @interlens = ();
+
+
+ for my $i (0 ... $#interposes){
+
+ my $interpos = $interposes[$i];
+ my $nexter = 0;
+ my $interruption = $interruptions[$i];
+ my $interlen = length($interruption);
+ push (@interlens, $interlen);
+
+
+ my $relativepos = (100 * $interpos) / $microsatlength;
+ print "relativepos = $relativepos ,interpos=$interpos, interruption=$interruption, interlen=$interlen \n" if $printer == 1;
+ $relativepos = (100 * ($interpos-$interlen)) / $microsatlength if $relativepos > 50;
+ print "--> = $relativepos\n" if $printer == 1;
+ $interruption = "IND" if length($interruption) < 1;
+
+ if ($turn == 1){
+ $fields[$microsatcord] = switch_micro($fields[$microsatcord]);
+ $interruption = switch_nucl($interruption) unless $interruption eq "IND";
+ $interpos = ($microsatlength - $interpos) - $interlen + 2;
+ print "turn interpos = $interpos for $fields[$microsatcord]\n" if $printer == 1;
+ $relativepos = (100 * $interpos) / $microsatlength;
+ $relativepos = (100 * ($interpos-$interlen)) / $microsatlength if $relativepos > 50;
+
+
+ $strand = '+' if $strand eq '-';
+ $strand = '-' if $strand eq '+';
+ }
+ print "final relativepos = $relativepos\n" if $printer == 1;
+ push(@relativeposes, $relativepos);
+ }
+ push(@endstatement,join("\t",($exacttags[$species_counter],$sp_chr, $sp_start, $sp_end, $firstmotif,length($firstmotif),$fields[$microsatcord],$strand,$microsatlength,join(",",@interposes),join(",",@relativeposes),join(",",@interruptions), join(",",@interlens))));
+ }
+
+ else{
+ push(@endstatement, join("\t",$exacttags[$species_counter],$sp_chr, $sp_start, $sp_end, $firstmotif,length($firstmotif),$fields[$microsatcord],$strand,$microsatlength,"NA","NA","NA", "NA"));
+ }
+
+ $species_counter++;
+ }
+
+ $locusmotif = $sameHash{$locusmotif} if exists $sameHash{$locusmotif};
+ $locusmotif = $revHash{$locusmotif} if exists $revHash{$locusmotif};
+
+ my $endst = join("\t", @endstatement, $orthHchr, $Hstart, $Hend);
+ print join("\t", @endstatement, $orthHchr, $Hstart, $Hend), "\n" if $printer == 1;
+
+
+ return (join("\t", @endstatement, $orthHchr, $Hstart, $Hend), $orthHchr, $Hstart, $Hend, $locusmotif, length($locusmotif));
+
+}
+
+sub switch_nucl{
+ my @strand = split(/\s*/,$_[0]);
+ for my $i (0 ... $#strand){
+ if ($strand[$i] =~ /c/i) {$strand[$i] = "G";next;}
+ if ($strand[$i] =~ /a/i) {$strand[$i] = "T";next;}
+ if ($strand[$i] =~ /t/i) { $strand[$i] = "A";next;}
+ if ($strand[$i] =~ /g/i) {$strand[$i] = "C";next;}
+ }
+ return join("",@strand);
+}
+
+
+sub switch_micro{
+ my $micro = reverse($_[0]);
+ my @strand = split(/\s*/,$micro);
+ for my $i (0 ... $#strand){
+ if ($strand[$i] =~ /c/i) {$strand[$i] = "G";next;}
+ if ($strand[$i] =~ /a/i) {$strand[$i] = "T";next;}
+ if ($strand[$i] =~ /t/i) { $strand[$i] = "A";next;}
+ if ($strand[$i] =~ /g/i) {$strand[$i] = "C";next;}
+ if ($strand[$i] =~ /\[/i) {$strand[$i] = "]";next;}
+ if ($strand[$i] =~ /\]/i) {$strand[$i] = "[";next;}
+ }
+ return join("",@strand);
+}
+sub decipher_history{
+ my $printer = 0;
+ my ($mutations_array, $tags_string, $nodes, $branches_hash, $tree_analysis, $confirmation_string, $alivehash) = @_;
+ my %mutations_hash=();
+ foreach my $mutation (@$mutations_array){
+ print "mutation = $mutation\n" if $printer == 1;
+ my %local = $mutation =~ /([\S ]+)=([\S ]+)/g;
+ push @{$mutations_hash{$local{"node"}}},$mutation;
+ print "just for confirmation: $local{node} pushed as: $mutation\n" if $printer == 1;
+ }
+ my @nodes;
+ my @birth_steps=();
+ my @death_steps=();
+
+ my @tags=split(/\s*/,$tags_string);
+ my @confirmation=split(/\s+/,$confirmation_string);
+ my %info=();
+
+ for my $i (0 ... $#tags){
+ $info{$tags[$i]}=$confirmation[$i];
+ print "feeding info: $tags[$i] = $info{$tags[$i]}\n" if $printer == 1;
+ }
+
+ for my $keys (@$nodes) {
+ foreach my $key (@$keys){
+# print "current key = $key\n";
+ my $copykey = $key;
+ $copykey =~ s/[\W ]+//g;
+ my @copykeys=split(/\s*/,$copykey);
+ my $states=();
+ foreach my $copy (@copykeys){
+ $states=$states.$info{$copy};
+ }
+
+ print "reduced key = $copykey and state = $states\n" if $printer == 1;
+
+ if (exists $mutations_hash{$key}) {
+
+ if ($states=~/\+/){
+ push @birth_steps, @{$mutations_hash{$key}};
+ $birth_steps[$#birth_steps] =~ s/\S+=//g;
+ delete $mutations_hash{$key};
+ }
+ else{
+ push @death_steps, @{$mutations_hash{$key}};
+ $death_steps[$#death_steps] =~ s/\S+=//g;
+ delete $mutations_hash{$key};
+ }
+ }
+ }
+ }
+ print "conformation = $confirmation_string\n" if $printer == 1;
+ push (@birth_steps, "NULL") if scalar(@birth_steps) == 0;
+ push (@death_steps, "NULL") if scalar(@death_steps) == 0;
+ print "birth steps = ",join("\n",@birth_steps)," and death steps = ",join("\n",@death_steps),"\n" if $printer == 1;
+ return \@birth_steps, \@death_steps;
+}
+
+sub fillAlignmentGaps{
+ my $printer = 0;
+ print "received: @_\n" if $printer == 1;
+ my ($tree, $sequences, $alignment, $tagarray, $microsathash, $nonmicrosathash, $motif, $tree_analysis, $threshold, $microsatstarts) = @_;
+ print "in fillAlignmentGaps.. tree = $tree \n" if $printer == 1;
+ my %sequence_hash=();
+
+ my @phases = ();
+ my $concat = $motif.$motif;
+ my $motifsize = length($motif);
+
+ for my $i (1 ... $motifsize){
+ push @phases, substr($concat, $i, $motifsize);
+ }
+
+ my $concatalignment = ();
+ foreach my $tag (@tags){
+ $concatalignment = $concatalignment.$alignment->{$tag};
+ }
+# print "returningg NULL","NULL","NULL", "NULL\n" if $concatalignment !~ /-/;
+ return 0, "NULL","NULL","NULL", "NULL","NULL" if $concatalignment !~ /-/;
+
+
+
+ my %node_sequences_temp=();
+ my %node_alignments_temp =(); #NEW, Nov 28 2008
+
+ my @tags=();
+ my @locus_sequences=();
+ my %alivehash=();
+
+# print "IN fillAlignmentGaps\n";# ;
+ my %fillrecord = ();
+
+ my $change = 0;
+ foreach my $tag (@$tagarray) {
+ #print "adding: $tag\n";
+ push(@tags, $tag);
+ if (exists $microsathash->{$tag}){
+ my $micro = $microsathash->{$tag};
+ my $orig_micro = $micro;
+ ($micro, $fillrecord{$tag}) = fillgaps($micro, \@phases);
+ $change = 1 if uc($micro) ne uc($orig_micro);
+ $node_sequences_temp{$tag}=$micro if $microsathash->{$tag} ne "NULL";
+ }
+ if (exists $nonmicrosathash->{$tag}){
+ my $micro = $nonmicrosathash->{$tag};
+ my $orig_micro = $micro;
+ ($micro, $fillrecord{$tag}) = fillgaps($micro, \@phases);
+ $change = 1 if uc($micro) ne uc($orig_micro);
+ $node_sequences_temp{$tag}=$micro if $nonmicrosathash->{$tag} ne "NULL";
+ }
+
+ if (exists $alignment->{$tag}){
+ my $micro = $alignment->{$tag};
+ my $orig_micro = $micro;
+ ($micro, $fillrecord{$tag}) = fillgaps($micro, \@phases);
+ $change = 1 if uc($micro) ne uc($orig_micro);
+ $node_alignments_temp{$tag}=$micro if $alignment->{$tag} ne "NULL";
+ }
+
+ #print "adding to node_sequences: $tag = ",$node_sequences_temp{$tag},"\n" if $printer == 1;
+ #print "adding to node_alignments: $tag = ",$node_alignments_temp{$tag},"\n" if $printer == 1;
+ }
+
+
+ my %node_sequences=();
+ my %node_alignments =(); #NEW, Nov 28 2008
+ foreach my $tag (@$tagarray) {
+ $node_sequences{$tag} = join ".",split(/\s*/,$node_sequences_temp{$tag});
+ $node_alignments{$tag} = join ".",split(/\s*/,$node_alignments_temp{$tag});
+ }
+
+ print "\n", "#" x 50, "\n" if $printer == 1;
+ foreach my $tag (@tags){
+ print "$tag: $alignment->{$tag} = $node_alignments{$tag}\n" if $printer == 1;
+ }
+ print "\n", "#" x 50, "\n" if $printer == 1;
+# print "change = $change\n";
+ # if $concatalignment=~/\-/;
+
+# if $printer == 1 && $concatalignment =~ /\-/;
+
+ return 0, "NULL","NULL","NULL", "NULL", "NULL" if $change == 0;
+
+ my ($nodes_arr, $branches_hash) = get_nodes($tree);
+ my @nodes=@$nodes_arr;
+ print "recieved nodes = @nodes\n" if $printer == 1;
+
+
+ #POPULATE branches_hash WITH INFORMATION ABOUT LIVESTATUS
+ foreach my $keys (@nodes){
+ my @pair = @$keys;
+ my $joint = "(".join(", ",@pair).")";
+ my $copykey = join "", @pair;
+ $copykey =~ s/[\W ]+//g;
+ print "for node: $keys, copykey = $copykey and joint = $joint\n" if $printer == 1;
+ my $livestatus = 1;
+ foreach my $copy (split(/\s*/,$copykey)){
+ $livestatus = 0 if !exists $alivehash{$copy};
+ }
+ $alivehash{$joint} = $joint if !exists $alivehash{$joint} && $livestatus == 1;
+ print "alivehash = $alivehash{$joint}\n" if exists $alivehash{$joint} && $printer == 1;
+ }
+
+
+
+ @nodes = reverse(@nodes); #1 THIS IS IN ORDER TO GO THROUGH THE TREE FROM LEAVES TO ROOT.
+
+ my @mutations_array=();
+
+ my $joint = ();
+ foreach my $node (@nodes){
+ my @pair = @$node;
+ print "now in the nodes for loop, pair = @pair\n and sequences=\n" if $printer == 1;
+ $joint = "(".join(", ",@pair).")";
+ print "joint = $joint \n" if $printer == 1;
+ my @pair_sequences=();
+
+ foreach my $tag (@pair){
+ print "tag = $tag: " if $printer == 1;
+ print $node_alignments{$tag},"\n" if $printer == 1;
+ push @pair_sequences, $node_alignments{$tag};
+ }
+# print "fillgap\n";
+ my ($compared, $substitutions_list) = base_by_base_simple($motif,\@pair_sequences, scalar(@pair_sequences), @pair, $joint);
+ $node_alignments{$joint}=$compared;
+ push( @mutations_array,split(/:/,$substitutions_list));
+ print "newly added to node_sequences: $node_alignments{$joint} and list of mutations = @mutations_array\n" if $printer == 1;
+ }
+ print "now sending for analyze_mutations: mutation_array=@mutations_array, nodes=@nodes, branches_hash=$branches_hash, alignment=$alignment, tags=@tags, alivehash=%alivehash, node_sequences=\%node_sequences, microsatstarts=$microsatstarts, motif=$motif\n" if $printer == 1;
+# if $printer == 1;
+
+ my $analayzed_mutations = analyze_mutations(\@mutations_array, \@nodes, $branches_hash, $alignment, \@tags, \%alivehash, \%node_sequences, $microsatstarts, $motif);
+
+# print "returningt: ", $analayzed_mutations, \@nodes,"\n" if scalar @mutations_array > 0;;
+# print "returningy: NULL, NULL, NULL " if scalar @mutations_array == 0 && $printer == 1;
+ print "final node alignment after filling for $joint= " if $printer == 1;
+ print "$node_alignments{$joint}\n" if $printer == 1;
+
+
+ return 1, $analayzed_mutations, \@nodes, $branches_hash, \%alivehash, $node_alignments{$joint} if scalar @mutations_array > 0 ;
+ return 1, "NULL","NULL","NULL", "NULL", "NULL" if scalar @mutations_array == 0;
+}
+
+
+
+sub add_mutation{
+ my $printer = 0;
+ print "IN SUBROUTUNE add_mutation.. information received = @_\n" if $printer == 1;
+ my ($i , $bite, $to, $from) = @_;
+ print "bite = $bite.. all received info = ",join("^", @_),"\n" if $printer == 1;
+ print "to=$to\n" if $printer == 1;
+ print "tis split = ",join(" and ",split(/!/,$to)),"\n" if $printer == 1;
+ my @toields = split "!",$to;
+ print "toilds = @toields\n" if $printer == 1;
+ my @mutations=();
+
+ foreach my $toield (@toields){
+ my @toinfo=split(":",$toield);
+ print " at toinfo=@toinfo \n" if $printer == 1;
+ next if $toinfo[1] =~ /$from/i;
+ my @mutation = @toinfo if $toinfo[1] !~ /$from/i;
+ print "adding to mutaton list: ", join(",", "node=$mutation[0]","type=substitution" ,"position=$i", "from=$from", "to=$mutation[1]", "insertion=", "deletion="),"\n" if $printer == 1;
+ push (@mutations, join("\t", "node=$mutation[0]","type=substitution" ,"position=$i", "from=$from", "to=$mutation[1]", "insertion=", "deletion="));
+ }
+ return @mutations;
+}
+
+
+sub add_bases{
+
+ my $printer = 0;
+ print "IN SUBROUTUNE add_bases.. information received = @_\n" if $printer == 1;
+ my ($optional0, $optional1, $pair0, $pair1,$joint) = @_;
+ my $total_list=();
+
+ my @total_list0=split(/!/,$optional0);
+ my @total_list1=split(/!/,$optional1);
+ my @all_list=();
+ my %total_hash0=();
+ foreach my $entry (@total_list0) {
+ $entry = uc $entry;
+ $entry =~ /(\S+):(\S+)/;
+ $total_hash0{$2}=$1;
+ push @all_list, $2;
+ }
+
+ my %total_hash1=();
+ foreach my $entry (@total_list1) {
+ $entry = uc $entry;
+ $entry =~ /(\S+):(\S+)/;
+ $total_hash1{$2}=$1;
+ push @all_list, $2;
+ }
+
+ my %alphabetical_hash=();
+ my @return_options=();
+
+ for my $i (0 ... $#all_list){
+ my $alph = $all_list[$i];
+ if (exists $total_hash0{$alph} && exists $total_hash1{$alph}){
+ push(@return_options, $joint.":".$alph);
+ delete $total_hash0{$alph}; delete $total_hash1{$alph};
+ }
+ if (exists $total_hash0{$alph} && !exists $total_hash1{$alph}){
+ push(@return_options, $pair0.":".$alph);
+ delete $total_hash0{$alph};
+ }
+ if (!exists $total_hash0{$alph} && exists $total_hash1{$alph}){
+ push(@return_options, $pair1.":".$alph);
+ delete $total_hash1{$alph};
+ }
+
+ }
+
+ print "returning ",join "!",@return_options,"\n" if $printer == 1;
+ return join "!",@return_options;
+
+}
+
+
+sub fillgaps{
+# print "IN fillgaps: @_\n";
+ my ($micro, $phasesinput) = @_;
+ #print "in microsathash ,,.. micro = $micro\n";
+ return $micro if $micro !~ /\-/;
+ my $orig_micro = $micro;
+ my @phases = @$phasesinput;
+
+ my %tested_patterns = ();
+
+ foreach my $phase (@phases){
+ # print "considering phase: $phase\n";
+ my @phase_prefixes = ();
+ my @prephase_left_contexts = ();
+ my @prephase_right_contexts = ();
+ my @pregapsize = ();
+ my @prepostfilins = ();
+
+ my @phase_suffixes;
+ my @suffphase_left_contexts;
+ my @suffphase_right_contexts;
+ my @suffgapsize;
+ my @suffpostfilins;
+
+ my @postfilins = ();
+ my $motifsize = length($phases[0]);
+
+ my $change = 0;
+
+ for my $u (0 ... $motifsize-1){
+ my $concat = $phase.$phase.$phase.$phase;
+ my @concatarr = split(/\s*/, $concat);
+ my $l = 0;
+ while ($l < $u){
+ shift @concatarr;
+ $l++;
+ }
+ $concat = join ("", @concatarr);
+
+ for my $t (0 ... $motifsize-1){
+ for my $k (1 ... $motifsize-1){
+ push @phase_prefixes, substr($concat, $motifsize+$t, $k);
+ push @prephase_left_contexts, substr ($concat, $t, $motifsize);
+ push @prephase_right_contexts, substr ($concat, $motifsize+$t+$k+($motifsize-$k), 1);
+ push @pregapsize, $k;
+ push @prepostfilins, substr($concat, $motifsize+$t+$k, ($motifsize-$k));
+ # print "reading: $concat, t=$t, k=$k prefix: $prephase_left_contexts[$#prephase_left_contexts] $phase_prefixes[$#phase_prefixes] -x$pregapsize[$#pregapsize] $prephase_right_contexts[$#prephase_right_contexts]\n";
+ # print "phase_prefixes = $phase_prefixes[$#phase_prefixes]\n";
+ # print "prephase_left_contexts = $prephase_left_contexts[$#prephase_left_contexts]\n";
+ # print "prephase_right_contexts = $prephase_right_contexts[$#prephase_right_contexts]\n";
+ # print "pregapsize = $pregapsize[$#pregapsize]\n";
+ # print "prepostfilins = $prepostfilins[$#prepostfilins]\n";
+ }
+ }
+ }
+
+ # print "looking if $micro =~ /($phase\-{$motifsize})/i || $micro =~ /^(\-{$motifsize,}$phase)/i\n";
+ if ($micro =~ /($phase\-{$motifsize,})$/i || $micro =~ /^(\-{$motifsize,}$phase)/i){
+ # print "micro: $micro needs further gap removal: $1\n";
+ while ($micro =~ /$phase(\-{$motifsize,})$/i || $micro =~ /^(\-{$motifsize,})$phase/i){
+ # print "micro: $micro needs further gap removal: $1\n";
+
+ # print "phase being considered = $phase\n";
+ my $num = ();
+ $num = $micro =~ s/$phase\-{$motifsize}/$phase$phase/gi if $micro =~ /$phase\-{$motifsize,}/i;
+ $num = $micro =~ s/\-{$motifsize}$phase/$phase$phase/gi if $micro =~ /\-{$motifsize,}$phase/i;
+ # print "num = $num\n";
+ $change = 1 if $num == 1;
+ }
+ }
+
+ elsif ($micro =~ /(($phase)+)\-{$motifsize,}(($phase)+)/i){
+ while ($micro =~ /(($phase)+)\-{$motifsize,}(($phase)+)/i){
+ # print "checking lengths of $1 and $3 for $micro... \n";
+ my $num = ();
+ if (length($1) >= length($3)){
+ # print "$micro matches (($phase)+)\-{$motifsize,}(($phase)+) = $1, >= , $3 \n";
+ $num = $micro =~ s/$phase\-{$motifsize}/$phase$phase/gi ;
+ }
+ if (length($1) < length($3)){
+ # print "$micro matches (($phase)+)\-{$motifsize,}(($phase)+) = $1, < , $3 \n";
+ $num = $micro =~ s/\-{$motifsize}$phase/$phase$phase/gi ;
+ }
+ # print "micro changed to $micro\n";
+ }
+ }
+ elsif ($micro =~ /([A-Z]+)\-{$motifsize,}(($phase)+)/i){
+ while ($micro =~ /([A-Z]+)\-{$motifsize,}(($phase)+)/i){
+ # print "$micro matches ([A-Z]+)\-{$motifsize}(($phase)+) = 1=$1, - , 3=$3 \n";
+ my $num = 0;
+ $num = $micro =~ s/\-{$motifsize}$phase/$phase$phase/gi ;
+ }
+ }
+ elsif ($micro =~ /(($phase)+)\-{$motifsize,}([A-Z]+)/i){
+ while ($micro =~ /(($phase)+)\-{$motifsize,}([A-Z]+)/i){
+ # print "$micro matches (($phase)+)\-{$motifsize,}([A-Z]+) = 1=$1, - , 3=$3 \n";
+ my $num = 0;
+ $num = $micro =~ s/$phase\-{$motifsize}/$phase$phase/gi ;
+ }
+ }
+
+ # print "$orig_micro to $micro\n";
+
+ #s ;
+
+ for my $h (0 ... $#phase_prefixes){
+ # print "searching using prefix : $prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h]\n";
+ my $pattern = $prephase_left_contexts[$h].$phase_prefixes[$h].$pregapsize[$h].$prephase_right_contexts[$h];
+ # print "returning orig_micro = $orig_micro, micro = $micro \n" if exists $tested_patterns{$pattern};
+ if ($micro =~ /$prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h]/i){
+ return $orig_micro if exists $tested_patterns{$pattern};
+ while ($micro =~ /($prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h])/i){
+ $tested_patterns{$pattern} = $pattern;
+ # print "micro: $micro needs further gap removal: $1\n";
+
+ # print "prefix being considered = $phase_prefixes[$h]\n";
+ my $num = ();
+ $num = ($micro =~ s/$prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h]/$prephase_left_contexts[$h]$phase_prefixes[$h]$prepostfilins[$h]$prephase_right_contexts[$h]/gi) ;
+ # print "num = $num, micro = $micro\n";
+ $change = 1 if $num == 1;
+
+ return $orig_micro if $num > 1;
+ }
+ }
+
+ }
+ }
+ return $orig_micro if length($micro) != length($orig_micro);
+ return $micro;
+}
+
+sub selectMutationArray{
+ my $printer =0;
+
+ my $oldmutspt = $_[0];
+ my $newmutspt = $_[1];
+ my $tagstringpt = $_[2];
+ my $alivehashpt = $_[3];
+ my $alignmentpt = $_[4];
+ my $motif = $_[5];
+
+ my @alivehasharr=();
+
+ my @tags = @$tagstringpt;
+ my $alignmentln = length($alignmentpt->{$tags[0]});
+
+ foreach my $key (keys %$alivehashpt) { push @alivehasharr, $key; print "we have alive: $key\n" if $printer == 1;}
+
+ my %newside = ();
+ my %oldside = ();
+ my %newmuts = ();
+
+ my %commons = ();
+ my %olds = ();
+ foreach my $old (@$oldmutspt){
+ $olds{$old} = 1;
+ }
+ foreach my $new (@$newmutspt){
+ $commons{$new} = 1 if exists $olds{$new};;
+ }
+
+
+ foreach my $pos ( 0 ... $alignmentln){
+ #print "pos = $pos\n" if $printer == 1;
+ my $newyes = 0;
+ foreach my $mut (@$newmutspt){
+ $newmuts{$mut} = 1;
+ chomp $mut;
+ $newyes++;
+ $mut =~ s/=\t/= \t/g;
+ $mut =~ s/=$/= /g;
+
+ $mut =~ /node=([A-Z\(\), ]+)\stype=([a-zA-Z ]+)\sposition=([0-9 ]+)\sfrom=([a-zA-Z\- ]+)\sto=([a-zA-Z\- ]+)\sinsertion=([a-zA-Z\- ]+)\sdeletion=([a-zA-Z\- ]+)/;
+ my $node = $1;
+ next if $3 != $pos;
+ print "new mut = $mut\n" if $printer == 1;
+ print "node = $node, pos = $3 ... and alivehasharr = >@alivehasharr<\n" if $printer == 1;
+ my $alivenode = 0;
+ foreach my $key (@alivehasharr){
+ $alivenode = 1 if $key =~ /$node/;
+ }
+ # next if $alivenode == 0;
+ my $indel_type = " ";
+ if ($2 eq "insertion" || $2 eq "deletion"){
+ my $thisindel = ();
+ $thisindel = $6 if $2 eq "insertion";
+ $thisindel = $7 if $2 eq "deletion";
+
+ $indel_type = "i".checkIndelType($node, $thisindel, $motif,$alignmentpt,$3, $2) if $2 eq "insertion";
+ $indel_type = "d".checkIndelType($node, $thisindel, $motif,$alignmentpt, $3, $2) if $2 eq "deletion";
+ $indel_type = $indel_type."f" if $indel_type =~ /mot/ && length($thisindel) >= length($motif);
+ }
+ print "indeltype = $indel_type\n" if $printer == 1;
+ my $added = 0;
+
+ if (exists $newside{$pos} && $indel_type =~ /[a-z]+/){
+ print "we have a preexisting one for $pos\n" if $printer == 1;
+ my @preexisting = @{$newside{$pos}};
+ foreach my $pre (@preexisting){
+ print "looking at $pre\n" if $printer == 1;
+ next if $pre !~ /node=$node/;
+ next if $pre !~ /indeltype=([a-z]+)/;
+ my $currtype = $1;
+
+ if ($currtype =~ /inon/ && $indel_type =~ /dmot/){
+ delete $newside{$pos};
+ push @{$newside{$pos}}, $pre;
+ $added = 1;
+ }
+ if ($currtype =~ /dnon/ && $indel_type =~ /imot/){
+ delete $newside{$pos};
+ push @{$newside{$pos}}, $pre;
+ $added = 1;
+ }
+ if ($currtype =~ /dmot/ && $indel_type =~ /inon/){
+ delete $newside{$pos};
+ push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
+ $added = 1;
+ }
+ if ($currtype =~ /imot/ && $indel_type =~ /dnon/){
+ delete $newside{$pos};
+ push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
+ $added = 1;
+ }
+ }
+ }
+
+ print "added = $added\n" if $printer == 1;
+ push @{$newside{$pos}}, $mut."\tindeltype=$indel_type" if $added == 0;
+
+ print "for new pos,: $pos we have: @{$newside{$pos}}\n " if $printer == 1;
+ }
+ }
+
+ foreach my $pos ( 0 ... $alignmentln){
+ my $oldyes = 0;
+ foreach my $mut (@$oldmutspt){
+ chomp $mut;
+ $oldyes++;
+ $mut =~ s/=\t/= \t/g;
+ $mut =~ s/=$/= /g;
+ $mut =~ /node=([A-Z\(\), ]+)\ttype=([a-zA-Z ]+)\tposition=([0-9 ]+)\tfrom=([a-zA-Z\- ]+)\tto=([a-zA-Z\- ]+)\tinsertion=([a-zA-Z\- ]+)\tdeletion=([a-zA-Z\- ]+)/;
+ my $node = $1;
+ next if $3 != $pos;
+ print "old mut = $mut\n" if $printer == 1;
+ my $alivenode = 0;
+ foreach my $key (@alivehasharr){
+ $alivenode = 1 if $key =~ /$node/;
+ }
+ #next if $alivenode == 0;
+ my $indel_type = " ";
+ if ($2 eq "insertion" || $2 eq "deletion"){
+ $indel_type = "i".checkIndelType($node, $6, $motif,$alignmentpt, $3, $2) if $2 eq "insertion";
+ $indel_type = "d".checkIndelType($node, $7, $motif,$alignmentpt, $3, $2) if $2 eq "deletion";
+ next if $indel_type =~/non/;
+ }
+ else{ next;}
+
+ my $imp=0;
+ $imp = 1 if $indel_type =~ /dmot/ && $alivenode == 0;
+ $imp = 1 if $indel_type =~ /imot/ && $alivenode == 1;
+
+
+ if (exists $newside{$pos} && $indel_type =~ /[a-z]+/){
+ my @preexisting = @{$newside{$pos}};
+ print "we have a preexisting one for $pos: @preexisting\n" if $printer == 1;
+ next if $imp == 0;
+
+ if (scalar(@preexisting) == 1){
+ my $foundmut = $preexisting[0];
+ $foundmut=~ /node=([A-Z, \(\)]+)/;
+ next if $1 eq $node;
+
+ if (exists $oldside{$pos} || exists $commons{$foundmut}){
+ print "not replacing, but just adding\n" if $printer == 1;
+ push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
+ push @{$oldside{$pos}}, $mut."\tindeltype=$indel_type";
+ next;
+ }
+
+ delete $newside{$pos};
+ push @{$oldside{$pos}}, $mut."\tindeltype=$indel_type";
+ push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
+ print "now new one is : @{$newside{$pos}}\n" if $printer == 1;
+ }
+
+ print "for pos: $pos: @{$newside{$pos}}\n" if $printer == 1;
+ next;
+ }
+
+
+ my @news = @{$newside{$pos}} if exists $newside{$pos};
+
+ print "mut = $mut and news = @news\n" if $printer == 1;
+ push @{$oldside{$pos}}, $mut."\tindeltype=$indel_type";
+ push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
+ }
+ }
+
+ print "in the end, our collected mutations = \n" if $printer == 1;
+ my @returnarr = ();
+ foreach my $key (keys %newside) {push @returnarr,@{$newside{$key}};}
+ print join("\n", @returnarr),"\n" if $printer == 1;
+ #;
+ return @returnarr;
+
+}
+
+
+sub checkIndelType{
+ my $printer = 0;
+ my $node = $_[0];
+ my $indel = $_[1];
+ my $motif = $_[2];
+ my $alignmentpt = $_[3];
+ my $posit = $_[4];
+ my $type = $_[5];
+ my @phases =();
+ my %prephases = ();
+ my %postphases = ();
+ #print "motif = $motif\n";
+ print "IN checkIndelType ... received: @_\n" if $printer == 1;
+ my $concat = $motif.$motif.$motif.$motif;
+ my $motiflength = length($motif);
+
+ if ($motiflength > length ($indel)){
+ return "non" if $motif !~ /$indel/i;
+ return checkIndelType_ComplexAnalysis($node, $indel, $motif, $alignmentpt, $posit, $type);
+ }
+
+ my $firstpass = 0;
+ for my $y (0 ... $motiflength-1){
+ my $phase = substr($concat, $motiflength+$y, $motiflength);
+ push @phases, $phase;
+ $firstpass = 1 if $indel =~ /$phase/i;
+ for my $k (0 ... length($motif)-1){
+ print "at: motiflength=$motiflength , y=$y , k=$k.. for pre: $motiflength+$y-$k and post: $motiflength+$y-$k+$motiflength in $concat\n" if $printer == 1;
+ my $pre = substr($concat, $motiflength+$y-$k, $k );
+ my $post = substr($concat, $motiflength+$y+$motiflength, $k);
+ print "adding to phases : $phase - $pre and $post\n" if $printer == 1;
+ push @{$prephases{$phase}} , $pre;
+ push @{$postphases{$phase}} , $post;
+ }
+
+ }
+ print "firstpass 1= $firstpass\n" if $printer == 1;
+ return "non" if $firstpass ==0;
+ $firstpass =0;
+
+ foreach my $phase (@phases){
+ my @pres = @{$prephases{$phase}};
+ my @posts = @{$postphases{$phase}};
+
+ foreach my $pre (@pres){
+ foreach my $post (@posts){
+
+ $firstpass = 1 if $indel =~ /($pre)?($phase)+($post)?/i && length($indel) > (3 * length($motif));
+ $firstpass = 1 if $indel =~ /^($pre)?($phase)+($post)?$/i && length($indel) < (3 * length($motif));
+ print "matched here : ($pre)?($phase)+($post)?\n" if $printer == 1;
+ last if $firstpass == 1;
+ }
+ last if $firstpass == 1;
+ }
+ last if $firstpass == 1;
+ }
+
+ print "firstpass 2= $firstpass\n" if $printer == 1;
+ return "non" if $firstpass ==0;
+ return "mot" if $firstpass ==1;
+}
+
+
+sub checkIndelType_ComplexAnalysis{
+ my $printer = 0;
+ my $node = $_[0];
+ my $indel = $_[1];
+ my $motif = $_[2];
+ my $alignmentpt = $_[3];
+ my $pos = $_[4];
+ my $type = $_[5];
+ my @speciesinvolved = $node =~ /[A-Z]+/g;
+
+ my @seqs = ();
+ my $residualseq = length($motif) - length($indel);
+ print "IN COMPLEX ANALYSIS ... received: @_ .... speciesinvolved = @speciesinvolved\n" if $printer == 1;
+ print "we have position = $pos, sseq = $alignmentpt->{$speciesinvolved[0]}\n" if $printer == 1;
+ print "residualseq = $residualseq\n" if $printer == 1;
+ print "pos=$pos... got: @_\n" if $printer == 1;
+ foreach my $sp (@speciesinvolved){
+ my $spseq = $alignmentpt->{$sp};
+ #print "orig spseq = $spseq\n";
+ my $subseq = ();
+
+ if ($type eq "deletion"){
+ my @indelparts = split(/\s*/,$indel);
+ my @seqparts = split(/\s*/,$spseq);
+
+ for my $p ($pos ... $pos+length($indel)-1){
+ $seqparts[$p] = shift @indelparts;
+ }
+ $spseq = join("",@seqparts);
+ }
+ #print "mod spseq = $spseq\n";
+ # $spseq=~ s/\-//g if $type !~ /deletion/;
+
+ print "substr($spseq, $pos-($residualseq), length($indel)+$residualseq+$residualseq)\n" if $pos > 0 && $pos < (length($spseq) - length($motif)) && $printer == 1;
+ print "substr($spseq, 0, length($indel)+$residualseq)\n" if $pos == 0 && $printer == 1;
+ print "substr($spseq, $pos - $residualseq, length($indel)+$residualseq)\n" if $pos >= (length($spseq) - length($motif)) && $printer == 1;
+
+ $subseq = substr($spseq, $pos-($residualseq), length($indel)+$residualseq+$residualseq) if $pos > 0 && $pos < (length($spseq) - length($motif)) ;
+ $subseq = substr($spseq, 0, length($indel)+$residualseq) if $pos == 0;
+ $subseq = substr($spseq, $pos - $residualseq, length($indel)+$residualseq) if $pos >= (length($spseq) - length($motif)) ;
+ print "spseq = $spseq . subseq=$subseq . type = $type\n" if $printer == 1;
+ # if $subseq !~ /[a-z\-]/i;
+ $subseq =~ s/\-/$indel/g if $type =~ /insertion/;
+ push @seqs, $subseq;
+ print "seqs = @seqs\n" if $printer == 1;
+ }
+ return "non" if checkIfSeqsIdentical(@seqs) eq "NO";
+
+ print "checking for $seqs[0] \n" if $printer == 1;
+
+ my @phases =();
+ my %prephases = ();
+ my %postphases = ();
+ my $concat = $motif.$motif.$motif.$motif;
+ my $motiflength = length($motif);
+
+ my $firstpass = 0;
+
+ for my $y (0 ... $motiflength-1){
+ my $phase = substr($concat, $motiflength+$y, $motiflength);
+ push @phases, $phase;
+ $firstpass = 1 if $seqs[0] =~ /$phase/i;
+ for my $k (0 ... length($motif)-1){
+ my $pre = substr($concat, $motiflength+$y-$k, $k );
+ my $post = substr($concat, $motiflength+$y+$motiflength, $k);
+ print "adding to phases : $phase - $pre and $post\n" if $printer == 1;
+ push @{$prephases{$phase}} , $pre;
+ push @{$postphases{$phase}} , $post;
+ }
+
+ }
+ print "firstpass 1= $firstpass.. also, res-d = ",(length($seqs[0]))%(length($motif)),"\n" if $printer == 1;
+ return "non" if $firstpass ==0;
+ $firstpass =0;
+ foreach my $phase (@phases){
+
+ $firstpass = 1 if $seqs[0] =~ /^($phase)+$/i && ((length($seqs[0]))%(length($motif))) == 0;
+
+ if (((length($seqs[0]))%(length($motif))) != 0){
+ my @pres = @{$prephases{$phase}};
+ my @posts = @{$postphases{$phase}};
+ foreach my $pre (@pres){
+ foreach my $post (@posts){
+ next if $pre !~ /\S/ && $post !~ /\S/;
+ $firstpass = 1 if ($seqs[0] =~ /^($pre)($phase)+($post)$/i || $seqs[0] =~ /^($pre)($phase)+$/i || $seqs[0] =~ /^($phase)+($post)$/i);
+ print "caught with $pre $phase $post\n" if $printer == 1;
+ last if $firstpass == 1;
+ }
+ last if $firstpass == 1;
+ }
+ }
+
+ last if $firstpass == 1;
+ }
+
+ #print "indel = $indel.. motif = $motif.. firstpass 2= mot\n" if $firstpass ==1;
+ #print "indel = $indel.. motif = $motif.. firstpass 2= non\n" if $firstpass ==0;
+ #;# if $firstpass ==1;
+ return "non" if $firstpass ==0;
+ return "mot" if $firstpass ==1;
+
+}
+
+sub checkIfSeqsIdentical{
+ my @seqs = @_;
+ my $identical = 1;
+
+ for my $j (1 ... $#seqs){
+ $identical = 0 if uc($seqs[0]) ne uc($seqs[$j]);
+ }
+ return "NO" if $identical == 0;
+ return "YES" if $identical == 1;
+
+}
+
+sub summarizeMutations{
+ my $mutspt = $_[0];
+ my @muts = @$mutspt;
+ my $tree = $_[1];
+
+ my @returnarr = ();
+
+ for (1 ... 38){
+ push @returnarr, "NA";
+ }
+ push @returnarr, "NULL";
+ return @returnarr if $tree eq "NULL" || scalar(@muts) < 1;
+
+
+ my @bspecies = ();
+ my @dspecies = ();
+ my $treecopy = $tree;
+ $treecopy =~ s/[\(\)]//g;
+ my @treeparts = split(/[\.,]+/, $treecopy);
+
+ for my $part (@treeparts){
+ if ($part =~ /\+/){
+ $part =~ s/\+//g;
+ #my @sp = split(/\s*/, $part);
+ #foreach my $p (@sp) {push @bspecies, $p;}
+ push @bspecies, $part;
+ }
+ if ($part =~ /\-/){
+ $part =~ s/\-//g;
+ #my @sp = split(/\s*/, $part);
+ #foreach my $p (@sp) {push @dspecies, $p;}
+ push @dspecies, $part;
+ }
+
+ }
+ #print "-------------------------------------------------------\n";
+
+ my ($insertions, $deletions, $motinsertions, $motinsertionsf, $motdeletions, $motdeletionsf, $noninsertions, $nondeletions) = (0,0,0,0,0,0,0,0);
+ my ($binsertions, $bdeletions, $bmotinsertions,$bmotinsertionsf, $bmotdeletions, $bmotdeletionsf, $bnoninsertions, $bnondeletions) = (0,0,0,0,0,0,0,0);
+ my ($dinsertions, $ddeletions, $dmotinsertions,$dmotinsertionsf, $dmotdeletions, $dmotdeletionsf, $dnoninsertions, $dnondeletions) = (0,0,0,0,0,0,0,0);
+ my ($ninsertions, $ndeletions, $nmotinsertions,$nmotinsertionsf, $nmotdeletions, $nmotdeletionsf, $nnoninsertions, $nnondeletions) = (0,0,0,0,0,0,0,0);
+ my ($substitutions, $bsubstitutions, $dsubstitutions, $nsubstitutions, $indels, $subs) = (0,0,0,0,"NA","NA");
+
+ my @insertionsarr = (" ");
+ my @deletionsarr = (" ");
+
+ my @substitutionsarr = (" ");
+
+
+ foreach my $mut (@muts){
+ # print "mut = $mut\n";
+ chomp $mut;
+ $mut =~ s/=\t/= /g;
+ $mut =~ s/=$/= /g;
+ my %mhash = ();
+ my @mields = split(/\t/,$mut);
+
+ foreach my $m (@mields){
+ my @fields = split(/=/,$m);
+ next if $fields[1] eq " ";
+ $mhash{$fields[0]} = $fields[1];
+ }
+
+ my $myutype = ();
+ my $decided = 0;
+
+ my $localnode = $mhash{"node"};
+ $localnode =~ s/[\(\)\. ,]//g;
+
+
+ foreach my $s (@bspecies){
+ if ($localnode eq $s) {
+ $decided = 1; $myutype = "b";
+ }
+ }
+
+ foreach my $s (@dspecies){
+ if ($localnode eq $s) {
+ $decided = 1; $myutype = "d";
+ }
+ }
+
+ $myutype = "n" if $decided != 1;
+
+
+ # print "tree=$tree, birth species=@bspecies, death species=@dspecies, node=$mhash{node} .. myutype=$myutype .. \n";
+ # if $mhash{"type"} eq "insertion" && $myutype eq "b";
+
+
+ if ($mhash{"type"} eq "substitution"){
+ $substitutions++;
+ $bsubstitutions++ if $myutype eq "b";
+ $dsubstitutions++ if $myutype eq "d";
+ $nsubstitutions++ if $myutype eq "n";
+ # print "substitution: from= $mhash{from}, to = $mhash{to}, and type = myutype\n";
+ push @substitutionsarr, "b:$mhash{position}:".$mhash{"from"}.">".$mhash{"to"} if $myutype eq "b";
+ push @substitutionsarr, "d:$mhash{position}:".$mhash{"from"}.">".$mhash{"to"} if $myutype eq "d";
+ push @substitutionsarr, "n:$mhash{position}:".$mhash{"from"}.">".$mhash{"to"} if $myutype eq "n";
+ # print "substitutionsarr = @substitutionsarr\n";
+ # ;
+ }
+ else{
+ #print "tree=$tree, birth species=@bspecies, death species=@dspecies, node=$mhash{node} .. myutype=$myutype .. indeltype=$mhash{indeltype}\n";
+ if ($mhash{"type"} eq "deletion"){
+ $deletions++;
+
+ $motdeletions++ if $mhash{"indeltype"} =~ /dmot/;
+ $motdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/;
+
+ $nondeletions++ if $mhash{"indeltype"} =~ /dnon/;
+
+ $bdeletions++ if $myutype eq "b";
+ $ddeletions++ if $myutype eq "d";
+ $ndeletions++ if $myutype eq "n";
+
+ $bmotdeletions++ if $mhash{"indeltype"} =~ /dmot/ && $myutype eq "b";
+ $bmotdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/ && $myutype eq "b";
+ $bnondeletions++ if $mhash{"indeltype"} =~ /dnon/ && $myutype eq "b";
+
+ $dmotdeletions++ if $mhash{"indeltype"} =~ /dmot/ && $myutype eq "d";
+ $dmotdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/ && $myutype eq "d";
+ $dnondeletions++ if $mhash{"indeltype"} =~ /dnon/ && $myutype eq "d";
+
+ $nmotdeletions++ if $mhash{"indeltype"} =~ /dmot/ && $myutype eq "n";
+ $nmotdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/ && $myutype eq "n";
+ $nnondeletions++ if $mhash{"indeltype"} =~ /dnon/ && $myutype eq "n";
+
+ push @deletionsarr, "b:$mhash{indeltype}:$mhash{position}:".$mhash{"deletion"} if $myutype eq "b";
+ push @deletionsarr, "d:$mhash{indeltype}:$mhash{position}:".$mhash{"deletion"} if $myutype eq "d";
+ push @deletionsarr, "n:$mhash{indeltype}:$mhash{position}:".$mhash{"deletion"} if $myutype eq "n";
+ }
+
+ if ($mhash{"type"} eq "insertion"){
+ $insertions++;
+
+ $motinsertions++ if $mhash{"indeltype"} =~ /imot/;
+ $motinsertionsf++ if $mhash{"indeltype"} =~ /imotf/;
+ $noninsertions++ if $mhash{"indeltype"} =~ /inon/;
+
+ $binsertions++ if $myutype eq "b";
+ $dinsertions++ if $myutype eq "d";
+ $ninsertions++ if $myutype eq "n";
+
+ $bmotinsertions++ if $mhash{"indeltype"} =~ /imot/ && $myutype eq "b";
+ $bmotinsertionsf++ if $mhash{"indeltype"} =~ /imotf/ && $myutype eq "b";
+ $bnoninsertions++ if $mhash{"indeltype"} =~ /inon/ && $myutype eq "b";
+
+ $dmotinsertions++ if $mhash{"indeltype"} =~ /imot/ && $myutype eq "d";
+ $dmotinsertionsf++ if $mhash{"indeltype"} =~ /imotf/ && $myutype eq "d";
+ $dnoninsertions++ if $mhash{"indeltype"} =~ /inon/ && $myutype eq "d";
+
+ $nmotinsertions++ if $mhash{"indeltype"} =~ /imot/ && $myutype eq "n";
+ $nmotinsertionsf++ if $mhash{"indeltype"} =~ /imotf/ && $myutype eq "n";
+ $nnoninsertions++ if $mhash{"indeltype"} =~ /inon/ && $myutype eq "n";
+
+ push @insertionsarr, "b:$mhash{indeltype}:$mhash{position}:".$mhash{"insertion"} if $myutype eq "b";
+ push @insertionsarr, "d:$mhash{indeltype}:$mhash{position}:".$mhash{"insertion"} if $myutype eq "d";
+ push @insertionsarr, "n:$mhash{indeltype}:$mhash{position}:".$mhash{"insertion"} if $myutype eq "n";
+
+ }
+ }
+ }
+
+
+
+ $indels = "ins=".join(",",@insertionsarr).";dels=".join(",",@deletionsarr) if scalar(@insertionsarr) > 1 || scalar(@deletionsarr) > 1 ;
+ $subs = join(",",@substitutionsarr) if scalar(@substitutionsarr) > 1;
+ $indels =~ s/ //g;
+ $subs =~ s/ //g ;
+
+ #print "indels = $indels, subs=$subs\n";
+ ## if $indels =~ /[a-zA-Z0-9]/ || $subs =~ /[a-zA-Z0-9]/ ;
+ #print "tree = $tree, indels = $indels, subs = $subs, bspecies = @bspecies, dspecies = @dspecies \n";
+ my @returnarray = ();
+
+ push (@returnarray, $insertions, $deletions, $motinsertions, $motinsertionsf, $motdeletions, $motdeletionsf, $noninsertions, $nondeletions) ;
+ push (@returnarray, $binsertions, $bdeletions, $bmotinsertions,$bmotinsertionsf, $bmotdeletions, $bmotdeletionsf, $bnoninsertions, $bnondeletions) ;
+ push (@returnarray, $dinsertions, $ddeletions, $dmotinsertions,$dmotinsertionsf, $dmotdeletions, $dmotdeletionsf, $dnoninsertions, $dnondeletions) ;
+ push (@returnarray, $ninsertions, $ndeletions, $nmotinsertions,$nmotinsertionsf, $nmotdeletions, $nmotdeletionsf, $nnoninsertions, $nnondeletions) ;
+ push (@returnarray, $substitutions, $bsubstitutions, $dsubstitutions, $nsubstitutions, $indels, $subs) ;
+
+ push @returnarray, $tree;
+
+ my @copy = @returnarray;
+ return (@returnarray);
+
+}
+
+sub selectBetterTree{
+ my $printer = 0;
+ my $treestudy = $_[0];
+ my $alt = $_[1];
+ my $mutspt = $_[2];
+ my @muts = @$mutspt;
+ my @trees = (); my @alternatetrees=();
+
+ @trees = split(/\|/,$treestudy) if $treestudy =~ /\|/;
+ @alternatetrees = split(/[\|;]/,$alt) if $alt =~ /[\|;\(\)]/;
+
+ $trees[0] = $treestudy if $treestudy !~ /\|/;
+ $alternatetrees[0] = $alt if $alt !~ /[\|;\(\)]/;
+
+ my @alltrees = (@trees, @alternatetrees);
+# push(@alltrees,@alternatetrees);
+
+ my %mutspecies = ();
+
+ print "IN selectBetterTree..treestudy=$treestudy. alt=$alt. for: @_. trees=@trees<. alternatetrees=@alternatetrees\n" if $printer == 1;
+ #;
+ foreach my $mut (@muts){
+ print colored ['green'],"mut = $mut\n" if $printer == 1;
+ $mut =~ /node=([A-Z,\(\) ]+)/;
+ my $node = $1;
+ $node =~s/[,\(\) ]+//g;
+ my @indivspecies = $node =~ /[A-Z]+/g;
+ #print "adding node: $node\n" if $printer == 1;
+ $mutspecies{$node} = $node;
+
+ #foreach (@indivspecies) {
+ #$mutspecies{$mut} = $_; #print "for $_ adding $mutspecies{$_}\n";
+ #}
+
+ }
+
+ my @treerecords = ();
+ my $treecount = -1;
+ foreach my $tree (@alltrees){
+ print "checking with tree $tree\n" if $printer == 1;
+ $treecount++;
+ $treerecords[$treecount] = 0;
+ my @indivspecies = ($tree =~ /[A-Z]+/g);
+ print "indivspecies=@indivspecies\n" if $printer == 1;
+ foreach my $species (@indivspecies){
+ print "checkin if exists species: $species\n" if $printer == 1;
+ $treerecords[$treecount]+=2 if exists $mutspecies{$species} && $mutspecies{$species} !~ /indeltype=[a-z]mot/;
+ $treerecords[$treecount]+=1.5 if exists $mutspecies{$species} && $mutspecies{$species} =~ /indeltype=[a-z]mot/;
+ $treerecords[$treecount]-- if !exists $mutspecies{$species};
+ }
+
+ print "for tree $tree, our treecount = $treerecords[$treecount]\n" if $printer == 1;
+ }
+
+ my @best_tree = array_largest_number_arrayPosition(@treerecords);
+ print "treerecords = @treerecords. hence, best tree = @best_tree\n" if $printer == 1;
+
+ return ($alltrees[$best_tree[0]], $treerecords[$best_tree[0]]) if scalar(@best_tree) == 1;
+ print "best_tree[0] = $best_tree[0], and treerecords = $treerecords[$best_tree[0]]\n" if $printer == 1;
+ return ("NULL", -1) if $treerecords[$best_tree[0]] < 1;
+ my $rando = int(rand($#trees));
+ return ($alltrees[$rando], $treerecords[$rando]) if scalar(@best_tree) > 1;
+
+}
+
+
+
+
+sub load_sameHash{
+ #my $g = %$_[0];
+ $sameHash{"CAGT"}="AGTC";
+ $sameHash{"ATGA"}="AATG";
+ $sameHash{"CAAC"}="AACC";
+ $sameHash{"GGAA"}="AAGG";
+ $sameHash{"TAAG"}="AAGT";
+ $sameHash{"CGAG"}="AGCG";
+ $sameHash{"TAGG"}="AGGT";
+ $sameHash{"GCAG"}="AGGC";
+ $sameHash{"TAGA"}="ATAG";
+ $sameHash{"TGA"}="ATG";
+ $sameHash{"CAAG"}="AAGC";
+ $sameHash{"CTAA"}="AACT";
+ $sameHash{"CAAT"}="AATC";
+ $sameHash{"GTAG"}="AGGT";
+ $sameHash{"GAAG"}="AAGG";
+ $sameHash{"CGA"}="ACG";
+ $sameHash{"GTAA"}="AAGT";
+ $sameHash{"ACAA"}="AAAC";
+ $sameHash{"GCGG"}="GGGC";
+ $sameHash{"ATCA"}="AATC";
+ $sameHash{"TAAC"}="AACT";
+ $sameHash{"GGCA"}="AGGC";
+ $sameHash{"TGAG"}="AGTG";
+ $sameHash{"AACA"}="AAAC";
+ $sameHash{"GAGC"}="AGCG";
+ $sameHash{"ACCA"}="AACC";
+ $sameHash{"TGAA"}="AATG";
+ $sameHash{"ACA"}="AAC";
+ $sameHash{"GAAC"}="AACG";
+ $sameHash{"GCA"}="AGC";
+ $sameHash{"CCAC"}="ACCC";
+ $sameHash{"CATA"}="ATAC";
+ $sameHash{"CAC"}="ACC";
+ $sameHash{"TACA"}="ATAC";
+ $sameHash{"GGAC"}="ACGG";
+ $sameHash{"AGA"}="AAG";
+ $sameHash{"ATAA"}="AAAT";
+ $sameHash{"CA"}="AC";
+ $sameHash{"CCCA"}="ACCC";
+ $sameHash{"TCAA"}="AATC";
+ $sameHash{"CAGA"}="AGAC";
+ $sameHash{"AATA"}="AAAT";
+ $sameHash{"CCA"}="ACC";
+ $sameHash{"AGAA"}="AAAG";
+ $sameHash{"AGTA"}="AAGT";
+ $sameHash{"GACG"}="ACGG";
+ $sameHash{"TCAG"}="AGTC";
+ $sameHash{"ACGA"}="AACG";
+ $sameHash{"CGCA"}="ACGC";
+ $sameHash{"GAGT"}="AGTG";
+ $sameHash{"GA"}="AG";
+ $sameHash{"TA"}="AT";
+ $sameHash{"TAA"}="AAT";
+ $sameHash{"CAG"}="AGC";
+ $sameHash{"GATA"}="ATAG";
+ $sameHash{"GTA"}="AGT";
+ $sameHash{"CCAA"}="AACC";
+ $sameHash{"TAG"}="AGT";
+ $sameHash{"CAAA"}="AAAC";
+ $sameHash{"AAGA"}="AAAG";
+ $sameHash{"CACG"}="ACGC";
+ $sameHash{"GTCA"}="AGTC";
+ $sameHash{"GGA"}="AGG";
+ $sameHash{"GGAT"}="ATGG";
+ $sameHash{"CGGG"}="GGGC";
+ $sameHash{"CGGA"}="ACGG";
+ $sameHash{"AGGA"}="AAGG";
+ $sameHash{"TAAA"}="AAAT";
+ $sameHash{"GAGA"}="AGAG";
+ $sameHash{"ACTA"}="AACT";
+ $sameHash{"GCGA"}="AGCG";
+ $sameHash{"CACA"}="ACAC";
+ $sameHash{"AGAT"}="ATAG";
+ $sameHash{"GAGG"}="AGGG";
+ $sameHash{"CGAC"}="ACCG";
+ $sameHash{"GGAG"}="AGGG";
+ $sameHash{"GCCA"}="AGCC";
+ $sameHash{"CCAG"}="AGCC";
+ $sameHash{"GAAA"}="AAAG";
+ $sameHash{"CAGG"}="AGGC";
+ $sameHash{"GAC"}="ACG";
+ $sameHash{"CAA"}="AAC";
+ $sameHash{"GACC"}="ACCG";
+ $sameHash{"GGCG"}="GGGC";
+ $sameHash{"GGTA"}="AGGT";
+ $sameHash{"AGCA"}="AAGC";
+ $sameHash{"GATG"}="ATGG";
+ $sameHash{"GTGA"}="AGTG";
+ $sameHash{"ACAG"}="AGAC";
+ $sameHash{"CGG"}="GGC";
+ $sameHash{"ATA"}="AAT";
+ $sameHash{"GACA"}="AGAC";
+ $sameHash{"GCAA"}="AAGC";
+ $sameHash{"CAGC"}="AGCC";
+ $sameHash{"GGGA"}="AGGG";
+ $sameHash{"GAG"}="AGG";
+ $sameHash{"ACAT"}="ATAC";
+ $sameHash{"GAAT"}="AATG";
+ $sameHash{"CACC"}="ACCC";
+ $sameHash{"GAT"}="ATG";
+ $sameHash{"GCG"}="GGC";
+ $sameHash{"GCAC"}="ACGC";
+ $sameHash{"GAA"}="AAG";
+ $sameHash{"TGGA"}="ATGG";
+ $sameHash{"CCGA"}="ACCG";
+ $sameHash{"CGAA"}="AACG";
+}
+
+
+
+sub load_revHash{
+ $revHash{"CTGA"}="AGTC";
+ $revHash{"TCTT"}="AAAG";
+ $revHash{"CTAG"}="AGCT";
+ $revHash{"GGTG"}="ACCC";
+ $revHash{"GCC"}="GGC";
+ $revHash{"GCTT"}="AAGC";
+ $revHash{"GCGT"}="ACGC";
+ $revHash{"GTTG"}="AACC";
+ $revHash{"CTCC"}="AGGG";
+ $revHash{"ATC"}="ATG";
+ $revHash{"CGAT"}="ATCG";
+ $revHash{"TTAA"}="AATT";
+ $revHash{"GTTC"}="AACG";
+ $revHash{"CTGC"}="AGGC";
+ $revHash{"TCGA"}="ATCG";
+ $revHash{"ATCT"}="ATAG";
+ $revHash{"GGTT"}="AACC";
+ $revHash{"CTTA"}="AAGT";
+ $revHash{"TGGC"}="AGCC";
+ $revHash{"CCG"}="GGC";
+ $revHash{"CGGC"}="GGCC";
+ $revHash{"TTAG"}="AACT";
+ $revHash{"GTG"}="ACC";
+ $revHash{"CTTT"}="AAAG";
+ $revHash{"TGCA"}="ATGC";
+ $revHash{"CGCT"}="AGCG";
+ $revHash{"TTCC"}="AAGG";
+ $revHash{"CT"}="AG";
+ $revHash{"C"}="G";
+ $revHash{"CTCT"}="AGAG";
+ $revHash{"ACTT"}="AAGT";
+ $revHash{"GGTC"}="ACCG";
+ $revHash{"ATTC"}="AATG";
+ $revHash{"GGGT"}="ACCC";
+ $revHash{"CCTA"}="AGGT";
+ $revHash{"CGCG"}="GCGC";
+ $revHash{"GTGT"}="ACAC";
+ $revHash{"GCCC"}="GGGC";
+ $revHash{"GTCG"}="ACCG";
+ $revHash{"TCCC"}="AGGG";
+ $revHash{"TTCA"}="AATG";
+ $revHash{"AGTT"}="AACT";
+ $revHash{"CCCT"}="AGGG";
+ $revHash{"CCGC"}="GGGC";
+ $revHash{"CTT"}="AAG";
+ $revHash{"TTGG"}="AACC";
+ $revHash{"ATT"}="AAT";
+ $revHash{"TAGC"}="AGCT";
+ $revHash{"ACTG"}="AGTC";
+ $revHash{"TCAC"}="AGTG";
+ $revHash{"CTGT"}="AGAC";
+ $revHash{"TGTG"}="ACAC";
+ $revHash{"ATCC"}="ATGG";
+ $revHash{"GTGG"}="ACCC";
+ $revHash{"TGGG"}="ACCC";
+ $revHash{"TCGG"}="ACCG";
+ $revHash{"CGGT"}="ACCG";
+ $revHash{"GCTC"}="AGCG";
+ $revHash{"TACG"}="ACGT";
+ $revHash{"GTTT"}="AAAC";
+ $revHash{"CAT"}="ATG";
+ $revHash{"CATG"}="ATGC";
+ $revHash{"GTTA"}="AACT";
+ $revHash{"CACT"}="AGTG";
+ $revHash{"TCAT"}="AATG";
+ $revHash{"TTA"}="AAT";
+ $revHash{"TGTA"}="ATAC";
+ $revHash{"TTTC"}="AAAG";
+ $revHash{"TACT"}="AAGT";
+ $revHash{"TGTT"}="AAAC";
+ $revHash{"CTA"}="AGT";
+ $revHash{"GACT"}="AGTC";
+ $revHash{"TTGC"}="AAGC";
+ $revHash{"TTC"}="AAG";
+ $revHash{"GCT"}="AGC";
+ $revHash{"GCAT"}="ATGC";
+ $revHash{"TGGT"}="AACC";
+ $revHash{"CCT"}="AGG";
+ $revHash{"CATC"}="ATGG";
+ $revHash{"CCAT"}="ATGG";
+ $revHash{"CCCG"}="GGGC";
+ $revHash{"TGCC"}="AGGC";
+ $revHash{"TG"}="AC";
+ $revHash{"TGCT"}="AAGC";
+ $revHash{"GCCG"}="GGCC";
+ $revHash{"TCTG"}="AGAC";
+ $revHash{"TGT"}="AAC";
+ $revHash{"TTAT"}="AAAT";
+ $revHash{"TAGT"}="AACT";
+ $revHash{"TATG"}="ATAC";
+ $revHash{"TTTA"}="AAAT";
+ $revHash{"CGTA"}="ACGT";
+ $revHash{"TA"}="AT";
+ $revHash{"TGTC"}="AGAC";
+ $revHash{"CTAT"}="ATAG";
+ $revHash{"TATA"}="ATAT";
+ $revHash{"TAC"}="AGT";
+ $revHash{"TC"}="AG";
+ $revHash{"CATT"}="AATG";
+ $revHash{"TCG"}="ACG";
+ $revHash{"ATTT"}="AAAT";
+ $revHash{"CGTG"}="ACGC";
+ $revHash{"CTG"}="AGC";
+ $revHash{"TCGT"}="AACG";
+ $revHash{"TCCG"}="ACGG";
+ $revHash{"GTT"}="AAC";
+ $revHash{"ATGT"}="ATAC";
+ $revHash{"CTTG"}="AAGC";
+ $revHash{"CCTT"}="AAGG";
+ $revHash{"GATC"}="ATCG";
+ $revHash{"CTGG"}="AGCC";
+ $revHash{"TTCT"}="AAAG";
+ $revHash{"CGTC"}="ACGG";
+ $revHash{"CG"}="GC";
+ $revHash{"TATT"}="AAAT";
+ $revHash{"CTCG"}="AGCG";
+ $revHash{"TCTC"}="AGAG";
+ $revHash{"TCCT"}="AAGG";
+ $revHash{"TGG"}="ACC";
+ $revHash{"ACTC"}="AGTG";
+ $revHash{"CTC"}="AGG";
+ $revHash{"CGC"}="GGC";
+ $revHash{"TTG"}="AAC";
+ $revHash{"ACCT"}="AGGT";
+ $revHash{"TCTA"}="ATAG";
+ $revHash{"GTAC"}="ACGT";
+ $revHash{"TTGA"}="AATC";
+ $revHash{"GTCC"}="ACGG";
+ $revHash{"GATT"}="AATC";
+ $revHash{"T"}="A";
+ $revHash{"CGTT"}="AACG";
+ $revHash{"GTC"}="ACG";
+ $revHash{"GCCT"}="AGGC";
+ $revHash{"TGC"}="AGC";
+ $revHash{"TTTG"}="AAAC";
+ $revHash{"GGCT"}="AGCC";
+ $revHash{"TCA"}="ATG";
+ $revHash{"GTGC"}="ACGC";
+ $revHash{"TGAT"}="AATC";
+ $revHash{"TAT"}="AAT";
+ $revHash{"CTAC"}="AGGT";
+ $revHash{"TGCG"}="ACGC";
+ $revHash{"CTCA"}="AGTG";
+ $revHash{"CTTC"}="AAGG";
+ $revHash{"GCTG"}="AGCC";
+ $revHash{"TATC"}="ATAG";
+ $revHash{"TAAT"}="AATT";
+ $revHash{"ACT"}="AGT";
+ $revHash{"TCGC"}="AGCG";
+ $revHash{"GGT"}="ACC";
+ $revHash{"TCC"}="AGG";
+ $revHash{"TTGT"}="AAAC";
+ $revHash{"TGAC"}="AGTC";
+ $revHash{"TTAC"}="AAGT";
+ $revHash{"CGT"}="ACG";
+ $revHash{"ATTA"}="AATT";
+ $revHash{"ATTG"}="AATC";
+ $revHash{"CCTC"}="AGGG";
+ $revHash{"CCGG"}="GGCC";
+ $revHash{"CCGT"}="ACGG";
+ $revHash{"TCCA"}="ATGG";
+ $revHash{"CGCC"}="GGGC";
+ $revHash{"GT"}="AC";
+ $revHash{"TTCG"}="AACG";
+ $revHash{"CCTG"}="AGGC";
+ $revHash{"TCT"}="AAG";
+ $revHash{"GTAT"}="ATAC";
+ $revHash{"GTCT"}="AGAC";
+ $revHash{"GCTA"}="AGCT";
+ $revHash{"TACC"}="AGGT";
+}
+
+
+sub allCaps{
+ my $motif = $_[0];
+ $motif =~ s/a/A/g;
+ $motif =~ s/c/C/g;
+ $motif =~ s/t/T/g;
+ $motif =~ s/g/G/g;
+ return $motif;
+}
+
+
+sub all_caps{
+ my @strand = split(/\s*/,$_[0]);
+ for my $i (0 ... $#strand){
+ if ($strand[$i] =~ /c/) {$strand[$i] = "C";next;}
+ if ($strand[$i] =~ /a/) {$strand[$i] = "A";next;}
+ if ($strand[$i] =~ /t/) { $strand[$i] = "T";next;}
+ if ($strand[$i] =~ /g/) {$strand[$i] = "G";next;}
+ }
+ return join("",@strand);
+}
+sub array_mean{
+ return "NA" if scalar(@_) == 0;
+ my $sum = 0;
+ foreach my $val (@_){
+ $sum = $sum + $val;
+ }
+ return ($sum/scalar(@_));
+}
+sub array_sum{
+ return "NA" if scalar(@_) == 0;
+ my $sum = 0;
+ foreach my $val (@_){
+ $sum = $sum + $val;
+ }
+ return ($sum);
+}
+
+sub variance{
+ return "NA" if scalar(@_) == 0;
+ return 0 if scalar(@_) == 1;
+ my $mean = array_mean(@_);
+ my $num = 0;
+ return 0 if scalar(@_) == 1;
+# print "mean = $mean .. array = >@_<\n";
+ foreach my $ele (@_){
+ # print "$num = $num + ($ele-$mean)*($ele-$mean)\n";
+ $num = $num + ($ele-$mean)*($ele-$mean);
+ }
+ my $var = $num / scalar(@_);
+ return $var;
+}
+
+sub array_95confIntervals{
+ return "NA" if scalar(@_) <= 0;
+ my @sorted = sort { $a <=> $b } @_;
+# print "@sorted=",scalar(@sorted), "\n";
+ my $aDeechNo = int((scalar(@sorted) * 2.5) / 100);
+ my $saaDeNo = int((scalar(@sorted) * 97.5) / 100);
+
+ return ($sorted[$aDeechNo], $sorted[$saaDeNo]);
+}
+
+sub array_median{
+ return "NA" if scalar(@_) == 0;
+ return $_[0] if scalar(@_) == 1;
+ my @sorted = sort { $a <=> $b } @_;
+ my $totalno = scalar(@sorted);
+
+ #print "sorted = @sorted\n";
+
+ my $pick = ();
+ if ($totalno % 2 == 1){
+ #print "odd set .. totalno = $totalno\n";
+ my $mid = $totalno / 2;
+ my $onehalfno = $mid - $mid % 1;
+ my $secondhalfno = $onehalfno + 1;
+ my $onehalf = $sorted[$onehalfno-1];
+ my $secondhalf = $sorted[$secondhalfno-1];
+ #print "onehalfno = $onehalfno and secondhalfno = $secondhalfno \n onehalf = $onehalf and secondhalf = $secondhalf\n";
+
+ $pick = $secondhalf;
+ }
+ else{
+ #print "even set .. totalno = $totalno\n";
+ my $mid = $totalno / 2;
+ my $onehalfno = $mid;
+ my $secondhalfno = $onehalfno + 1;
+ my $onehalf = $sorted[$onehalfno-1];
+ my $secondhalf = $sorted[$secondhalfno-1];
+ #print "onehalfno = $onehalfno and secondhalfno = $secondhalfno \n onehalf = $onehalf and secondhalf = $secondhalf\n";
+ $pick = ($onehalf + $secondhalf )/2;
+
+ }
+ #print "pick = $pick..\n";
+ return $pick;
+
+}
+
+
+sub array_numerical_sort{
+ return "NA" if scalar(@_) == 0;
+ my @sorted = sort { $a <=> $b } @_;
+ return (@sorted);
+}
+
+sub array_smallest_number{
+ return "NA" if scalar(@_) == 0;
+ return $_[0] if scalar(@_) == 1;
+ my @sorted = sort { $a <=> $b } @_;
+ return $sorted[0];
+}
+
+
+sub array_largest_number{
+ return "NA" if scalar(@_) == 0;
+ return $_[0] if scalar(@_) == 1;
+ my @sorted = sort { $a <=> $b } @_;
+ return $sorted[$#sorted];
+}
+
+
+sub array_largest_number_arrayPosition{
+ return "NA" if scalar(@_) == 0;
+ return 0 if scalar(@_) == 1;
+ my $maxpos = 0;
+ my @maxposes = ();
+ my @maxvals = ();
+ my $maxval = array_smallest_number(@_);
+ for my $i (0 ... $#_){
+ if ($_[$i] > $maxval){
+ $maxval = $_[$i];
+ $maxpos = $i;
+ }
+ if ($_[$i] == $maxval){
+ $maxval = $_[$i];
+ if (scalar(@maxposes) == 0){
+ push @maxposes, $i;
+ push @maxvals, $_[$i];
+
+ }
+ elsif ($maxvals[0] == $maxval){
+ push @maxposes, $i;
+ push @maxvals, $_[$i];
+ }
+ else{
+ @maxposes = (); @maxvals = ();
+ push @maxposes, $i;
+ push @maxvals, $_[$i];
+ }
+
+ }
+
+ }
+ return $maxpos if scalar(@maxposes) < 2;
+ return (@maxposes);
+}
+
+sub array_smallest_number_arrayPosition{
+ return "NA" if scalar(@_) == 0;
+ return 0 if scalar(@_) == 1;
+ my $minpos = 0;
+ my @minposes = ();
+ my @minvals = ();
+ my $minval = array_largest_number(@_);
+ my $maxval = array_smallest_number(@_);
+ #print "starting with $maxval, ending with $minval\n";
+ for my $i (0 ... $#_){
+ if ($_[$i] < $minval){
+ $minval = $_[$i];
+ $minpos = $i;
+ }
+ if ($_[$i] == $minval){
+ $minval = $_[$i];
+ if (scalar(@minposes) == 0){
+ push @minposes, $i;
+ push @minvals, $_[$i];
+
+ }
+ elsif ($minvals[0] == $minval){
+ push @minposes, $i;
+ push @minvals, $_[$i];
+ }
+ else{
+ @minposes = (); @minvals = ();
+ push @minposes, $i;
+ push @minvals, $_[$i];
+ }
+
+ }
+
+ }
+ #print "minposes=@minposes\n";
+
+ return $minpos if scalar(@minposes) < 2;
+ return (@minposes);
+}
+
+sub basic_stats{
+ my @arr = @_;
+# print " array_smallest_number= ", array_smallest_number(@arr)," array_largest_number= ", array_largest_number(@arr), " array_mean= ",array_mean(@arr),"\n";
+ return ":";
+}
+#xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx
+
+sub maftoAxt_multispecies {
+ my $printer = 0;
+# print "in maftoAxt_multispecies : got @_\n";
+ my $fname=$_[0];
+ open(IN,"<$_[0]") or die "Cannot open $_[0]: $! \n";
+ my $treedefinition = $_[1];
+ open(OUT,">$_[2]") or die "Cannot open $_[2]: $! \n";
+ my $counter = 0;
+ my $exactspeciesset = $_[3];
+ my @exactspeciesset_unarranged = split(/,/,$exactspeciesset);
+
+ $treedefinition=~s/[\)\(, ]/\t/g;
+ my @species=split(/\t+/,$treedefinition);
+ my @exactspecies=();
+
+ foreach my $spec (@species){
+ foreach my $espec (@exactspeciesset_unarranged){
+ push @exactspecies, $spec if $spec eq $espec;
+ }
+ }
+# print "exactspecies=@exactspecies\n";
+
+ ###########
+ my $select = 2;
+ #select = 1 if all species need sequences to be present for each block otherwise, it is 0
+ #select = 2 only the allowed set make up the alignment. use the removeset
+ # information to detect alignmenets that have other important genomes aligned.
+ ###########
+ my @allowedset = ();
+ @allowedset = split(/;/,allowedSetOfSpecies(join("_",@species))) if $select == 0;
+ @allowedset = join("_",0,@species) if $select == 1;
+ #print "species = @species , allowedset =",join("\n", @allowedset) ," \n";
+ @allowedset = join("_",0,@exactspecies) if $select == 2;
+ #print "allowedset = @allowedset and exactspecies = @exactspecies\n";
+
+ my $start = 0;
+ my @sequences = ();
+ my @titles = ();
+ my $species_counter = "0";
+ my $countermatch = 0;
+ my $outsideSpecies=0;
+
+ while(my $line = ){
+ next if $line =~ /^#/;
+ next if $line =~ /^i/;
+ chomp $line;
+ #print "$line";
+ my @fields = split(/\s+/,$line);
+ chomp $line;
+ if ($line =~ /^a /){
+ $start = 1;
+ }
+
+ if ($line =~ /^s /){
+ # print "fields1 = $fields[1] , start = $start\n";
+
+ foreach my $sp (@species){
+ if ($fields[1] =~ /$sp/){
+ $species_counter = $species_counter."_".$sp;
+ push(@sequences, $fields[6]);
+ my @sp_info = split(/\./,$fields[1]);
+ my $title = join(" ",@sp_info, $fields[2], ($fields[2]+$fields[3]), $fields[4]);
+ push(@titles, $title);
+
+ }
+ }
+ }
+
+ if (($line !~ /^a/) && ($line !~ /^s/) && ($line !~ /^#/) && ($line !~ /^i/) && ($start = 1)){
+
+ my $arranged = reorderSpecies($species_counter, @species);
+ my $stopper = 1;
+ my $arrno = 0;
+ foreach my $set (@allowedset){
+ if ($arranged eq $set){
+ # print "$arranged == $set\n";
+ $stopper = 0; last;
+ }
+ $arrno++;
+ }
+
+ if ($stopper == 0) {
+ # print " accepted\n";
+ @titles = split ";", orderInfo(join(";", @titles), $species_counter, $arranged) if $species_counter ne $arranged;
+
+ @sequences = split ";", orderInfo(join(";", @sequences), $species_counter, $arranged) if $species_counter ne $arranged;
+ my $filteredseq = filter_gaps(@sequences);
+
+ if ($filteredseq ne "SHORT"){
+ $counter++;
+ print OUT join (" ",$counter, @titles), "\n";
+ print OUT $filteredseq, "\n";
+ print OUT "\n";
+ $countermatch++;
+ }
+ # my @filtered_seq = split(/\t/,filter_gaps(@sequences) );
+ }
+ else{#print "\n";
+ }
+
+ @sequences = (); @titles = (); $start = 0;$species_counter = "0";
+ next;
+
+ }
+ }
+# print "countermatch = $countermatch\n";
+}
+
+sub reorderSpecies{
+ my @inarr=@_;
+ my $currSpecies = shift (@inarr);
+ my $ordered_species = 0;
+ my @species=@inarr;
+ foreach my $order (@species){
+ $ordered_species = $ordered_species."_".$order if $currSpecies=~ /$order/;
+ }
+ return $ordered_species;
+
+}
+
+sub filter_gaps{
+ my @sequences = @_;
+# print "sequences sent are @sequences\n";
+ my $seq_length = length($sequences[0]);
+ my $seq_no = scalar(@sequences);
+ my $allgaps = ();
+ for (1 ... $seq_no){
+ $allgaps = $allgaps."-";
+ }
+
+ my @seq_array = ();
+ my $seq_counter = 0;
+ foreach my $seq (@sequences){
+# my @sequence = split(/\s*/,$seq);
+ $seq_array[$seq_counter] = [split(/\s*/,$seq)];
+# push @seq_array, [@sequence];
+ $seq_counter++;
+ }
+ my $g = 0;
+ while ( $g < $seq_length){
+ last if (!exists $seq_array[0][$g]);
+ my $bases = ();
+ for my $u (0 ... $#seq_array){
+ $bases = $bases.$seq_array[$u][$g];
+ }
+# print $bases, "\n";
+ if ($bases eq $allgaps){
+# print "bases are $bases, position is $g \n";
+ for my $seq (@seq_array){
+ splice(@$seq , $g, 1);
+ }
+ }
+ else {
+ $g++;
+ }
+ }
+
+ my @outs = ();
+
+ foreach my $seq (@seq_array){
+ push(@outs, join("",@$seq));
+ }
+ return "SHORT" if length($outs[0]) <=100;
+ return (join("\n", @outs));
+}
+
+
+sub allowedSetOfSpecies{
+ my @allowed_species = split(/_/,$_[0]);
+ unshift @allowed_species, 0;
+# print "allowed set = @allowed_species \n";
+ my @output = ();
+ for (0 ... scalar(@allowed_species) - 4){
+ push(@output, join("_",@allowed_species));
+ pop @allowed_species;
+ }
+ return join(";",reverse(@output));
+
+}
+
+
+sub orderInfo{
+ my @info = split(/;/,$_[0]);
+# print "info = @info";
+ my @old = split(/_/,$_[1]);
+ my @new = split(/_/,$_[2]);
+ shift @old; shift @new;
+ my @outinfo = ();
+ foreach my $spe (@new){
+ for my $no (0 ... $#old){
+ if ($spe eq $old[$no]){
+ push(@outinfo, $info[$no]);
+ }
+ }
+ }
+# print "outinfo = @outinfo \n";
+ return join(";", @outinfo);
+}
+
+#xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx
+
+sub printarr {
+ print ">::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\n";
+ foreach my $line (@_) {print "$line\n";}
+ print "::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::<\n";
+}
diff --git a/tools/regVariation/microsatellite_birthdeath.xml b/tools/regVariation/microsatellite_birthdeath.xml
new file mode 100644
index 00000000000..87f2dc12d76
--- /dev/null
+++ b/tools/regVariation/microsatellite_birthdeath.xml
@@ -0,0 +1,64 @@
+
+ and causal mutational mechanisms from previously identified orthologous microsatellite sets
+
+ microsatellite_birthdeath.pl
+ $alignment
+ $orthfile
+ $outfile
+ ${alignment.metadata.species}
+ "$tree_definition"
+ $thresholds
+ $separation
+ $simthresh
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+.. class:: infomark
+
+**What it does**
+
+This tool uses raw orthologous microsatellite clusters (identified by the tool "Extract orthologous microsatellites") to identify microsatellite births and deaths along individual lineages of a phylogenetic tree.
+
+
+
+
+
diff --git a/tools/regVariation/multispecies_MicrosatDataGenerator_interrupted_GALAXY.xml b/tools/regVariation/multispecies_MicrosatDataGenerator_interrupted_GALAXY.xml
index 1053340882c..1f8b0b7e192 100644
--- a/tools/regVariation/multispecies_MicrosatDataGenerator_interrupted_GALAXY.xml
+++ b/tools/regVariation/multispecies_MicrosatDataGenerator_interrupted_GALAXY.xml
@@ -51,13 +51,6 @@
This tool finds ortholgous microsatellite blocks between aligned species
------
-
-.. class:: warningmark
-
-**Note**
-
-NA