|
0
|
1 #!/usr/bin/perl -w
|
|
|
2 use strict;
|
|
|
3 use warnings;
|
|
|
4 use Term::ANSIColor;
|
|
|
5 use Pod::Checker;
|
|
|
6 use File::Basename;
|
|
|
7 use IO::Handle;
|
|
|
8 use Cwd;
|
|
|
9 use File::Path qw(make_path remove_tree);
|
|
|
10 use File::Temp qw/ tempfile tempdir /;
|
|
|
11 my $tdir = tempdir( CLEANUP => 1 );
|
|
|
12 chdir $tdir;
|
|
|
13 my $dir = getcwd;
|
|
|
14 # print "current dit=$dir\n";
|
|
|
15 #<STDIN>;
|
|
|
16 use vars qw (%treesToReject %template $printer $interr_poscord $interrcord $no_of_interruptionscord $stringfile @tags
|
|
|
17 $infocord $typecord $startcord $strandcord $endcord $microsatcord $motifcord $sequencepos $no_of_species
|
|
|
18 $gapcord %thresholdhash $tree_decipherer @sp_ident %revHash %sameHash %treesToIgnore %alternate @exactspecies_orig @exactspecies @exacttags
|
|
|
19 $mono_flanksimplicityRepno $di_flanksimplicityRepno $prop_of_seq_allowedtoAT $prop_of_seq_allowedtoCG);
|
|
|
20 use FileHandle;
|
|
|
21 use IO::Handle; # 5.004 or higher
|
|
|
22
|
|
|
23
|
|
|
24
|
|
|
25 #my @ar = ("/Users/ydk/work/rhesus_microsat/results/galay/chr22_5sp.maf.txt", "/Users/ydk/work/rhesus_microsat/results/galay/dataset_11.dat",
|
|
|
26 #"/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",
|
|
|
27 #"10","0.8");
|
|
|
28 my @ar = @ARGV;
|
|
|
29 my ($maf, $orth, $summout, $species_set, $tree_definition, $thresholds, $FLANK_SUPPORT, $SIMILARITY_THRESH) = @ar;
|
|
|
30 $SIMILARITY_THRESH=$SIMILARITY_THRESH/100;
|
|
|
31 #########################
|
|
|
32 $SIMILARITY_THRESH = $SIMILARITY_THRESH/100;
|
|
|
33 my $EDGE_DISTANCE = 10;
|
|
|
34 my $COMPLEXITY_SUPPORT = 20;
|
|
|
35 load_thresholds("9_10_12_12");
|
|
|
36 my $FLANKINDEL_MAXTHRESH = 0.3;
|
|
|
37
|
|
|
38 my $mono_flanksimplicityRepno=6;
|
|
|
39 my $di_flanksimplicityRepno=10;
|
|
|
40 my $prop_of_seq_allowedtoAT=0.5;
|
|
|
41 my $prop_of_seq_allowedtoCG=0.66;
|
|
|
42
|
|
|
43 #########################
|
|
|
44 my $tspecies_set = $species_set;
|
|
|
45
|
|
|
46 my %speciesReplacement = ();
|
|
|
47 my %speciesReplacementTag = ();
|
|
|
48 my %replacementArr= ();
|
|
|
49 my %replacementArrTag= ();
|
|
|
50 my %backReplacementArr= ();
|
|
|
51 my %backReplacementArrTag= ();
|
|
|
52 $tree_definition=~s/\s+//g;
|
|
|
53
|
|
|
54 my $tree_definition_split = $tree_definition;
|
|
|
55 $tree_definition_split =~ s/[\(\)]//g;
|
|
|
56 my @gotSpecies = ($tree_definition_split =~ /(,)/g);
|
|
|
57 # print "gotSpecies = @gotSpecies\n";
|
|
|
58
|
|
|
59 if (scalar(@gotSpecies)+1 ==5){
|
|
|
60
|
|
|
61 $speciesReplacement{1}="calJac1";
|
|
|
62 $speciesReplacement{2}="rheMac2";
|
|
|
63 $speciesReplacement{3}="ponAbe2";
|
|
|
64 $speciesReplacement{4}="panTro2";
|
|
|
65 $speciesReplacement{5}="hg18";
|
|
|
66
|
|
|
67 $speciesReplacementTag{1}="M";
|
|
|
68 $speciesReplacementTag{2}="R";
|
|
|
69 $speciesReplacementTag{3}="O";
|
|
|
70 $speciesReplacementTag{4}="C";
|
|
|
71 $speciesReplacementTag{5}="H";
|
|
|
72 $species_set="hg18,panTro2,ponAbe2,rheMac2,calJac1";
|
|
|
73
|
|
|
74 }
|
|
|
75 if (scalar(@gotSpecies)+1 ==4){
|
|
|
76
|
|
|
77 $speciesReplacement{1}="rheMac2";
|
|
|
78 $speciesReplacement{2}="ponAbe2";
|
|
|
79 $speciesReplacement{3}="panTro2";
|
|
|
80 $speciesReplacement{4}="hg18";
|
|
|
81
|
|
|
82 $speciesReplacementTag{1}="R";
|
|
|
83 $speciesReplacementTag{2}="O";
|
|
|
84 $speciesReplacementTag{3}="C";
|
|
|
85 $speciesReplacementTag{4}="H";
|
|
|
86 $species_set="hg18,panTro2,ponAbe2,rheMac2";
|
|
|
87
|
|
|
88 }
|
|
|
89
|
|
|
90 # $tree_definition = "((((hg18,panTro2),ponAbe2),rheMac2),calJac1)";
|
|
|
91 my $tree_definition_copy = $tree_definition;
|
|
|
92 my $tree_definition_orig = $tree_definition;
|
|
|
93 my $brackets = 0;
|
|
|
94
|
|
|
95 while (1){
|
|
|
96 #last if $tree_definition_copy !~ /\(/;
|
|
|
97 $brackets++;
|
|
|
98 # print "brackets = $brackets\n";
|
|
|
99 last if $brackets > 6;
|
|
|
100 $tree_definition_copy =~ s/^\(//g;
|
|
|
101 $tree_definition_copy =~ s/\)$//g;
|
|
|
102 # print "tree_definition_copy = $tree_definition_copy\n";
|
|
|
103 my @arr = ();
|
|
|
104
|
|
|
105 if ($tree_definition_copy =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_\(\),]+)\)$/){
|
|
|
106 @arr = $tree_definition_copy =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_\(\),]+)$/;
|
|
|
107 # print "arr = @arr\n";
|
|
|
108 $tree_definition_copy = $2;
|
|
|
109 $replacementArr{$1} = $speciesReplacement{$brackets};
|
|
|
110 $backReplacementArr{$speciesReplacement{$brackets}}=$1;
|
|
|
111
|
|
|
112 $replacementArrTag{$1} = $speciesReplacementTag{$brackets};
|
|
|
113 $backReplacementArrTag{$speciesReplacementTag{$brackets}}=$1;
|
|
|
114 # print "replacing $1 with $replacementArr{$1}\n";
|
|
|
115
|
|
|
116 $sp_ident[$brackets-1] = $1;
|
|
|
117
|
|
|
118 }
|
|
|
119 elsif ($tree_definition_copy =~ /^\(([a-zA-Z0-9_\(\),]+),([a-zA-Z0-9_]+)$/){
|
|
|
120 @arr = $tree_definition_copy =~ /^([a-zA-Z0-9_\(\),]+),([a-zA-Z0-9_]+)$/;
|
|
|
121 # print "arr = @arr\n";
|
|
|
122 $tree_definition_copy = $1;
|
|
|
123 $replacementArr{$2} = $speciesReplacement{$brackets};
|
|
|
124 $backReplacementArr{$speciesReplacement{$brackets}}=$2;
|
|
|
125
|
|
|
126 $replacementArrTag{$2} = $speciesReplacementTag{$brackets};
|
|
|
127 $backReplacementArrTag{$speciesReplacementTag{$brackets}}=$2;
|
|
|
128 # print "replacing $2 with $replacementArr{$2}\n";
|
|
|
129
|
|
|
130 $sp_ident[$brackets-1] = $2;
|
|
|
131 }
|
|
|
132 elsif ($tree_definition_copy =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_]+)$/){
|
|
|
133 @arr = $tree_definition_copy =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_]+)$/;
|
|
|
134 # print "arr = @arr .. TERMINAL\n";
|
|
|
135 $tree_definition_copy = $1;
|
|
|
136 $replacementArr{$2} = $speciesReplacement{$brackets};
|
|
|
137 $replacementArr{$1} = $speciesReplacement{$brackets+1};
|
|
|
138 $backReplacementArr{$speciesReplacement{$brackets}}=$2;
|
|
|
139 $backReplacementArr{$speciesReplacement{$brackets+1}}=$1;
|
|
|
140
|
|
|
141 $replacementArrTag{$1} = $speciesReplacementTag{$brackets+1};
|
|
|
142 $backReplacementArrTag{$speciesReplacementTag{$brackets+1}}=$1;
|
|
|
143
|
|
|
144 $replacementArrTag{$2} = $speciesReplacementTag{$brackets};
|
|
|
145 $backReplacementArrTag{$speciesReplacementTag{$brackets}}=$2;
|
|
|
146 # print "replacing $1 with $replacementArr{$1}\n";
|
|
|
147 # print "replacing $2 with $replacementArr{$2}\n";
|
|
|
148 # print "replacing $1 with $replacementArrTag{$1}\n";
|
|
|
149 # print "replacing $2 with $replacementArrTag{$2}\n";
|
|
|
150
|
|
|
151 $sp_ident[$brackets-1] = $2;
|
|
|
152 $sp_ident[$brackets] = $1;
|
|
|
153
|
|
|
154
|
|
|
155 last;
|
|
|
156
|
|
|
157 }
|
|
|
158 elsif ($tree_definition_copy =~ /^\(([a-zA-Z0-9_\(\),]+),([a-zA-Z0-9_\(\),]+)\)$/){
|
|
|
159 $tree_definition_copy =~ s/^\(//g;
|
|
|
160 $tree_definition_copy =~ s/\)$//g;
|
|
|
161 $brackets--;
|
|
|
162 }
|
|
|
163 }
|
|
|
164
|
|
|
165 foreach my $key (keys %replacementArr){
|
|
|
166 my $replacement = $replacementArr{$key};
|
|
|
167 $tree_definition =~ s/$key/$replacement/g;
|
|
|
168 }
|
|
|
169 @sp_ident = reverse(@sp_ident);
|
|
|
170 # print "modified tree_definition = $tree_definition\n";
|
|
|
171 # print "done .. tree_definition = $tree_definition\n";
|
|
|
172 # print "sp_ident = @sp_ident\n";
|
|
|
173 #<STDIN>;
|
|
|
174
|
|
|
175
|
|
|
176 my $complexity=int($COMPLEXITY_SUPPORT * (1/40));
|
|
|
177
|
|
|
178 #print "complexity=$complexity\n";
|
|
|
179 #<STDIN>;
|
|
|
180
|
|
|
181 $printer = 1;
|
|
|
182
|
|
|
183 my $rando = int(rand(1000));
|
|
|
184 my $localdate = `date`;
|
|
|
185 $localdate =~ /([0-9]+):([0-9]+):([0-9]+)/;
|
|
|
186 my $info = $rando.$1.$2.$3;
|
|
|
187
|
|
|
188 #---------------------------------------------------------------------------
|
|
|
189 # GETTING INPUT INFORMATION AND OPENING INPUT AND OUTPUT FILES
|
|
|
190
|
|
|
191
|
|
|
192 my @thresharr = (0, split(/,/,$thresholds));
|
|
|
193 my $randno=int(rand(100000));
|
|
|
194 my $megamatch = $randno.".megamatch.net.axt"; #"/gpfs/home/ydk104/work/rhesus_microsat/axtNet/hg18.panTro2.ponAbe2.rheMac2.calJac1/chr1.hg18.panTro2.ponAbe2.rheMac2.calJac1.net.axt";
|
|
|
195 my $megamatchlck = $megamatch.".lck";
|
|
|
196 unlink $megamatchlck;
|
|
|
197
|
|
|
198 my $selected= $orth;
|
|
|
199 #my $eventfile = $orth;
|
|
|
200 $selected = $selected."_SELECTED";
|
|
|
201 #$selected = $selected."_".$SIMILARITY_THRESH;
|
|
|
202 #my $runtime = $selected.".runtime";
|
|
|
203
|
|
|
204 my $inputtags = "H:C:O:R:M";
|
|
|
205 $inputtags = $ARGV[3] if exists $ARGV[3] && $ARGV[3] =~ /[A-Z]:[A-Z]/;
|
|
|
206
|
|
|
207 my @all_tags = split(/:/, $inputtags);
|
|
|
208 my $inputsp = "hg18:panTro2:ponAbe2:rheMac2:calJac1";
|
|
|
209 $inputsp = $ARGV[4] if exists $ARGV[4] && $ARGV[3] =~ /[0-9]+:/;
|
|
|
210 #@sp_ident = split(/:/,$inputsp);
|
|
|
211 my $junkfile = $orth."_junk";
|
|
|
212
|
|
|
213 my $sh = load_sameHash(1);
|
|
|
214 my $rh = load_revHash(1);
|
|
|
215
|
|
|
216 #print "inputs are : \n"; foreach(@ARGV){print $_,"\n";}
|
|
|
217 #open (SELECT, ">$selected") or die "Cannot open selected file: $selected: $!";
|
|
|
218 open (SUMMARY, ">$summout") or die "Cannot open summout file: $summout: $!";
|
|
|
219 #open (RUN, ">$runtime") or die "Cannot open orth file: $runtime: $!";
|
|
|
220 #my $ctlfile = "baseml\.ctl"; #$ARGV[4];
|
|
|
221 #my $treefile = "/gpfs/home/ydk104/work/rhesus_microsat/codes/lib/"; #1 THIS IS THE THE TREE UNDER CONSIDERATION, IN NEWICK
|
|
|
222 my %registeredTrees = ();
|
|
|
223 my @removalReasons =
|
|
|
224 ("microsatellite is compound",
|
|
|
225 "complex structure",
|
|
|
226 "if no. if micros is more than no. of species",
|
|
|
227 "if more than one micro per species ",
|
|
|
228 "if microsat contains N",
|
|
|
229 "different motif than required ",
|
|
|
230 "more than zero interruptions",
|
|
|
231 "microsat could not form key ",
|
|
|
232 "orthologous microsats of different motif size ",
|
|
|
233 "orthologous microsats of different motifs ",
|
|
|
234 "microsats belong to different alignment blocks altogether",
|
|
|
235 "microsat near edge",
|
|
|
236 "microsat in low complexity region",
|
|
|
237 "microsat flanks dont align well",
|
|
|
238 "phylogeny not informative");
|
|
|
239 my %allowedhash=();
|
|
|
240 #---------------------------------------------------------------------------
|
|
|
241 # WORKING ON MAKING THE MEGAMATCH FILE
|
|
|
242 my $chromt=int(rand(10000));
|
|
|
243 my $p_chr=$chromt;
|
|
|
244
|
|
|
245 my $tree_definition_orig_copy = $tree_definition_orig;
|
|
|
246
|
|
|
247 $tree_definition=~s/,/, /g;
|
|
|
248 $tree_definition =~ s/, +/, /g;
|
|
|
249 $tree_definition_orig=~s/,/, /g;
|
|
|
250 $tree_definition_orig =~ s/, +/, /g;
|
|
|
251 my @exactspeciesset_unarranged = split(/,/,$species_set);
|
|
|
252 my @exactspeciesset_unarranged_orig = split(/,/,$tspecies_set);
|
|
|
253 my $largesttree = "$tree_definition;";
|
|
|
254 my $largesttree_orig = "$tree_definition_orig;";
|
|
|
255 # print "largesttree = $largesttree\n";
|
|
|
256 $tree_definition =~ s/\(//g;
|
|
|
257 $tree_definition =~ s/\)//g;
|
|
|
258 $tree_definition=~s/[\)\(, ]/\t/g;
|
|
|
259 $tree_definition =~ s/\t+/\t/g;
|
|
|
260
|
|
|
261 $tree_definition_orig =~ s/\(//g;
|
|
|
262 $tree_definition_orig =~ s/\)//g;
|
|
|
263 $tree_definition_orig =~s/[\)\(, ]/\t/g;
|
|
|
264 $tree_definition_orig =~ s/\t+/\t/g;
|
|
|
265 # print "tree_definition = $tree_definition tree_definition_orig = $tree_definition_orig\n";
|
|
|
266
|
|
|
267 my @treespecies=split(/\t+/,$tree_definition);
|
|
|
268 my @treespecies_orig=split(/\t+/,$tree_definition_orig);
|
|
|
269 # print "tree_definition = $tree_definition .. treespecies=@treespecies ... treespecies_orig=@treespecies_orig\n";
|
|
|
270 #<STDIN>;
|
|
|
271
|
|
|
272 foreach my $spec (@treespecies){
|
|
|
273 foreach my $espec (@exactspeciesset_unarranged){
|
|
|
274 # print "spec=$spec and espec=$espec\n";
|
|
|
275 push @exactspecies, $spec if $spec eq $espec;
|
|
|
276 }
|
|
|
277 }
|
|
|
278
|
|
|
279 foreach my $spec (@treespecies_orig){
|
|
|
280 foreach my $espec (@exactspeciesset_unarranged_orig){
|
|
|
281 # print "spec=$spec and espec=$espec\n";
|
|
|
282 push @exactspecies_orig, $spec if $spec eq $espec;
|
|
|
283 }
|
|
|
284 }
|
|
|
285
|
|
|
286 my $focalspec = $exactspecies[0];
|
|
|
287 my $focalspec_orig = $exactspecies_orig[0];
|
|
|
288 # print "exactspecies=@exactspecies ... focalspec=$focalspec\n";
|
|
|
289 # print "texactspecies=@exactspecies_orig ... focalspec_orig=$focalspec_orig\n";
|
|
|
290 #<STDIN>;
|
|
|
291 my $arranged_species_set = join(".",@exactspecies);
|
|
|
292 my $arranged_species_set_orig = join(".",@exactspecies_orig);
|
|
|
293
|
|
|
294 @exacttags=@exactspecies;
|
|
|
295 my @exacttags_orig=@exactspecies_orig;
|
|
|
296
|
|
|
297 foreach my $extag (@exacttags){
|
|
|
298 $extag =~ s/hg18/H/g;
|
|
|
299 $extag =~ s/panTro2/C/g;
|
|
|
300 $extag =~ s/ponAbe2/O/g;
|
|
|
301 $extag =~ s/rheMac2/R/g;
|
|
|
302 $extag =~ s/calJac1/M/g;
|
|
|
303 }
|
|
|
304
|
|
|
305 foreach my $extag (@exacttags_orig){
|
|
|
306 $extag =~ s/hg18/H/g;
|
|
|
307 $extag =~ s/panTro2/C/g;
|
|
|
308 $extag =~ s/ponAbe2/O/g;
|
|
|
309 $extag =~ s/rheMac2/R/g;
|
|
|
310 $extag =~ s/calJac1/M/g;
|
|
|
311 }
|
|
|
312
|
|
|
313 my $chr_name = join(".",("chr".$p_chr),$arranged_species_set, "net", "axt");
|
|
|
314 #print "sending to maftoAxt_multispecies: $maf, $tree_definition, $chr_name, $species_set .. focalspec=$focalspec \n";
|
|
|
315
|
|
|
316 maftoAxt_multispecies($maf, $tree_definition_orig_copy, $chr_name, $tspecies_set);
|
|
|
317 #print "made files\n";
|
|
|
318 my @filterseqfiles= ($chr_name);
|
|
|
319 $largesttree =~ s/hg18/H/g;
|
|
|
320 $largesttree =~ s/panTro2/C/g;
|
|
|
321 $largesttree =~ s/ponAbe2/O/g;
|
|
|
322 $largesttree =~ s/rheMac2/R/g;
|
|
|
323 $largesttree =~ s/calJac1/M/g;
|
|
|
324 #<STDIN>;
|
|
|
325 #---------------------------------------------------------------------------
|
|
|
326
|
|
|
327 my ($lagestnodes, $largestbranches) = get_nodes($largesttree);
|
|
|
328 shift (@$lagestnodes);
|
|
|
329 my @extendedtitle=();
|
|
|
330
|
|
|
331 my $title = ();
|
|
|
332 my $parttitle = ();
|
|
|
333 my @titlearr = ();
|
|
|
334 my @firsttitle=($focalspec_orig."chrom", $focalspec_orig."start", $focalspec_orig."end", $focalspec_orig."motif", $focalspec_orig."motifsize", $focalspec_orig."threshold");
|
|
|
335
|
|
|
336 my @finames= qw(chr start end motif motifsize microsat mutation mutation.position mutation.from mutation.to insertion.details deletion.details);
|
|
|
337
|
|
|
338 my @fititle=();
|
|
|
339
|
|
|
340 foreach my $spec (split(",",$tspecies_set)){
|
|
|
341 push @fititle, $spec;
|
|
|
342 foreach my $name (@finames){
|
|
|
343 push @fititle, $spec.".".$name;
|
|
|
344 }
|
|
|
345 }
|
|
|
346
|
|
|
347
|
|
|
348 my @othertitle=qw(somechr somestart somened event source);
|
|
|
349
|
|
|
350 my @fnames = ();
|
|
|
351 push @fnames, qw(insertions_num deletions_num motinsertions_num motinsertionsf_num motdeletions_num motdeletionsf_num noninsertions_num nondeletions_num) ;
|
|
|
352 push @fnames, qw(binsertions_num bdeletions_num bmotinsertions_num bmotinsertionsf_num bmotdeletions_num bmotdeletionsf_num bnoninsertions_num bnondeletions_num) ;
|
|
|
353 push @fnames, qw(dinsertions_num ddeletions_num dmotinsertions_num dmotinsertionsf_num dmotdeletions_num dmotdeletionsf_num dnoninsertions_num dnondeletions_num) ;
|
|
|
354 push @fnames, qw(ninsertions_num ndeletions_num nmotinsertions_num nmotinsertionsf_num nmotdeletions_num nmotdeletionsf_num nnoninsertions_num nnondeletions_num) ;
|
|
|
355 push @fnames, qw(substitutions_num bsubstitutions_num dsubstitutions_num nsubstitutions_num indels_num subs_num);
|
|
|
356
|
|
|
357 my @fullnames = ();
|
|
|
358 # print "revising\n";
|
|
|
359 # print "H = $backReplacementArrTag{H}\n";
|
|
|
360 # print "C = $backReplacementArrTag{C}\n";
|
|
|
361 # print "O = $backReplacementArrTag{O}\n";
|
|
|
362 # print "R = $backReplacementArrTag{R}\n";
|
|
|
363 # print "M = $backReplacementArrTag{M}\n";
|
|
|
364
|
|
|
365 foreach my $lnode (@$lagestnodes){
|
|
|
366 my @pair = @$lnode;
|
|
|
367 my @nodemutarr = ();
|
|
|
368 for my $p (@pair){
|
|
|
369 # print "p = $p\n";
|
|
|
370 $p =~ s/[\(\), ]+//g;
|
|
|
371 $p =~ s/([A-Z])/$1./g;
|
|
|
372 $p =~ s/\.$//g;
|
|
|
373
|
|
|
374 $p =~ s/H/$backReplacementArrTag{H}/g;
|
|
|
375 $p =~ s/C/$backReplacementArrTag{C}/g;
|
|
|
376 $p =~ s/O/$backReplacementArrTag{O}/g;
|
|
|
377 $p =~ s/R/$backReplacementArrTag{R}/g;
|
|
|
378 $p =~ s/M/$backReplacementArrTag{M}/g;
|
|
|
379 foreach my $n (@fnames) { push @fullnames, $p.".".$n;}
|
|
|
380 }
|
|
|
381 }
|
|
|
382
|
|
|
383 #print SUMMARY "#",join("\t", @firsttitle, @fititle, @othertitle);
|
|
|
384
|
|
|
385 #print SUMMARY "\t",join("\t", @fullnames);
|
|
|
386 my $header = join("\t",@firsttitle, @fititle, @othertitle, @fullnames, @fnames, "tree", "cleancase");
|
|
|
387 # print "header= $header\n";
|
|
|
388 #<STDIN>;
|
|
|
389
|
|
|
390 #print SUMMARY "\t",join("\t", @fnames);
|
|
|
391 #$title= $title."\t".join("\t", @fnames);
|
|
|
392
|
|
|
393 #print SUMMARY "\t","tree","\t", "cleancase", "\n";
|
|
|
394 #$title= $title."\t"."tree"."\t"."cleancase". "\n";
|
|
|
395
|
|
|
396 #print $title; #<STDIN>;
|
|
|
397
|
|
|
398 #print "all_tags = @all_tags\n";
|
|
|
399
|
|
|
400 for my $no (3 ... $#all_tags+1){
|
|
|
401 # print "no=$no\n"; #<STDIN>;
|
|
|
402 @tags = @all_tags[0 ... $no-1];
|
|
|
403 # print "all_tags=>@all_tags< , tags = >@tags<\n" if $printer == 1; #<STDIN>;
|
|
|
404 %template=();
|
|
|
405 my @nextcounter = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
|
|
|
406 #next if scalar(@tags) < 4;
|
|
|
407 # print "now doing tags = @tags, no = $no\n";
|
|
|
408 open (ORTH, "<$orth") or die "Cannot open orth file: $orth: $!";
|
|
|
409
|
|
|
410 # print SUMMARY join "\t", qw (species chr start end branch motif microsat mutation position from to insertion deletion);
|
|
|
411
|
|
|
412
|
|
|
413 ##################### T E M P O R A R Y #####################
|
|
|
414 my @finaltitle=();
|
|
|
415 my @singletitle = qw (species chr start end motif motifsize microsat strand microsatsize col10 col11 col12 col13);
|
|
|
416 my $endtitle = ();
|
|
|
417 foreach my $tag (@tags){
|
|
|
418 my @tempsingle = ();
|
|
|
419
|
|
|
420 foreach my $single (@singletitle){
|
|
|
421 push @tempsingle, $tag.$single;
|
|
|
422 }
|
|
|
423 @finaltitle = (@finaltitle, @tempsingle);
|
|
|
424 }
|
|
|
425
|
|
|
426 # print SUMMARY join("\t",@finaltitle),"\n";
|
|
|
427
|
|
|
428 #############################################################
|
|
|
429
|
|
|
430 #---------------------------------------------------------------------------
|
|
|
431 # GET THE TREE FROM TREE FILE
|
|
|
432 my $tree = ();
|
|
|
433 $tree = "((H, C), O)" if $no == 3;
|
|
|
434 $tree = "(((H, C), O), R)" if $no == 4;
|
|
|
435 $tree = "((((H, C), O), R), M)" if $no == 5;
|
|
|
436 # $tree=~s/;$//g;
|
|
|
437 # print "our tree = $tree\n";
|
|
|
438 #---------------------------------------------------------------------------
|
|
|
439 # LOADING HASH CONTAINING ALL POSSIBLE TREES:
|
|
|
440 $tree_decipherer = "/gpfs/home/ydk104/work/rhesus_microsat/codes/lib/tree_analysis_".join("",@tags).".txt";
|
|
|
441 %template=();
|
|
|
442 %alternate=();
|
|
|
443 load_allPossibleTrees($tree_decipherer, \%template, \%alternate);
|
|
|
444
|
|
|
445 #---------------------------------------------------------------------------
|
|
|
446 # LOADING THE TREES TO REJECT FOR BIRTH ANALYSIS
|
|
|
447 %treesToReject=();
|
|
|
448 %treesToIgnore=();
|
|
|
449 load_treesToReject(@tags);
|
|
|
450 load_treesToIgnore(@tags);
|
|
|
451 #---------------------------------------------------------------------------
|
|
|
452 # LOADING INPUT DATA INTO HASHES AND ARRAYS
|
|
|
453
|
|
|
454
|
|
|
455 #1 THIS IS THE POINT WHERE WE CAN FILTER OUT LARGE MICROSAT CLUSTERS
|
|
|
456 #2 AS WELL AS MULTIPLE-ALIGNMENT-BLOCKS-SPANNING MICROSATS (KIND OF
|
|
|
457 #3 IMPLICIT IN THE FIRST PART OF THE SENTENCE ITSELF IN MOST CASES).
|
|
|
458
|
|
|
459 my %orths=();
|
|
|
460 my $counterm = 0;
|
|
|
461 my $loaded = 0;
|
|
|
462 my %seen = ();
|
|
|
463 my @allowedchrs = ();
|
|
|
464 # print "no = $no\n"; #<STDIN>;
|
|
|
465
|
|
|
466 while (my $line = <ORTH>){
|
|
|
467 # print "line=$line\n";
|
|
|
468 my $register1 = $line =~ s/>$exactspecies_orig[0]/>$replacementArrTag{$exactspecies_orig[0]}/g;
|
|
|
469 my $register2 = $line =~ s/>$exactspecies_orig[1]/>$replacementArrTag{$exactspecies_orig[1]}/g;
|
|
|
470 my $register3 = $line =~ s/>$exactspecies_orig[2]/>$replacementArrTag{$exactspecies_orig[2]}/g;
|
|
|
471 my $register4 = $line =~ s/>$exactspecies_orig[3]/>$replacementArrTag{$exactspecies_orig[3]}/g;
|
|
|
472 my $register5 = $line =~ s/>$exactspecies_orig[4]/>$replacementArrTag{$exactspecies_orig[4]}/g if exists $exactspecies_orig[4];
|
|
|
473
|
|
|
474 # print "line = $line\n"; <STDIN>;
|
|
|
475
|
|
|
476
|
|
|
477 # next if $register1 + $register2 + $register3 + $register4 + $register5 > scalar(@tags);
|
|
|
478 my @micros = split(/>/,$line); # LOADING ALL THE MICROSAT ENTRIES FROM THE CLUSTER INTO @micros
|
|
|
479 #print "micros=",printarr(@micros),"\n"; #<STDIN>;
|
|
|
480 shift @micros; # EMPTYING THE FIRST, EMTPY ELEMENT OF THE ARRAY
|
|
|
481
|
|
|
482 $no_of_species = adjustCoordinates($micros[0]);
|
|
|
483 # print "A: $no_of_species\n";
|
|
|
484 next if $no_of_species != $no;
|
|
|
485 # print "no = $no ... no_of_species=$no_of_species\n";#<STDIN>;
|
|
|
486 $counterm++;
|
|
|
487 #------------------------------------------------
|
|
|
488 $nextcounter[0]++ if $line =~ /compound/;
|
|
|
489 next if $line =~ /compound/; # GETTING RID OF COMPOUND MICROSATS
|
|
|
490 #------------------------------------------------
|
|
|
491 #next if $line =~ /[A-Za-z]>[a-zA-Z]/;
|
|
|
492 #------------------------------------------------
|
|
|
493 chomp $line;
|
|
|
494 my $match_count = ($line =~ s/>/>/g); # COUNTING THE NUMBER OF MICROSAT ENTRIES IN THE CLUSTER
|
|
|
495 #print "number of species = $match_count\n";
|
|
|
496 my $stopper = 0;
|
|
|
497 foreach my $mic (@micros){
|
|
|
498 my @local = split(/\t/,$mic);
|
|
|
499 if ($local[$typecord] =~ /\./ || exists($local[$no_of_interruptionscord+2])) {$stopper = 1; $nextcounter[1]++;
|
|
|
500 last; }
|
|
|
501 # REMOVING CLUSTERS WITH THE CYRPTIC, (UNRESOLVABLY COMPLEX) MICROSAT ENTRIES IN THEM
|
|
|
502 }
|
|
|
503 next if $stopper ==1;
|
|
|
504 #------------------------------------------------
|
|
|
505 $nextcounter[2]++ if (scalar(@micros) >$no_of_species);
|
|
|
506
|
|
|
507 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.
|
|
|
508 #2 THIS IS SO BECAUSE SUCH CLUSTERS IMPLY THAT IN AT LEAST ONE SPECIES, THERE IS MORE THAN ONE MICROSAT ENTRY
|
|
|
509 #3 IN THE CLUSTER. THUS, HERE WE ARE GETTING RID OF MICROSATS CLUSTERS THAT INCLUDE MULTUPLE, NEIGHBORING
|
|
|
510 #4 MICROSATS, AND STICK TO CLEAN MICROSATS THAT DO NOT HAVE ANY MICROSATS IN NEIGHBORHOOD.
|
|
|
511 #5 THIS 'NEIGHBORHOOD-RANGE' HAD BEEN DECIDED PREVIOUSLY IN OUR CODE multiSpecies_orthFinder4.pl
|
|
|
512 my $nexter = 0;
|
|
|
513 foreach my $tag (@tags){
|
|
|
514 my $tagcount = ($line =~ s/>$tag\t/>$tag\t/g);
|
|
|
515 if ($tagcount > 1) { $nexter =1; #print colored ['red'],"multiple entires per species : $tagcount of $tag\n" if $printer == 1;
|
|
|
516 next;
|
|
|
517 }
|
|
|
518 }
|
|
|
519
|
|
|
520 if ($nexter == 1){
|
|
|
521 $nextcounter[3]++;
|
|
|
522 next;
|
|
|
523 }
|
|
|
524 #------------------------------------------------
|
|
|
525 foreach my $mic (@micros){ #1 REMOVING MICROSATELLITES WITH ANY 'N's IN THEM
|
|
|
526 my @local = split(/\t/,$mic);
|
|
|
527 if ($local[$microsatcord] =~ /N/) {$stopper =1; $nextcounter[4]++;
|
|
|
528 last;}
|
|
|
529 }
|
|
|
530 next if $stopper ==1;
|
|
|
531 #print "till here 1\n"; #<STDIN>;
|
|
|
532 #------------------------------------------------
|
|
|
533 my @micros_copy = @micros;
|
|
|
534
|
|
|
535 my $tempmicro = shift(@micros_copy); #1 CURRENTLY OBTAINING INFORMATION FOR THE FIRST
|
|
|
536 #2 MICROSAT IN THE CLUSTER.
|
|
|
537 my @tempfields = split(/\t/,$tempmicro);
|
|
|
538 my $prevtype = $tempfields[$typecord];
|
|
|
539 my $tempmotif = $tempfields[$motifcord];
|
|
|
540
|
|
|
541 my $tempfirstmotif = ();
|
|
|
542 if (scalar(@tempfields) > $microsatcord + 2){
|
|
|
543 if ($tempfields[$no_of_interruptionscord] >= 1) { #1 DISCARDING MICROSATS WITH MORE THAN ZERO INTERRUPTIONS
|
|
|
544 #2 IN THE FIRST MICROSAT OF THE CLUSTER
|
|
|
545 $nexter =1; #print colored ['blue'],"more than one interruptions \n" if $printer == 1;
|
|
|
546 }
|
|
|
547 }
|
|
|
548 if ($nexter == 1){
|
|
|
549 $nextcounter[6]++;
|
|
|
550 next;
|
|
|
551 } #1 DONE OBTAINING INFORMATION REGARDING
|
|
|
552 #2 THE FIRST MICROSAT FROM THE CLUSTER
|
|
|
553
|
|
|
554 if ($tempmotif =~ /^\[/){
|
|
|
555 $tempmotif =~ s/^\[//g;
|
|
|
556 $tempmotif =~ /([a-zA-Z]+)\].*/;
|
|
|
557 $tempfirstmotif = $1; #1 OBTAINING THE FIRTS MOTIF OF MICROSAT
|
|
|
558 }
|
|
|
559 else {$tempfirstmotif = $tempmotif;}
|
|
|
560 my $prevmotif = $tempfirstmotif;
|
|
|
561
|
|
|
562 my $key = ();
|
|
|
563 # print "searching temp micro for 0-9 $focalspec chr0-9a-zA-Z 0-9 0-9 \n";
|
|
|
564 # print "tempmicro = $tempmicro .. looking for ([0-9]+)\s+($focalspec_orig)\s(chr[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)\n"; <STDIN>;
|
|
|
565 if ($tempmicro =~ /([0-9]+)\s+($focalspec_orig)\s(chr[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)/ ) {
|
|
|
566 # print "B: $no_of_species\n";
|
|
|
567 $key = join("_",$2, $3, $4, $5);
|
|
|
568 }
|
|
|
569 else{
|
|
|
570 # print "counld not form a key for temp\n"; # if $printer == 1;
|
|
|
571 $nextcounter[7]++;
|
|
|
572 next;
|
|
|
573 }
|
|
|
574 #----------------- #1 NOW, AFTER OBTAINING INFORMATION ABOUT
|
|
|
575 #2 THE FIRST MICROSAT IN THE CLUSTER, THE
|
|
|
576 #3 FOLLOWING LOOP GOES THROUGH THE OTHER MICROSATS
|
|
|
577 #4 TO SEE IF THEY SHARE THE REQUIRED FEATURES (BELOW)
|
|
|
578
|
|
|
579 foreach my $micro (@micros_copy){
|
|
|
580 my @fields = split(/\t/,$micro);
|
|
|
581 #-----------------
|
|
|
582 if (scalar(@fields) > $microsatcord + 2){ #1 DISCARDING MICROSATS WITH MORE THAN ONE INTERRUPTIONS
|
|
|
583 if ($fields[$no_of_interruptionscord] >= 1) {$nexter =1; #print colored ['blue'],"more than one interruptions \n" if $printer == 1;
|
|
|
584 $nextcounter[6]++;
|
|
|
585 last; }
|
|
|
586 }
|
|
|
587 #-----------------
|
|
|
588 if (($prevtype ne "0") && ($prevtype ne $fields[$typecord])) {
|
|
|
589 $nexter =1; #print colored ['yellow'],"microsat of different type \n" if $printer == 1;
|
|
|
590 $nextcounter[8]++;
|
|
|
591 last; } #1 DISCARDING MICROSAT CLUSTERS WHERE MICROSATS BELONG
|
|
|
592 #----------------- #2 TO DIFFERENT TYPES (MONOS, DIS, TRIS ETC.)
|
|
|
593 $prevtype = $fields[$typecord];
|
|
|
594
|
|
|
595 my $motif = $fields[$motifcord];
|
|
|
596 my $firstmotif = ();
|
|
|
597
|
|
|
598 if ($motif =~ /^\[/){
|
|
|
599 $motif =~ s/^\[//g;
|
|
|
600 $motif =~ /([a-zA-Z]+)\].*/;
|
|
|
601 $firstmotif = $1;
|
|
|
602 }
|
|
|
603 else {$firstmotif = $motif;}
|
|
|
604
|
|
|
605 my $motifpattern = $firstmotif.$firstmotif;
|
|
|
606 my $prevmotifpattern = $prevmotif.$prevmotif;
|
|
|
607
|
|
|
608 if (($prevmotif ne "0")&&(($motifpattern !~ /$prevmotif/i)||($prevmotifpattern !~ /$firstmotif/i)) ) {
|
|
|
609 $nexter =1; #print colored ['green'],"different motifs used \n$line\n" if $printer == 1;
|
|
|
610 $nextcounter[9]++;
|
|
|
611 last;
|
|
|
612 } #1 DISCARDING MICROSAT CLUSTERS WHERE MICROSATS BELONG
|
|
|
613 #2 TO DIFFERENT MOTIFS
|
|
|
614 my $prevmotif = $firstmotif;
|
|
|
615 #-----------------
|
|
|
616
|
|
|
617 for my $t (0 ... $#tags){ #1 DISCARDING MICROSAT CLUSTERS WHERE MICROSAT ENTRIES BELONG
|
|
|
618 #2 DIFFERENT ALIGNMENT BLOCKS
|
|
|
619 if ($micro =~ /([0-9]+)\s+($focalspec_orig)\s([_0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)/ ) {
|
|
|
620 my $key2 = join("_",$2, $3, $4, $5);
|
|
|
621 # print "key = $key .. key2 = $key2\n"; #<STDIN>;
|
|
|
622 if ($key2 ne $key){
|
|
|
623 # print "microsats belong to diffferent alignment blocks altogether\n" if $printer == 1;
|
|
|
624 $nextcounter[10]++;
|
|
|
625 $nexter = 1; last;
|
|
|
626 }
|
|
|
627 }
|
|
|
628 else{
|
|
|
629 # print "counld not form a key for $line\n"; # if $printer == 1;
|
|
|
630 #<STDIN>;
|
|
|
631 $nexter = 1; last;
|
|
|
632 }
|
|
|
633 }
|
|
|
634 }
|
|
|
635 #print "D2: $no_of_species\n";
|
|
|
636
|
|
|
637 #####################
|
|
|
638 if ($nexter == 1){
|
|
|
639 # print "nexting\n"; # if $printer == 1;
|
|
|
640 next;
|
|
|
641 }
|
|
|
642 else{
|
|
|
643 # print "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n$key:\n$line\nvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv\n" if $printer == 1;
|
|
|
644 push (@{$orths{$key}},$line);
|
|
|
645 $loaded++;
|
|
|
646 if ($line =~ /($focalspec_orig)\s([_a-zA-Z0-9]+)\s([0-9]+)\s([0-9]+)/ ) {
|
|
|
647
|
|
|
648 # print "$line\n" if $printer == 1; #if $line =~ /Contig/;
|
|
|
649 # print "################ ################\n" if $printer == 1;
|
|
|
650 push @allowedchrs, $2 if !exists $allowedhash{$2};
|
|
|
651 $allowedhash{$2} = 1;
|
|
|
652 my $key = join("\t",$1, $2, $3, $4);
|
|
|
653 # print "C: $no_of_species .. key = $key\n";#<STDIN>;
|
|
|
654 # print "print the shit: $key\n" ; #if $printer == 1;
|
|
|
655 $seen{$key} = 1;
|
|
|
656 }
|
|
|
657 else { #print "Key could not be formed in SPUT for ($focalspec_orig) (chrom) ([0-9]+) ([0-9]+)\n";
|
|
|
658 }
|
|
|
659 }
|
|
|
660 }
|
|
|
661 close ORTH;
|
|
|
662 # print "now studying where we lost microsatellites: @nextcounter\n";
|
|
|
663 for my $reason (0 ... $#nextcounter){
|
|
|
664 #print $removalReasons[$reason]."\t".$nextcounter[$reason],"\n";
|
|
|
665 }
|
|
|
666 # print "\ntotal number of keys formed = ", scalar(keys %orths), " = \n";
|
|
|
667 # print "done filtering .. counterm = $counterm and loaded = $loaded\n";
|
|
|
668 #----------------------------------------------------------------------------------------------------------------
|
|
|
669 # NOW GENERATING THE ALIGNMENT FILE WITH RELELEVENT ALIGNMENTS STORED ONLY.
|
|
|
670 # print "adding files @filterseqfiles \n";
|
|
|
671 #<STDIN>;
|
|
|
672
|
|
|
673 while (1){
|
|
|
674 if (-e $megamatchlck){
|
|
|
675 # print "waiting to write into $megamatchlck\n";
|
|
|
676 sleep 10;
|
|
|
677 }
|
|
|
678 else{
|
|
|
679 open (MEGAMLCK, ">$megamatchlck") or die "Cannot open megamatchlck file $megamatchlck: $!";
|
|
|
680 open (MEGAM, ">$megamatch") or die "Cannot open megamatch file $megamatch: $!";
|
|
|
681 last;
|
|
|
682 }
|
|
|
683 }
|
|
|
684 foreach my $seqfile (@filterseqfiles){
|
|
|
685 my $fullpath = $seqfile;
|
|
|
686 open (MATCH, "<$fullpath") or die "Cannot open MATCH file $fullpath: $!";
|
|
|
687 my $matchlines = 0;
|
|
|
688
|
|
|
689 while (my $line = <MATCH>) {
|
|
|
690 #print "checking $line";
|
|
|
691 if ($line =~ /($focalspec_orig)\s([a-zA-Z0-9]+)\s([0-9]+)\s([0-9]+)/ ) {
|
|
|
692 my $key = join("\t",$1, $2, $3, $4);
|
|
|
693 # print "key = $key\n";
|
|
|
694 #print "------------------------------------------------------\n";
|
|
|
695 #print "asking $line\n";
|
|
|
696 if (exists $seen{$key}){
|
|
|
697
|
|
|
698 #print "seen $line \n"; <STDIN>;
|
|
|
699 while (1){
|
|
|
700 $matchlines++;
|
|
|
701 print MEGAM $line;
|
|
|
702 $line = <MATCH>;
|
|
|
703 print MEGAM "\n" if $line !~ /[0-9a-zA-Z]/;
|
|
|
704 last if $line !~/[0-9a-zA-Z]/;
|
|
|
705 }
|
|
|
706 }
|
|
|
707 else{
|
|
|
708 # print "not seen\n";
|
|
|
709 }
|
|
|
710 }
|
|
|
711 }
|
|
|
712 # print "matchlines = $matchlines\n";
|
|
|
713 close MATCH;
|
|
|
714 }
|
|
|
715 close MEGAMLCK;
|
|
|
716
|
|
|
717 unlink $megamatchlck;
|
|
|
718 close MEGAM;
|
|
|
719 undef %seen;
|
|
|
720 # print "done writitn to $megamatch\n";#<STDIN>;
|
|
|
721 #----------------------------------------------------------------------------------------------------------------
|
|
|
722 #<STDIN>;
|
|
|
723 #---------------------------------------------------------------------------
|
|
|
724 # NOW, AFTER FILTERING MANY MICROSATS, AND LOADING THE FILTERED ONES INTO
|
|
|
725 # THE HASH %orths , WE GO THROUGH THE ALIGNMENT FILE, AND STUDY THE
|
|
|
726 # FLANKING SEQUENCES OF ALL THESE MICROSATS, TO FILTER THEM FURTHER
|
|
|
727 #$printer = 1;
|
|
|
728
|
|
|
729 my $microreadcounter=0;
|
|
|
730 my $contigsentered=0;
|
|
|
731 my $contignotrightcounter=0;
|
|
|
732 my $keynotformedcounter=0;
|
|
|
733 my $keynotfoundcounter= 0;
|
|
|
734 my $dotcounter = 0;
|
|
|
735
|
|
|
736 # print "opening $megamatch\n";
|
|
|
737
|
|
|
738 open (BO, "<$megamatch") or die "Cannot open alignment file: $megamatch: $!";
|
|
|
739 # print "doing $megamatch\n " ;
|
|
|
740
|
|
|
741 #<STDIN>;
|
|
|
742
|
|
|
743 while (my $line = <BO>){
|
|
|
744 # print $line; #<STDIN>;
|
|
|
745 # print "." if $dotcounter % 100 ==0;
|
|
|
746 # print "\n" if $dotcounter % 5000 ==0;
|
|
|
747 # print "dotcounter = $dotcounter\n " if $printer == 1;
|
|
|
748 next if $line !~ /^[0-9]+/;
|
|
|
749 $dotcounter++;
|
|
|
750 # print colored ['green'], "~" x 60, "\n" if $printer == 1;
|
|
|
751 # print colored ['green'], $line;# if $printer == 1;
|
|
|
752 chomp $line;
|
|
|
753 my @fields2 = split(/\t/,$line);
|
|
|
754 my $key2 = ();
|
|
|
755 my $alignment_no = (); #1 TEMPORARY
|
|
|
756 if ($line =~ /([0-9]+)\s+($focalspec_orig)\s([_\-s0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)/ ) {
|
|
|
757 # $key2 = join("\t",$1, $2, $4, $5);
|
|
|
758 $key2 = join("_",$2, $3, $4, $5);
|
|
|
759 # print "key = $key2\n";
|
|
|
760
|
|
|
761 # print "key = $key2\n";
|
|
|
762 $alignment_no=$1;
|
|
|
763 }
|
|
|
764 else {print "seq line $line incompatible\n"; $keynotformedcounter++; next;}
|
|
|
765
|
|
|
766 $no_of_species = adjustCoordinates($line);
|
|
|
767
|
|
|
768 $contignotrightcounter++ if $no_of_species != $no;
|
|
|
769 # print "contignotrightcounter=$contignotrightcounter\n";
|
|
|
770 # print "no_of_species=$no_of_species\n";
|
|
|
771 # print "no=$no\n";
|
|
|
772
|
|
|
773 next if $no_of_species != $no;
|
|
|
774 # print "D: $no_of_species\n";
|
|
|
775 # print "E: $no_of_species\n";
|
|
|
776 #<STDIN>;
|
|
|
777 # print "key = $key2\n" if $printer == 1;
|
|
|
778 my @clusters = (); #1 EXTRACTING MICROSATS CORRESPONDING TO THIS
|
|
|
779 #2 ALIGNMENT BLOCK
|
|
|
780 if (exists($orths{$key2})){
|
|
|
781 @clusters = @{$orths{$key2}};
|
|
|
782 $contigsentered++;
|
|
|
783 delete $orths{$key2};
|
|
|
784 }
|
|
|
785 else{
|
|
|
786 # print "orth does not exist\n";
|
|
|
787 $keynotfoundcounter++;
|
|
|
788 next;
|
|
|
789 }
|
|
|
790
|
|
|
791 my %sequences=(); #1 WILL STORE SEQUENCES IN THE CURRENT ALIGNMENT BLOCK
|
|
|
792 my $humseq = ();
|
|
|
793 foreach my $tag (@tags){ #1 READING THE ALIGNMENT FILE AND CAPTURING SEQUENCES
|
|
|
794 my $seq = <BO>; #2 OF ALL SPECIES.
|
|
|
795 chomp $seq;
|
|
|
796 $sequences{$tag} = " ".$seq;
|
|
|
797 #print "sequences = $sequences{$tag}\n" if $printer == 1;
|
|
|
798 $humseq = $seq if $tag =~ /H/;
|
|
|
799 }
|
|
|
800
|
|
|
801
|
|
|
802 foreach my $cluster (@clusters){ #1 NOW, GOING THROUGH THE CLUSTER OF MICROSATS
|
|
|
803 # print "x" x 60, "\n" if $printer == 1;
|
|
|
804 # print colored ['red'],"cluster = $cluster\n";
|
|
|
805 $largesttree =~ s/hg18/H/g;
|
|
|
806 $largesttree =~ s/panTro2/C/g;
|
|
|
807 $largesttree =~ s/ponAbe2/O/g;
|
|
|
808 $largesttree =~ s/rheMac2/R/g;
|
|
|
809 $largesttree =~ s/calJac1/M/g;
|
|
|
810
|
|
|
811 $microreadcounter++;
|
|
|
812 my @micros = split(/>/,$cluster);
|
|
|
813
|
|
|
814
|
|
|
815
|
|
|
816
|
|
|
817
|
|
|
818
|
|
|
819 shift @micros;
|
|
|
820
|
|
|
821 my $edge_microsat=0; #1 THIS WILL HAVE VALUE "1" IF MICROSAT IS FOUND
|
|
|
822 #2 TO BE TOO CLOSE TO THE EDGES OF ALIGNMENT BLOCK
|
|
|
823
|
|
|
824 my @starts= (); my %start_hash=(); #1 STORES THE START AND END COORDINATES OF MICROSATELLITES
|
|
|
825 my @ends = (); my %end_hash=(); #2 SO THAT LATER, WE WILL BE ABLE TO FIND THE EXTREME
|
|
|
826 #3 COORDINATE VALUES OF THE ORTHOLOGOUS MIROSATELLITES.
|
|
|
827
|
|
|
828 my %microhash=();
|
|
|
829 my %microsathash=();
|
|
|
830 my %nonmicrosathash=();
|
|
|
831 my $motif=(); #1 BASIC MOTIF OF THE MICROSATELLITE.. THERE'S ONLY 1
|
|
|
832 # print "tags=@tags\n";
|
|
|
833 for my $i (0 ... $#tags){ #1 FINDING THE MICROSAT, AND THE ALIGNMENT SEQUENCE
|
|
|
834 #2 CORRESPONDING TO THE PARTICULAR SPECIES (AS PER
|
|
|
835 #3 THE VARIABLE $TAG;
|
|
|
836 my $tag = $tags[$i];
|
|
|
837 # print $seq;
|
|
|
838 my $locus="NULL"; #1 THIS WILL STORE THE MICROSAT OF THIS SPECIES.
|
|
|
839 #2 IF THERE IS NO MICROSAT, IT WILL REMAIN "NULL"
|
|
|
840
|
|
|
841 foreach my $micro (@micros){
|
|
|
842 # print "micro=$micro, tag=$tag\n";
|
|
|
843 if ($micro =~ /^$tag/){ #1 MICROSAT OF THIS SPECIES FOUND..
|
|
|
844 $locus = $micro;
|
|
|
845 my @fields = split(/\t/,$micro);
|
|
|
846 $motif = $fields[$motifcord];
|
|
|
847 $microsathash{$tag}=$fields[$microsatcord];
|
|
|
848 # print "fields=@fields, and startcord=$startcord = $fields[$startcord]\n";
|
|
|
849 push(@starts, $fields[$startcord]);
|
|
|
850 push(@ends, $fields[$endcord]);
|
|
|
851 $start_hash{$tag}=$fields[$startcord];
|
|
|
852 $end_hash{$tag}=$fields[$endcord];
|
|
|
853 last;
|
|
|
854 }
|
|
|
855 else{$microsathash{$tag}="NULL"}
|
|
|
856 }
|
|
|
857 $microhash{$tag}=$locus;
|
|
|
858
|
|
|
859 }
|
|
|
860
|
|
|
861
|
|
|
862
|
|
|
863 my $extreme_start = smallest_number(@starts); #1 THESE TWO ARE THE EXTREME COORDINATES OF THE
|
|
|
864 my $extreme_end = largest_number(@ends); #2 MICROSAT CLUSTER ACCROSS ALL THE SPECIES IN
|
|
|
865 #3 WHOM IT IS FOUND TO BE ORTHOLOGOUS.
|
|
|
866
|
|
|
867 # print "starts=@starts... ends=@ends\n";
|
|
|
868
|
|
|
869 my %up_flanks = (); #1 CONTAINS UPSTEAM FLANKING REGIONS FOR EACH SPECIES
|
|
|
870 my %down_flanks = (); #1 CONTAINS DOWNDTREAM FLANKING REGIONS FOR EACH SPECIES
|
|
|
871
|
|
|
872 my %up_largeflanks = ();
|
|
|
873 my %down_largeflanks = ();
|
|
|
874
|
|
|
875 my %locusandflanks = ();
|
|
|
876 my %locusandlargeflanks = ();
|
|
|
877
|
|
|
878 my %up_internal_flanks=(); #1 CONTAINS SEQUENCE BETWEEN THE $extreme_start and the
|
|
|
879 #2 ACTUAL START OF MICROSATELLITE IN THE SPECIES
|
|
|
880 my %down_internal_flanks=(); #1 CONTAINS SEQUENCE BETWEEN THE $extreme_end and the
|
|
|
881 #2 ACTUAL end OF MICROSATELLITE IN THE SPECIES
|
|
|
882
|
|
|
883 my %alignment=(); #1 CONTAINS ACTUAL ALIGNMENT SEQUENCE BETWEEN THE TWO
|
|
|
884 #2 EXTEME VALUES.
|
|
|
885
|
|
|
886 my %microsatstarts=(); #1 WITHIN EACH ALIGNMENT, IF THERE EXISTS A MICROSATELLITE
|
|
|
887 #2 THIS HASH CONTAINS THE START SITE OF THE MICROSATELLITE
|
|
|
888 #3 WIHIN THE ALIGNMENT
|
|
|
889 next if !defined $extreme_start;
|
|
|
890 next if !defined $extreme_end;
|
|
|
891 next if $extreme_start > length($sequences{$tags[0]});
|
|
|
892 next if $extreme_start < 0;
|
|
|
893 next if $extreme_end > length($sequences{$tags[0]});
|
|
|
894
|
|
|
895 for my $i (0 ... $#tags){ #1 NOW THAT WE HAVE GATHERED INFORMATION REGARDING
|
|
|
896 #2 SEQUENCE ALIGNMENT AND MICROSATELLITE COORDINATES
|
|
|
897 #3 AS WELL AS THE EXTREME COORDINATES OF THE
|
|
|
898 #4 MICROSAT CLUSTER, WE WILL PROCEED TO EXTRACT THE
|
|
|
899 #5 FLANKING SEQUENCE OF ALL ORGS, AND STUDY IT IN
|
|
|
900 #6 MORE DETAIL.
|
|
|
901 my $tag = $tags[$i];
|
|
|
902 # print "tag=$tag.. seqlength = ",length($sequences{$tag})," extreme_start=$extreme_start and extreme_end=$extreme_end\n";
|
|
|
903 my $upstream_gaps = (substr($sequences{$tag}, 0, $extreme_start) =~ s/\-/-/g); #1 NOW MEASURING THE NUMBER OF GAPS IN THE UPSTEAM
|
|
|
904 #2 AND DOWNSTREAM SEQUENCES OF THE MICROSATs IN THIS
|
|
|
905 #3 CLUSTER.
|
|
|
906 # print "seq length $tag = $sequences{$tag} = ",length($sequences{$tag})," extreme_end=$extreme_end\n" ;
|
|
|
907 my $downstream_gaps = (substr($sequences{$tag}, $extreme_end) =~ s/\-/-/g);
|
|
|
908 if (($extreme_start - $upstream_gaps )< $EDGE_DISTANCE || (length($sequences{$tag}) - $extreme_end - $downstream_gaps) < $EDGE_DISTANCE){
|
|
|
909 $edge_microsat=1;
|
|
|
910
|
|
|
911 last;
|
|
|
912 }
|
|
|
913 else{
|
|
|
914 $up_flanks{$tag} = substr($sequences{$tag}, $extreme_start - $FLANK_SUPPORT, $FLANK_SUPPORT);
|
|
|
915 $down_flanks{$tag} = substr($sequences{$tag}, $extreme_end+1, $FLANK_SUPPORT);
|
|
|
916
|
|
|
917 $up_largeflanks{$tag} = substr($sequences{$tag}, $extreme_start - $COMPLEXITY_SUPPORT, $COMPLEXITY_SUPPORT);
|
|
|
918 $down_largeflanks{$tag} = substr($sequences{$tag}, $extreme_end+1, $COMPLEXITY_SUPPORT);
|
|
|
919
|
|
|
920
|
|
|
921 $alignment{$tag} = substr($sequences{$tag}, $extreme_start, $extreme_end-$extreme_start+1);
|
|
|
922 $locusandflanks{$tag} = $up_flanks{$tag}."[".$alignment{$tag}."]".$down_flanks{$tag};
|
|
|
923 $locusandlargeflanks{$tag} = $up_largeflanks{$tag}."[".$alignment{$tag}."]".$down_largeflanks{$tag};
|
|
|
924
|
|
|
925 if ($microhash{$tag} ne "NULL"){
|
|
|
926 $up_internal_flanks{$tag} = substr($sequences{$tag}, $extreme_start , $start_hash{$tag}-$extreme_start);
|
|
|
927 $down_internal_flanks{$tag} = substr($sequences{$tag}, $end_hash{$tag} , $extreme_end-$end_hash{$tag});
|
|
|
928 $microsatstarts{$tag}=$start_hash{$tag}-$extreme_start;
|
|
|
929 # print "tag = $tag, internal flanks = $up_internal_flanks{$tag} and $down_internal_flanks{$tag} and start = $microsatstarts{$tag}\n" if $printer == 1;
|
|
|
930 }
|
|
|
931 else{
|
|
|
932 $nonmicrosathash{$tag}=substr($sequences{$tag}, $extreme_start, $extreme_end-$extreme_start+1);
|
|
|
933
|
|
|
934 }
|
|
|
935 # print "up flank for species $tag = $up_flanks{$tag} \ndown flank for species $tag = $down_flanks{$tag} \n" if $printer == 1;
|
|
|
936
|
|
|
937 }
|
|
|
938
|
|
|
939 }
|
|
|
940 $nextcounter[11]++ if $edge_microsat==1;
|
|
|
941 next if $edge_microsat==1;
|
|
|
942
|
|
|
943
|
|
|
944 my $low_complexity = 0; #1 VALUE WILL BE 1 IF ANY OF THE FLANKING REGIONS
|
|
|
945 #2 IS FOUND TO BE OF LOW COMPLEXITY, BY USING THE
|
|
|
946 #3 FUNCTION sub test_complexity
|
|
|
947
|
|
|
948
|
|
|
949 for my $i (0 ... $#tags){
|
|
|
950 # print "i = $tags[$i]\n" if $printer == 1;
|
|
|
951 if (test_complexity($up_largeflanks{$tags[$i]}, $COMPLEXITY_SUPPORT) eq "LOW" || test_complexity($down_largeflanks{$tags[$i]}, $COMPLEXITY_SUPPORT) eq "LOW"){
|
|
|
952 # 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;
|
|
|
953 $low_complexity =1; last;
|
|
|
954 }
|
|
|
955 }
|
|
|
956
|
|
|
957 $nextcounter[12]++ if $low_complexity==1;
|
|
|
958 next if $low_complexity == 1;
|
|
|
959
|
|
|
960
|
|
|
961 my $sequence_dissimilarity = 0; #1 THIS VALYE WILL BE 1 IF THE SEQUENCE SIMILARITY
|
|
|
962 #2 BETWEEN ANY OF THE SPECIES AGAINST THE HUMAN
|
|
|
963 #3 FLANKING SEQUENCES IS BELOW A CERTAIN THRESHOLD
|
|
|
964 #4 AS DESCRIBED IN FUNCTION sub sequence_similarity
|
|
|
965 my %donepair = ();
|
|
|
966 for my $i (0 ... $#tags){
|
|
|
967 # print "i = $tags[$i]\n" if $printer == 1;
|
|
|
968 # next if $i == 0;
|
|
|
969 # print colored ['magenta'],"THIS IS UP\n" if $printer == 1;
|
|
|
970
|
|
|
971 for my $b (0 ... $#tags){
|
|
|
972 next if $b == $i;
|
|
|
973 my $pair = ();
|
|
|
974 $pair = $i."_".$b if $i < $b;
|
|
|
975 $pair = $b."_".$i if $b < $i;
|
|
|
976 next if exists $donepair{$pair};
|
|
|
977 my ($up_similarity,$upnucdiffs, $upindeldiffs) = sequence_similarity($up_flanks{$tags[$i]}, $up_flanks{$tags[$b]}, $SIMILARITY_THRESH, $info);
|
|
|
978 my ($down_similarity,$downnucdiffs, $downindeldiffs) = sequence_similarity($down_flanks{$tags[$i]}, $down_flanks{$tags[$b]}, $SIMILARITY_THRESH, $info);
|
|
|
979 $donepair{$pair} = $up_similarity."_".$down_similarity;
|
|
|
980
|
|
|
981 # print RUN "$up_similarity $upnucdiffs $upindeldiffs $down_similarity $downnucdiffs $downindeldiffs\n";
|
|
|
982
|
|
|
983 if ( $up_similarity < $SIMILARITY_THRESH || $down_similarity < $SIMILARITY_THRESH){
|
|
|
984 $sequence_dissimilarity =1;
|
|
|
985 last;
|
|
|
986 }
|
|
|
987 }
|
|
|
988 }
|
|
|
989 $nextcounter[13]++ if $sequence_dissimilarity==1;
|
|
|
990
|
|
|
991 next if $sequence_dissimilarity == 1;
|
|
|
992 my ($simplified_microsat, $Hchrom, $Hstart, $Hend, $locusmotif, $locusmotifsize) = summarize_microsat($cluster, $humseq);
|
|
|
993 # print "simplified_microsat=$simplified_microsat\n";
|
|
|
994 my ($tree_analysis, $conformation) = treeStudy($simplified_microsat);
|
|
|
995 # print "tree_analysis = $tree_analysis .. conformation=$conformation\n";
|
|
|
996 #<STDIN>;
|
|
|
997
|
|
|
998 # print SELECT "\"$conformation\"\t$tree_analysis\n";
|
|
|
999
|
|
|
1000 next if $tree_analysis =~ /DISCARD/;
|
|
|
1001 if (exists $treesToReject{$tree_analysis}){
|
|
|
1002 $nextcounter[14]++;
|
|
|
1003 next;
|
|
|
1004 }
|
|
|
1005
|
|
|
1006 # print "F: $no_of_species\n";
|
|
|
1007
|
|
|
1008 # my $adjuster=();
|
|
|
1009 # if ($no_of_species == 4){
|
|
|
1010 # my @sields = split(/\t/,$simplified_microsat);
|
|
|
1011 # my $somend = pop(@sields);
|
|
|
1012 # my $somestart = pop(@sields);
|
|
|
1013 # my $somechr = pop(@sields);
|
|
|
1014 # $adjuster = "NA\t" x 13 ;
|
|
|
1015 # $simplified_microsat = join ("\t", @sields, $adjuster).$somechr."\t".$somestart."\t".$somend;
|
|
|
1016 # }
|
|
|
1017 # if ($no_of_species == 3){
|
|
|
1018 # my @sields = split(/\t/,$simplified_microsat);
|
|
|
1019 # my $somend = pop(@sields);
|
|
|
1020 # my $somestart = pop(@sields);
|
|
|
1021 # my $somechr = pop(@sields);
|
|
|
1022 # $adjuster = "NA\t" x 26 ;
|
|
|
1023 # $simplified_microsat = join ("\t", @sields, $adjuster).$somechr."\t".$somestart."\t".$somend;
|
|
|
1024 # }
|
|
|
1025 #
|
|
|
1026 $registeredTrees{$tree_analysis} = 1 if !exists $registeredTrees{$tree_analysis};
|
|
|
1027 $registeredTrees{$tree_analysis}++ if exists $registeredTrees{$tree_analysis};
|
|
|
1028
|
|
|
1029 if (exists $treesToIgnore{$tree_analysis}){
|
|
|
1030 my @appendarr = ();
|
|
|
1031
|
|
|
1032 # print SUMMARY $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize], "\t", $simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
|
|
|
1033 #print "SUMMARY ",$Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t", $thresharr[$locusmotifsize], "\t", $simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t";
|
|
|
1034 # print SELECT $Hchrom,"\t",$Hstart,"\t",$Hend,"\t","NOEVENT", "\t\t", $cluster,"\n";
|
|
|
1035
|
|
|
1036 foreach my $lnode (@$lagestnodes){
|
|
|
1037 my @pair = @$lnode;
|
|
|
1038 my @nodemutarr = ();
|
|
|
1039 for my $p (@pair){
|
|
|
1040 my @mutinfoarray1 = ();
|
|
|
1041 for (1 ... 38){
|
|
|
1042 # push (@mutinfoarray1, "NA")
|
|
|
1043 }
|
|
|
1044 # print SUMMARY join ("\t", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
|
|
|
1045 # print join ("\t", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
|
|
|
1046 }
|
|
|
1047
|
|
|
1048 }
|
|
|
1049 for (1 ... 38){
|
|
|
1050 push (@appendarr, "NA")
|
|
|
1051 }
|
|
|
1052 # print SUMMARY join ("\t", @appendarr,"NULL", "NULL"),"\n";
|
|
|
1053 # print join ("\t", @appendarr,"NULL", "NULL"),"\n";
|
|
|
1054 # print "SUMMARY ",join ("\t", @appendarr,"NULL", "NULL"),"\n"; #<STDIN>;
|
|
|
1055 next;
|
|
|
1056 }
|
|
|
1057 # print colored ['blue'],"cluster = $cluster\n";
|
|
|
1058
|
|
|
1059 my ($mutations_array, $nodes, $branches_hash, $alivehash, $primaryalignment) = peel_onion($tree, \%sequences, \%alignment, \@tags, \%microsathash, \%nonmicrosathash, $motif, $tree_analysis, $thresholdhash{length($motif)}, \%microsatstarts);
|
|
|
1060
|
|
|
1061 if ($mutations_array eq "NULL"){
|
|
|
1062 # print "cluster = $cluster \n"; <STDIN>;
|
|
|
1063 my @appendarr = ();
|
|
|
1064
|
|
|
1065 # print SUMMARY $Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize, "\t";
|
|
|
1066
|
|
|
1067 # foreach my $lnode (@$lagestnodes){
|
|
|
1068 # my @pair = @$lnode;
|
|
|
1069 # my @nodemutarr = ();
|
|
|
1070 # for my $p (@pair){
|
|
|
1071 # my @mutinfoarray1 = ();
|
|
|
1072 # for (1 ... 38){
|
|
|
1073 # push (@mutinfoarray1, "NA")
|
|
|
1074 # }
|
|
|
1075 # print SUMMARY join ("\t", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
|
|
|
1076 # print join ("\t", @mutinfoarray1[0...($#mutinfoarray1)] ),"\t";
|
|
|
1077 # }
|
|
|
1078 # }
|
|
|
1079 # for (1 ... 38){
|
|
|
1080 # push (@appendarr, "NA")
|
|
|
1081 # }
|
|
|
1082 # print SUMMARY join ("\t", @appendarr,"NULL", "NULL"),"\n";
|
|
|
1083 # print join ("\t", @appendarr,"NULL", "NULL"),"\n";
|
|
|
1084 # print join ("\t","SUMMARY", @appendarr,"NULL", "NULL"),"\n"; #<STDIN>;
|
|
|
1085 next;
|
|
|
1086 }
|
|
|
1087
|
|
|
1088
|
|
|
1089 # print "sent: \n" if $printer == 1;
|
|
|
1090 # print "nodes = @$nodes, branches array:\n" if $mutations_array ne "NULL" && $printer == 1;
|
|
|
1091
|
|
|
1092 my ($newstatus, $newmutations_array, $newnodes, $newbranches_hash, $newalivehash, $finalalignment) = fillAlignmentGaps($tree, \%sequences, \%alignment, \@tags, \%microsathash, \%nonmicrosathash, $motif, $tree_analysis, $thresholdhash{length($motif)}, \%microsatstarts);
|
|
|
1093 # print "newmutations_array returned = \n",join("\n",@$newmutations_array),"\n" if $newmutations_array ne "NULL" && $printer == 1;
|
|
|
1094 my @finalmutations_array= ();
|
|
|
1095 @finalmutations_array = selectMutationArray($mutations_array, $newmutations_array, \@tags, $alivehash, \%alignment, $motif) if $newmutations_array ne "NULL";
|
|
|
1096 @finalmutations_array = selectMutationArray($mutations_array, $mutations_array, \@tags, $alivehash, \%alignment, $motif) if $newmutations_array eq "NULL";
|
|
|
1097 # print "alt = $alternate{$conformation}\n";
|
|
|
1098
|
|
|
1099 my ($besttree, $treescore) = selectBetterTree($tree_analysis, $alternate{$conformation}, \@finalmutations_array);
|
|
|
1100 my $cleancase = "UNCLEAN";
|
|
|
1101 $cleancase = checkCleanCase($besttree, $finalalignment) if $treescore > 0 && $finalalignment ne "NULL" && $finalalignment =~ /\!/;
|
|
|
1102 $cleancase = checkCleanCase($besttree, $primaryalignment) if $treescore > 0 && $finalalignment eq "NULL" && $primaryalignment =~ /\!/ && $primaryalignment ne "NULL";
|
|
|
1103 $cleancase = "CLEAN" if $finalalignment eq "NULL" && $primaryalignment !~ /\!/ && $primaryalignment ne "NULL";
|
|
|
1104 $cleancase = "CLEAN" if $finalalignment ne "NULL" && $finalalignment !~ /\!/ ;
|
|
|
1105 # print "besttree = $besttree ... cleancase=$cleancase\n"; #<STDIN>;
|
|
|
1106
|
|
|
1107 my @selects = ("-C","+C","-H","+H","-HC","+HC","-O","+O","-H.-C","-H.-O","-HC,+C","-HC,+H","-HC.-O","-HCO,+HC","-HCO,+O","-O.-C","-O.-H",
|
|
|
1108 "+C.+O","+H.+C","+H.+O","+HC,-C","+HC,-H","+HC.+O","+HCO,-C","+HCO,-H","+HCO,-HC","+HCO,-O","+O.+C","+O.+H","+H.+C.+O","-H.-C.-O","+HCO","-HCO");
|
|
|
1109 next if (oneOf(@selects, $besttree) == 0);
|
|
|
1110 if ( ($besttree =~ /,/ || $besttree =~ /\./) && $cleancase eq "UNCLEAN"){
|
|
|
1111 $besttree = "$besttree / $tree_analysis";
|
|
|
1112 }
|
|
|
1113
|
|
|
1114 $besttree = "NULL" if $treescore <= 0;
|
|
|
1115
|
|
|
1116 while ($besttree =~ /[A-Z][A-Z]/){
|
|
|
1117 $besttree =~ s/([A-Z])([A-Z])/$1:$2/g;
|
|
|
1118 }
|
|
|
1119
|
|
|
1120 if ($besttree !~ /NULL/){
|
|
|
1121 my @elements = ($besttree =~ /([A-Z])/g);
|
|
|
1122
|
|
|
1123 foreach my $ele (@elements){
|
|
|
1124 # print "replacing $ele with $backReplacementArrTag{$ele}\n";
|
|
|
1125 $besttree =~ s/$ele/$backReplacementArrTag{$ele}/g if exists $backReplacementArrTag{$ele};
|
|
|
1126 }
|
|
|
1127 }
|
|
|
1128 my $endendstate = $focalspec_orig.".".$Hchrom."\t".$Hstart."\t".$Hend."\t".$locusmotif."\t".$locusmotifsize."\t".$tree_analysis."\t";
|
|
|
1129 next if $endendstate =~ /NA\tNA\tNA/;
|
|
|
1130
|
|
|
1131 print SUMMARY $focalspec_orig,".",$Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t";
|
|
|
1132 # print "SUMMARY\t", $focalspec_orig,".",$Hchrom,"\t",$Hstart,"\t",$Hend,"\t",$locusmotif,"\t",$locusmotifsize,"\t",$tree_analysis,"\t" ;
|
|
|
1133
|
|
|
1134 my @mutinfoarray =();
|
|
|
1135
|
|
|
1136 foreach my $lnode (@$lagestnodes){
|
|
|
1137 my @pair = @$lnode;
|
|
|
1138 my $joint = "(".join(", ",@pair).")";
|
|
|
1139 my @nodemutarr = ();
|
|
|
1140
|
|
|
1141 for my $p (@pair){
|
|
|
1142 foreach my $mut (@finalmutations_array){
|
|
|
1143 $mut =~ /node=([A-Z, \(\)]+)/;
|
|
|
1144 push @nodemutarr, $mut if $p eq $1;
|
|
|
1145 }
|
|
|
1146 @mutinfoarray = summarizeMutations(\@nodemutarr, $besttree);
|
|
|
1147
|
|
|
1148 # print SUMMARY join ("\t", @mutinfoarray[0...($#mutinfoarray-1)] ),"\t";
|
|
|
1149 # print join ("\t", @mutinfoarray[0...($#mutinfoarray-1)] ),"\t";
|
|
|
1150 }
|
|
|
1151 }
|
|
|
1152
|
|
|
1153 # print "G: $no_of_species\n";
|
|
|
1154
|
|
|
1155 my @alignmentarr = ();
|
|
|
1156
|
|
|
1157 foreach my $key (keys %alignment){
|
|
|
1158 push @alignmentarr, $backReplacementArrTag{$key}.":".$alignment{$key};
|
|
|
1159
|
|
|
1160 }
|
|
|
1161 # print "alignmentarr = @alignmentarr"; <STDIN>;
|
|
|
1162
|
|
|
1163 @mutinfoarray = summarizeMutations(\@finalmutations_array, $besttree);
|
|
|
1164 print SUMMARY join ("\t", @mutinfoarray ),"\t";
|
|
|
1165 print SUMMARY join(",",@alignmentarr),"\n";
|
|
|
1166 # print join("\t","--------------","\n",$besttree, join("",@tags)),"\n" if scalar(@tags) < 5;
|
|
|
1167 # <STDIN> if scalar(@tags) < 5;
|
|
|
1168 # print $cleancase, "\n";
|
|
|
1169 # print join ("\t", @mutinfoarray,$cleancase,join(",",@alignmentarr)),"\n"; #<STDIN>;
|
|
|
1170 # print "summarized\n"; #<STDIN>;
|
|
|
1171
|
|
|
1172
|
|
|
1173
|
|
|
1174 my %indelcatch = ();
|
|
|
1175 my %substcatch = ();
|
|
|
1176 my %typecatch = ();
|
|
|
1177 my %nodescatch = ();
|
|
|
1178 my $mutconcat = join("\t", @finalmutations_array)."\n";
|
|
|
1179 my %indelposcatch = ();
|
|
|
1180 my %subsposcatch = ();
|
|
|
1181
|
|
|
1182 foreach my $fmut ( @finalmutations_array){
|
|
|
1183 # next if $fmut !~ /indeltype=[a-zA-Z]+/;
|
|
|
1184 #print RUN $fmut, "\n";
|
|
|
1185 $fmut =~ /node=([a-zA-Z, \(\)]+)/;
|
|
|
1186 my $lnode = $1;
|
|
|
1187 $nodescatch{$1}=1;
|
|
|
1188
|
|
|
1189 if ($fmut =~ /type=substitution/){
|
|
|
1190 # print "fmut=$fmut\n";
|
|
|
1191 $fmut =~ /from=([a-zA-Z\-]+)\tto=([a-zA-Z\-]+)/;
|
|
|
1192 my $from=$1;
|
|
|
1193 # print "from=$from\n";
|
|
|
1194 my $to=$2;
|
|
|
1195 # print "to=$to\n";
|
|
|
1196 push @{$substcatch{$lnode}} , ("from:".$from." to:".$to);
|
|
|
1197 $fmut =~ /position=([0-9]+)/;
|
|
|
1198 push @{$subsposcatch{$lnode}}, $1;
|
|
|
1199 }
|
|
|
1200
|
|
|
1201 if ($fmut =~ /insertion=[a-zA-Z\-]+/){
|
|
|
1202 $fmut =~ /insertion=([a-zA-Z\-]+)/;
|
|
|
1203 push @{$indelcatch{$lnode}} , $1;
|
|
|
1204 $fmut =~ /indeltype=([a-zA-Z]+)/;
|
|
|
1205 push @{$typecatch{$lnode}}, $1;
|
|
|
1206 $fmut =~ /position=([0-9]+)/;
|
|
|
1207 push @{$indelposcatch{$lnode}}, $1;
|
|
|
1208 }
|
|
|
1209 if ($fmut =~ /deletion=[a-zA-Z\-]+/){
|
|
|
1210 $fmut =~ /deletion=([a-zA-Z\-]+)/;
|
|
|
1211 push @{$indelcatch{$lnode}} , $1;
|
|
|
1212 $fmut =~ /indeltype=([a-zA-Z]+)/;
|
|
|
1213 push @{$typecatch{$lnode}}, $1;
|
|
|
1214 $fmut =~ /position=([0-9]+)/;
|
|
|
1215 push @{$indelposcatch{$lnode}}, $1;
|
|
|
1216 }
|
|
|
1217 }
|
|
|
1218
|
|
|
1219 # print $simplified_microsat,"\t", $tree_analysis,"\t", join("",@tags), "\t" if $printer == 1;
|
|
|
1220 # print join ("<\t>", @mutinfoarray),"\n" if $printer == 1;
|
|
|
1221 # print "where mutinfoarray = @mutinfoarray\n" if $printer == 1;
|
|
|
1222 # #print RUN ".";
|
|
|
1223
|
|
|
1224 # print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
|
|
|
1225 # print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
|
|
|
1226
|
|
|
1227 # print colored ['red'],"finalmutations_array=\n" if $printer == 1;
|
|
|
1228 # foreach (@finalmutations_array) {
|
|
|
1229 # print colored ['red'], "$_\n" if $_ =~ /type=substitution/ && $printer == 1 ;
|
|
|
1230 # print colored ['yellow'], "$_\n" if $_ !~ /type=substitution/ && $printer == 1 ;
|
|
|
1231
|
|
|
1232 # }# if $line =~ /cal/;# && $line =~ /chr4/;
|
|
|
1233
|
|
|
1234 # print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
|
|
|
1235 # print colored ['red'], "-------------------------------------------------------------\n" if $printer == 1;
|
|
|
1236 # print "tree analysis = $tree_analysis\n" if $printer == 1;
|
|
|
1237
|
|
|
1238 # my $mutations = "@$mutations_array";
|
|
|
1239
|
|
|
1240
|
|
|
1241 next;
|
|
|
1242 for my $keys (@$nodes) {foreach my $key (@$keys){
|
|
|
1243 #print "key = $key, => $branches_hash->{$key}\n";
|
|
|
1244 }
|
|
|
1245 # print "x" x 50, "\n";
|
|
|
1246 }
|
|
|
1247 my ($birth_steps, $death_steps) = decipher_history($mutations_array,join("",@tags),$nodes,$branches_hash,$tree_analysis,$conformation, $alivehash, $simplified_microsat);
|
|
|
1248 }
|
|
|
1249 }
|
|
|
1250 close BO;
|
|
|
1251 # print "now studying where we lost microsatellites:\n";
|
|
|
1252 # print "x" x 60,"\n";
|
|
|
1253 for my $reason (0 ... $#nextcounter){
|
|
|
1254 # print $removalReasons[$reason]."\t".$nextcounter[$reason],"\n";
|
|
|
1255 }
|
|
|
1256 # print "x" x 60,"\n";
|
|
|
1257 # print "In total we read $microreadcounter microsatellites after reading through $contigsentered contigs\n";
|
|
|
1258 # print " we lost $keynotformedcounter contigs as they did not form the key, \n";
|
|
|
1259 # print "$contignotrightcounter contigs as they were not of the right species configuration\n";
|
|
|
1260 # print "$keynotfoundcounter contigs as they did not contain the microsats\n";
|
|
|
1261 # print "... In total we went through a file that had $dotcounter contigs...\n";
|
|
|
1262 # print join ("\n","remaining orth keys = ", (keys %orths),"");
|
|
|
1263 # print "------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ \n";
|
|
|
1264 # print "now printing counted trees: \n";
|
|
|
1265 if (scalar(keys %registeredTrees) > 0){
|
|
|
1266 foreach my $keyb ( sort (keys %registeredTrees) )
|
|
|
1267 {
|
|
|
1268 # print "$keyb : $registeredTrees{$keyb}\n";
|
|
|
1269 }
|
|
|
1270 }
|
|
|
1271
|
|
|
1272
|
|
|
1273 }
|
|
|
1274 close SUMMARY;
|
|
|
1275
|
|
|
1276 my @summarizarr = ("+C=+C +R.+C -HCOR,+C",
|
|
|
1277 "+H=+H +R.+H -HCOR,+H",
|
|
|
1278 "-C=-C -R.-C +HCOR,-C",
|
|
|
1279 "-H=-H -R.-H +HCOR,-H",
|
|
|
1280 "+HC=+HC",
|
|
|
1281 "-HC=-HC",
|
|
|
1282 "+O=+O -HCOR,+O",
|
|
|
1283 "-O=-O +HCOR,-O",
|
|
|
1284 "+HCO=+HCO",
|
|
|
1285 "-HCO=-HCO",
|
|
|
1286 "+R=+R +R.+C +R.+H",
|
|
|
1287 "-R=-R -R.-C -R.-H");
|
|
|
1288
|
|
|
1289 foreach my $line (@summarizarr){
|
|
|
1290 next if $line !~ /[A-Za-z0-9]/;
|
|
|
1291 # print $line;
|
|
|
1292 chomp $line;
|
|
|
1293 my @fields = split(/=/,$line);
|
|
|
1294 # print "title = $fields[0]\n";
|
|
|
1295 my @parts=split(/ +/, $fields[1]);
|
|
|
1296 my %partshash = ();
|
|
|
1297 foreach my $part (@parts){$partshash{$part}=1;}
|
|
|
1298 my $count=0;
|
|
|
1299 foreach my $key ( sort keys %registeredTrees ){
|
|
|
1300 next if !exists $partshash{$key};
|
|
|
1301 # print "now adding $registeredTrees{$key} from $key\n";
|
|
|
1302 $count+=$registeredTrees{$key};
|
|
|
1303 }
|
|
|
1304 # print "$fields[0] : $count\n";
|
|
|
1305 }
|
|
|
1306 my $rootdir = $dir;
|
|
|
1307 $rootdir =~ s/\/[A-Za-z0-9\-_]+$//;
|
|
|
1308 chdir $rootdir;
|
|
|
1309 remove_tree($dir);
|
|
|
1310
|
|
|
1311
|
|
|
1312 #--------------------------------------------------------------------------------------------------------
|
|
|
1313 #--------------------------------------------------------------------------------------------------------
|
|
|
1314 #--------------------------------------------------------------------------------------------------------
|
|
|
1315 #--------------------------------------------------------------------------------------------------------
|
|
|
1316 #--------------------------------------------------------------------------------------------------------
|
|
|
1317 sub largest_number{
|
|
|
1318 my $counter = 0;
|
|
|
1319 my($max) = shift(@_);
|
|
|
1320 foreach my $temp (@_) {
|
|
|
1321
|
|
|
1322 #print "finding largest array: $maxcounter \n";
|
|
|
1323 if($temp > $max){
|
|
|
1324 $max = $temp;
|
|
|
1325 }
|
|
|
1326 }
|
|
|
1327 return($max);
|
|
|
1328 }
|
|
|
1329
|
|
|
1330 sub smallest_number{
|
|
|
1331 my $counter = 0;
|
|
|
1332 my($min) = shift(@_);
|
|
|
1333 foreach my $temp (@_) {
|
|
|
1334 #print "finding largest array: $maxcounter \n";
|
|
|
1335 if($temp < $min){
|
|
|
1336 $min = $temp;
|
|
|
1337 }
|
|
|
1338 }
|
|
|
1339 return($min);
|
|
|
1340 }
|
|
|
1341 #--------------------------------------------------------------------------------------------------------
|
|
|
1342 #--------------------------------------------------------------------------------------------------------
|
|
|
1343 #--------------------------------------------------------------------------------------------------------
|
|
|
1344 #--------------------------------------------------------------------------------------------------------
|
|
|
1345 #--------------------------------------------------------------------------------------------------------
|
|
|
1346 sub baseml_parser{
|
|
|
1347 my $outputfile = $_[0];
|
|
|
1348 open(BOUT,"<$outputfile") or die "Cannot open output of upstream baseml $outputfile: $!";
|
|
|
1349 my @info = ();
|
|
|
1350 my @branchields = ();
|
|
|
1351 my @distanceields = ();
|
|
|
1352 my @bout = <BOUT>;
|
|
|
1353 #print colored ['red'], @bout ,"\n";
|
|
|
1354 for my $b (0 ... $#bout){
|
|
|
1355 my $bine=$bout[$b];
|
|
|
1356 #print colored ['yellow'], "sentence = ",$bine;
|
|
|
1357 if ($bine =~ /TREE/){
|
|
|
1358 $bine=$bout[$b++];
|
|
|
1359 $bine=$bout[$b++];
|
|
|
1360 $bine=$bout[$b++];
|
|
|
1361 #print "FOUND",$bine;
|
|
|
1362 chomp $bine;
|
|
|
1363 $bine =~ s/^\s+//g;
|
|
|
1364 @branchields = split(/\s+/,$bine);
|
|
|
1365 $bine=$bout[$b++];
|
|
|
1366 chomp $bine;
|
|
|
1367 $bine =~ s/^\s+//g;
|
|
|
1368 @distanceields = split(/\s+/,$bine);
|
|
|
1369 #print "LASTING..............\n";
|
|
|
1370 last;
|
|
|
1371 }
|
|
|
1372 else{
|
|
|
1373 }
|
|
|
1374 }
|
|
|
1375
|
|
|
1376 close BOUT;
|
|
|
1377 # print "branchfields = @branchields and distanceields = @distanceields\n" if $printer == 1;
|
|
|
1378 my %distance_hash=();
|
|
|
1379 for my $d (0 ... $#branchields){
|
|
|
1380 $distance_hash{$branchields[$d]} = $distanceields[$d];
|
|
|
1381 }
|
|
|
1382
|
|
|
1383 $info[0] = $distance_hash{"9..1"} + $distance_hash{"9..2"};
|
|
|
1384 $info[1] = $distance_hash{"9..1"} + $distance_hash{"8..9"}+ $distance_hash{"8..3"};
|
|
|
1385 $info[2] = $distance_hash{"9..1"} + $distance_hash{"8..9"}+$distance_hash{"7..8"}+$distance_hash{"7..4"};
|
|
|
1386 $info[3] = $distance_hash{"9..1"} + $distance_hash{"8..9"}+$distance_hash{"7..8"}+$distance_hash{"6..7"}+$distance_hash{"6..5"};
|
|
|
1387
|
|
|
1388 # print "\nsending back: @info\n" if $printer == 1;
|
|
|
1389
|
|
|
1390 return join("\t",@info);
|
|
|
1391
|
|
|
1392 }
|
|
|
1393
|
|
|
1394
|
|
|
1395 #--------------------------------------------------------------------------------------------------------
|
|
|
1396 #--------------------------------------------------------------------------------------------------------
|
|
|
1397 #--------------------------------------------------------------------------------------------------------
|
|
|
1398 #--------------------------------------------------------------------------------------------------------
|
|
|
1399 #--------------------------------------------------------------------------------------------------------
|
|
|
1400 sub test_complexity{
|
|
|
1401 my $printer = 0;
|
|
|
1402 my $sequence = $_[0];
|
|
|
1403 #print "sequence = $sequence\n";
|
|
|
1404 my $COMPLEXITY_SUPPORT = $_[1];
|
|
|
1405 my $complexity=int($COMPLEXITY_SUPPORT * (1/40)); #1 THIS IS AN ARBITRARY THRESHOLD SET FOR LOW COMPLEXITY.
|
|
|
1406 #2 THE INSPIRATION WAS WEB MILLER'S MAIL SENT ON
|
|
|
1407 #3 19 Apr 2008 WHERE HE CLASSED AS HIGH COMPLEXITY
|
|
|
1408 #4 REGION, IF 40 BP OF SEQUENCE HAS AT LEAST 3 OF
|
|
|
1409 #5 EACH NUCLEOTIDE. HENCE, I NORMALIZE THIS PARAMETER
|
|
|
1410 #6 FOR THE ACTUAL LENGTH OF $FLANK_SUPPORT SET BY
|
|
|
1411 #7 THE USER.
|
|
|
1412 #8 WEB MILLER SENT THE MAIL TO YDK104@PSU.EDU
|
|
|
1413
|
|
|
1414
|
|
|
1415
|
|
|
1416 my $As = ($sequence=~ s/A/A/gi);
|
|
|
1417 my $Ts = ($sequence=~ s/T/T/gi);
|
|
|
1418 my $Gs = ($sequence=~ s/G/G/gi);
|
|
|
1419 my $Cs = ($sequence=~ s/C/C/gi);
|
|
|
1420 my $dashes = ($sequence=~ s/\-/-/gi);
|
|
|
1421 $dashes = 0 if $sequence !~ /\-/;
|
|
|
1422 # print "seq = $sequence, As=$As, Ts=$Ts, Gs=$Gs, Cs=$Cs, dashes=$dashes\n";
|
|
|
1423 return "LOW" if $dashes > length($sequence)/2;
|
|
|
1424
|
|
|
1425 my $ans = ();
|
|
|
1426
|
|
|
1427 return "HIGH" if $As >= $complexity && $Ts >= $complexity && $Cs >= $complexity && $Gs >= $complexity;
|
|
|
1428
|
|
|
1429 my @nts = ("A","T","G","C","-");
|
|
|
1430
|
|
|
1431 my $lowcomplex = 0;
|
|
|
1432
|
|
|
1433 foreach my $nt (@nts){
|
|
|
1434 $lowcomplex =1 if $sequence =~ /(($nt\-*){$mono_flanksimplicityRepno,})/i;
|
|
|
1435 $lowcomplex =1 if $sequence =~ /(($nt[A-Za-z]){$di_flanksimplicityRepno,})/i;
|
|
|
1436 $lowcomplex =1 if $sequence =~ /(([A-Za-z]$nt){$di_flanksimplicityRepno,})/i;
|
|
|
1437 my $nont = ($sequence=~ s/$nt/$nt/gi);
|
|
|
1438 $lowcomplex = 1 if $nont > (length($sequence) * $prop_of_seq_allowedtoAT) && ($nt =~ /[AT\-]/);
|
|
|
1439 $lowcomplex = 1 if $nont > (length($sequence) * $prop_of_seq_allowedtoCG) && ($nt =~ /[CG]/);
|
|
|
1440 }
|
|
|
1441 # print "leaving for now.. $sequence\n" if $printer == 1 && $lowcomplex == 0;
|
|
|
1442 #<STDIN>;
|
|
|
1443 return "HIGH" if $lowcomplex == 0;
|
|
|
1444 return "LOW" ;
|
|
|
1445 }
|
|
|
1446 #--------------------------------------------------------------------------------------------------------
|
|
|
1447 #--------------------------------------------------------------------------------------------------------
|
|
|
1448 #--------------------------------------------------------------------------------------------------------
|
|
|
1449 #--------------------------------------------------------------------------------------------------------
|
|
|
1450 #--------------------------------------------------------------------------------------------------------
|
|
|
1451 sub sequence_similarity{
|
|
|
1452 my $printer = 0;
|
|
|
1453 my @seq1 = split(/\s*/, $_[0]);
|
|
|
1454 my @seq2 = split(/\s*/, $_[1]);
|
|
|
1455 my $similarity_thresh = $_[2];
|
|
|
1456 my $info = $_[3];
|
|
|
1457 # print "input = @_\n" if $printer == 1;
|
|
|
1458 my $seq1str = $_[0];
|
|
|
1459 my $seq2str = $_[1];
|
|
|
1460 $seq1str=~s/\-//g; $seq2str=~s/\-//g;
|
|
|
1461 my $similarity=0;
|
|
|
1462
|
|
|
1463 my $nucdiffs=0;
|
|
|
1464 my $nucsims=0;
|
|
|
1465 my $indeldiffs=0;
|
|
|
1466
|
|
|
1467 for my $i (0...$#seq1){
|
|
|
1468 $similarity++ if $seq1[$i] =~ /$seq2[$i]/i ; #|| $seq1[$i] =~ /\-/i || $seq2[$i] =~ /\-/i ;
|
|
|
1469 $nucsims++ if $seq1[$i] =~ /$seq2[$i]/i && ($seq1[$i] =~ /[a-zA-Z]/i && $seq2[$i] =~ /[a-zA-Z]/i);
|
|
|
1470 $nucdiffs++ if $seq1[$i] !~ /$seq2[$i]/i && ($seq1[$i] =~ /[a-zA-Z]/i && $seq2[$i] =~ /[a-zA-Z]/i);
|
|
|
1471 $indeldiffs++ if $seq1[$i] !~ /$seq2[$i]/i && $seq1[$i] =~ /\-/i || $seq2[$i] =~ /\-/i;
|
|
|
1472 }
|
|
|
1473 my $sim = $similarity/length($_[0]);
|
|
|
1474 return ( $sim, $nucdiffs, $indeldiffs ); #<= $similarity_thresh;
|
|
|
1475 }
|
|
|
1476 #--------------------------------------------------------------------------------------------------------
|
|
|
1477 #--------------------------------------------------------------------------------------------------------
|
|
|
1478 #--------------------------------------------------------------------------------------------------------
|
|
|
1479
|
|
|
1480 sub load_treesToReject{
|
|
|
1481 my @rejectlist = ();
|
|
|
1482 my $alltags = join("",@_);
|
|
|
1483 @rejectlist = qw (-HCOR +HCOR) if $alltags eq "HCORM";
|
|
|
1484 @rejectlist = qw ( -HCO|+R +HCO|-R) if $alltags eq "HCOR";
|
|
|
1485 @rejectlist = qw ( -HC|+O +HC|-O) if $alltags eq "HCO";
|
|
|
1486
|
|
|
1487 %treesToReject=();
|
|
|
1488 $treesToReject{$_} = $_ foreach (@rejectlist);
|
|
|
1489 #print "loaded to reject for $alltags; ", $treesToReject{$_},"\n" foreach (@rejectlist); #<STDIN>;
|
|
|
1490 }
|
|
|
1491 #--------------------------------------------------------------------------------------------------------
|
|
|
1492 sub load_treesToIgnore{
|
|
|
1493 my @rejectlist = ();
|
|
|
1494 my $alltags = join("",@_);
|
|
|
1495 @rejectlist = qw (-HCOR +HCOR +HCORM -HCORM) if $alltags eq "HCORM";
|
|
|
1496 @rejectlist = qw ( -HCO|+R +HCO|-R +HCOR -HCOR) if $alltags eq "HCOR";
|
|
|
1497 @rejectlist = qw ( -HC|+O +HC|-O +HCO -HCO) if $alltags eq "HCO";
|
|
|
1498
|
|
|
1499 %treesToIgnore=();
|
|
|
1500 $treesToIgnore{$_} = $_ foreach (@rejectlist);
|
|
|
1501 #print "loaded ", $treesToIgnore{$_},"\n" foreach (@rejectlist);
|
|
|
1502 }
|
|
|
1503 #--------------------------------------------------------------------------------------------------------
|
|
|
1504 sub load_thresholds{
|
|
|
1505 my @threshold_array=split(/[,_]/,$_[0]);
|
|
|
1506 unshift @threshold_array, "0";
|
|
|
1507 for my $size (1 ... 4){
|
|
|
1508 $thresholdhash{$size}=$threshold_array[$size];
|
|
|
1509 }
|
|
|
1510 }
|
|
|
1511 #--------------------------------------------------------------------------------------------------------
|
|
|
1512 sub load_allPossibleTrees{
|
|
|
1513 #1 THIS FILE STORES ALL POSSIBLE SCENARIOS OF MICROSATELLITE
|
|
|
1514 #2 BIRTH AND DEATH EVENTS ON A 5-PRIMATE TREE OF H,C,O,R,M
|
|
|
1515 #3 IN FORM OF A TEXT FILE. THIS WILL BE USED AS A TEMPLET
|
|
|
1516 #4 TO COMPARE EACH MICROSATELLITE CLUSTER TO UNDERSTAND THE
|
|
|
1517 #5 EVOLUTION OF EACH LOCUS. WE WILL THEN DISCARD SOME
|
|
|
1518 #6 MICROSATS ACCRODING TO THEIR EVOLUTIONARY BEHAVIOUR ON
|
|
|
1519 #7 THE TREE. MOST PROBABLY WE WILL REMOVE THOSE MICROSATS
|
|
|
1520 #8 THAT ARE NOT SUFFICIENTLY INFORMATIVE, LIKE IN CASE OF
|
|
|
1521 #9 AN OUTGROUP MICROSATELLITE BEING DIFFERENT FRON ALL OTHER
|
|
|
1522 #10 SPECIES IN THE TREE.
|
|
|
1523 my $tree_list = $_[0];
|
|
|
1524 # print "file to be loaded: $tree_list\n";
|
|
|
1525
|
|
|
1526 my @trarr = ();
|
|
|
1527 @trarr = ("#H C O CONCLUSION ALTERNATE",
|
|
|
1528 "+ + + +HCO NA",
|
|
|
1529 "+ _ _ +H NA",
|
|
|
1530 "_ + _ +C NA",
|
|
|
1531 "_ _ + -HC|+O NA",
|
|
|
1532 "+ _ + -C +H",
|
|
|
1533 "_ + + -H +C",
|
|
|
1534 "+ + _ +HC|-O NA",
|
|
|
1535 "_ _ _ -HCO NA") if $tree_list =~ /_HCO\.txt/;
|
|
|
1536 @trarr = ("#H C O R CONCLUSION ALTERNATE",
|
|
|
1537 "_ _ _ _ -HCOR NA",
|
|
|
1538 "+ + + + +HCOR NA",
|
|
|
1539 "+ + + _ +HCO|-R +H.+C.+O",
|
|
|
1540 "+ + _ _ +HC +H.+C;-O",
|
|
|
1541 "+ _ _ _ +H +HC,-C;+HC,-C",
|
|
|
1542 "_ + _ _ +C +HC,-H;+HC,-H",
|
|
|
1543 "_ _ + _ +O -HC|-H.-C",
|
|
|
1544 "_ _ + + -HC -H.-C",
|
|
|
1545 "+ _ _ + +H|-C.-O +HC,-C",
|
|
|
1546 "_ + _ + +C -H.-O",
|
|
|
1547 "_ + + _ -H +C.+O",
|
|
|
1548 "_ _ _ + -HCO|+R NA",
|
|
|
1549 "+ _ + _ +H.+O|-C NA",
|
|
|
1550 "_ + + + -H -HC,+C",
|
|
|
1551 "+ _ + + -C -HC,+H",
|
|
|
1552 "+ + _ + -O +HC") if $tree_list =~ /_HCOR\.txt/;
|
|
|
1553
|
|
|
1554 @trarr = ("#H C O R M CONCLUSION ALTERNATE",
|
|
|
1555 "+ + _ + + -O -HCO,+HC|-HCO,+HC;-HCO,(+H.+C)",
|
|
|
1556 "+ _ + + + -C -HC,+H;+HCO,(+H.+O)",
|
|
|
1557 "_ + + + + -H -HC,+C;-HCO,(+C.+O)",
|
|
|
1558 "_ _ + _ _ +O +HCO,-HC;+HCO,(-H.-C)",
|
|
|
1559 "_ + _ _ _ +C +HC,-H;+HCO,(-H.-O)",
|
|
|
1560 "+ _ _ _ _ +H +HC,-C;+HCO,(-C.-O)",
|
|
|
1561 "+ + + _ _ +HCO +H.+C.+O",
|
|
|
1562 "_ _ _ + + -HCO -HC.-O;-H.-C.-O",
|
|
|
1563 "+ _ _ + + -O.-C|-HCO,+H +R.+H;-HCO,(+R.+H)",
|
|
|
1564 "_ + _ + + -O.-H|-HCO,+C +R.+C;-HCO,(+R.+C)",
|
|
|
1565 "_ + + _ _ +HCO,-H|+O.+C NA",
|
|
|
1566 "+ _ + _ _ +HCO,-C|+O.+H NA",
|
|
|
1567 "_ _ + + + -HC -H.-C|-HCO,+O",
|
|
|
1568 "+ + _ _ _ +HC +H.+C|+HCO,-O|-HCO,+HC;-HCO,(+H.+C)",
|
|
|
1569 "+ + + + + +HCORM NA",
|
|
|
1570 "_ _ + _ + DISCARD +O;+HCO,-HC;+HCO,(-H.-C)",
|
|
|
1571 "_ + _ _ + +C +HC,-H;+HCO,(-H.-O)",
|
|
|
1572 "+ _ _ _ + +H +HC,-C;+HCO,(-C.-O)",
|
|
|
1573 "+ + _ _ + +HC -R.-O|+HCO,-O|+H.+C;-HCO,+HC;-HCO,(+H.+C)",
|
|
|
1574 "+ _ + _ + DISCARD -R.-C|+HCO,-C|+H.+O NA",
|
|
|
1575 "_ + + _ + DISCARD -R.-H|+HCO,-H|+C.+O NA",
|
|
|
1576 "_ _ _ _ + DISCARD -HCOR NA",
|
|
|
1577 "_ _ _ + _ DISCARD +R;-HC.-O;-H.-C.-O",
|
|
|
1578 "+ + _ + _ -O +R.+HC|-HCO,+HC;+H.+C.+R|-HCO,(+H.+C)",
|
|
|
1579 "+ + + + _ +HCOR NA",
|
|
|
1580 "+ + + _ + DISCARD -R;+HCO;+HC.+O;+H.+C.+O",
|
|
|
1581 "+ _ + + _ -C -HC,+H;+H.+O.+R|-HCO,(+H.+O)",
|
|
|
1582 "_ + + + _ -H -HC,+C;+C.+O.+R|-HCO,(+C.+O)",
|
|
|
1583 "_ _ + + _ -HC +R.+O|-HCO,+O|+HCO,-HC",
|
|
|
1584 "_ + _ + _ +C +R.+C|-HCO,+C|-HC,+C +HCO,(-H.-O)",
|
|
|
1585 "+ _ _ + _ +H +R.+H|-C.-O +HCO,(-C.-O)"
|
|
|
1586 ) if $tree_list =~ /_HCORM\.txt/;
|
|
|
1587
|
|
|
1588
|
|
|
1589 my $template_p = $_[1];
|
|
|
1590 my $alternate_p = $_[2];
|
|
|
1591 #1 THIS IS THE HASH IN WHICH INFORMATION FROM THE ABOVE FILE
|
|
|
1592 #2 GETS STORED, USING THE WHILE LOOP BELOW. HERE, THE KEY
|
|
|
1593 #3 OF EACH ROW IS THE EVOLUTIONARY CONFIGURATION OF A LOCUS
|
|
|
1594 #4 ON THE PRIMATE TREE, BASED ON PRESENCE/ABSENCE OF A MICROSAT
|
|
|
1595 #5 AT THAT LOCUS, LIKE SAY "+ + + _ _" .. EACH COLUMN BELONGS
|
|
|
1596 #6 TO ONE SPECIES; HERE THE COLUMN NAMES ARE "H C O R M".
|
|
|
1597 #7 THE VALUE FOR EACH ENTRY IS THE MEANING OF THE ABOVE
|
|
|
1598 #8 CONFIGURATION (I.E., CONFIGURAION OF THE KEY. HERE, THE
|
|
|
1599 #9 VALUE WILL BE +HCO, SIGNIFYING A BIRTH IN HUMAN-CHIMP-ORANG
|
|
|
1600 #10 COMMON ANCESTOR. THIS HASH HAS BEEN LOADED HERE TO BE USED
|
|
|
1601 #11 LATER BY THE SUBROUTINE sub treeStudy{} THAT STUDIES
|
|
|
1602 #12 EVOLUTIONARY CONFIGURAION OF EACH MICROSAT LOCUS, AS
|
|
|
1603 #13 MENTIONED ABOVE.
|
|
|
1604 my @keys_array=();
|
|
|
1605 foreach my $line (@trarr){
|
|
|
1606 # print $line,"\n";
|
|
|
1607 next if $line =~ /^#/;
|
|
|
1608 chomp $line;
|
|
|
1609 my @fields = split("\t", $line);
|
|
|
1610 push @keys_array, $fields[0];
|
|
|
1611 # print "loading: $fields[0]\n";
|
|
|
1612 $template_p->{$fields[0]}[0] = $fields[1];
|
|
|
1613 $template_p->{$fields[0]}[1] = 0;
|
|
|
1614 $alternate_p->{$fields[0]} = $fields[2];
|
|
|
1615 # $alternate_p->{$fields[1]} = $fields[2];
|
|
|
1616 # print "loading alternate_p $fields[1] $fields[2]\n"; #<STDIN> if $fields[1] eq "+H";
|
|
|
1617 }
|
|
|
1618 # print "loaded the trees with keys: @keys_array\n";
|
|
|
1619 return $template_p, \@keys_array, $alternate_p;
|
|
|
1620 }
|
|
|
1621
|
|
|
1622 #--------------------------------------------------------------------------------------------------------
|
|
|
1623 #--------------------------------------------------------------------------------------------------------
|
|
|
1624 #--------------------------------------------------------------------------------------------------------
|
|
|
1625 #--------------------------------------------------------------------------------------------------------
|
|
|
1626 #--------------------------------------------------------------------------------------------------------
|
|
|
1627 sub checkCleanCase{
|
|
|
1628 my $printer = 0;
|
|
|
1629 my $tree = $_[0];
|
|
|
1630 my $finalalignment = $_[1];
|
|
|
1631
|
|
|
1632 #print "IN checkCleanCase: @_\n";
|
|
|
1633 #<STDIN>;
|
|
|
1634 my @indivspecies = $tree =~ /[A-Z]/g;
|
|
|
1635 $finalalignment =~ s/\./_/g;
|
|
|
1636 my @captured = $finalalignment =~ /[A-Za-z, \(\):]+\![:A-Za-z, \(\)]/g;
|
|
|
1637
|
|
|
1638 my $unclean = 0;
|
|
|
1639
|
|
|
1640 foreach my $sp (@indivspecies){
|
|
|
1641 foreach my $cap (@captured){
|
|
|
1642 $cap =~ s/:[A-Za-z\-]+//g;
|
|
|
1643 my @sps = $cap =~ /[A-Z]+/g;
|
|
|
1644 my $spsc = join("", @sps);
|
|
|
1645 # print "checking whether imp species $sp is present in $cap i.e, in $spsc\n " if $printer == 1;
|
|
|
1646 if ($spsc =~ /$sp/){
|
|
|
1647 # print "foind : $sp\n";
|
|
|
1648 $unclean = 1; last;
|
|
|
1649 }
|
|
|
1650 }
|
|
|
1651 last if $unclean == 1;
|
|
|
1652 }
|
|
|
1653 #<STDIN>;
|
|
|
1654 return "CLEAN" if $unclean == 0;
|
|
|
1655 return "UNCLEAN";
|
|
|
1656 }
|
|
|
1657
|
|
|
1658 #--------------------------------------------------------------------------------------------------------
|
|
|
1659 #--------------------------------------------------------------------------------------------------------
|
|
|
1660 #--------------------------------------------------------------------------------------------------------
|
|
|
1661 #--------------------------------------------------------------------------------------------------------
|
|
|
1662 #--------------------------------------------------------------------------------------------------------
|
|
|
1663
|
|
|
1664
|
|
|
1665 sub adjustCoordinates{
|
|
|
1666 my $line = $_[0];
|
|
|
1667 return 0 if !defined $line;
|
|
|
1668 #print "------x------x------x------x------x------x------x------x------\n";
|
|
|
1669 #print $line,"\n\n";
|
|
|
1670 my $no_of_species = $line =~ s/(chr[0-9a-zA-Z]+)|(Contig[0-9a-zA-Z\._\-]+)|(scaffold[0-9a-zA-Z\._\-]+)|(supercontig[0-9a-zA-Z\._\-]+)/x/ig;
|
|
|
1671 #print $line,"\n";
|
|
|
1672 #print "------x------x------x------x------x------x------x------x------\n\n\n";
|
|
|
1673 # my @got = ($line =~ s/(chr[0-9a-zA-Z]+)|(Contig[0-9a-zA-Z\._\-]+)/x/g);
|
|
|
1674 # print "line = $line\n";
|
|
|
1675 $infocord = 2 + (4*$no_of_species) - 1;
|
|
|
1676 $typecord = 2 + (4*$no_of_species) + 1 - 1;
|
|
|
1677 $motifcord = 2 + (4*$no_of_species) + 2 - 1;
|
|
|
1678 $gapcord = $motifcord+1;
|
|
|
1679 $startcord = $gapcord+1;
|
|
|
1680 $strandcord = $startcord+1;
|
|
|
1681 $endcord = $strandcord + 1;
|
|
|
1682 $microsatcord = $endcord + 1;
|
|
|
1683 $sequencepos = 2 + (5*$no_of_species) + 1 -1 ;
|
|
|
1684 $interr_poscord = $microsatcord + 3;
|
|
|
1685 $no_of_interruptionscord = $microsatcord + 4;
|
|
|
1686 $interrcord = $microsatcord + 2;
|
|
|
1687 # print "$line\n startcord = $startcord, and endcord = $endcord and no_of_species = $no_of_species\n" if $line !~ /calJac/i;
|
|
|
1688
|
|
|
1689 return $no_of_species;
|
|
|
1690 }
|
|
|
1691
|
|
|
1692
|
|
|
1693 sub printhash{
|
|
|
1694 my $alivehash = $_[0];
|
|
|
1695 my @tags = @$_[1];
|
|
|
1696 # print "print hash\n";
|
|
|
1697 foreach my $tag (@tags){
|
|
|
1698 # print "$tag=",$alivehash->{$tag},"\n" if exists $alivehash->{$tag};
|
|
|
1699 }
|
|
|
1700
|
|
|
1701 return "\n"
|
|
|
1702 }
|
|
|
1703 sub peel_onion{
|
|
|
1704 my $printer = 0;
|
|
|
1705 # print "received: @_\n" ; #<STDIN>;
|
|
|
1706 $printer = 0;
|
|
|
1707 my ($tree, $sequences, $alignment, $tagarray, $microsathash, $nonmicrosathash, $motif, $tree_analysis, $threshold, $microsatstarts) = @_;
|
|
|
1708 # print "in peel onion.. tree = $tree \n" if $printer == 1;
|
|
|
1709 my %sequence_hash=();
|
|
|
1710
|
|
|
1711
|
|
|
1712 # for my $i (0 ... $#sequences){ $sequence_hash{$species[$i]}=$sequences->[$i]; }
|
|
|
1713
|
|
|
1714
|
|
|
1715 my %node_sequences=();
|
|
|
1716
|
|
|
1717 my %node_alignments = (); #NEW, Nov 28 2008
|
|
|
1718 my @tags=();
|
|
|
1719 my @locus_sequences=();
|
|
|
1720 my %alivehash=();
|
|
|
1721 foreach my $tag (@$tagarray) {
|
|
|
1722 #print "adding: $tag\n";
|
|
|
1723 push(@tags, $tag);
|
|
|
1724 $node_sequences{$tag}=join ".",split(/\s*/,$microsathash->{$tag}) if $microsathash->{$tag} ne "NULL";
|
|
|
1725 $alivehash{$tag}= $tag if $microsathash->{$tag} ne "NULL";
|
|
|
1726 $node_sequences{$tag}=join ".",split(/\s*/,$nonmicrosathash->{$tag}) if $microsathash->{$tag} eq "NULL";
|
|
|
1727 $node_alignments{$tag}=join ".",split(/\s*/,$alignment->{$tag}) ;
|
|
|
1728 push @locus_sequences, $node_sequences{$tag};
|
|
|
1729 # print "adding to node_seq: $tag = ",$node_alignments{$tag},"\n";
|
|
|
1730 }
|
|
|
1731
|
|
|
1732 #<STDIN>;
|
|
|
1733
|
|
|
1734 my ($nodes_arr, $branches_hash) = get_nodes($tree);
|
|
|
1735 my @nodes=@$nodes_arr;
|
|
|
1736 # print "recieved nodes = " if $printer == 1;
|
|
|
1737 # foreach my $key (@nodes) {print "@$key " if $printer == 1;}
|
|
|
1738
|
|
|
1739 # print "\n" if $printer == 1;
|
|
|
1740
|
|
|
1741 #POPULATE branches_hash WITH INFORMATION ABOUT LIVESTATUS
|
|
|
1742 foreach my $keys (@nodes){
|
|
|
1743 my @pair = @$keys;
|
|
|
1744 my $joint = "(".join(", ",@pair).")";
|
|
|
1745 my $copykey = join "", @pair;
|
|
|
1746 $copykey =~ s/[\W ]+//g;
|
|
|
1747 # print "for node: $keys, copykey = $copykey and joint = $joint\n" if $printer == 1;
|
|
|
1748 my $livestatus = 1;
|
|
|
1749 foreach my $copy (split(/\s*/,$copykey)){
|
|
|
1750 $livestatus = 0 if !exists $alivehash{$copy};
|
|
|
1751 }
|
|
|
1752 $alivehash{$joint} = $joint if !exists $alivehash{$joint} && $livestatus == 1;
|
|
|
1753 # print "alivehash = $alivehash{$joint}\n" if exists $alivehash{$joint} && $printer == 1;
|
|
|
1754 }
|
|
|
1755
|
|
|
1756 @nodes = reverse(@nodes); #1 THIS IS IN ORDER TO GO THROUGH THE TREE FROM LEAVES TO ROOT.
|
|
|
1757
|
|
|
1758 my @mutations_array=();
|
|
|
1759
|
|
|
1760 my $joint = ();
|
|
|
1761 foreach my $node (@nodes){
|
|
|
1762 my @pair = @$node;
|
|
|
1763 # print "now in the nodes for loop, pair = @pair\n and sequences=\n" if $printer == 1;
|
|
|
1764 $joint = "(".join(", ",@pair).")";
|
|
|
1765 my @pair_sequences=();
|
|
|
1766
|
|
|
1767 foreach my $tag (@pair){
|
|
|
1768 # print "$tag: $node_alignments{$tag}\n" if $printer == 1;
|
|
|
1769 # print $node_alignments{$tag},"\n" if $printer == 1;
|
|
|
1770 push @pair_sequences, $node_alignments{$tag};
|
|
|
1771 }
|
|
|
1772 # print "ppeel onion joint = $joint , pair_sequences=>@pair_sequences< , pair=>@pair<\n" if $printer == 1;
|
|
|
1773
|
|
|
1774 my ($compared, $substitutions_list) = base_by_base_simple($motif,\@pair_sequences, scalar(@pair_sequences), @pair, $joint);
|
|
|
1775 $node_alignments{$joint}=$compared;
|
|
|
1776 push( @mutations_array,split(/:/,$substitutions_list));
|
|
|
1777 # print "newly added to node_sequences: $node_alignments{$joint} and list of mutations =\n", join("\n",@mutations_array),"\n" if $printer == 1;
|
|
|
1778 }
|
|
|
1779
|
|
|
1780
|
|
|
1781 my $analayzed_mutations = analyze_mutations(\@mutations_array, \@nodes, $branches_hash, $alignment, \@tags, \%alivehash, \%node_sequences, $microsatstarts, $motif);
|
|
|
1782
|
|
|
1783 return ($analayzed_mutations, \@nodes, $branches_hash, \%alivehash, $node_alignments{$joint}) if scalar @mutations_array > 0;
|
|
|
1784 return ("NULL",\@nodes,$branches_hash, \%alivehash, "NULL") if scalar @mutations_array == 0;
|
|
|
1785 }
|
|
|
1786
|
|
|
1787 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1788 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1789
|
|
|
1790 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1791 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1792
|
|
|
1793 sub get_nodes{
|
|
|
1794 my $printer = 0;
|
|
|
1795
|
|
|
1796 my $tree=$_[0];
|
|
|
1797 #$tree =~ s/ +//g;
|
|
|
1798 $tree =~ s/\t+//g;
|
|
|
1799 $tree=~s/;//g;
|
|
|
1800 # print "tree=$tree\n" if $printer == 1;
|
|
|
1801 my @nodes = ();
|
|
|
1802 my @onions=($tree);
|
|
|
1803 my %branches=();
|
|
|
1804 foreach my $bite (@onions){
|
|
|
1805 $bite=~ s/^\(|\)$//g;
|
|
|
1806 chomp $bite;
|
|
|
1807 # print "tree = $bite \n";
|
|
|
1808 # <STDIN>;
|
|
|
1809 $bite=~ /([ ,\(\)A-Z]+)\,\s*([ ,\(\)A-Z]+)/;
|
|
|
1810 #$tree =~ /(\(\(\(H, C\), O\), R\))\, (M)/;
|
|
|
1811 my @raw_nodes = ($1, $2);
|
|
|
1812 # print "raw nodes = $1 and $2\n" if $printer == 1;
|
|
|
1813 push(@nodes, [@raw_nodes]);
|
|
|
1814 foreach my $node (@raw_nodes) {push (@onions, $node) if $node =~ /,/;}
|
|
|
1815 foreach my $node (@raw_nodes) {$branches{$node}="(".$bite.")"; }
|
|
|
1816 # print "onions = @onions\n" if $printer == 1;<STDIN> if $printer == 1;
|
|
|
1817 }
|
|
|
1818 $printer = 0;
|
|
|
1819 return \@nodes, \%branches;
|
|
|
1820 }
|
|
|
1821
|
|
|
1822
|
|
|
1823 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1824 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1825 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1826 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1827 sub analyze_mutations{
|
|
|
1828 my ($mutations_array, $nodes, $branches_hash, $alignment, $tags, $alivehash, $node_sequences, $microsatstarts, $motif) = @_;
|
|
|
1829 my $locuslength = length($alignment->{$tags->[0]});
|
|
|
1830 my $printer = 0;
|
|
|
1831
|
|
|
1832
|
|
|
1833 # print " IN analyzed_mutations....\n" if $printer == 1; # \n mutations array = @$mutations_array, \nAND locuslength = $locuslength\n" if $printer == 1;
|
|
|
1834 my %mutation_hash=();
|
|
|
1835 my %froms_megahash=();
|
|
|
1836 my %tos_megahash=();
|
|
|
1837 my %position_hash=();
|
|
|
1838 my @solutions_array=();
|
|
|
1839 foreach my $mutation (@$mutations_array){
|
|
|
1840 # print "loadin mutation: $mutation\n" if $printer == 1;
|
|
|
1841 my %localhash= $mutation =~ /([\S ]+)=([\S ]+)/g;
|
|
|
1842 $mutation_hash{$localhash{"position"}} = {%localhash};
|
|
|
1843 push @{$position_hash{$localhash{"position"}}},$localhash{"node"};
|
|
|
1844 # print "feeding position hash with $localhash{position}: $position_hash{$localhash{position}}[0]\n" if $printer == 1;
|
|
|
1845 $froms_megahash{$localhash{"position"}}{$localhash{"node"}}=$localhash{"from"};
|
|
|
1846 $tos_megahash{$localhash{"position"}}{$localhash{"node"}}=$localhash{"to"};
|
|
|
1847 # print "just a trial: $mutation_hash{$localhash{position}}{position}\n" if $printer == 1;
|
|
|
1848 # print "loadin in tos_megahash: $localhash{position} {$localhash{node} = $localhash{to}\n" if $printer == 1;
|
|
|
1849 # print "loadin in from: $localhash{position} {$localhash{node} = $localhash{from}\n" if $printer == 1;
|
|
|
1850 }
|
|
|
1851
|
|
|
1852 # print "now going through each position in loculength:\n" if $printer == 1;
|
|
|
1853 ## <STDIN> if $printer == 1;
|
|
|
1854
|
|
|
1855 for my $pos (0 ... $locuslength-1){
|
|
|
1856 # print "at position: $pos\n" if $printer == 1;
|
|
|
1857
|
|
|
1858 if (exists($mutation_hash{$pos})){
|
|
|
1859 my @local_nodes=@{$position_hash{$pos}};
|
|
|
1860 # print "found mutation: @{$position_hash{$pos}} : @local_nodes\n" if $printer == 1;
|
|
|
1861
|
|
|
1862 foreach my $local_node (@local_nodes){
|
|
|
1863 # print "at local node: $local_node ... from state = $froms_megahash{$pos}{$local_node}\n" if $printer == 1;
|
|
|
1864 my $open_insertion=();
|
|
|
1865 my $open_deletion=();
|
|
|
1866 my $open_to_substitution=();
|
|
|
1867 my $open_from_substitution=();
|
|
|
1868 if ($froms_megahash{$pos}{$local_node} eq "-"){
|
|
|
1869 # 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}};;
|
|
|
1870 # 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;
|
|
|
1871 #next if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
|
|
|
1872 $open_insertion=$tos_megahash{$pos}{$local_node};
|
|
|
1873 for my $posnext ($pos+1 ... $locuslength-1){
|
|
|
1874 # print "in first if .... studying posnext: $posnext\n" if $printer == 1;
|
|
|
1875 last if !exists ($froms_megahash{$posnext}{$local_node});
|
|
|
1876 # 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;
|
|
|
1877 $open_insertion = $open_insertion.$tos_megahash{$posnext}{$local_node} if $froms_megahash{$posnext}{$local_node} eq "-";
|
|
|
1878 # print "now open_insertion=$open_insertion\n" if $printer == 1;
|
|
|
1879 delete $mutation_hash{$posnext} if $froms_megahash{$posnext}{$local_node} eq "-";
|
|
|
1880 }
|
|
|
1881 # print "1 Feeding in: ", join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_insertion", "deletion="),"\n" if $printer == 1;
|
|
|
1882 push (@solutions_array, join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_insertion", "deletion="));
|
|
|
1883 }
|
|
|
1884 elsif ($tos_megahash{$pos}{$local_node} eq "-"){
|
|
|
1885 # 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}};;
|
|
|
1886 # 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}};
|
|
|
1887 #next if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
|
|
|
1888 $open_deletion=$froms_megahash{$pos}{$local_node};
|
|
|
1889 for my $posnext ($pos+1 ... $locuslength-1){
|
|
|
1890 # print "in 1st elsif studying posnext: $posnext\n" if $printer == 1;
|
|
|
1891 # print "nexting as nextpos does not exist\n" if !exists ($tos_megahash{$posnext}{$local_node}) && $printer == 1;
|
|
|
1892 last if !exists ($tos_megahash{$posnext}{$local_node});
|
|
|
1893 # print "for posnext: $posnext, there exists $tos_megahash{$posnext}{$local_node}\n" if $printer == 1;
|
|
|
1894 $open_deletion = $open_deletion.$froms_megahash{$posnext}{$local_node} if $tos_megahash{$posnext}{$local_node} eq "-";
|
|
|
1895 delete $mutation_hash{$posnext} if $tos_megahash{$posnext}{$local_node} eq "-";
|
|
|
1896 }
|
|
|
1897 # print "2 Feeding in:", join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_deletion"), "\n" if $printer == 1;
|
|
|
1898 push (@solutions_array, join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_deletion"));
|
|
|
1899 }
|
|
|
1900 elsif ($tos_megahash{$pos}{$local_node} ne "-"){
|
|
|
1901 # 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}};;
|
|
|
1902 # 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}};
|
|
|
1903 #next if exists $alivehash->{$local_node} && exists $alivehash->{$branches_hash->{$local_node}};
|
|
|
1904 # print "microsatstart = $microsatstarts->{$local_node} \n" if exists $microsatstarts->{$local_node} && $pos < $microsatstarts->{$local_node} && $printer == 1;
|
|
|
1905 next if exists $microsatstarts->{$local_node} && $pos < $microsatstarts->{$local_node};
|
|
|
1906 $open_to_substitution=$tos_megahash{$pos}{$local_node};
|
|
|
1907 $open_from_substitution=$froms_megahash{$pos}{$local_node};
|
|
|
1908 # print "open from substitution: $open_from_substitution \n" if $printer == 1;
|
|
|
1909 for my $posnext ($pos+1 ... $locuslength-1){
|
|
|
1910 #print "in last elsif studying posnext: $posnext\n";
|
|
|
1911 last if !exists ($tos_megahash{$posnext}{$local_node});
|
|
|
1912 # print "for posnext: $posnext, there exists $tos_megahash{$posnext}{$local_node}\n" if $printer == 1;
|
|
|
1913 $open_to_substitution = $open_to_substitution.$tos_megahash{$posnext}{$local_node} if $tos_megahash{$posnext}{$local_node} ne "-";
|
|
|
1914 $open_from_substitution = $open_from_substitution.$froms_megahash{$posnext}{$local_node} if $tos_megahash{$posnext}{$local_node} ne "-";
|
|
|
1915 delete $mutation_hash{$posnext} if $tos_megahash{$posnext}{$local_node} ne "-" && $froms_megahash{$posnext}{$local_node} ;
|
|
|
1916 }
|
|
|
1917 # print "open from substitution: $open_from_substitution \n" if $printer == 1;
|
|
|
1918
|
|
|
1919 #IS THE STRETCH OF SUBSTITUTION MICROSATELLITE-LIKE?
|
|
|
1920 my @motif_parts=split(/\s*/,$motif);
|
|
|
1921 #GENERATING THE FLEXIBLE LEFT END
|
|
|
1922 my $left_query=();
|
|
|
1923 for my $k (1 ... $#motif_parts) {
|
|
|
1924 $left_query= $motif_parts[$k]."|)";
|
|
|
1925 $left_query="(".$left_query;
|
|
|
1926 }
|
|
|
1927 $left_query=$left_query."?";
|
|
|
1928 # print "left_quewry = $left_query\n" if $printer == 1;
|
|
|
1929 #GENERATING THE FLEXIBLE RIGHT END
|
|
|
1930 my $right_query=();
|
|
|
1931 for my $k (0 ... ($#motif_parts-1)) {
|
|
|
1932 $right_query= "(|".$motif_parts[$k];
|
|
|
1933 $right_query=$right_query.")";
|
|
|
1934 }
|
|
|
1935 $right_query=$right_query."?";
|
|
|
1936 # print "right_query = $right_query\n" if $printer == 1;
|
|
|
1937 # print "Hence, searching for: ^$left_query($motif)+$right_query\$\n" if $printer == 1;
|
|
|
1938
|
|
|
1939 my $motifcomb=$motif x 50;
|
|
|
1940 # print "motifcomb = $motifcomb\n" if $printer == 1;
|
|
|
1941 if ( ($motifcomb =~/$open_to_substitution/i) && (length ($open_to_substitution) >= length($motif)) ){
|
|
|
1942 # print "sequence microsat-like\n" if $printer == 1;
|
|
|
1943 my $all_microsat_like = 0;
|
|
|
1944 # print "3 feeding in: ", join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_from_substitution"), "\n" if $printer == 1;
|
|
|
1945 push (@solutions_array, join("\t", "node=$local_node","type=deletion" ,"position=$pos", "from=", "to=", "insertion=", "deletion=$open_from_substitution"));
|
|
|
1946 # print "4 feeding in: ", join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_to_substitution", "deletion="), "\n" if $printer == 1;
|
|
|
1947 push (@solutions_array, join("\t", "node=$local_node","type=insertion" ,"position=$pos", "from=", "to=", "insertion=$open_to_substitution", "deletion="));
|
|
|
1948
|
|
|
1949 }
|
|
|
1950 else{
|
|
|
1951 # 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;
|
|
|
1952 push (@solutions_array, join("\t", "node=$local_node","type=substitution" ,"position=$pos", "from=$open_from_substitution", "to=$open_to_substitution", "insertion=", "deletion="));
|
|
|
1953 }
|
|
|
1954 #IS THE FROM-SEQUENCE MICROSATELLITE-LIKE?
|
|
|
1955
|
|
|
1956 }
|
|
|
1957 #<STDIN> if $printer ==1;
|
|
|
1958 }
|
|
|
1959 #<STDIN> if $printer ==1;
|
|
|
1960 }
|
|
|
1961 }
|
|
|
1962 # print "\n", "#" x 50, "\n" if $printer == 1;
|
|
|
1963 foreach my $tag (@$tags){
|
|
|
1964 # print "$tag: $alignment->{$tag}\n" if $printer == 1;
|
|
|
1965 }
|
|
|
1966 # print "\n", "#" x 50, "\n" if $printer == 1;
|
|
|
1967 # print "returning SOLUTIONS ARRAY : \n",join("\n", @solutions_array),"\n" if $printer == 1;
|
|
|
1968 #print "end\n";
|
|
|
1969 #<STDIN> if
|
|
|
1970 return \@solutions_array;
|
|
|
1971 }
|
|
|
1972 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1973 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1974 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1975 #+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#+++++++++++#
|
|
|
1976
|
|
|
1977 sub base_by_base_simple{
|
|
|
1978 my $printer = 0;
|
|
|
1979 my ($motif, $locus, $no, $pair0, $pair1, $joint) = @_;
|
|
|
1980 my @seq_array=();
|
|
|
1981 # print "IN SUBROUTUNE base_by_base_simple.. information received = @_\n" if $printer == 1;
|
|
|
1982 # print "pair0 = $pair0 and pair1 = $pair1\n" if $printer == 1;
|
|
|
1983
|
|
|
1984 my @example=split(/\./,$locus->[0]);
|
|
|
1985 # print "example, for length = @example\n" if $printer == 1;
|
|
|
1986 for my $i (0...$no-1){push(@seq_array, [split(/\./,$locus->[$i])]); }
|
|
|
1987
|
|
|
1988 my @compared_sequence=();
|
|
|
1989 my @substitutions_list;
|
|
|
1990 for my $i (0...scalar(@example)-1){
|
|
|
1991
|
|
|
1992 #print "i = $i\n" if $printer == 1;
|
|
|
1993 #print "comparing $seq_array[0][$i] and $seq_array[1][$i] \n" ;#if $printer == 1;
|
|
|
1994 if ($seq_array[0][$i] =~ /!/ && $seq_array[1][$i] !~ /!/){
|
|
|
1995
|
|
|
1996 my $resolution= resolve_base($seq_array[0][$i],$seq_array[1][$i], $pair1 ,"keep" );
|
|
|
1997 # print "ancestral = $resolution\n" if $printer == 1;
|
|
|
1998
|
|
|
1999 if ($resolution =~ /$seq_array[1][$i]/i && $resolution !~ /!/){
|
|
|
2000 push @substitutions_list, add_mutation($i, $pair0, $seq_array[0][$i], $resolution );
|
|
|
2001 }
|
|
|
2002 elsif ( $resolution !~ /!/){
|
|
|
2003 push @substitutions_list, add_mutation($i, $pair1, $seq_array[1][$i], $resolution);
|
|
|
2004 }
|
|
|
2005 push @compared_sequence,$resolution;
|
|
|
2006 }
|
|
|
2007 elsif ($seq_array[0][$i] !~ /!/ && $seq_array[1][$i] =~ /!/){
|
|
|
2008
|
|
|
2009 my $resolution= resolve_base($seq_array[1][$i],$seq_array[0][$i], $pair0, "invert" );
|
|
|
2010 # print "ancestral = $resolution\n" if $printer == 1;
|
|
|
2011
|
|
|
2012 if ($resolution =~ /$seq_array[0][$i]/i && $resolution !~ /!/){
|
|
|
2013 push @substitutions_list, add_mutation($i, $pair1, $seq_array[1][$i], $resolution);
|
|
|
2014 }
|
|
|
2015 elsif ( $resolution !~ /!/){
|
|
|
2016 push @substitutions_list, add_mutation($i, $pair0, $seq_array[0][$i], $resolution);
|
|
|
2017 }
|
|
|
2018 push @compared_sequence,$resolution;
|
|
|
2019 }
|
|
|
2020 elsif($seq_array[0][$i] =~ /!/ && $seq_array[1][$i] =~ /!/){
|
|
|
2021 push @compared_sequence, add_bases($seq_array[0][$i],$seq_array[1][$i], $pair0, $pair1, $joint );
|
|
|
2022 }
|
|
|
2023 else{
|
|
|
2024 if($seq_array[0][$i] !~ /^$seq_array[1][$i]$/i){
|
|
|
2025 push @compared_sequence, $pair0.":".$seq_array[0][$i]."!".$pair1.":".$seq_array[1][$i];
|
|
|
2026 }
|
|
|
2027 else{
|
|
|
2028 # print "perfect match\n" if $printer == 1;
|
|
|
2029 push @compared_sequence, $seq_array[0][$i];
|
|
|
2030 }
|
|
|
2031 }
|
|
|
2032 }
|
|
|
2033 # print "returning: comared = @compared_sequence \nand substitutions list =\n", join("\n",@substitutions_list),"\n" if $printer == 1;
|
|
|
2034 return join(".",@compared_sequence), join(":", @substitutions_list) if scalar (@substitutions_list) > 0;
|
|
|
2035 return join(".",@compared_sequence), "" if scalar (@substitutions_list) == 0;
|
|
|
2036 }
|
|
|
2037
|
|
|
2038
|
|
|
2039 sub resolve_base{
|
|
|
2040 my $printer = 0;
|
|
|
2041 # print "IN SUBROUTUNE resolve_base.. information received = @_\n" if $printer == 1;
|
|
|
2042 my ($optional, $single, $singlesp, $arg) = @_;
|
|
|
2043 my @options=split(/!/,$optional);
|
|
|
2044 foreach my $option(@options) {
|
|
|
2045 $option=~s/[A-Z\(\) ,]+://g;
|
|
|
2046 if ($option =~ /$single/i){
|
|
|
2047 # print "option = $option , returning single: $single\n" if $printer == 1;
|
|
|
2048 return $single;
|
|
|
2049 }
|
|
|
2050 }
|
|
|
2051 # print "returning ",$optional."!".$singlesp.":".$single. "\n" if $arg eq "keep" && $printer == 1;
|
|
|
2052 # print "returning ",$singlesp.":".$single."!".$optional. "\n" if $arg eq "invert" && $printer == 1;
|
|
|
2053 return $optional."!".$singlesp.":".$single if $arg eq "keep";
|
|
|
2054 return $singlesp.":".$single."!".$optional if $arg eq "invert";
|
|
|
2055
|
|
|
2056 }
|
|
|
2057
|
|
|
2058 sub same_length{
|
|
|
2059 my $printer = 0;
|
|
|
2060 my @locus = @_;
|
|
|
2061 my $temp = shift @locus;
|
|
|
2062 $temp=~s/-|,//g;
|
|
|
2063 foreach my $l (@locus){
|
|
|
2064 $l=~s/-|,//g;
|
|
|
2065 return 0 if length($l) != length($temp);
|
|
|
2066 $temp = $l;
|
|
|
2067 }
|
|
|
2068 return 1;
|
|
|
2069 }
|
|
|
2070 sub treeStudy{
|
|
|
2071 my $printer = 1;
|
|
|
2072 # print "template DEFINED.. received: @_\n" if defined %template;
|
|
|
2073 # print "only received = @_" if !defined %template;
|
|
|
2074 my $stopper = 0;
|
|
|
2075 # if (!defined %template){ TEMP MASKED OCT 18 2012
|
|
|
2076 $stopper = 1;
|
|
|
2077 %template=();
|
|
|
2078 # print "tree decipherer = $tree_decipherer\n" if $printer == 1;
|
|
|
2079 my ( $template_ref, $keys_array)=load_allPossibleTrees($tree_decipherer, \%template);
|
|
|
2080 # print "return = $template_ref and @{$keys_array}\n" if $printer == 1;
|
|
|
2081 foreach my $key (@$keys_array){
|
|
|
2082 # print "addding : $template_ref->{$key} for $key\n" if $printer == 1;
|
|
|
2083 $template{$key} = $template_ref->{$key};
|
|
|
2084 }
|
|
|
2085 # } TEMP MASK OCT 18 2012 END
|
|
|
2086 # <STDIN>;
|
|
|
2087 for my $templet ( keys %template ) {
|
|
|
2088 # print "$templet => @{$template{$templet}}\n";
|
|
|
2089 }
|
|
|
2090 # <STDIN> if !defined %template;
|
|
|
2091
|
|
|
2092 my $strict = 0;
|
|
|
2093
|
|
|
2094 my $H = 0;
|
|
|
2095 my $Hchr = 1;
|
|
|
2096 my $Hstart = 2;
|
|
|
2097 my $Hend = 3;
|
|
|
2098 my $Hmotif = 4;
|
|
|
2099 my $Hmotiflen = 5;
|
|
|
2100 my $Hmicro = 6;
|
|
|
2101 my $Hstrand = 7;
|
|
|
2102 my $Hmicrolen = 8;
|
|
|
2103 my $Hinterpos = 9;
|
|
|
2104 my $Hrelativepos = 10;
|
|
|
2105 my $Hinter = 11;
|
|
|
2106 my $Hinterlen = 12;
|
|
|
2107
|
|
|
2108 my $C = 13;
|
|
|
2109 my $Cchr = 14;
|
|
|
2110 my $Cstart = 15;
|
|
|
2111 my $Cend = 16;
|
|
|
2112 my $Cmotif = 17;
|
|
|
2113 my $Cmotiflen = 18;
|
|
|
2114 my $Cmicro = 19;
|
|
|
2115 my $Cstrand = 20;
|
|
|
2116 my $Cmicrolen = 21;
|
|
|
2117 my $Cinterpos = 22;
|
|
|
2118 my $Crelativepos = 23;
|
|
|
2119 my $Cinter = 24;
|
|
|
2120 my $Cinterlen = 25;
|
|
|
2121
|
|
|
2122 my $O = 26;
|
|
|
2123 my $Ochr = 27;
|
|
|
2124 my $Ostart = 28;
|
|
|
2125 my $Oend = 29;
|
|
|
2126 my $Omotif = 30;
|
|
|
2127 my $Omotiflen = 31;
|
|
|
2128 my $Omicro = 32;
|
|
|
2129 my $Ostrand = 33;
|
|
|
2130 my $Omicrolen = 34;
|
|
|
2131 my $Ointerpos = 35;
|
|
|
2132 my $Orelativepos = 36;
|
|
|
2133 my $Ointer = 37;
|
|
|
2134 my $Ointerlen = 38;
|
|
|
2135
|
|
|
2136 my $R = 39;
|
|
|
2137 my $Rchr = 40;
|
|
|
2138 my $Rstart = 41;
|
|
|
2139 my $Rend = 42;
|
|
|
2140 my $Rmotif = 43;
|
|
|
2141 my $Rmotiflen = 44;
|
|
|
2142 my $Rmicro = 45;
|
|
|
2143 my $Rstrand = 46;
|
|
|
2144 my $Rmicrolen = 47;
|
|
|
2145 my $Rinterpos = 48;
|
|
|
2146 my $Rrelativepos = 49;
|
|
|
2147 my $Rinter = 50;
|
|
|
2148 my $Rinterlen = 51;
|
|
|
2149
|
|
|
2150 my $Mchr = 52;
|
|
|
2151 my $Mstart = 53;
|
|
|
2152 my $Mend = 54;
|
|
|
2153 my $M = 55;
|
|
|
2154 my $Mmotif = 56;
|
|
|
2155 my $Mmotiflen = 57;
|
|
|
2156 my $Mmicro = 58;
|
|
|
2157 my $Mstrand = 59;
|
|
|
2158 my $Mmicrolen = 60;
|
|
|
2159 my $Minterpos = 61;
|
|
|
2160 my $Mrelativepos = 62;
|
|
|
2161 my $Minter = 63;
|
|
|
2162 my $Minterlen = 64;
|
|
|
2163
|
|
|
2164 #-------------------------------------------------------------------------------#
|
|
|
2165 my @analysis=();
|
|
|
2166
|
|
|
2167
|
|
|
2168 my %speciesOrder = ();
|
|
|
2169 $speciesOrder{"H"} = 0;
|
|
|
2170 $speciesOrder{"C"} = 1;
|
|
|
2171 $speciesOrder{"O"} = 2;
|
|
|
2172 $speciesOrder{"R"} = 3;
|
|
|
2173 $speciesOrder{"M"} = 4;
|
|
|
2174 #-------------------------------------------------------------------------------#
|
|
|
2175
|
|
|
2176 my $line = $_[0];
|
|
|
2177 chomp $line;
|
|
|
2178
|
|
|
2179 my @f = split(/\t/,$line);
|
|
|
2180 # print "received array : @f.. recieved tags = @tags\n" if $printer == 1;
|
|
|
2181
|
|
|
2182 # collect all motifs
|
|
|
2183 my @motifs=();
|
|
|
2184 @motifs = ($f[$Hmotif], $f[$Cmotif], $f[$Omotif], $f[$Rmotif], $f[$Mmotif]) if $tags[$#tags] =~ /M/;
|
|
|
2185 @motifs = ($f[$Hmotif], $f[$Cmotif], $f[$Omotif], $f[$Rmotif]) if $tags[$#tags] =~ /R/;
|
|
|
2186 @motifs = ($f[$Hmotif], $f[$Cmotif], $f[$Omotif]) if $tags[$#tags] =~ /O/;
|
|
|
2187 # print "motifs in the array = $f[$Hmotif], $f[$Cmotif], $f[$Omotif], $f[$Rmotif]\n" if $tags[$#tags] =~ /R/;;
|
|
|
2188 # print "motifs = @motifs\n" if $printer == 1;
|
|
|
2189 my @translation = ();
|
|
|
2190 foreach my $motif (@motifs){
|
|
|
2191 push(@translation, "_") if $motif eq "NA";
|
|
|
2192 push(@translation, "+") if $motif ne "NA";
|
|
|
2193 }
|
|
|
2194 my $translate = join(" ", @translation);
|
|
|
2195 # print "translate = >$translate< and analysis = $template{$translate}[0].. on the other hand, ",$template{"- - +"}[0],"\n";
|
|
|
2196 my @analyses = split(/\|/,$template{$translate}[0]);
|
|
|
2197 # print "motifs = @motifs, analyses = @analyses\n" if $printer == 1;
|
|
|
2198
|
|
|
2199 if (scalar(@analyses) == 1) {
|
|
|
2200 #print "analysis = $analyses[0]\n";
|
|
|
2201 if ($analyses[0] !~ /,|\./ ){
|
|
|
2202 if ($analyses[0] =~ /\+/){
|
|
|
2203 my $analysis = $analyses[0];
|
|
|
2204 $analysis =~ s/\+|\-//g;
|
|
|
2205 my @species = split(/\s*/,$analysis);
|
|
|
2206 my @currentMotifs = ();
|
|
|
2207 foreach my $specie (@species){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); #print "pushing into currentMotifs: $speciesOrder{$specie}: $motifs[$speciesOrder{$specie}]\n" if $printer == 1;
|
|
|
2208 }
|
|
|
2209 # print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
|
|
|
2210 $template{$translate}[1]++ if $strict == 1 && consistency(@currentMotifs) ne "NULL";
|
|
|
2211 $template{$translate}[1]++ if $strict == 0;
|
|
|
2212 # print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
|
|
|
2213 }
|
|
|
2214 else{
|
|
|
2215 my $analysis = $analyses[0];
|
|
|
2216 $analysis =~ s/\+|\-//g;
|
|
|
2217 my @species = split(/\s*/,$analysis);
|
|
|
2218 my @currentMotifs = ();
|
|
|
2219 my @complementarySpecies = ();
|
|
|
2220 my $allSpecies = join("",@tags);
|
|
|
2221 foreach my $specie (@species){ $allSpecies =~ s/$specie//g; }
|
|
|
2222 foreach my $specie (split(/\s*/,$allSpecies)){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); #print "pushing into currentMotifs: $speciesOrder{$specie}: $motifs[$speciesOrder{$specie}]\n" if $printer == 1;;
|
|
|
2223 }
|
|
|
2224 # print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
|
|
|
2225 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
|
|
|
2226 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
|
|
|
2227 # print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
|
|
|
2228 }
|
|
|
2229 }
|
|
|
2230
|
|
|
2231 elsif ($analyses[0] =~ /,/) {
|
|
|
2232 my @events = split(/,/,$analyses[0]);
|
|
|
2233 # print "events = @events \n " if $printer == 1;
|
|
|
2234 if ($events[0] =~ /\+/){
|
|
|
2235 my $analysis1 = $events[0];
|
|
|
2236 $analysis1 =~ s/\+|\-//g;
|
|
|
2237 my $analysis2 = $events[1];
|
|
|
2238 $analysis2 =~ s/\+|\-//g;
|
|
|
2239 my @nSpecies = split(/\s*/,$analysis2);
|
|
|
2240 # print "original anslysis = $analysis1 " if $printer == 1;
|
|
|
2241 foreach my $specie (@nSpecies){ $analysis1=~ s/$specie//g;}
|
|
|
2242 # print "processed anslysis = $analysis1 \n" if $printer == 1;
|
|
|
2243 my @currentMotifs = ();
|
|
|
2244 foreach my $specie (split(/\s*/,$analysis1)){push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
|
|
|
2245 # print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
|
|
|
2246 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
|
|
|
2247 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
|
|
|
2248 # print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
|
|
|
2249 }
|
|
|
2250 else{
|
|
|
2251 my $analysis1 = $events[0];
|
|
|
2252 $analysis1 =~ s/\+|\-//g;
|
|
|
2253 my $analysis2 = $events[1];
|
|
|
2254 $analysis2 =~ s/\+|\-//g;
|
|
|
2255 my @pSpecies = split(/\s*/,$analysis2);
|
|
|
2256 my @currentMotifs = ();
|
|
|
2257 foreach my $specie (@pSpecies){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
|
|
|
2258 # print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
|
|
|
2259 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
|
|
|
2260 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
|
|
|
2261 # print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
|
|
|
2262
|
|
|
2263 }
|
|
|
2264
|
|
|
2265 }
|
|
|
2266 elsif ($analyses[0] =~ /\./) {
|
|
|
2267 my @events = split(/\./,$analyses[0]);
|
|
|
2268 foreach my $event (@events){
|
|
|
2269 # print "event = $event \n" if $printer == 1;
|
|
|
2270 if ($event =~ /\+/){
|
|
|
2271 my $analysis = $event;
|
|
|
2272 $analysis =~ s/\+|\-//g;
|
|
|
2273 my @species = split(/\s*/,$analysis);
|
|
|
2274 my @currentMotifs = ();
|
|
|
2275 foreach my $specie (@species){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
|
|
|
2276 #print consistency(@currentMotifs),"<- \n";
|
|
|
2277 # print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
|
|
|
2278 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
|
|
|
2279 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
|
|
|
2280 # print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
|
|
|
2281 }
|
|
|
2282 else{
|
|
|
2283 my $analysis = $event;
|
|
|
2284 $analysis =~ s/\+|\-//g;
|
|
|
2285 my @species = split(/\s*/,$analysis);
|
|
|
2286 my @currentMotifs = ();
|
|
|
2287 my @complementarySpecies = ();
|
|
|
2288 my $allSpecies = join("",@tags);
|
|
|
2289 foreach my $specie (@species){ $allSpecies =~ s/$specie//g; }
|
|
|
2290 foreach my $specie (split(/\s*/,$allSpecies)){ push(@currentMotifs, $motifs[$speciesOrder{$specie}]); }
|
|
|
2291 #print consistency(@currentMotifs),"<- \n";
|
|
|
2292 # print "current motifs = @currentMotifs and consistency? ", (consistency(@currentMotifs))," \n" if $printer == 1;
|
|
|
2293 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 1 && consistency(@currentMotifs) ne "NULL";
|
|
|
2294 $template{$translate}[1]=$template{$translate}[1]+1 if $strict == 0;
|
|
|
2295 # print "adding to template $translate: $template{$translate}[1]\n" if $printer == 1;
|
|
|
2296 }
|
|
|
2297 }
|
|
|
2298
|
|
|
2299 }
|
|
|
2300 }
|
|
|
2301 else{
|
|
|
2302 my $finalanalysis = ();
|
|
|
2303 $template{$translate}[1]++;
|
|
|
2304 foreach my $analysis (@analyses){ ;}
|
|
|
2305 }
|
|
|
2306 # test if motifs where microsats are present, as indeed of same the motif composition
|
|
|
2307
|
|
|
2308
|
|
|
2309
|
|
|
2310 for my $templet ( keys %template ) {
|
|
|
2311 if (@{ $template{$templet} }[1] > 0){
|
|
|
2312
|
|
|
2313 $template{$templet}[1] = 0;
|
|
|
2314 # print "now returning: @{$template{$templet}}[0], $templet\n";
|
|
|
2315 return (@{$template{$templet}}[0], $templet);
|
|
|
2316 }
|
|
|
2317 }
|
|
|
2318 undef %template;
|
|
|
2319 # print "sending NULL\n" if $printer == 1;
|
|
|
2320 return ("NULL", "NULL");
|
|
|
2321
|
|
|
2322 }
|
|
|
2323
|
|
|
2324
|
|
|
2325 sub consistency{
|
|
|
2326 my @motifs = @_;
|
|
|
2327 # print "in consistency \n" if $printer == 1;
|
|
|
2328 # print "motifs sent = >",join("|",@motifs),"< \n" if $printer == 1;
|
|
|
2329 return $motifs[0] if scalar(@motifs) == 1;
|
|
|
2330 my $prevmotif = shift(@motifs);
|
|
|
2331 my $stopper = 0;
|
|
|
2332 for my $i (0 ... $#motifs){
|
|
|
2333 next if $motifs[$i] eq "NA";
|
|
|
2334 my $templet = $motifs[$i].$motifs[$i];
|
|
|
2335 if ($templet !~ /$prevmotif/i){
|
|
|
2336 $stopper = 1; last;
|
|
|
2337 }
|
|
|
2338 }
|
|
|
2339 return $prevmotif if $stopper == 0;
|
|
|
2340 return "NULL" if $stopper == 1;
|
|
|
2341 }
|
|
|
2342 sub summarize_microsat{
|
|
|
2343 my $printer = 1;
|
|
|
2344 my $line = $_[0];
|
|
|
2345 my $humseq = $_[1];
|
|
|
2346
|
|
|
2347 my @gaps = $line =~ /[0-9]+\t[0-9]+\t[\+\-]/g;
|
|
|
2348 my @starts = $line =~ /[0-9]+\t[\+\-]/g;
|
|
|
2349 my @ends = $line =~ /[\+\-]\t[0-9]+/g;
|
|
|
2350 # print "starts = @starts\tends = @ends\n" if $printer == 1;
|
|
|
2351 for my $i (0 ... $#gaps) {$gaps[$i] =~ s/\t[0-9]+\t[\+\-]//g;}
|
|
|
2352 for my $i (0 ... $#starts) {$starts[$i] =~ s/\t[\+\-]//g;}
|
|
|
2353 for my $i (0 ... $#ends) {$ends[$i] =~ s/[\+\-]\t//g;}
|
|
|
2354
|
|
|
2355 my $minstart = array_smallest_number(@starts);
|
|
|
2356 my $maxend = array_largest_number(@ends);
|
|
|
2357
|
|
|
2358 my $humupstream_st = substr($humseq, 0, $minstart);
|
|
|
2359 my $humupstream_en = substr($humseq, 0, $maxend);
|
|
|
2360 my $no_of_gaps_to_start = 0;
|
|
|
2361 my $no_of_gaps_to_end = 0;
|
|
|
2362 $no_of_gaps_to_start = ($humupstream_st =~ s/\-/x/g) if $humupstream_st=~/\-/;
|
|
|
2363 $no_of_gaps_to_end = ($humupstream_en =~ s/\-/x/g) if $humupstream_en=~/\-/;
|
|
|
2364
|
|
|
2365 my $locusmotif = ();
|
|
|
2366 # print "IN SUB SUMMARIZE_MICROSAT $line\n" if $printer == 1;
|
|
|
2367 #return "NULL" if $line =~ /compound/;
|
|
|
2368 my $Hstart = "NA";
|
|
|
2369 my $Hend = "NA";
|
|
|
2370 chomp $line;
|
|
|
2371 my $match_count = ($line =~ s/>/>/g);
|
|
|
2372 #print "number of species = $match_count\n";
|
|
|
2373 my @micros = split(/>/,$line);
|
|
|
2374 shift @micros;
|
|
|
2375 my $stopper = 0;
|
|
|
2376
|
|
|
2377
|
|
|
2378 foreach my $mic (@micros){
|
|
|
2379 my @local = split(/\t/,$mic);
|
|
|
2380 if ($local[$microsatcord] =~ /N/) {$stopper =1; last;}
|
|
|
2381 }
|
|
|
2382 return "NULL" if $stopper ==1;
|
|
|
2383
|
|
|
2384 #------------------------------------------------------
|
|
|
2385
|
|
|
2386 my @arranged = ();
|
|
|
2387 for my $arr (0 ... $#exacttags) {$arranged[$arr] = '0';}
|
|
|
2388
|
|
|
2389 foreach my $micro (@micros){
|
|
|
2390 for my $i (0 ... $#exacttags){
|
|
|
2391 if ($micro =~ /^$exacttags[$i]/){
|
|
|
2392 $arranged[$i] = $micro;
|
|
|
2393 last;
|
|
|
2394 }
|
|
|
2395 }
|
|
|
2396 }
|
|
|
2397 # print "arranged = @arranged \n" ; <STDIN>;;
|
|
|
2398
|
|
|
2399 my @endstatement = ();
|
|
|
2400 my $turn = 0;
|
|
|
2401 my $species_counter = 0;
|
|
|
2402 # print scalar(@arranged),"\n";
|
|
|
2403
|
|
|
2404 my $species_no=0;
|
|
|
2405
|
|
|
2406 my $orthHchr = 0;
|
|
|
2407
|
|
|
2408 foreach my $micro (@arranged) {
|
|
|
2409 $micro =~ s/\t\t/\t \t/g;
|
|
|
2410 $micro =~ s/\t,/\t ,/g;
|
|
|
2411 $micro =~ s/,\t/, \t/g;
|
|
|
2412 # print "------------------------------------------------------------------------------------------\n" if $printer == 1;
|
|
|
2413 chomp $micro;
|
|
|
2414 if ($micro eq '0'){
|
|
|
2415 push(@endstatement, join("\t",$exacttags[$species_counter],"NA","NA","NA","NA",0 ,"NA", "NA", 0,"NA","NA","NA", "NA" ));
|
|
|
2416 $species_counter++;
|
|
|
2417 # print join("|","ENDSTATEMENT:",@endstatement),"\n" if $printer == 1;
|
|
|
2418 next;
|
|
|
2419 }
|
|
|
2420 # print $micro,"\n";
|
|
|
2421 # print "micro = $micro \n" if $printer == 1;
|
|
|
2422 my @fields = split(/\t/,$micro);
|
|
|
2423 my $microcopy = $fields[$microsatcord];
|
|
|
2424 $microcopy =~ s/\[|\]|-//g;
|
|
|
2425 my $microsatlength = length($microcopy);
|
|
|
2426 # print "microsat = $fields[$microsatcord] and microsatlength = $microsatlength\n" if $printer == 1;
|
|
|
2427 # print "sp_ident = @sp_ident.. species_no=$species_no\n";
|
|
|
2428 $micro =~ /$sp_ident[$species_no]\s([0-9a-zA-Z_]+)\s([0-9]+)\s([0-9]+)/;
|
|
|
2429 # print "$micro =~ /$sp_ident[$species_no] ([0-9a-zA-Z_]+) ([0-9]+) ([0-9]+)/\n";
|
|
|
2430 my $sp_chr=$1;
|
|
|
2431 my $sp_start=$2 + $fields[$startcord] - $fields[$gapcord];
|
|
|
2432 my $sp_end= $sp_start + $microsatlength - 1;
|
|
|
2433
|
|
|
2434 $species_no++;
|
|
|
2435
|
|
|
2436 $micro =~ /$focalspec_orig\s(\S+)\s([0-9]+)\s([0-9]+)/;
|
|
|
2437 $orthHchr=$1;
|
|
|
2438 $Hstart=$2+$minstart-$no_of_gaps_to_start;
|
|
|
2439 $Hend=$2+$maxend-$no_of_gaps_to_end;
|
|
|
2440 # print "Hstart = $Hstart = $fields[4] + $fields[$startcord] - $fields[$gapcord]\n" if $printer == 1;
|
|
|
2441
|
|
|
2442 my $motif = $fields[$motifcord];
|
|
|
2443 my $firstmotif = ();
|
|
|
2444 my $strand = $fields[$strandcord];
|
|
|
2445 # print "strand = $strand\n";
|
|
|
2446
|
|
|
2447
|
|
|
2448 if ($motif =~ /^\[/){
|
|
|
2449 $motif =~ s/^\[//g;
|
|
|
2450 $motif =~ /([a-zA-Z]+)\].*/;
|
|
|
2451 $firstmotif = $1;
|
|
|
2452 }
|
|
|
2453
|
|
|
2454 else {$firstmotif = $motif;}
|
|
|
2455 # print "firstmotif =$firstmotif : \n" if $printer == 1;
|
|
|
2456 $firstmotif = allCaps($firstmotif);
|
|
|
2457
|
|
|
2458 if (exists $revHash{$firstmotif} && $turn == 0) {
|
|
|
2459 $turn=1 if $species_counter==0;
|
|
|
2460 $firstmotif = $revHash{$firstmotif};
|
|
|
2461 }
|
|
|
2462
|
|
|
2463 elsif (exists $revHash{$firstmotif} && $turn == 1) {$firstmotif = $revHash{$firstmotif}; $turn = 1;}
|
|
|
2464 # print "changed firstmotif =$firstmotif\n" if $printer == 1;
|
|
|
2465 # <STDIN>;
|
|
|
2466 $locusmotif = $firstmotif;
|
|
|
2467
|
|
|
2468 if (scalar(@fields) > $microsatcord + 2){
|
|
|
2469 # print "fields = @fields ... interr_poscord=$interr_poscord=$fields[$interr_poscord] .. interrcord=$interrcord=$fields[$interrcord]\n" if $printer == 1;
|
|
|
2470
|
|
|
2471 my @interposes = ();
|
|
|
2472 @interposes = split(",",$fields[$interr_poscord]) if $fields[$interr_poscord] =~ /,/;
|
|
|
2473 $interposes[0] = $fields[$interr_poscord] if $fields[$interr_poscord] !~ /,/ ;
|
|
|
2474 # print "interposes=@interposes\n" if $printer == 1;
|
|
|
2475 my @relativeposes = ();
|
|
|
2476 my @interruptions = ();
|
|
|
2477 @interruptions = split(",",$fields[$interrcord]) if $fields[$interrcord] =~ /,/;
|
|
|
2478 $interruptions[0] = $fields[$interrcord] if $fields[$interrcord] !~ /,/;
|
|
|
2479 my @interlens = ();
|
|
|
2480
|
|
|
2481
|
|
|
2482 for my $i (0 ... $#interposes){
|
|
|
2483
|
|
|
2484 my $interpos = $interposes[$i];
|
|
|
2485 my $nexter = 0;
|
|
|
2486 my $interruption = $interruptions[$i];
|
|
|
2487 my $interlen = length($interruption);
|
|
|
2488 push (@interlens, $interlen);
|
|
|
2489
|
|
|
2490
|
|
|
2491 my $relativepos = (100 * $interpos) / $microsatlength;
|
|
|
2492 # print "relativepos = $relativepos ,interpos=$interpos, interruption=$interruption, interlen=$interlen \n" if $printer == 1;
|
|
|
2493 $relativepos = (100 * ($interpos-$interlen)) / $microsatlength if $relativepos > 50;
|
|
|
2494 # print "--> = $relativepos\n" if $printer == 1;
|
|
|
2495 $interruption = "IND" if length($interruption) < 1;
|
|
|
2496
|
|
|
2497 if ($turn == 1){
|
|
|
2498 $fields[$microsatcord] = switch_micro($fields[$microsatcord]);
|
|
|
2499 $interruption = switch_nucl($interruption) unless $interruption eq "IND";
|
|
|
2500 $interpos = ($microsatlength - $interpos) - $interlen + 2;
|
|
|
2501 # print "turn interpos = $interpos for $fields[$microsatcord]\n" if $printer == 1;
|
|
|
2502 $relativepos = (100 * $interpos) / $microsatlength;
|
|
|
2503 $relativepos = (100 * ($interpos-$interlen)) / $microsatlength if $relativepos > 50;
|
|
|
2504
|
|
|
2505
|
|
|
2506 $strand = '+' if $strand eq '-';
|
|
|
2507 $strand = '-' if $strand eq '+';
|
|
|
2508 }
|
|
|
2509 # print "final relativepos = $relativepos\n" if $printer == 1;
|
|
|
2510 push(@relativeposes, $relativepos);
|
|
|
2511 }
|
|
|
2512 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))));
|
|
|
2513 }
|
|
|
2514
|
|
|
2515 else{
|
|
|
2516 push(@endstatement, join("\t",$exacttags[$species_counter],$sp_chr, $sp_start, $sp_end, $firstmotif,length($firstmotif),$fields[$microsatcord],$strand,$microsatlength,"NA","NA","NA", "NA"));
|
|
|
2517 }
|
|
|
2518
|
|
|
2519 $species_counter++;
|
|
|
2520 }
|
|
|
2521
|
|
|
2522 $locusmotif = $sameHash{$locusmotif} if exists $sameHash{$locusmotif};
|
|
|
2523 $locusmotif = $revHash{$locusmotif} if exists $revHash{$locusmotif};
|
|
|
2524
|
|
|
2525 my $endst = join("\t", @endstatement, $orthHchr, $Hstart, $Hend);
|
|
|
2526 # print join("\t", @endstatement, $orthHchr, $Hstart, $Hend), "\n" if $printer == 1;
|
|
|
2527
|
|
|
2528
|
|
|
2529 return (join("\t", @endstatement, $orthHchr, $Hstart, $Hend), $orthHchr, $Hstart, $Hend, $locusmotif, length($locusmotif));
|
|
|
2530
|
|
|
2531 }
|
|
|
2532
|
|
|
2533 sub switch_nucl{
|
|
|
2534 my @strand = split(/\s*/,$_[0]);
|
|
|
2535 for my $i (0 ... $#strand){
|
|
|
2536 if ($strand[$i] =~ /c/i) {$strand[$i] = "G";next;}
|
|
|
2537 if ($strand[$i] =~ /a/i) {$strand[$i] = "T";next;}
|
|
|
2538 if ($strand[$i] =~ /t/i) { $strand[$i] = "A";next;}
|
|
|
2539 if ($strand[$i] =~ /g/i) {$strand[$i] = "C";next;}
|
|
|
2540 }
|
|
|
2541 return join("",@strand);
|
|
|
2542 }
|
|
|
2543
|
|
|
2544
|
|
|
2545 sub switch_micro{
|
|
|
2546 my $micro = reverse($_[0]);
|
|
|
2547 my @strand = split(/\s*/,$micro);
|
|
|
2548 for my $i (0 ... $#strand){
|
|
|
2549 if ($strand[$i] =~ /c/i) {$strand[$i] = "G";next;}
|
|
|
2550 if ($strand[$i] =~ /a/i) {$strand[$i] = "T";next;}
|
|
|
2551 if ($strand[$i] =~ /t/i) { $strand[$i] = "A";next;}
|
|
|
2552 if ($strand[$i] =~ /g/i) {$strand[$i] = "C";next;}
|
|
|
2553 if ($strand[$i] =~ /\[/i) {$strand[$i] = "]";next;}
|
|
|
2554 if ($strand[$i] =~ /\]/i) {$strand[$i] = "[";next;}
|
|
|
2555 }
|
|
|
2556 return join("",@strand);
|
|
|
2557 }
|
|
|
2558 sub decipher_history{
|
|
|
2559 my $printer = 0;
|
|
|
2560 my ($mutations_array, $tags_string, $nodes, $branches_hash, $tree_analysis, $confirmation_string, $alivehash) = @_;
|
|
|
2561 my %mutations_hash=();
|
|
|
2562 foreach my $mutation (@$mutations_array){
|
|
|
2563 # print "mutation = $mutation\n" if $printer == 1;
|
|
|
2564 my %local = $mutation =~ /([\S ]+)=([\S ]+)/g;
|
|
|
2565 push @{$mutations_hash{$local{"node"}}},$mutation;
|
|
|
2566 # print "just for confirmation: $local{node} pushed as: $mutation\n" if $printer == 1;
|
|
|
2567 }
|
|
|
2568 my @nodes;
|
|
|
2569 my @birth_steps=();
|
|
|
2570 my @death_steps=();
|
|
|
2571
|
|
|
2572 my @tags=split(/\s*/,$tags_string);
|
|
|
2573 my @confirmation=split(/\s+/,$confirmation_string);
|
|
|
2574 my %info=();
|
|
|
2575
|
|
|
2576 for my $i (0 ... $#tags){
|
|
|
2577 $info{$tags[$i]}=$confirmation[$i];
|
|
|
2578 # print "feeding info: $tags[$i] = $info{$tags[$i]}\n" if $printer == 1;
|
|
|
2579 }
|
|
|
2580
|
|
|
2581 for my $keys (@$nodes) {
|
|
|
2582 foreach my $key (@$keys){
|
|
|
2583 # print "current key = $key\n";
|
|
|
2584 my $copykey = $key;
|
|
|
2585 $copykey =~ s/[\W ]+//g;
|
|
|
2586 my @copykeys=split(/\s*/,$copykey);
|
|
|
2587 my $states=();
|
|
|
2588 foreach my $copy (@copykeys){
|
|
|
2589 $states=$states.$info{$copy};
|
|
|
2590 }
|
|
|
2591 # print "reduced key = $copykey and state = $states\n" if $printer == 1;
|
|
|
2592
|
|
|
2593 if (exists $mutations_hash{$key}) {
|
|
|
2594
|
|
|
2595 if ($states=~/\+/){
|
|
|
2596 push @birth_steps, @{$mutations_hash{$key}};
|
|
|
2597 $birth_steps[$#birth_steps] =~ s/\S+=//g;
|
|
|
2598 delete $mutations_hash{$key};
|
|
|
2599 }
|
|
|
2600 else{
|
|
|
2601 push @death_steps, @{$mutations_hash{$key}};
|
|
|
2602 $death_steps[$#death_steps] =~ s/\S+=//g;
|
|
|
2603 delete $mutations_hash{$key};
|
|
|
2604 }
|
|
|
2605 }
|
|
|
2606 }
|
|
|
2607 }
|
|
|
2608 # print "conformation = $confirmation_string\n" if $printer == 1;
|
|
|
2609 push (@birth_steps, "NULL") if scalar(@birth_steps) == 0;
|
|
|
2610 push (@death_steps, "NULL") if scalar(@death_steps) == 0;
|
|
|
2611 # print "birth steps = ",join("\n",@birth_steps)," and death steps = ",join("\n",@death_steps),"\n" if $printer == 1;
|
|
|
2612 return \@birth_steps, \@death_steps;
|
|
|
2613 }
|
|
|
2614
|
|
|
2615 sub fillAlignmentGaps{
|
|
|
2616 my $printer = 0;
|
|
|
2617 # print "received: @_\n" if $printer == 1;
|
|
|
2618 my ($tree, $sequences, $alignment, $tagarray, $microsathash, $nonmicrosathash, $motif, $tree_analysis, $threshold, $microsatstarts) = @_;
|
|
|
2619 # print "in fillAlignmentGaps.. tree = $tree \n" if $printer == 1;
|
|
|
2620 my %sequence_hash=();
|
|
|
2621
|
|
|
2622 my @phases = ();
|
|
|
2623 my $concat = $motif.$motif;
|
|
|
2624 my $motifsize = length($motif);
|
|
|
2625
|
|
|
2626 for my $i (1 ... $motifsize){
|
|
|
2627 push @phases, substr($concat, $i, $motifsize);
|
|
|
2628 }
|
|
|
2629
|
|
|
2630 my $concatalignment = ();
|
|
|
2631 foreach my $tag (@tags){
|
|
|
2632 $concatalignment = $concatalignment.$alignment->{$tag};
|
|
|
2633 }
|
|
|
2634 # print "returningg NULL","NULL","NULL", "NULL\n" if $concatalignment !~ /-/;
|
|
|
2635 return 0, "NULL","NULL","NULL", "NULL","NULL" if $concatalignment !~ /-/;
|
|
|
2636
|
|
|
2637
|
|
|
2638
|
|
|
2639 my %node_sequences_temp=();
|
|
|
2640 my %node_alignments_temp =(); #NEW, Nov 28 2008
|
|
|
2641
|
|
|
2642 my @tags=();
|
|
|
2643 my @locus_sequences=();
|
|
|
2644 my %alivehash=();
|
|
|
2645
|
|
|
2646 # print "IN fillAlignmentGaps\n";# <STDIN>;
|
|
|
2647 my %fillrecord = ();
|
|
|
2648
|
|
|
2649 my $change = 0;
|
|
|
2650 foreach my $tag (@$tagarray) {
|
|
|
2651 #print "adding: $tag\n";
|
|
|
2652 push(@tags, $tag);
|
|
|
2653 if (exists $microsathash->{$tag}){
|
|
|
2654 my $micro = $microsathash->{$tag};
|
|
|
2655 my $orig_micro = $micro;
|
|
|
2656 ($micro, $fillrecord{$tag}) = fillgaps($micro, \@phases);
|
|
|
2657 $change = 1 if uc($micro) ne uc($orig_micro);
|
|
|
2658 $node_sequences_temp{$tag}=$micro if $microsathash->{$tag} ne "NULL";
|
|
|
2659 }
|
|
|
2660 if (exists $nonmicrosathash->{$tag}){
|
|
|
2661 my $micro = $nonmicrosathash->{$tag};
|
|
|
2662 my $orig_micro = $micro;
|
|
|
2663 ($micro, $fillrecord{$tag}) = fillgaps($micro, \@phases);
|
|
|
2664 $change = 1 if uc($micro) ne uc($orig_micro);
|
|
|
2665 $node_sequences_temp{$tag}=$micro if $nonmicrosathash->{$tag} ne "NULL";
|
|
|
2666 }
|
|
|
2667
|
|
|
2668 if (exists $alignment->{$tag}){
|
|
|
2669 my $micro = $alignment->{$tag};
|
|
|
2670 my $orig_micro = $micro;
|
|
|
2671 ($micro, $fillrecord{$tag}) = fillgaps($micro, \@phases);
|
|
|
2672 $change = 1 if uc($micro) ne uc($orig_micro);
|
|
|
2673 $node_alignments_temp{$tag}=$micro if $alignment->{$tag} ne "NULL";
|
|
|
2674 }
|
|
|
2675
|
|
|
2676 #print "adding to node_sequences: $tag = ",$node_sequences_temp{$tag},"\n" if $printer == 1;
|
|
|
2677 #print "adding to node_alignments: $tag = ",$node_alignments_temp{$tag},"\n" if $printer == 1;
|
|
|
2678 }
|
|
|
2679
|
|
|
2680
|
|
|
2681 my %node_sequences=();
|
|
|
2682 my %node_alignments =(); #NEW, Nov 28 2008
|
|
|
2683 foreach my $tag (@$tagarray) {
|
|
|
2684 $node_sequences{$tag} = join ".",split(/\s*/,$node_sequences_temp{$tag});
|
|
|
2685 $node_alignments{$tag} = join ".",split(/\s*/,$node_alignments_temp{$tag});
|
|
|
2686 }
|
|
|
2687 # print "\n", "#" x 50, "\n" if $printer == 1;
|
|
|
2688 foreach my $tag (@tags){
|
|
|
2689 # print "$tag: $alignment->{$tag} = $node_alignments{$tag}\n" if $printer == 1;
|
|
|
2690 }
|
|
|
2691 # print "\n", "#" x 50, "\n" if $printer == 1;
|
|
|
2692 # print "change = $change\n";
|
|
|
2693 #<STDIN> if $concatalignment=~/\-/;
|
|
|
2694
|
|
|
2695 # <STDIN> if $printer == 1 && $concatalignment =~ /\-/;
|
|
|
2696
|
|
|
2697 return 0, "NULL","NULL","NULL", "NULL", "NULL" if $change == 0;
|
|
|
2698
|
|
|
2699 my ($nodes_arr, $branches_hash) = get_nodes($tree);
|
|
|
2700 my @nodes=@$nodes_arr;
|
|
|
2701 # print "recieved nodes = @nodes\n" if $printer == 1;
|
|
|
2702
|
|
|
2703
|
|
|
2704 #POPULATE branches_hash WITH INFORMATION ABOUT LIVESTATUS
|
|
|
2705 foreach my $keys (@nodes){
|
|
|
2706 my @pair = @$keys;
|
|
|
2707 my $joint = "(".join(", ",@pair).")";
|
|
|
2708 my $copykey = join "", @pair;
|
|
|
2709 $copykey =~ s/[\W ]+//g;
|
|
|
2710 # print "for node: $keys, copykey = $copykey and joint = $joint\n" if $printer == 1;
|
|
|
2711 my $livestatus = 1;
|
|
|
2712 foreach my $copy (split(/\s*/,$copykey)){
|
|
|
2713 $livestatus = 0 if !exists $alivehash{$copy};
|
|
|
2714 }
|
|
|
2715 $alivehash{$joint} = $joint if !exists $alivehash{$joint} && $livestatus == 1;
|
|
|
2716 # print "alivehash = $alivehash{$joint}\n" if exists $alivehash{$joint} && $printer == 1;
|
|
|
2717 }
|
|
|
2718
|
|
|
2719
|
|
|
2720
|
|
|
2721 @nodes = reverse(@nodes); #1 THIS IS IN ORDER TO GO THROUGH THE TREE FROM LEAVES TO ROOT.
|
|
|
2722
|
|
|
2723 my @mutations_array=();
|
|
|
2724
|
|
|
2725 my $joint = ();
|
|
|
2726 foreach my $node (@nodes){
|
|
|
2727 my @pair = @$node;
|
|
|
2728 # print "now in the nodes for loop, pair = @pair\n and sequences=\n" if $printer == 1;
|
|
|
2729 $joint = "(".join(", ",@pair).")";
|
|
|
2730 # print "joint = $joint \n" if $printer == 1;
|
|
|
2731 my @pair_sequences=();
|
|
|
2732
|
|
|
2733 foreach my $tag (@pair){
|
|
|
2734 # print "tag = $tag: " if $printer == 1;
|
|
|
2735 # print $node_alignments{$tag},"\n" if $printer == 1;
|
|
|
2736 push @pair_sequences, $node_alignments{$tag};
|
|
|
2737 }
|
|
|
2738 # print "fillgap\n";
|
|
|
2739 my ($compared, $substitutions_list) = base_by_base_simple($motif,\@pair_sequences, scalar(@pair_sequences), @pair, $joint);
|
|
|
2740 $node_alignments{$joint}=$compared;
|
|
|
2741 push( @mutations_array,split(/:/,$substitutions_list));
|
|
|
2742 # print "newly added to node_sequences: $node_alignments{$joint} and list of mutations = @mutations_array\n" if $printer == 1;
|
|
|
2743 }
|
|
|
2744 # 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;
|
|
|
2745 # <STDIN> if $printer == 1;
|
|
|
2746
|
|
|
2747 my $analayzed_mutations = analyze_mutations(\@mutations_array, \@nodes, $branches_hash, $alignment, \@tags, \%alivehash, \%node_sequences, $microsatstarts, $motif);
|
|
|
2748
|
|
|
2749 # print "returningt: ", $analayzed_mutations, \@nodes,"\n" if scalar @mutations_array > 0;;
|
|
|
2750 # print "returningy: NULL, NULL, NULL " if scalar @mutations_array == 0 && $printer == 1;
|
|
|
2751 # print "final node alignment after filling for $joint= " if $printer == 1;
|
|
|
2752 # print "$node_alignments{$joint}\n" if $printer == 1;
|
|
|
2753
|
|
|
2754
|
|
|
2755 return 1, $analayzed_mutations, \@nodes, $branches_hash, \%alivehash, $node_alignments{$joint} if scalar @mutations_array > 0 ;
|
|
|
2756 return 1, "NULL","NULL","NULL", "NULL", "NULL" if scalar @mutations_array == 0;
|
|
|
2757 }
|
|
|
2758
|
|
|
2759
|
|
|
2760
|
|
|
2761 sub add_mutation{
|
|
|
2762 my $printer = 0;
|
|
|
2763 # print "IN SUBROUTUNE add_mutation.. information received = @_\n" if $printer == 1;
|
|
|
2764 my ($i , $bite, $to, $from) = @_;
|
|
|
2765 # print "bite = $bite.. all received info = ",join("^", @_),"\n" if $printer == 1;
|
|
|
2766 # print "to=$to\n" if $printer == 1;
|
|
|
2767 # print "tis split = ",join(" and ",split(/!/,$to)),"\n" if $printer == 1;
|
|
|
2768 my @toields = split "!",$to;
|
|
|
2769 # print "toilds = @toields\n" if $printer == 1;
|
|
|
2770 my @mutations=();
|
|
|
2771
|
|
|
2772 foreach my $toield (@toields){
|
|
|
2773 my @toinfo=split(":",$toield);
|
|
|
2774 # print " at toinfo=@toinfo \n" if $printer == 1;
|
|
|
2775 next if $toinfo[1] =~ /$from/i;
|
|
|
2776 my @mutation = @toinfo if $toinfo[1] !~ /$from/i;
|
|
|
2777 # print "adding to mutaton list: ", join(",", "node=$mutation[0]","type=substitution" ,"position=$i", "from=$from", "to=$mutation[1]", "insertion=", "deletion="),"\n" if $printer == 1;
|
|
|
2778 push (@mutations, join("\t", "node=$mutation[0]","type=substitution" ,"position=$i", "from=$from", "to=$mutation[1]", "insertion=", "deletion="));
|
|
|
2779 }
|
|
|
2780 return @mutations;
|
|
|
2781 }
|
|
|
2782
|
|
|
2783
|
|
|
2784 sub add_bases{
|
|
|
2785
|
|
|
2786 my $printer = 0;
|
|
|
2787 # print "IN SUBROUTUNE add_bases.. information received = @_\n" if $printer == 1;
|
|
|
2788 my ($optional0, $optional1, $pair0, $pair1,$joint) = @_;
|
|
|
2789 my $total_list=();
|
|
|
2790
|
|
|
2791 my @total_list0=split(/!/,$optional0);
|
|
|
2792 my @total_list1=split(/!/,$optional1);
|
|
|
2793 my @all_list=();
|
|
|
2794 my %total_hash0=();
|
|
|
2795 foreach my $entry (@total_list0) {
|
|
|
2796 $entry = uc $entry;
|
|
|
2797 $entry =~ /(\S+):(\S+)/;
|
|
|
2798 $total_hash0{$2}=$1;
|
|
|
2799 push @all_list, $2;
|
|
|
2800 }
|
|
|
2801
|
|
|
2802 my %total_hash1=();
|
|
|
2803 foreach my $entry (@total_list1) {
|
|
|
2804 $entry = uc $entry;
|
|
|
2805 $entry =~ /(\S+):(\S+)/;
|
|
|
2806 $total_hash1{$2}=$1;
|
|
|
2807 push @all_list, $2;
|
|
|
2808 }
|
|
|
2809
|
|
|
2810 my %alphabetical_hash=();
|
|
|
2811 my @return_options=();
|
|
|
2812
|
|
|
2813 for my $i (0 ... $#all_list){
|
|
|
2814 my $alph = $all_list[$i];
|
|
|
2815 if (exists $total_hash0{$alph} && exists $total_hash1{$alph}){
|
|
|
2816 push(@return_options, $joint.":".$alph);
|
|
|
2817 delete $total_hash0{$alph}; delete $total_hash1{$alph};
|
|
|
2818 }
|
|
|
2819 if (exists $total_hash0{$alph} && !exists $total_hash1{$alph}){
|
|
|
2820 push(@return_options, $pair0.":".$alph);
|
|
|
2821 delete $total_hash0{$alph};
|
|
|
2822 }
|
|
|
2823 if (!exists $total_hash0{$alph} && exists $total_hash1{$alph}){
|
|
|
2824 push(@return_options, $pair1.":".$alph);
|
|
|
2825 delete $total_hash1{$alph};
|
|
|
2826 }
|
|
|
2827
|
|
|
2828 }
|
|
|
2829 # print "returning ",join "!",@return_options,"\n" if $printer == 1;
|
|
|
2830 return join "!",@return_options;
|
|
|
2831
|
|
|
2832 }
|
|
|
2833
|
|
|
2834
|
|
|
2835 sub fillgaps{
|
|
|
2836 # print "IN fillgaps: @_\n";
|
|
|
2837 my ($micro, $phasesinput) = @_;
|
|
|
2838 #print "in microsathash ,,.. micro = $micro\n";
|
|
|
2839 return $micro if $micro !~ /\-/;
|
|
|
2840 my $orig_micro = $micro;
|
|
|
2841 my @phases = @$phasesinput;
|
|
|
2842
|
|
|
2843 my %tested_patterns = ();
|
|
|
2844
|
|
|
2845 foreach my $phase (@phases){
|
|
|
2846 # print "considering phase: $phase\n";
|
|
|
2847 my @phase_prefixes = ();
|
|
|
2848 my @prephase_left_contexts = ();
|
|
|
2849 my @prephase_right_contexts = ();
|
|
|
2850 my @pregapsize = ();
|
|
|
2851 my @prepostfilins = ();
|
|
|
2852
|
|
|
2853 my @phase_suffixes;
|
|
|
2854 my @suffphase_left_contexts;
|
|
|
2855 my @suffphase_right_contexts;
|
|
|
2856 my @suffgapsize;
|
|
|
2857 my @suffpostfilins;
|
|
|
2858
|
|
|
2859 my @postfilins = ();
|
|
|
2860 my $motifsize = length($phases[0]);
|
|
|
2861
|
|
|
2862 my $change = 0;
|
|
|
2863
|
|
|
2864 for my $u (0 ... $motifsize-1){
|
|
|
2865 my $concat = $phase.$phase.$phase.$phase;
|
|
|
2866 my @concatarr = split(/\s*/, $concat);
|
|
|
2867 my $l = 0;
|
|
|
2868 while ($l < $u){
|
|
|
2869 shift @concatarr;
|
|
|
2870 $l++;
|
|
|
2871 }
|
|
|
2872 $concat = join ("", @concatarr);
|
|
|
2873
|
|
|
2874 for my $t (0 ... $motifsize-1){
|
|
|
2875 for my $k (1 ... $motifsize-1){
|
|
|
2876 push @phase_prefixes, substr($concat, $motifsize+$t, $k);
|
|
|
2877 push @prephase_left_contexts, substr ($concat, $t, $motifsize);
|
|
|
2878 push @prephase_right_contexts, substr ($concat, $motifsize+$t+$k+($motifsize-$k), 1);
|
|
|
2879 push @pregapsize, $k;
|
|
|
2880 push @prepostfilins, substr($concat, $motifsize+$t+$k, ($motifsize-$k));
|
|
|
2881 # 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";
|
|
|
2882 # print "phase_prefixes = $phase_prefixes[$#phase_prefixes]\n";
|
|
|
2883 # print "prephase_left_contexts = $prephase_left_contexts[$#prephase_left_contexts]\n";
|
|
|
2884 # print "prephase_right_contexts = $prephase_right_contexts[$#prephase_right_contexts]\n";
|
|
|
2885 # print "pregapsize = $pregapsize[$#pregapsize]\n";
|
|
|
2886 # print "prepostfilins = $prepostfilins[$#prepostfilins]\n";
|
|
|
2887 }
|
|
|
2888 }
|
|
|
2889 }
|
|
|
2890
|
|
|
2891 # print "looking if $micro =~ /($phase\-{$motifsize})/i || $micro =~ /^(\-{$motifsize,}$phase)/i\n";
|
|
|
2892 if ($micro =~ /($phase\-{$motifsize,})$/i || $micro =~ /^(\-{$motifsize,}$phase)/i){
|
|
|
2893 # print "micro: $micro needs further gap removal: $1\n";
|
|
|
2894 while ($micro =~ /$phase(\-{$motifsize,})$/i || $micro =~ /^(\-{$motifsize,})$phase/i){
|
|
|
2895 # print "micro: $micro needs further gap removal: $1\n";
|
|
|
2896
|
|
|
2897 # print "phase being considered = $phase\n";
|
|
|
2898 my $num = ();
|
|
|
2899 $num = $micro =~ s/$phase\-{$motifsize}/$phase$phase/gi if $micro =~ /$phase\-{$motifsize,}/i;
|
|
|
2900 $num = $micro =~ s/\-{$motifsize}$phase/$phase$phase/gi if $micro =~ /\-{$motifsize,}$phase/i;
|
|
|
2901 # print "num = $num\n";
|
|
|
2902 $change = 1 if $num == 1;
|
|
|
2903 }
|
|
|
2904 }
|
|
|
2905
|
|
|
2906 elsif ($micro =~ /(($phase)+)\-{$motifsize,}(($phase)+)/i){
|
|
|
2907 while ($micro =~ /(($phase)+)\-{$motifsize,}(($phase)+)/i){
|
|
|
2908 # print "checking lengths of $1 and $3 for $micro... \n";
|
|
|
2909 my $num = ();
|
|
|
2910 if (length($1) >= length($3)){
|
|
|
2911 # print "$micro matches (($phase)+)\-{$motifsize,}(($phase)+) = $1, >= , $3 \n";
|
|
|
2912 $num = $micro =~ s/$phase\-{$motifsize}/$phase$phase/gi ;
|
|
|
2913 }
|
|
|
2914 if (length($1) < length($3)){
|
|
|
2915 # print "$micro matches (($phase)+)\-{$motifsize,}(($phase)+) = $1, < , $3 \n";
|
|
|
2916 $num = $micro =~ s/\-{$motifsize}$phase/$phase$phase/gi ;
|
|
|
2917 }
|
|
|
2918 # print "micro changed to $micro\n";
|
|
|
2919 }
|
|
|
2920 }
|
|
|
2921 elsif ($micro =~ /([A-Z]+)\-{$motifsize,}(($phase)+)/i){
|
|
|
2922 while ($micro =~ /([A-Z]+)\-{$motifsize,}(($phase)+)/i){
|
|
|
2923 # print "$micro matches ([A-Z]+)\-{$motifsize}(($phase)+) = 1=$1, - , 3=$3 \n";
|
|
|
2924 my $num = 0;
|
|
|
2925 $num = $micro =~ s/\-{$motifsize}$phase/$phase$phase/gi ;
|
|
|
2926 }
|
|
|
2927 }
|
|
|
2928 elsif ($micro =~ /(($phase)+)\-{$motifsize,}([A-Z]+)/i){
|
|
|
2929 while ($micro =~ /(($phase)+)\-{$motifsize,}([A-Z]+)/i){
|
|
|
2930 # print "$micro matches (($phase)+)\-{$motifsize,}([A-Z]+) = 1=$1, - , 3=$3 \n";
|
|
|
2931 my $num = 0;
|
|
|
2932 $num = $micro =~ s/$phase\-{$motifsize}/$phase$phase/gi ;
|
|
|
2933 }
|
|
|
2934 }
|
|
|
2935
|
|
|
2936 # print "$orig_micro to $micro\n";
|
|
|
2937
|
|
|
2938 #s <STDIN>;
|
|
|
2939
|
|
|
2940 for my $h (0 ... $#phase_prefixes){
|
|
|
2941 # print "searching using prefix : $prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h]\n";
|
|
|
2942 my $pattern = $prephase_left_contexts[$h].$phase_prefixes[$h].$pregapsize[$h].$prephase_right_contexts[$h];
|
|
|
2943 # print "returning orig_micro = $orig_micro, micro = $micro \n" if exists $tested_patterns{$pattern};
|
|
|
2944 if ($micro =~ /$prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h]/i){
|
|
|
2945 return $orig_micro if exists $tested_patterns{$pattern};
|
|
|
2946 while ($micro =~ /($prephase_left_contexts[$h]$phase_prefixes[$h]\-{$pregapsize[$h]}$prephase_right_contexts[$h])/i){
|
|
|
2947 $tested_patterns{$pattern} = $pattern;
|
|
|
2948 # print "micro: $micro needs further gap removal: $1\n";
|
|
|
2949
|
|
|
2950 # print "prefix being considered = $phase_prefixes[$h]\n";
|
|
|
2951 my $num = ();
|
|
|
2952 $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) ;
|
|
|
2953 # print "num = $num, micro = $micro\n";
|
|
|
2954 $change = 1 if $num == 1;
|
|
|
2955
|
|
|
2956 return $orig_micro if $num > 1;
|
|
|
2957 }
|
|
|
2958 }
|
|
|
2959
|
|
|
2960 }
|
|
|
2961 }
|
|
|
2962 return $orig_micro if length($micro) != length($orig_micro);
|
|
|
2963 return $micro;
|
|
|
2964 }
|
|
|
2965
|
|
|
2966 sub selectMutationArray{
|
|
|
2967 my $printer =0;
|
|
|
2968
|
|
|
2969 my $oldmutspt = $_[0];
|
|
|
2970 my $newmutspt = $_[1];
|
|
|
2971 my $tagstringpt = $_[2];
|
|
|
2972 my $alivehashpt = $_[3];
|
|
|
2973 my $alignmentpt = $_[4];
|
|
|
2974 my $motif = $_[5];
|
|
|
2975
|
|
|
2976 my @alivehasharr=();
|
|
|
2977
|
|
|
2978 my @tags = @$tagstringpt;
|
|
|
2979 my $alignmentln = length($alignmentpt->{$tags[0]});
|
|
|
2980
|
|
|
2981 foreach my $key (keys %$alivehashpt) { push @alivehasharr, $key; }
|
|
|
2982
|
|
|
2983 my %newside = ();
|
|
|
2984 my %oldside = ();
|
|
|
2985 my %newmuts = ();
|
|
|
2986
|
|
|
2987 my %commons = ();
|
|
|
2988 my %olds = ();
|
|
|
2989 foreach my $old (@$oldmutspt){
|
|
|
2990 $olds{$old} = 1;
|
|
|
2991 }
|
|
|
2992 foreach my $new (@$newmutspt){
|
|
|
2993 $commons{$new} = 1 if exists $olds{$new};;
|
|
|
2994 }
|
|
|
2995
|
|
|
2996
|
|
|
2997 foreach my $pos ( 0 ... $alignmentln){
|
|
|
2998 #print "pos = $pos\n" if $printer == 1;
|
|
|
2999 my $newyes = 0;
|
|
|
3000 foreach my $mut (@$newmutspt){
|
|
|
3001 $newmuts{$mut} = 1;
|
|
|
3002 chomp $mut;
|
|
|
3003 $newyes++;
|
|
|
3004 $mut =~ s/=\t/= \t/g;
|
|
|
3005 $mut =~ s/=$/= /g;
|
|
|
3006
|
|
|
3007 $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\- ]+)/;
|
|
|
3008 my $node = $1;
|
|
|
3009 next if $3 != $pos;
|
|
|
3010 # print "new mut = $mut\n" if $printer == 1;
|
|
|
3011 # print "node = $node, pos = $3 ... and alivehasharr = >@alivehasharr<\n" if $printer == 1;
|
|
|
3012 my $alivenode = 0;
|
|
|
3013 foreach my $key (@alivehasharr){
|
|
|
3014 $alivenode = 1 if $key =~ /$node/;
|
|
|
3015 }
|
|
|
3016 # next if $alivenode == 0;
|
|
|
3017 my $indel_type = " ";
|
|
|
3018 if ($2 eq "insertion" || $2 eq "deletion"){
|
|
|
3019 my $thisindel = ();
|
|
|
3020 $thisindel = $6 if $2 eq "insertion";
|
|
|
3021 $thisindel = $7 if $2 eq "deletion";
|
|
|
3022
|
|
|
3023 $indel_type = "i".checkIndelType($node, $thisindel, $motif,$alignmentpt,$3, $2) if $2 eq "insertion";
|
|
|
3024 $indel_type = "d".checkIndelType($node, $thisindel, $motif,$alignmentpt, $3, $2) if $2 eq "deletion";
|
|
|
3025 $indel_type = $indel_type."f" if $indel_type =~ /mot/ && length($thisindel) >= length($motif);
|
|
|
3026 }
|
|
|
3027 # print "indeltype = $indel_type\n" if $printer == 1;
|
|
|
3028 my $added = 0;
|
|
|
3029
|
|
|
3030 if (exists $newside{$pos} && $indel_type =~ /[a-z]+/){
|
|
|
3031 # print "we have a preexisting one for $pos\n" if $printer == 1;
|
|
|
3032 my @preexisting = @{$newside{$pos}};
|
|
|
3033 foreach my $pre (@preexisting){
|
|
|
3034 # print "looking at $pre\n" if $printer == 1;
|
|
|
3035 next if $pre !~ /node=$node/;
|
|
|
3036 next if $pre !~ /indeltype=([a-z]+)/;
|
|
|
3037 my $currtype = $1;
|
|
|
3038
|
|
|
3039 if ($currtype =~ /inon/ && $indel_type =~ /dmot/){
|
|
|
3040 delete $newside{$pos};
|
|
|
3041 push @{$newside{$pos}}, $pre;
|
|
|
3042 $added = 1;
|
|
|
3043 }
|
|
|
3044 if ($currtype =~ /dnon/ && $indel_type =~ /imot/){
|
|
|
3045 delete $newside{$pos};
|
|
|
3046 push @{$newside{$pos}}, $pre;
|
|
|
3047 $added = 1;
|
|
|
3048 }
|
|
|
3049 if ($currtype =~ /dmot/ && $indel_type =~ /inon/){
|
|
|
3050 delete $newside{$pos};
|
|
|
3051 push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3052 $added = 1;
|
|
|
3053 }
|
|
|
3054 if ($currtype =~ /imot/ && $indel_type =~ /dnon/){
|
|
|
3055 delete $newside{$pos};
|
|
|
3056 push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3057 $added = 1;
|
|
|
3058 }
|
|
|
3059 }
|
|
|
3060 }
|
|
|
3061 # print "added = $added\n" if $printer == 1;
|
|
|
3062 push @{$newside{$pos}}, $mut."\tindeltype=$indel_type" if $added == 0;
|
|
|
3063 # print "for new pos,: $pos we have: @{$newside{$pos}}\n " if $printer == 1;
|
|
|
3064 }
|
|
|
3065 }
|
|
|
3066
|
|
|
3067 foreach my $pos ( 0 ... $alignmentln){
|
|
|
3068 my $oldyes = 0;
|
|
|
3069 foreach my $mut (@$oldmutspt){
|
|
|
3070 chomp $mut;
|
|
|
3071 $oldyes++;
|
|
|
3072 $mut =~ s/=\t/= \t/g;
|
|
|
3073 $mut =~ s/=$/= /g;
|
|
|
3074 $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\- ]+)/;
|
|
|
3075 my $node = $1;
|
|
|
3076 next if $3 != $pos;
|
|
|
3077 # print "old mut = $mut\n" if $printer == 1;
|
|
|
3078 my $alivenode = 0;
|
|
|
3079 foreach my $key (@alivehasharr){
|
|
|
3080 $alivenode = 1 if $key =~ /$node/;
|
|
|
3081 }
|
|
|
3082 #next if $alivenode == 0;
|
|
|
3083 my $indel_type = " ";
|
|
|
3084 if ($2 eq "insertion" || $2 eq "deletion"){
|
|
|
3085 $indel_type = "i".checkIndelType($node, $6, $motif,$alignmentpt, $3, $2) if $2 eq "insertion";
|
|
|
3086 $indel_type = "d".checkIndelType($node, $7, $motif,$alignmentpt, $3, $2) if $2 eq "deletion";
|
|
|
3087 next if $indel_type =~/non/;
|
|
|
3088 }
|
|
|
3089 else{ next;}
|
|
|
3090
|
|
|
3091 my $imp=0;
|
|
|
3092 $imp = 1 if $indel_type =~ /dmot/ && $alivenode == 0;
|
|
|
3093 $imp = 1 if $indel_type =~ /imot/ && $alivenode == 1;
|
|
|
3094
|
|
|
3095
|
|
|
3096 if (exists $newside{$pos} && $indel_type =~ /[a-z]+/){
|
|
|
3097 my @preexisting = @{$newside{$pos}};
|
|
|
3098 # print "we have a preexisting one for $pos: @preexisting\n" if $printer == 1;
|
|
|
3099 next if $imp == 0;
|
|
|
3100
|
|
|
3101 if (scalar(@preexisting) == 1){
|
|
|
3102 my $foundmut = $preexisting[0];
|
|
|
3103 $foundmut=~ /node=([A-Z, \(\)]+)/;
|
|
|
3104 next if $1 eq $node;
|
|
|
3105
|
|
|
3106 if (exists $oldside{$pos} || exists $commons{$foundmut}){
|
|
|
3107 # print "not replacing, but just adding\n" if $printer == 1;
|
|
|
3108 push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3109 push @{$oldside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3110 next;
|
|
|
3111 }
|
|
|
3112
|
|
|
3113 delete $newside{$pos};
|
|
|
3114 push @{$oldside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3115 push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3116 # print "now new one is : @{$newside{$pos}}\n" if $printer == 1;
|
|
|
3117 }
|
|
|
3118 # print "for pos: $pos: @{$newside{$pos}}\n" if $printer == 1;
|
|
|
3119 next;
|
|
|
3120 }
|
|
|
3121
|
|
|
3122
|
|
|
3123 my @news = @{$newside{$pos}} if exists $newside{$pos};
|
|
|
3124 # print "mut = $mut and news = @news\n" if $printer == 1;
|
|
|
3125 push @{$oldside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3126 push @{$newside{$pos}}, $mut."\tindeltype=$indel_type";
|
|
|
3127 }
|
|
|
3128 }
|
|
|
3129 # print "in the end, our collected mutations = \n" if $printer == 1;
|
|
|
3130 my @returnarr = ();
|
|
|
3131 foreach my $key (keys %newside) {push @returnarr,@{$newside{$key}};}
|
|
|
3132 # print join("\n", @returnarr),"\n" if $printer == 1;
|
|
|
3133 #<STDIN>;
|
|
|
3134 return @returnarr;
|
|
|
3135
|
|
|
3136 }
|
|
|
3137
|
|
|
3138
|
|
|
3139 sub checkIndelType{
|
|
|
3140 my $printer = 0;
|
|
|
3141 my $node = $_[0];
|
|
|
3142 my $indel = $_[1];
|
|
|
3143 my $motif = $_[2];
|
|
|
3144 my $alignmentpt = $_[3];
|
|
|
3145 my $posit = $_[4];
|
|
|
3146 my $type = $_[5];
|
|
|
3147 my @phases =();
|
|
|
3148 my %prephases = ();
|
|
|
3149 my %postphases = ();
|
|
|
3150 #print "motif = $motif\n";
|
|
|
3151 # print "IN checkIndelType ... received: @_\n" if $printer == 1;
|
|
|
3152 my $concat = $motif.$motif.$motif.$motif;
|
|
|
3153 my $motiflength = length($motif);
|
|
|
3154
|
|
|
3155 if ($motiflength > length ($indel)){
|
|
|
3156 return "non" if $motif !~ /$indel/i;
|
|
|
3157 return checkIndelType_ComplexAnalysis($node, $indel, $motif, $alignmentpt, $posit, $type);
|
|
|
3158 }
|
|
|
3159
|
|
|
3160 my $firstpass = 0;
|
|
|
3161 for my $y (0 ... $motiflength-1){
|
|
|
3162 my $phase = substr($concat, $motiflength+$y, $motiflength);
|
|
|
3163 push @phases, $phase;
|
|
|
3164 $firstpass = 1 if $indel =~ /$phase/i;
|
|
|
3165 for my $k (0 ... length($motif)-1){
|
|
|
3166 # 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;
|
|
|
3167 my $pre = substr($concat, $motiflength+$y-$k, $k );
|
|
|
3168 my $post = substr($concat, $motiflength+$y+$motiflength, $k);
|
|
|
3169 # print "adding to phases : $phase - $pre and $post\n" if $printer == 1;
|
|
|
3170 push @{$prephases{$phase}} , $pre;
|
|
|
3171 push @{$postphases{$phase}} , $post;
|
|
|
3172 }
|
|
|
3173
|
|
|
3174 }
|
|
|
3175 # print "firstpass 1= $firstpass\n" if $printer == 1;
|
|
|
3176 return "non" if $firstpass ==0;
|
|
|
3177 $firstpass =0;
|
|
|
3178
|
|
|
3179 foreach my $phase (@phases){
|
|
|
3180 my @pres = @{$prephases{$phase}};
|
|
|
3181 my @posts = @{$postphases{$phase}};
|
|
|
3182
|
|
|
3183 foreach my $pre (@pres){
|
|
|
3184 foreach my $post (@posts){
|
|
|
3185
|
|
|
3186 $firstpass = 1 if $indel =~ /($pre)?($phase)+($post)?/i && length($indel) > (3 * length($motif));
|
|
|
3187 $firstpass = 1 if $indel =~ /^($pre)?($phase)+($post)?$/i && length($indel) < (3 * length($motif));
|
|
|
3188 # print "matched here : ($pre)?($phase)+($post)?\n" if $printer == 1;
|
|
|
3189 last if $firstpass == 1;
|
|
|
3190 }
|
|
|
3191 last if $firstpass == 1;
|
|
|
3192 }
|
|
|
3193 last if $firstpass == 1;
|
|
|
3194 }
|
|
|
3195 # print "firstpass 2= $firstpass\n" if $printer == 1;
|
|
|
3196 return "non" if $firstpass ==0;
|
|
|
3197 return "mot" if $firstpass ==1;
|
|
|
3198 }
|
|
|
3199
|
|
|
3200
|
|
|
3201 sub checkIndelType_ComplexAnalysis{
|
|
|
3202 my $printer = 0;
|
|
|
3203 my $node = $_[0];
|
|
|
3204 my $indel = $_[1];
|
|
|
3205 my $motif = $_[2];
|
|
|
3206 my $alignmentpt = $_[3];
|
|
|
3207 my $pos = $_[4];
|
|
|
3208 my $type = $_[5];
|
|
|
3209 my @speciesinvolved = $node =~ /[A-Z]+/g;
|
|
|
3210
|
|
|
3211 my @seqs = ();
|
|
|
3212 my $residualseq = length($motif) - length($indel);
|
|
|
3213 # print "IN COMPLEX ANALYSIS ... received: @_ .... speciesinvolved = @speciesinvolved\n" if $printer == 1;
|
|
|
3214 # print "we have position = $pos, sseq = $alignmentpt->{$speciesinvolved[0]}\n" if $printer == 1;
|
|
|
3215 # print "residualseq = $residualseq\n" if $printer == 1;
|
|
|
3216 # print "pos=$pos... got: @_\n" if $printer == 1;
|
|
|
3217 foreach my $sp (@speciesinvolved){
|
|
|
3218 my $spseq = $alignmentpt->{$sp};
|
|
|
3219 #print "orig spseq = $spseq\n";
|
|
|
3220 my $subseq = ();
|
|
|
3221
|
|
|
3222 if ($type eq "deletion"){
|
|
|
3223 my @indelparts = split(/\s*/,$indel);
|
|
|
3224 my @seqparts = split(/\s*/,$spseq);
|
|
|
3225
|
|
|
3226 for my $p ($pos ... $pos+length($indel)-1){
|
|
|
3227 $seqparts[$p] = shift @indelparts;
|
|
|
3228 }
|
|
|
3229 $spseq = join("",@seqparts);
|
|
|
3230 }
|
|
|
3231 #print "mod spseq = $spseq\n";
|
|
|
3232 # $spseq=~ s/\-//g if $type !~ /deletion/;
|
|
|
3233 # print "substr($spseq, $pos-($residualseq), length($indel)+$residualseq+$residualseq)\n" if $pos > 0 && $pos < (length($spseq) - length($motif)) && $printer == 1;
|
|
|
3234 # print "substr($spseq, 0, length($indel)+$residualseq)\n" if $pos == 0 && $printer == 1;
|
|
|
3235 # print "substr($spseq, $pos - $residualseq, length($indel)+$residualseq)\n" if $pos >= (length($spseq) - length($motif)) && $printer == 1;
|
|
|
3236
|
|
|
3237 $subseq = substr($spseq, $pos-($residualseq), length($indel)+$residualseq+$residualseq) if $pos > 0 && $pos < (length($spseq) - length($motif)) ;
|
|
|
3238 $subseq = substr($spseq, 0, length($indel)+$residualseq) if $pos == 0;
|
|
|
3239 $subseq = substr($spseq, $pos - $residualseq, length($indel)+$residualseq) if $pos >= (length($spseq) - length($motif)) ;
|
|
|
3240 # print "spseq = $spseq . subseq=$subseq . type = $type\n" if $printer == 1;
|
|
|
3241 #<STDIN> if $subseq !~ /[a-z\-]/i;
|
|
|
3242 $subseq =~ s/\-/$indel/g if $type =~ /insertion/;
|
|
|
3243 push @seqs, $subseq;
|
|
|
3244 # print "seqs = @seqs\n" if $printer == 1;
|
|
|
3245 }
|
|
|
3246 return "non" if checkIfSeqsIdentical(@seqs) eq "NO";
|
|
|
3247 # print "checking for $seqs[0] \n" if $printer == 1;
|
|
|
3248
|
|
|
3249 my @phases =();
|
|
|
3250 my %prephases = ();
|
|
|
3251 my %postphases = ();
|
|
|
3252 my $concat = $motif.$motif.$motif.$motif;
|
|
|
3253 my $motiflength = length($motif);
|
|
|
3254
|
|
|
3255 my $firstpass = 0;
|
|
|
3256
|
|
|
3257 for my $y (0 ... $motiflength-1){
|
|
|
3258 my $phase = substr($concat, $motiflength+$y, $motiflength);
|
|
|
3259 push @phases, $phase;
|
|
|
3260 $firstpass = 1 if $seqs[0] =~ /$phase/i;
|
|
|
3261 for my $k (0 ... length($motif)-1){
|
|
|
3262 my $pre = substr($concat, $motiflength+$y-$k, $k );
|
|
|
3263 my $post = substr($concat, $motiflength+$y+$motiflength, $k);
|
|
|
3264 # print "adding to phases : $phase - $pre and $post\n" if $printer == 1;
|
|
|
3265 push @{$prephases{$phase}} , $pre;
|
|
|
3266 push @{$postphases{$phase}} , $post;
|
|
|
3267 }
|
|
|
3268
|
|
|
3269 }
|
|
|
3270 # print "firstpass 1= $firstpass.. also, res-d = ",(length($seqs[0]))%(length($motif)),"\n" if $printer == 1;
|
|
|
3271 return "non" if $firstpass ==0;
|
|
|
3272 $firstpass =0;
|
|
|
3273 foreach my $phase (@phases){
|
|
|
3274
|
|
|
3275 $firstpass = 1 if $seqs[0] =~ /^($phase)+$/i && ((length($seqs[0]))%(length($motif))) == 0;
|
|
|
3276
|
|
|
3277 if (((length($seqs[0]))%(length($motif))) != 0){
|
|
|
3278 my @pres = @{$prephases{$phase}};
|
|
|
3279 my @posts = @{$postphases{$phase}};
|
|
|
3280 foreach my $pre (@pres){
|
|
|
3281 foreach my $post (@posts){
|
|
|
3282 next if $pre !~ /\S/ && $post !~ /\S/;
|
|
|
3283 $firstpass = 1 if ($seqs[0] =~ /^($pre)($phase)+($post)$/i || $seqs[0] =~ /^($pre)($phase)+$/i || $seqs[0] =~ /^($phase)+($post)$/i);
|
|
|
3284 # print "caught with $pre $phase $post\n" if $printer == 1;
|
|
|
3285 last if $firstpass == 1;
|
|
|
3286 }
|
|
|
3287 last if $firstpass == 1;
|
|
|
3288 }
|
|
|
3289 }
|
|
|
3290
|
|
|
3291 last if $firstpass == 1;
|
|
|
3292 }
|
|
|
3293
|
|
|
3294 #print "indel = $indel.. motif = $motif.. firstpass 2= mot\n" if $firstpass ==1;
|
|
|
3295 #print "indel = $indel.. motif = $motif.. firstpass 2= non\n" if $firstpass ==0;
|
|
|
3296 #<STDIN>;# if $firstpass ==1;
|
|
|
3297 return "non" if $firstpass ==0;
|
|
|
3298 return "mot" if $firstpass ==1;
|
|
|
3299
|
|
|
3300 }
|
|
|
3301
|
|
|
3302 sub checkIfSeqsIdentical{
|
|
|
3303 my @seqs = @_;
|
|
|
3304 my $identical = 1;
|
|
|
3305
|
|
|
3306 for my $j (1 ... $#seqs){
|
|
|
3307 $identical = 0 if uc($seqs[0]) ne uc($seqs[$j]);
|
|
|
3308 }
|
|
|
3309 return "NO" if $identical == 0;
|
|
|
3310 return "YES" if $identical == 1;
|
|
|
3311
|
|
|
3312 }
|
|
|
3313
|
|
|
3314 sub summarizeMutations{
|
|
|
3315 my $mutspt = $_[0];
|
|
|
3316 my @muts = @$mutspt;
|
|
|
3317 my $tree = $_[1];
|
|
|
3318
|
|
|
3319 my @returnarr = ();
|
|
|
3320
|
|
|
3321 for (1 ... 38){
|
|
|
3322 push @returnarr, "NA";
|
|
|
3323 }
|
|
|
3324 push @returnarr, "NULL";
|
|
|
3325 return @returnarr if $tree eq "NULL" || scalar(@muts) < 1;
|
|
|
3326
|
|
|
3327
|
|
|
3328 my @bspecies = ();
|
|
|
3329 my @dspecies = ();
|
|
|
3330 my $treecopy = $tree;
|
|
|
3331 $treecopy =~ s/[\(\)]//g;
|
|
|
3332 my @treeparts = split(/[\.,]+/, $treecopy);
|
|
|
3333
|
|
|
3334 for my $part (@treeparts){
|
|
|
3335 if ($part =~ /\+/){
|
|
|
3336 $part =~ s/\+//g;
|
|
|
3337 #my @sp = split(/\s*/, $part);
|
|
|
3338 #foreach my $p (@sp) {push @bspecies, $p;}
|
|
|
3339 push @bspecies, $part;
|
|
|
3340 }
|
|
|
3341 if ($part =~ /\-/){
|
|
|
3342 $part =~ s/\-//g;
|
|
|
3343 #my @sp = split(/\s*/, $part);
|
|
|
3344 #foreach my $p (@sp) {push @dspecies, $p;}
|
|
|
3345 push @dspecies, $part;
|
|
|
3346 }
|
|
|
3347
|
|
|
3348 }
|
|
|
3349 #print "-------------------------------------------------------\n";
|
|
|
3350
|
|
|
3351 my ($insertions, $deletions, $motinsertions, $motinsertionsf, $motdeletions, $motdeletionsf, $noninsertions, $nondeletions) = (0,0,0,0,0,0,0,0);
|
|
|
3352 my ($binsertions, $bdeletions, $bmotinsertions,$bmotinsertionsf, $bmotdeletions, $bmotdeletionsf, $bnoninsertions, $bnondeletions) = (0,0,0,0,0,0,0,0);
|
|
|
3353 my ($dinsertions, $ddeletions, $dmotinsertions,$dmotinsertionsf, $dmotdeletions, $dmotdeletionsf, $dnoninsertions, $dnondeletions) = (0,0,0,0,0,0,0,0);
|
|
|
3354 my ($ninsertions, $ndeletions, $nmotinsertions,$nmotinsertionsf, $nmotdeletions, $nmotdeletionsf, $nnoninsertions, $nnondeletions) = (0,0,0,0,0,0,0,0);
|
|
|
3355 my ($substitutions, $bsubstitutions, $dsubstitutions, $nsubstitutions, $indels, $subs) = (0,0,0,0,"NA","NA");
|
|
|
3356
|
|
|
3357 my @insertionsarr = (" ");
|
|
|
3358 my @deletionsarr = (" ");
|
|
|
3359
|
|
|
3360 my @substitutionsarr = (" ");
|
|
|
3361
|
|
|
3362
|
|
|
3363 foreach my $mut (@muts){
|
|
|
3364 # print "mut = $mut\n";
|
|
|
3365 chomp $mut;
|
|
|
3366 $mut =~ s/=\t/= /g;
|
|
|
3367 $mut =~ s/=$/= /g;
|
|
|
3368 my %mhash = ();
|
|
|
3369 my @mields = split(/\t/,$mut);
|
|
|
3370
|
|
|
3371 foreach my $m (@mields){
|
|
|
3372 my @fields = split(/=/,$m);
|
|
|
3373 next if $fields[1] eq " ";
|
|
|
3374 $mhash{$fields[0]} = $fields[1];
|
|
|
3375 }
|
|
|
3376
|
|
|
3377 my $myutype = ();
|
|
|
3378 my $decided = 0;
|
|
|
3379
|
|
|
3380 my $localnode = $mhash{"node"};
|
|
|
3381 $localnode =~ s/[\(\)\. ,]//g;
|
|
|
3382
|
|
|
3383
|
|
|
3384 foreach my $s (@bspecies){
|
|
|
3385 if ($localnode eq $s) {
|
|
|
3386 $decided = 1; $myutype = "b";
|
|
|
3387 }
|
|
|
3388 }
|
|
|
3389
|
|
|
3390 foreach my $s (@dspecies){
|
|
|
3391 if ($localnode eq $s) {
|
|
|
3392 $decided = 1; $myutype = "d";
|
|
|
3393 }
|
|
|
3394 }
|
|
|
3395
|
|
|
3396 $myutype = "n" if $decided != 1;
|
|
|
3397
|
|
|
3398
|
|
|
3399 # print "tree=$tree, birth species=@bspecies, death species=@dspecies, node=$mhash{node} .. myutype=$myutype .. \n";
|
|
|
3400 # <STDIN> if $mhash{"type"} eq "insertion" && $myutype eq "b";
|
|
|
3401
|
|
|
3402
|
|
|
3403 if ($mhash{"type"} eq "substitution"){
|
|
|
3404 $substitutions++;
|
|
|
3405 $bsubstitutions++ if $myutype eq "b";
|
|
|
3406 $dsubstitutions++ if $myutype eq "d";
|
|
|
3407 $nsubstitutions++ if $myutype eq "n";
|
|
|
3408 # print "substitution: from= $mhash{from}, to = $mhash{to}, and type = myutype\n";
|
|
|
3409 push @substitutionsarr, "$mhash{position}:".$mhash{"from"}.">".$mhash{"to"} if $myutype eq "b";
|
|
|
3410 push @substitutionsarr, "$mhash{position}:".$mhash{"from"}.">".$mhash{"to"} if $myutype eq "d";
|
|
|
3411 push @substitutionsarr, "$mhash{position}:".$mhash{"from"}.">".$mhash{"to"} if $myutype eq "n";
|
|
|
3412 # print "substitutionsarr = @substitutionsarr\n";
|
|
|
3413 # <STDIN>;
|
|
|
3414 }
|
|
|
3415 else{
|
|
|
3416 #print "tree=$tree, birth species=@bspecies, death species=@dspecies, node=$mhash{node} .. myutype=$myutype .. indeltype=$mhash{indeltype}\n";
|
|
|
3417 if ($mhash{"type"} eq "deletion"){
|
|
|
3418 $deletions++;
|
|
|
3419
|
|
|
3420 $motdeletions++ if $mhash{"indeltype"} =~ /dmot/;
|
|
|
3421 $motdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/;
|
|
|
3422
|
|
|
3423 $nondeletions++ if $mhash{"indeltype"} =~ /dnon/;
|
|
|
3424
|
|
|
3425 $bdeletions++ if $myutype eq "b";
|
|
|
3426 $ddeletions++ if $myutype eq "d";
|
|
|
3427 $ndeletions++ if $myutype eq "n";
|
|
|
3428
|
|
|
3429 $bmotdeletions++ if $mhash{"indeltype"} =~ /dmot/ && $myutype eq "b";
|
|
|
3430 $bmotdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/ && $myutype eq "b";
|
|
|
3431 $bnondeletions++ if $mhash{"indeltype"} =~ /dnon/ && $myutype eq "b";
|
|
|
3432
|
|
|
3433 $dmotdeletions++ if $mhash{"indeltype"} =~ /dmot/ && $myutype eq "d";
|
|
|
3434 $dmotdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/ && $myutype eq "d";
|
|
|
3435 $dnondeletions++ if $mhash{"indeltype"} =~ /dnon/ && $myutype eq "d";
|
|
|
3436
|
|
|
3437 $nmotdeletions++ if $mhash{"indeltype"} =~ /dmot/ && $myutype eq "n";
|
|
|
3438 $nmotdeletionsf++ if $mhash{"indeltype"} =~ /dmotf/ && $myutype eq "n";
|
|
|
3439 $nnondeletions++ if $mhash{"indeltype"} =~ /dnon/ && $myutype eq "n";
|
|
|
3440
|
|
|
3441 push @deletionsarr, "$mhash{indeltype}:$mhash{position}:".$mhash{"deletion"} if $myutype eq "b";
|
|
|
3442 push @deletionsarr, "$mhash{indeltype}:$mhash{position}:".$mhash{"deletion"} if $myutype eq "d";
|
|
|
3443 push @deletionsarr, "$mhash{indeltype}:$mhash{position}:".$mhash{"deletion"} if $myutype eq "n";
|
|
|
3444 }
|
|
|
3445
|
|
|
3446 if ($mhash{"type"} eq "insertion"){
|
|
|
3447 $insertions++;
|
|
|
3448
|
|
|
3449 $motinsertions++ if $mhash{"indeltype"} =~ /imot/;
|
|
|
3450 $motinsertionsf++ if $mhash{"indeltype"} =~ /imotf/;
|
|
|
3451 $noninsertions++ if $mhash{"indeltype"} =~ /inon/;
|
|
|
3452
|
|
|
3453 $binsertions++ if $myutype eq "b";
|
|
|
3454 $dinsertions++ if $myutype eq "d";
|
|
|
3455 $ninsertions++ if $myutype eq "n";
|
|
|
3456
|
|
|
3457 $bmotinsertions++ if $mhash{"indeltype"} =~ /imot/ && $myutype eq "b";
|
|
|
3458 $bmotinsertionsf++ if $mhash{"indeltype"} =~ /imotf/ && $myutype eq "b";
|
|
|
3459 $bnoninsertions++ if $mhash{"indeltype"} =~ /inon/ && $myutype eq "b";
|
|
|
3460
|
|
|
3461 $dmotinsertions++ if $mhash{"indeltype"} =~ /imot/ && $myutype eq "d";
|
|
|
3462 $dmotinsertionsf++ if $mhash{"indeltype"} =~ /imotf/ && $myutype eq "d";
|
|
|
3463 $dnoninsertions++ if $mhash{"indeltype"} =~ /inon/ && $myutype eq "d";
|
|
|
3464
|
|
|
3465 $nmotinsertions++ if $mhash{"indeltype"} =~ /imot/ && $myutype eq "n";
|
|
|
3466 $nmotinsertionsf++ if $mhash{"indeltype"} =~ /imotf/ && $myutype eq "n";
|
|
|
3467 $nnoninsertions++ if $mhash{"indeltype"} =~ /inon/ && $myutype eq "n";
|
|
|
3468
|
|
|
3469 push @insertionsarr, "$mhash{indeltype}:$mhash{position}:".$mhash{"insertion"} if $myutype eq "b";
|
|
|
3470 push @insertionsarr, "$mhash{indeltype}:$mhash{position}:".$mhash{"insertion"} if $myutype eq "d";
|
|
|
3471 push @insertionsarr, "$mhash{indeltype}:$mhash{position}:".$mhash{"insertion"} if $myutype eq "n";
|
|
|
3472
|
|
|
3473 }
|
|
|
3474 }
|
|
|
3475 }
|
|
|
3476
|
|
|
3477
|
|
|
3478
|
|
|
3479 $indels = "ins=".join(",",@insertionsarr).";dels=".join(",",@deletionsarr) if scalar(@insertionsarr) > 1 || scalar(@deletionsarr) > 1 ;
|
|
|
3480 $subs = join(",",@substitutionsarr) if scalar(@substitutionsarr) > 1;
|
|
|
3481 $indels =~ s/ //g;
|
|
|
3482 $subs =~ s/ //g ;
|
|
|
3483
|
|
|
3484 #print "indels = $indels, subs=$subs\n";
|
|
|
3485 ##<STDIN> if $indels =~ /[a-zA-Z0-9]/ || $subs =~ /[a-zA-Z0-9]/ ;
|
|
|
3486 #print "tree = $tree, indels = $indels, subs = $subs, bspecies = @bspecies, dspecies = @dspecies \n";
|
|
|
3487 my @returnarray = ();
|
|
|
3488
|
|
|
3489 push (@returnarray, $indels, $subs) ;
|
|
|
3490
|
|
|
3491 push @returnarray, $tree;
|
|
|
3492
|
|
|
3493 my @copy = @returnarray;
|
|
|
3494 #print "\n\nreturnarray = @returnarray ... binsertions=$binsertions dinsertions=$dinsertions bsubstitutions=$bsubstitutions dsubstitutions=$dsubstitutions\n";
|
|
|
3495 #<STDIN>;
|
|
|
3496 return (@returnarray);
|
|
|
3497
|
|
|
3498 }
|
|
|
3499
|
|
|
3500 sub selectBetterTree{
|
|
|
3501 my $printer = 1;
|
|
|
3502 my $treestudy = $_[0];
|
|
|
3503 my $alt = $_[1];
|
|
|
3504 my $mutspt = $_[2];
|
|
|
3505 my @muts = @$mutspt;
|
|
|
3506 my @trees = (); my @alternatetrees=();
|
|
|
3507
|
|
|
3508 @trees = split(/\|/,$treestudy) if $treestudy =~ /\|/;
|
|
|
3509 @alternatetrees = split(/[\|;]/,$alt) if $alt =~ /[\|;\(\)]/;
|
|
|
3510
|
|
|
3511 $trees[0] = $treestudy if $treestudy !~ /\|/;
|
|
|
3512 $alternatetrees[0] = $alt if $alt !~ /[\|;\(\)]/;
|
|
|
3513
|
|
|
3514 my @alltrees = (@trees, @alternatetrees);
|
|
|
3515 # push(@alltrees,@alternatetrees);
|
|
|
3516
|
|
|
3517 my %mutspecies = ();
|
|
|
3518 # print "IN selectBetterTree..treestudy=$treestudy. alt=$alt. for: @_. trees=@trees<. alternatetrees=@alternatetrees\n" if $printer == 1;
|
|
|
3519 #<STDIN>;
|
|
|
3520 foreach my $mut (@muts){
|
|
|
3521 # print colored ['green'],"mut = $mut\n" if $printer == 1;
|
|
|
3522 $mut =~ /node=([A-Z,\(\) ]+)/;
|
|
|
3523 my $node = $1;
|
|
|
3524 $node =~s/[,\(\) ]+//g;
|
|
|
3525 my @indivspecies = $node =~ /[A-Z]+/g;
|
|
|
3526 #print "adding node: $node\n" if $printer == 1;
|
|
|
3527 $mutspecies{$node} = $node;
|
|
|
3528
|
|
|
3529 }
|
|
|
3530
|
|
|
3531 my @treerecords = ();
|
|
|
3532 my $treecount = -1;
|
|
|
3533 foreach my $tree (@alltrees){
|
|
|
3534 # print "checking with tree $tree\n" if $printer == 1;
|
|
|
3535 $treecount++;
|
|
|
3536 $treerecords[$treecount] = 0;
|
|
|
3537 my @indivspecies = ($tree =~ /[A-Z]+/g);
|
|
|
3538 # print "indivspecies=@indivspecies\n" if $printer == 1;
|
|
|
3539 foreach my $species (@indivspecies){
|
|
|
3540 # print "checkin if exists species: $species\n" if $printer == 1;
|
|
|
3541 $treerecords[$treecount]+=2 if exists $mutspecies{$species} && $mutspecies{$species} !~ /indeltype=[a-z]mot/;
|
|
|
3542 $treerecords[$treecount]+=1.5 if exists $mutspecies{$species} && $mutspecies{$species} =~ /indeltype=[a-z]mot/;
|
|
|
3543 $treerecords[$treecount]-- if !exists $mutspecies{$species};
|
|
|
3544 }
|
|
|
3545 # print "for tree $tree, our treecount = $treerecords[$treecount]\n" if $printer == 1;
|
|
|
3546 }
|
|
|
3547
|
|
|
3548 my @best_tree = array_largest_number_arrayPosition(@treerecords);
|
|
|
3549 # print "treerecords = @treerecords. hence, best tree = @best_tree = $alltrees[$best_tree[0]], $treerecords[$best_tree[0]]\n" if $printer == 1;
|
|
|
3550 #<STDIN>;
|
|
|
3551 return ($alltrees[$best_tree[0]], $treerecords[$best_tree[0]]) if scalar(@best_tree) == 1;
|
|
|
3552 # print "best_tree[0] = $best_tree[0], and treerecords = $treerecords[$best_tree[0]]\n" if $printer == 1;
|
|
|
3553 return ("NULL", -1) if $treerecords[$best_tree[0]] < 1;
|
|
|
3554 my $rando = int(rand($#trees));
|
|
|
3555 return ($alltrees[$rando], $treerecords[$rando]) if scalar(@best_tree) > 1;
|
|
|
3556
|
|
|
3557 }
|
|
|
3558
|
|
|
3559
|
|
|
3560
|
|
|
3561
|
|
|
3562 sub load_sameHash{
|
|
|
3563 #my $g = %$_[0];
|
|
|
3564 $sameHash{"CAGT"}="AGTC";
|
|
|
3565 $sameHash{"ATGA"}="AATG";
|
|
|
3566 $sameHash{"CAAC"}="AACC";
|
|
|
3567 $sameHash{"GGAA"}="AAGG";
|
|
|
3568 $sameHash{"TAAG"}="AAGT";
|
|
|
3569 $sameHash{"CGAG"}="AGCG";
|
|
|
3570 $sameHash{"TAGG"}="AGGT";
|
|
|
3571 $sameHash{"GCAG"}="AGGC";
|
|
|
3572 $sameHash{"TAGA"}="ATAG";
|
|
|
3573 $sameHash{"TGA"}="ATG";
|
|
|
3574 $sameHash{"CAAG"}="AAGC";
|
|
|
3575 $sameHash{"CTAA"}="AACT";
|
|
|
3576 $sameHash{"CAAT"}="AATC";
|
|
|
3577 $sameHash{"GTAG"}="AGGT";
|
|
|
3578 $sameHash{"GAAG"}="AAGG";
|
|
|
3579 $sameHash{"CGA"}="ACG";
|
|
|
3580 $sameHash{"GTAA"}="AAGT";
|
|
|
3581 $sameHash{"ACAA"}="AAAC";
|
|
|
3582 $sameHash{"GCGG"}="GGGC";
|
|
|
3583 $sameHash{"ATCA"}="AATC";
|
|
|
3584 $sameHash{"TAAC"}="AACT";
|
|
|
3585 $sameHash{"GGCA"}="AGGC";
|
|
|
3586 $sameHash{"TGAG"}="AGTG";
|
|
|
3587 $sameHash{"AACA"}="AAAC";
|
|
|
3588 $sameHash{"GAGC"}="AGCG";
|
|
|
3589 $sameHash{"ACCA"}="AACC";
|
|
|
3590 $sameHash{"TGAA"}="AATG";
|
|
|
3591 $sameHash{"ACA"}="AAC";
|
|
|
3592 $sameHash{"GAAC"}="AACG";
|
|
|
3593 $sameHash{"GCA"}="AGC";
|
|
|
3594 $sameHash{"CCAC"}="ACCC";
|
|
|
3595 $sameHash{"CATA"}="ATAC";
|
|
|
3596 $sameHash{"CAC"}="ACC";
|
|
|
3597 $sameHash{"TACA"}="ATAC";
|
|
|
3598 $sameHash{"GGAC"}="ACGG";
|
|
|
3599 $sameHash{"AGA"}="AAG";
|
|
|
3600 $sameHash{"ATAA"}="AAAT";
|
|
|
3601 $sameHash{"CA"}="AC";
|
|
|
3602 $sameHash{"CCCA"}="ACCC";
|
|
|
3603 $sameHash{"TCAA"}="AATC";
|
|
|
3604 $sameHash{"CAGA"}="AGAC";
|
|
|
3605 $sameHash{"AATA"}="AAAT";
|
|
|
3606 $sameHash{"CCA"}="ACC";
|
|
|
3607 $sameHash{"AGAA"}="AAAG";
|
|
|
3608 $sameHash{"AGTA"}="AAGT";
|
|
|
3609 $sameHash{"GACG"}="ACGG";
|
|
|
3610 $sameHash{"TCAG"}="AGTC";
|
|
|
3611 $sameHash{"ACGA"}="AACG";
|
|
|
3612 $sameHash{"CGCA"}="ACGC";
|
|
|
3613 $sameHash{"GAGT"}="AGTG";
|
|
|
3614 $sameHash{"GA"}="AG";
|
|
|
3615 $sameHash{"TA"}="AT";
|
|
|
3616 $sameHash{"TAA"}="AAT";
|
|
|
3617 $sameHash{"CAG"}="AGC";
|
|
|
3618 $sameHash{"GATA"}="ATAG";
|
|
|
3619 $sameHash{"GTA"}="AGT";
|
|
|
3620 $sameHash{"CCAA"}="AACC";
|
|
|
3621 $sameHash{"TAG"}="AGT";
|
|
|
3622 $sameHash{"CAAA"}="AAAC";
|
|
|
3623 $sameHash{"AAGA"}="AAAG";
|
|
|
3624 $sameHash{"CACG"}="ACGC";
|
|
|
3625 $sameHash{"GTCA"}="AGTC";
|
|
|
3626 $sameHash{"GGA"}="AGG";
|
|
|
3627 $sameHash{"GGAT"}="ATGG";
|
|
|
3628 $sameHash{"CGGG"}="GGGC";
|
|
|
3629 $sameHash{"CGGA"}="ACGG";
|
|
|
3630 $sameHash{"AGGA"}="AAGG";
|
|
|
3631 $sameHash{"TAAA"}="AAAT";
|
|
|
3632 $sameHash{"GAGA"}="AGAG";
|
|
|
3633 $sameHash{"ACTA"}="AACT";
|
|
|
3634 $sameHash{"GCGA"}="AGCG";
|
|
|
3635 $sameHash{"CACA"}="ACAC";
|
|
|
3636 $sameHash{"AGAT"}="ATAG";
|
|
|
3637 $sameHash{"GAGG"}="AGGG";
|
|
|
3638 $sameHash{"CGAC"}="ACCG";
|
|
|
3639 $sameHash{"GGAG"}="AGGG";
|
|
|
3640 $sameHash{"GCCA"}="AGCC";
|
|
|
3641 $sameHash{"CCAG"}="AGCC";
|
|
|
3642 $sameHash{"GAAA"}="AAAG";
|
|
|
3643 $sameHash{"CAGG"}="AGGC";
|
|
|
3644 $sameHash{"GAC"}="ACG";
|
|
|
3645 $sameHash{"CAA"}="AAC";
|
|
|
3646 $sameHash{"GACC"}="ACCG";
|
|
|
3647 $sameHash{"GGCG"}="GGGC";
|
|
|
3648 $sameHash{"GGTA"}="AGGT";
|
|
|
3649 $sameHash{"AGCA"}="AAGC";
|
|
|
3650 $sameHash{"GATG"}="ATGG";
|
|
|
3651 $sameHash{"GTGA"}="AGTG";
|
|
|
3652 $sameHash{"ACAG"}="AGAC";
|
|
|
3653 $sameHash{"CGG"}="GGC";
|
|
|
3654 $sameHash{"ATA"}="AAT";
|
|
|
3655 $sameHash{"GACA"}="AGAC";
|
|
|
3656 $sameHash{"GCAA"}="AAGC";
|
|
|
3657 $sameHash{"CAGC"}="AGCC";
|
|
|
3658 $sameHash{"GGGA"}="AGGG";
|
|
|
3659 $sameHash{"GAG"}="AGG";
|
|
|
3660 $sameHash{"ACAT"}="ATAC";
|
|
|
3661 $sameHash{"GAAT"}="AATG";
|
|
|
3662 $sameHash{"CACC"}="ACCC";
|
|
|
3663 $sameHash{"GAT"}="ATG";
|
|
|
3664 $sameHash{"GCG"}="GGC";
|
|
|
3665 $sameHash{"GCAC"}="ACGC";
|
|
|
3666 $sameHash{"GAA"}="AAG";
|
|
|
3667 $sameHash{"TGGA"}="ATGG";
|
|
|
3668 $sameHash{"CCGA"}="ACCG";
|
|
|
3669 $sameHash{"CGAA"}="AACG";
|
|
|
3670 }
|
|
|
3671
|
|
|
3672
|
|
|
3673
|
|
|
3674 sub load_revHash{
|
|
|
3675 $revHash{"CTGA"}="AGTC";
|
|
|
3676 $revHash{"TCTT"}="AAAG";
|
|
|
3677 $revHash{"CTAG"}="AGCT";
|
|
|
3678 $revHash{"GGTG"}="ACCC";
|
|
|
3679 $revHash{"GCC"}="GGC";
|
|
|
3680 $revHash{"GCTT"}="AAGC";
|
|
|
3681 $revHash{"GCGT"}="ACGC";
|
|
|
3682 $revHash{"GTTG"}="AACC";
|
|
|
3683 $revHash{"CTCC"}="AGGG";
|
|
|
3684 $revHash{"ATC"}="ATG";
|
|
|
3685 $revHash{"CGAT"}="ATCG";
|
|
|
3686 $revHash{"TTAA"}="AATT";
|
|
|
3687 $revHash{"GTTC"}="AACG";
|
|
|
3688 $revHash{"CTGC"}="AGGC";
|
|
|
3689 $revHash{"TCGA"}="ATCG";
|
|
|
3690 $revHash{"ATCT"}="ATAG";
|
|
|
3691 $revHash{"GGTT"}="AACC";
|
|
|
3692 $revHash{"CTTA"}="AAGT";
|
|
|
3693 $revHash{"TGGC"}="AGCC";
|
|
|
3694 $revHash{"CCG"}="GGC";
|
|
|
3695 $revHash{"CGGC"}="GGCC";
|
|
|
3696 $revHash{"TTAG"}="AACT";
|
|
|
3697 $revHash{"GTG"}="ACC";
|
|
|
3698 $revHash{"CTTT"}="AAAG";
|
|
|
3699 $revHash{"TGCA"}="ATGC";
|
|
|
3700 $revHash{"CGCT"}="AGCG";
|
|
|
3701 $revHash{"TTCC"}="AAGG";
|
|
|
3702 $revHash{"CT"}="AG";
|
|
|
3703 $revHash{"C"}="G";
|
|
|
3704 $revHash{"CTCT"}="AGAG";
|
|
|
3705 $revHash{"ACTT"}="AAGT";
|
|
|
3706 $revHash{"GGTC"}="ACCG";
|
|
|
3707 $revHash{"ATTC"}="AATG";
|
|
|
3708 $revHash{"GGGT"}="ACCC";
|
|
|
3709 $revHash{"CCTA"}="AGGT";
|
|
|
3710 $revHash{"CGCG"}="GCGC";
|
|
|
3711 $revHash{"GTGT"}="ACAC";
|
|
|
3712 $revHash{"GCCC"}="GGGC";
|
|
|
3713 $revHash{"GTCG"}="ACCG";
|
|
|
3714 $revHash{"TCCC"}="AGGG";
|
|
|
3715 $revHash{"TTCA"}="AATG";
|
|
|
3716 $revHash{"AGTT"}="AACT";
|
|
|
3717 $revHash{"CCCT"}="AGGG";
|
|
|
3718 $revHash{"CCGC"}="GGGC";
|
|
|
3719 $revHash{"CTT"}="AAG";
|
|
|
3720 $revHash{"TTGG"}="AACC";
|
|
|
3721 $revHash{"ATT"}="AAT";
|
|
|
3722 $revHash{"TAGC"}="AGCT";
|
|
|
3723 $revHash{"ACTG"}="AGTC";
|
|
|
3724 $revHash{"TCAC"}="AGTG";
|
|
|
3725 $revHash{"CTGT"}="AGAC";
|
|
|
3726 $revHash{"TGTG"}="ACAC";
|
|
|
3727 $revHash{"ATCC"}="ATGG";
|
|
|
3728 $revHash{"GTGG"}="ACCC";
|
|
|
3729 $revHash{"TGGG"}="ACCC";
|
|
|
3730 $revHash{"TCGG"}="ACCG";
|
|
|
3731 $revHash{"CGGT"}="ACCG";
|
|
|
3732 $revHash{"GCTC"}="AGCG";
|
|
|
3733 $revHash{"TACG"}="ACGT";
|
|
|
3734 $revHash{"GTTT"}="AAAC";
|
|
|
3735 $revHash{"CAT"}="ATG";
|
|
|
3736 $revHash{"CATG"}="ATGC";
|
|
|
3737 $revHash{"GTTA"}="AACT";
|
|
|
3738 $revHash{"CACT"}="AGTG";
|
|
|
3739 $revHash{"TCAT"}="AATG";
|
|
|
3740 $revHash{"TTA"}="AAT";
|
|
|
3741 $revHash{"TGTA"}="ATAC";
|
|
|
3742 $revHash{"TTTC"}="AAAG";
|
|
|
3743 $revHash{"TACT"}="AAGT";
|
|
|
3744 $revHash{"TGTT"}="AAAC";
|
|
|
3745 $revHash{"CTA"}="AGT";
|
|
|
3746 $revHash{"GACT"}="AGTC";
|
|
|
3747 $revHash{"TTGC"}="AAGC";
|
|
|
3748 $revHash{"TTC"}="AAG";
|
|
|
3749 $revHash{"GCT"}="AGC";
|
|
|
3750 $revHash{"GCAT"}="ATGC";
|
|
|
3751 $revHash{"TGGT"}="AACC";
|
|
|
3752 $revHash{"CCT"}="AGG";
|
|
|
3753 $revHash{"CATC"}="ATGG";
|
|
|
3754 $revHash{"CCAT"}="ATGG";
|
|
|
3755 $revHash{"CCCG"}="GGGC";
|
|
|
3756 $revHash{"TGCC"}="AGGC";
|
|
|
3757 $revHash{"TG"}="AC";
|
|
|
3758 $revHash{"TGCT"}="AAGC";
|
|
|
3759 $revHash{"GCCG"}="GGCC";
|
|
|
3760 $revHash{"TCTG"}="AGAC";
|
|
|
3761 $revHash{"TGT"}="AAC";
|
|
|
3762 $revHash{"TTAT"}="AAAT";
|
|
|
3763 $revHash{"TAGT"}="AACT";
|
|
|
3764 $revHash{"TATG"}="ATAC";
|
|
|
3765 $revHash{"TTTA"}="AAAT";
|
|
|
3766 $revHash{"CGTA"}="ACGT";
|
|
|
3767 $revHash{"TA"}="AT";
|
|
|
3768 $revHash{"TGTC"}="AGAC";
|
|
|
3769 $revHash{"CTAT"}="ATAG";
|
|
|
3770 $revHash{"TATA"}="ATAT";
|
|
|
3771 $revHash{"TAC"}="AGT";
|
|
|
3772 $revHash{"TC"}="AG";
|
|
|
3773 $revHash{"CATT"}="AATG";
|
|
|
3774 $revHash{"TCG"}="ACG";
|
|
|
3775 $revHash{"ATTT"}="AAAT";
|
|
|
3776 $revHash{"CGTG"}="ACGC";
|
|
|
3777 $revHash{"CTG"}="AGC";
|
|
|
3778 $revHash{"TCGT"}="AACG";
|
|
|
3779 $revHash{"TCCG"}="ACGG";
|
|
|
3780 $revHash{"GTT"}="AAC";
|
|
|
3781 $revHash{"ATGT"}="ATAC";
|
|
|
3782 $revHash{"CTTG"}="AAGC";
|
|
|
3783 $revHash{"CCTT"}="AAGG";
|
|
|
3784 $revHash{"GATC"}="ATCG";
|
|
|
3785 $revHash{"CTGG"}="AGCC";
|
|
|
3786 $revHash{"TTCT"}="AAAG";
|
|
|
3787 $revHash{"CGTC"}="ACGG";
|
|
|
3788 $revHash{"CG"}="GC";
|
|
|
3789 $revHash{"TATT"}="AAAT";
|
|
|
3790 $revHash{"CTCG"}="AGCG";
|
|
|
3791 $revHash{"TCTC"}="AGAG";
|
|
|
3792 $revHash{"TCCT"}="AAGG";
|
|
|
3793 $revHash{"TGG"}="ACC";
|
|
|
3794 $revHash{"ACTC"}="AGTG";
|
|
|
3795 $revHash{"CTC"}="AGG";
|
|
|
3796 $revHash{"CGC"}="GGC";
|
|
|
3797 $revHash{"TTG"}="AAC";
|
|
|
3798 $revHash{"ACCT"}="AGGT";
|
|
|
3799 $revHash{"TCTA"}="ATAG";
|
|
|
3800 $revHash{"GTAC"}="ACGT";
|
|
|
3801 $revHash{"TTGA"}="AATC";
|
|
|
3802 $revHash{"GTCC"}="ACGG";
|
|
|
3803 $revHash{"GATT"}="AATC";
|
|
|
3804 $revHash{"T"}="A";
|
|
|
3805 $revHash{"CGTT"}="AACG";
|
|
|
3806 $revHash{"GTC"}="ACG";
|
|
|
3807 $revHash{"GCCT"}="AGGC";
|
|
|
3808 $revHash{"TGC"}="AGC";
|
|
|
3809 $revHash{"TTTG"}="AAAC";
|
|
|
3810 $revHash{"GGCT"}="AGCC";
|
|
|
3811 $revHash{"TCA"}="ATG";
|
|
|
3812 $revHash{"GTGC"}="ACGC";
|
|
|
3813 $revHash{"TGAT"}="AATC";
|
|
|
3814 $revHash{"TAT"}="AAT";
|
|
|
3815 $revHash{"CTAC"}="AGGT";
|
|
|
3816 $revHash{"TGCG"}="ACGC";
|
|
|
3817 $revHash{"CTCA"}="AGTG";
|
|
|
3818 $revHash{"CTTC"}="AAGG";
|
|
|
3819 $revHash{"GCTG"}="AGCC";
|
|
|
3820 $revHash{"TATC"}="ATAG";
|
|
|
3821 $revHash{"TAAT"}="AATT";
|
|
|
3822 $revHash{"ACT"}="AGT";
|
|
|
3823 $revHash{"TCGC"}="AGCG";
|
|
|
3824 $revHash{"GGT"}="ACC";
|
|
|
3825 $revHash{"TCC"}="AGG";
|
|
|
3826 $revHash{"TTGT"}="AAAC";
|
|
|
3827 $revHash{"TGAC"}="AGTC";
|
|
|
3828 $revHash{"TTAC"}="AAGT";
|
|
|
3829 $revHash{"CGT"}="ACG";
|
|
|
3830 $revHash{"ATTA"}="AATT";
|
|
|
3831 $revHash{"ATTG"}="AATC";
|
|
|
3832 $revHash{"CCTC"}="AGGG";
|
|
|
3833 $revHash{"CCGG"}="GGCC";
|
|
|
3834 $revHash{"CCGT"}="ACGG";
|
|
|
3835 $revHash{"TCCA"}="ATGG";
|
|
|
3836 $revHash{"CGCC"}="GGGC";
|
|
|
3837 $revHash{"GT"}="AC";
|
|
|
3838 $revHash{"TTCG"}="AACG";
|
|
|
3839 $revHash{"CCTG"}="AGGC";
|
|
|
3840 $revHash{"TCT"}="AAG";
|
|
|
3841 $revHash{"GTAT"}="ATAC";
|
|
|
3842 $revHash{"GTCT"}="AGAC";
|
|
|
3843 $revHash{"GCTA"}="AGCT";
|
|
|
3844 $revHash{"TACC"}="AGGT";
|
|
|
3845 }
|
|
|
3846
|
|
|
3847
|
|
|
3848 sub allCaps{
|
|
|
3849 my $motif = $_[0];
|
|
|
3850 $motif =~ s/a/A/g;
|
|
|
3851 $motif =~ s/c/C/g;
|
|
|
3852 $motif =~ s/t/T/g;
|
|
|
3853 $motif =~ s/g/G/g;
|
|
|
3854 return $motif;
|
|
|
3855 }
|
|
|
3856
|
|
|
3857
|
|
|
3858 sub all_caps{
|
|
|
3859 my @strand = split(/\s*/,$_[0]);
|
|
|
3860 for my $i (0 ... $#strand){
|
|
|
3861 if ($strand[$i] =~ /c/) {$strand[$i] = "C";next;}
|
|
|
3862 if ($strand[$i] =~ /a/) {$strand[$i] = "A";next;}
|
|
|
3863 if ($strand[$i] =~ /t/) { $strand[$i] = "T";next;}
|
|
|
3864 if ($strand[$i] =~ /g/) {$strand[$i] = "G";next;}
|
|
|
3865 }
|
|
|
3866 return join("",@strand);
|
|
|
3867 }
|
|
|
3868 sub array_mean{
|
|
|
3869 return "NA" if scalar(@_) == 0;
|
|
|
3870 my $sum = 0;
|
|
|
3871 foreach my $val (@_){
|
|
|
3872 $sum = $sum + $val;
|
|
|
3873 }
|
|
|
3874 return ($sum/scalar(@_));
|
|
|
3875 }
|
|
|
3876 sub array_sum{
|
|
|
3877 return "NA" if scalar(@_) == 0;
|
|
|
3878 my $sum = 0;
|
|
|
3879 foreach my $val (@_){
|
|
|
3880 $sum = $sum + $val;
|
|
|
3881 }
|
|
|
3882 return ($sum);
|
|
|
3883 }
|
|
|
3884
|
|
|
3885 sub variance{
|
|
|
3886 return "NA" if scalar(@_) == 0;
|
|
|
3887 return 0 if scalar(@_) == 1;
|
|
|
3888 my $mean = array_mean(@_);
|
|
|
3889 my $num = 0;
|
|
|
3890 return 0 if scalar(@_) == 1;
|
|
|
3891 # print "mean = $mean .. array = >@_<\n";
|
|
|
3892 foreach my $ele (@_){
|
|
|
3893 # print "$num = $num + ($ele-$mean)*($ele-$mean)\n";
|
|
|
3894 $num = $num + ($ele-$mean)*($ele-$mean);
|
|
|
3895 }
|
|
|
3896 my $var = $num / scalar(@_);
|
|
|
3897 return $var;
|
|
|
3898 }
|
|
|
3899
|
|
|
3900 sub array_95confIntervals{
|
|
|
3901 return "NA" if scalar(@_) <= 0;
|
|
|
3902 my @sorted = sort { $a <=> $b } @_;
|
|
|
3903 # print "@sorted=",scalar(@sorted), "\n";
|
|
|
3904 my $aDeechNo = int((scalar(@sorted) * 2.5) / 100);
|
|
|
3905 my $saaDeNo = int((scalar(@sorted) * 97.5) / 100);
|
|
|
3906
|
|
|
3907 return ($sorted[$aDeechNo], $sorted[$saaDeNo]);
|
|
|
3908 }
|
|
|
3909
|
|
|
3910 sub array_median{
|
|
|
3911 return "NA" if scalar(@_) == 0;
|
|
|
3912 return $_[0] if scalar(@_) == 1;
|
|
|
3913 my @sorted = sort { $a <=> $b } @_;
|
|
|
3914 my $totalno = scalar(@sorted);
|
|
|
3915
|
|
|
3916 #print "sorted = @sorted\n";
|
|
|
3917
|
|
|
3918 my $pick = ();
|
|
|
3919 if ($totalno % 2 == 1){
|
|
|
3920 #print "odd set .. totalno = $totalno\n";
|
|
|
3921 my $mid = $totalno / 2;
|
|
|
3922 my $onehalfno = $mid - $mid % 1;
|
|
|
3923 my $secondhalfno = $onehalfno + 1;
|
|
|
3924 my $onehalf = $sorted[$onehalfno-1];
|
|
|
3925 my $secondhalf = $sorted[$secondhalfno-1];
|
|
|
3926 #print "onehalfno = $onehalfno and secondhalfno = $secondhalfno \n onehalf = $onehalf and secondhalf = $secondhalf\n";
|
|
|
3927
|
|
|
3928 $pick = $secondhalf;
|
|
|
3929 }
|
|
|
3930 else{
|
|
|
3931 #print "even set .. totalno = $totalno\n";
|
|
|
3932 my $mid = $totalno / 2;
|
|
|
3933 my $onehalfno = $mid;
|
|
|
3934 my $secondhalfno = $onehalfno + 1;
|
|
|
3935 my $onehalf = $sorted[$onehalfno-1];
|
|
|
3936 my $secondhalf = $sorted[$secondhalfno-1];
|
|
|
3937 #print "onehalfno = $onehalfno and secondhalfno = $secondhalfno \n onehalf = $onehalf and secondhalf = $secondhalf\n";
|
|
|
3938 $pick = ($onehalf + $secondhalf )/2;
|
|
|
3939
|
|
|
3940 }
|
|
|
3941 #print "pick = $pick..\n";
|
|
|
3942 return $pick;
|
|
|
3943
|
|
|
3944 }
|
|
|
3945
|
|
|
3946
|
|
|
3947 sub array_numerical_sort{
|
|
|
3948 return "NA" if scalar(@_) == 0;
|
|
|
3949 my @sorted = sort { $a <=> $b } @_;
|
|
|
3950 return (@sorted);
|
|
|
3951 }
|
|
|
3952
|
|
|
3953 sub array_smallest_number{
|
|
|
3954 return "NA" if scalar(@_) == 0;
|
|
|
3955 return $_[0] if scalar(@_) == 1;
|
|
|
3956 my @sorted = sort { $a <=> $b } @_;
|
|
|
3957 return $sorted[0];
|
|
|
3958 }
|
|
|
3959
|
|
|
3960
|
|
|
3961 sub array_largest_number{
|
|
|
3962 return "NA" if scalar(@_) == 0;
|
|
|
3963 return $_[0] if scalar(@_) == 1;
|
|
|
3964 my @sorted = sort { $a <=> $b } @_;
|
|
|
3965 return $sorted[$#sorted];
|
|
|
3966 }
|
|
|
3967
|
|
|
3968
|
|
|
3969 sub array_largest_number_arrayPosition{
|
|
|
3970 return "NA" if scalar(@_) == 0;
|
|
|
3971 return 0 if scalar(@_) == 1;
|
|
|
3972 my $maxpos = 0;
|
|
|
3973 my @maxposes = ();
|
|
|
3974 my @maxvals = ();
|
|
|
3975 my $maxval = array_smallest_number(@_);
|
|
|
3976 for my $i (0 ... $#_){
|
|
|
3977 if ($_[$i] > $maxval){
|
|
|
3978 $maxval = $_[$i];
|
|
|
3979 $maxpos = $i;
|
|
|
3980 }
|
|
|
3981 if ($_[$i] == $maxval){
|
|
|
3982 $maxval = $_[$i];
|
|
|
3983 if (scalar(@maxposes) == 0){
|
|
|
3984 push @maxposes, $i;
|
|
|
3985 push @maxvals, $_[$i];
|
|
|
3986
|
|
|
3987 }
|
|
|
3988 elsif ($maxvals[0] == $maxval){
|
|
|
3989 push @maxposes, $i;
|
|
|
3990 push @maxvals, $_[$i];
|
|
|
3991 }
|
|
|
3992 else{
|
|
|
3993 @maxposes = (); @maxvals = ();
|
|
|
3994 push @maxposes, $i;
|
|
|
3995 push @maxvals, $_[$i];
|
|
|
3996 }
|
|
|
3997
|
|
|
3998 }
|
|
|
3999
|
|
|
4000 }
|
|
|
4001 return $maxpos if scalar(@maxposes) < 2;
|
|
|
4002 return (@maxposes);
|
|
|
4003 }
|
|
|
4004
|
|
|
4005 sub array_smallest_number_arrayPosition{
|
|
|
4006 return "NA" if scalar(@_) == 0;
|
|
|
4007 return 0 if scalar(@_) == 1;
|
|
|
4008 my $minpos = 0;
|
|
|
4009 my @minposes = ();
|
|
|
4010 my @minvals = ();
|
|
|
4011 my $minval = array_largest_number(@_);
|
|
|
4012 my $maxval = array_smallest_number(@_);
|
|
|
4013 #print "starting with $maxval, ending with $minval\n";
|
|
|
4014 for my $i (0 ... $#_){
|
|
|
4015 if ($_[$i] < $minval){
|
|
|
4016 $minval = $_[$i];
|
|
|
4017 $minpos = $i;
|
|
|
4018 }
|
|
|
4019 if ($_[$i] == $minval){
|
|
|
4020 $minval = $_[$i];
|
|
|
4021 if (scalar(@minposes) == 0){
|
|
|
4022 push @minposes, $i;
|
|
|
4023 push @minvals, $_[$i];
|
|
|
4024
|
|
|
4025 }
|
|
|
4026 elsif ($minvals[0] == $minval){
|
|
|
4027 push @minposes, $i;
|
|
|
4028 push @minvals, $_[$i];
|
|
|
4029 }
|
|
|
4030 else{
|
|
|
4031 @minposes = (); @minvals = ();
|
|
|
4032 push @minposes, $i;
|
|
|
4033 push @minvals, $_[$i];
|
|
|
4034 }
|
|
|
4035
|
|
|
4036 }
|
|
|
4037
|
|
|
4038 }
|
|
|
4039 #print "minposes=@minposes\n";
|
|
|
4040
|
|
|
4041 return $minpos if scalar(@minposes) < 2;
|
|
|
4042 return (@minposes);
|
|
|
4043 }
|
|
|
4044
|
|
|
4045 sub basic_stats{
|
|
|
4046 my @arr = @_;
|
|
|
4047 # print " array_smallest_number= ", array_smallest_number(@arr)," array_largest_number= ", array_largest_number(@arr), " array_mean= ",array_mean(@arr),"\n";
|
|
|
4048 return ":";
|
|
|
4049 }
|
|
|
4050 #xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx
|
|
|
4051
|
|
|
4052 sub maftoAxt_multispecies {
|
|
|
4053 my $printer = 0;
|
|
|
4054 #print "in maftoAxt_multispecies : got @_\n";
|
|
|
4055 my $fname=$_[0];
|
|
|
4056 open(IN,"<$_[0]") or die "Cannot open $_[0]: $! \n";
|
|
|
4057 my $treedefinition = $_[1];
|
|
|
4058 #print "treedefinition= $treedefinition\n";
|
|
|
4059
|
|
|
4060 my @treedefinitions = MakeTrees($treedefinition);
|
|
|
4061
|
|
|
4062
|
|
|
4063 open(OUT,">$_[2]") or die "Cannot open $_[2]: $! \n";
|
|
|
4064 my $counter = 0;
|
|
|
4065 my $exactspeciesset = $_[3];
|
|
|
4066 my @exactspeciesset_unarranged = split(/,/,$exactspeciesset);
|
|
|
4067
|
|
|
4068 $treedefinition=~s/[\)\(, ]/\t/g;
|
|
|
4069 my @species=split(/\t+/,$treedefinition);
|
|
|
4070 @exactspeciesset_unarranged = @species;
|
|
|
4071 # print "species=@species\n";
|
|
|
4072
|
|
|
4073 my @exactspeciesarr=();
|
|
|
4074
|
|
|
4075
|
|
|
4076 foreach my $def (@treedefinitions){
|
|
|
4077 $def=~s/[\)\(, ]/\t/g;
|
|
|
4078 my @specs = split(/\t+/,$def);
|
|
|
4079 my @exactspecies=();
|
|
|
4080 foreach my $spec (@specs){
|
|
|
4081 foreach my $espec (@exactspeciesset_unarranged){
|
|
|
4082 # print "pushing >$spec< nd >$espec<\n" if $spec eq $espec && $espec =~ /[a-zA-Z0-9]/;
|
|
|
4083 push @exactspecies, $spec if $spec eq $espec && $espec =~ /[a-zA-Z0-9]/;
|
|
|
4084 }
|
|
|
4085
|
|
|
4086 }
|
|
|
4087 #print "exactspecies = >@exactspecies<\n";
|
|
|
4088 push @exactspeciesarr, [@exactspecies];
|
|
|
4089 }
|
|
|
4090 #<STDIN>;
|
|
|
4091 #print "exactspeciesarr=@exactspeciesarr\n";
|
|
|
4092
|
|
|
4093 ###########
|
|
|
4094 my $select = 2;
|
|
|
4095 #select = 1 if all species need sequences to be present for each block otherwise, it is 0
|
|
|
4096 #select = 2 only the allowed set make up the alignment. use the removeset
|
|
|
4097 # information to detect alignmenets that have other important genomes aligned.
|
|
|
4098 ###########
|
|
|
4099 my @allowedset = ();
|
|
|
4100 @allowedset = split(/;/,allowedSetOfSpecies(join("_",@species))) if $select == 0;
|
|
|
4101 @allowedset = join("_",0,@species) if $select == 1;
|
|
|
4102 #print "species = @species , allowedset =",join("\n", @allowedset) ," \n";
|
|
|
4103
|
|
|
4104
|
|
|
4105 foreach my $set (@exactspeciesarr){
|
|
|
4106 my @openset = @$set;
|
|
|
4107 push @allowedset, join("_",0,@openset) if $select == 2;
|
|
|
4108 # print "openset = >@openset<, allowedset = @allowedset and exactspecies = @exactspecies\n";
|
|
|
4109 }
|
|
|
4110 # <STDIN>;
|
|
|
4111 my $start = 0;
|
|
|
4112 my @sequences = ();
|
|
|
4113 my @titles = ();
|
|
|
4114 my $species_counter = "0";
|
|
|
4115 my $countermatch = 0;
|
|
|
4116 my $outsideSpecies=0;
|
|
|
4117
|
|
|
4118 while(my $line = <IN>){
|
|
|
4119 next if $line =~ /^#/;
|
|
|
4120 next if $line =~ /^i/;
|
|
|
4121 # print "$line .. species = @species\n";
|
|
|
4122 chomp $line;
|
|
|
4123 my @fields = split(/\s+/,$line);
|
|
|
4124 chomp $line;
|
|
|
4125 if ($line =~ /^a /){
|
|
|
4126 $start = 1;
|
|
|
4127 }
|
|
|
4128
|
|
|
4129 if ($line =~ /^s /){
|
|
|
4130 # print "fields1 = $fields[1] , start = $start\n";
|
|
|
4131
|
|
|
4132 foreach my $sp (@species){
|
|
|
4133 if ($fields[1] =~ /$sp/){
|
|
|
4134 $species_counter = $species_counter."_".$sp;
|
|
|
4135 push(@sequences, $fields[6]);
|
|
|
4136 my @sp_info = split(/\./,$fields[1]);
|
|
|
4137 my $title = join(" ",@sp_info, $fields[2], ($fields[2]+$fields[3]), $fields[4]);
|
|
|
4138 push(@titles, $title);
|
|
|
4139
|
|
|
4140 }
|
|
|
4141 }
|
|
|
4142 }
|
|
|
4143
|
|
|
4144 if (($line !~ /^a/) && ($line !~ /^s/) && ($line !~ /^#/) && ($line !~ /^i/) && ($start = 1)){
|
|
|
4145
|
|
|
4146 my $arranged = reorderSpecies($species_counter, @species);
|
|
|
4147 # print "species_counter=$species_counter .. arranged = $arranged\n";
|
|
|
4148
|
|
|
4149 my $stopper = 1;
|
|
|
4150 my $arrno = 0;
|
|
|
4151 foreach my $set (@allowedset){
|
|
|
4152 # print "checking $set with $arranged\n";
|
|
|
4153 if ($arranged eq $set){
|
|
|
4154 # print "checked $set with $arranged\n";
|
|
|
4155 $stopper = 0; last;
|
|
|
4156 }
|
|
|
4157 $arrno++;
|
|
|
4158 }
|
|
|
4159
|
|
|
4160 if ($stopper == 0) {
|
|
|
4161 # print " accepted\n";
|
|
|
4162 @titles = split ";", orderInfo(join(";", @titles), $species_counter, $arranged) if $species_counter ne $arranged;
|
|
|
4163
|
|
|
4164 @sequences = split ";", orderInfo(join(";", @sequences), $species_counter, $arranged) if $species_counter ne $arranged;
|
|
|
4165 my $filteredseq = filter_gaps(@sequences);
|
|
|
4166
|
|
|
4167 if ($filteredseq ne "SHORT"){
|
|
|
4168 $counter++;
|
|
|
4169 print OUT join (" ",$counter, @titles), "\n";
|
|
|
4170 print OUT $filteredseq, "\n";
|
|
|
4171 print OUT "\n";
|
|
|
4172 $countermatch++;
|
|
|
4173 }
|
|
|
4174 # my @filtered_seq = split(/\t/,filter_gaps(@sequences) );
|
|
|
4175 }
|
|
|
4176 else{#print "\n";
|
|
|
4177 }
|
|
|
4178
|
|
|
4179 @sequences = (); @titles = (); $start = 0;$species_counter = "0";
|
|
|
4180 next;
|
|
|
4181
|
|
|
4182 }
|
|
|
4183 }
|
|
|
4184 close IN;
|
|
|
4185 close OUT;
|
|
|
4186 # print "countermatch = $countermatch\n";
|
|
|
4187 # <STDIN>;
|
|
|
4188 }
|
|
|
4189
|
|
|
4190 sub reorderSpecies{
|
|
|
4191 my @inarr=@_;
|
|
|
4192 my $currSpecies = shift (@inarr);
|
|
|
4193 my $ordered_species = 0;
|
|
|
4194 my @species=@inarr;
|
|
|
4195 foreach my $order (@species){
|
|
|
4196 $ordered_species = $ordered_species."_".$order if $currSpecies=~ /$order/;
|
|
|
4197 }
|
|
|
4198 return $ordered_species;
|
|
|
4199
|
|
|
4200 }
|
|
|
4201
|
|
|
4202 sub filter_gaps{
|
|
|
4203 my @sequences = @_;
|
|
|
4204 # print "sequences sent are @sequences\n";
|
|
|
4205 my $seq_length = length($sequences[0]);
|
|
|
4206 my $seq_no = scalar(@sequences);
|
|
|
4207 my $allgaps = ();
|
|
|
4208 for (1 ... $seq_no){
|
|
|
4209 $allgaps = $allgaps."-";
|
|
|
4210 }
|
|
|
4211
|
|
|
4212 my @seq_array = ();
|
|
|
4213 my $seq_counter = 0;
|
|
|
4214 foreach my $seq (@sequences){
|
|
|
4215 # my @sequence = split(/\s*/,$seq);
|
|
|
4216 $seq_array[$seq_counter] = [split(/\s*/,$seq)];
|
|
|
4217 # push @seq_array, [@sequence];
|
|
|
4218 $seq_counter++;
|
|
|
4219 }
|
|
|
4220 my $g = 0;
|
|
|
4221 while ( $g < $seq_length){
|
|
|
4222 last if (!exists $seq_array[0][$g]);
|
|
|
4223 my $bases = ();
|
|
|
4224 for my $u (0 ... $#seq_array){
|
|
|
4225 $bases = $bases.$seq_array[$u][$g];
|
|
|
4226 }
|
|
|
4227 # print $bases, "\n";
|
|
|
4228 if ($bases eq $allgaps){
|
|
|
4229 # print "bases are $bases, position is $g \n";
|
|
|
4230 for my $seq (@seq_array){
|
|
|
4231 splice(@$seq , $g, 1);
|
|
|
4232 }
|
|
|
4233 }
|
|
|
4234 else {
|
|
|
4235 $g++;
|
|
|
4236 }
|
|
|
4237 }
|
|
|
4238
|
|
|
4239 my @outs = ();
|
|
|
4240
|
|
|
4241 foreach my $seq (@seq_array){
|
|
|
4242 push(@outs, join("",@$seq));
|
|
|
4243 }
|
|
|
4244 return "SHORT" if length($outs[0]) <=100;
|
|
|
4245 return (join("\n", @outs));
|
|
|
4246 }
|
|
|
4247
|
|
|
4248
|
|
|
4249 sub allowedSetOfSpecies{
|
|
|
4250 my @allowed_species = split(/_/,$_[0]);
|
|
|
4251 unshift @allowed_species, 0;
|
|
|
4252 # print "allowed set = @allowed_species \n";
|
|
|
4253 my @output = ();
|
|
|
4254 for (0 ... scalar(@allowed_species) - 4){
|
|
|
4255 push(@output, join("_",@allowed_species));
|
|
|
4256 pop @allowed_species;
|
|
|
4257 }
|
|
|
4258 return join(";",reverse(@output));
|
|
|
4259
|
|
|
4260 }
|
|
|
4261
|
|
|
4262
|
|
|
4263 sub orderInfo{
|
|
|
4264 my @info = split(/;/,$_[0]);
|
|
|
4265 # print "info = @info";
|
|
|
4266 my @old = split(/_/,$_[1]);
|
|
|
4267 my @new = split(/_/,$_[2]);
|
|
|
4268 shift @old; shift @new;
|
|
|
4269 my @outinfo = ();
|
|
|
4270 foreach my $spe (@new){
|
|
|
4271 for my $no (0 ... $#old){
|
|
|
4272 if ($spe eq $old[$no]){
|
|
|
4273 push(@outinfo, $info[$no]);
|
|
|
4274 }
|
|
|
4275 }
|
|
|
4276 }
|
|
|
4277 # print "outinfo = @outinfo \n";
|
|
|
4278 return join(";", @outinfo);
|
|
|
4279 }
|
|
|
4280
|
|
|
4281 #xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx
|
|
|
4282
|
|
|
4283 sub printarr {
|
|
|
4284 # print ">::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\n";
|
|
|
4285 # foreach my $line (@_) {print "$line\n";}
|
|
|
4286 # print "::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::<\n";
|
|
|
4287 }
|
|
|
4288
|
|
|
4289 sub oneOf{
|
|
|
4290 my @arr = @_;
|
|
|
4291 my $element = $arr[$#arr];
|
|
|
4292 @arr = @arr[0 ... $#arr-1];
|
|
|
4293 my $present = 0;
|
|
|
4294
|
|
|
4295 foreach my $el (@arr){
|
|
|
4296 $present = 1 if $el eq $element;
|
|
|
4297 }
|
|
|
4298 return $present;
|
|
|
4299 }
|
|
|
4300
|
|
|
4301 #xxxxxxxxxxxxxx MakeTrees xxxxxxxxxxxxxxxxxxxxxxxxxxxx MakeTrees xxxxxxxxxxxxxxxxxxxxxxxxxxxx MakeTrees xxxxxxxxxxxxxxxxxxxxxxxxxxxx
|
|
|
4302
|
|
|
4303 sub MakeTrees{
|
|
|
4304 my $tree = $_[0];
|
|
|
4305 # my @parts=($tree);
|
|
|
4306 my @parts=();
|
|
|
4307
|
|
|
4308 # print "parts=@parts\n";
|
|
|
4309
|
|
|
4310 while (1){
|
|
|
4311 $tree =~ s/^\(//g;
|
|
|
4312 $tree =~ s/\)$//g;
|
|
|
4313 my @arr = ();
|
|
|
4314
|
|
|
4315 if ($tree =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_\(\),]+)\)$/){
|
|
|
4316 @arr = $tree =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_\(\),]+)$/;
|
|
|
4317 push @parts, "(".$tree.")";
|
|
|
4318 $tree = $2;
|
|
|
4319 }
|
|
|
4320 elsif ($tree =~ /^\(([a-zA-Z0-9_\(\),]+),([a-zA-Z0-9_]+)$/){
|
|
|
4321 @arr = $tree =~ /^([a-zA-Z0-9_\(\),]+),([a-zA-Z0-9_]+)$/;
|
|
|
4322 push @parts, "(".$tree.")";
|
|
|
4323 $tree = $1;
|
|
|
4324 }
|
|
|
4325 elsif ($tree =~ /^([a-zA-Z0-9_]+),([a-zA-Z0-9_]+)$/){
|
|
|
4326 last;
|
|
|
4327 }
|
|
|
4328 }
|
|
|
4329 #print "parts=@parts\n";
|
|
|
4330 return @parts;
|
|
|
4331 }
|
|
|
4332
|
|
|
4333
|
