Mercurial > repos > xuebing > sharplabtool

comparison tools/regVariation/parseMAF_smallIndels.pl @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
1 #!/usr/bin/perl -w
2 # a program to get indels
3 # input is a MAF format 3-way alignment file
4 # from 3-way blocks only at this time
5 # translate seq2, seq3, etc coordinates to + if align orient is reverse complement
6
7 use strict;
8 use warnings;
9
10 # declare and initialize variables
11 my $fh; # variable to store filehandle
12 my $record;
13 my $offset;
14 my $library = $ARGV[0];
15 my $count = 0;
16 my $count2 = 0;
17 my $count3 = 0;
18 my $count4 = 0;
19 my $start1 = my $start2 = my $start3 = my $start4 = my $start5 = my $start6 = 0;
20 my $orient = "";
21 my $outgroup = $ARGV[2];
22 my $ingroup1 = my $ingroup2 = "";
23 my $count_seq1insert = my $count_seq1delete = 0;
24 my $count_seq2insert = my $count_seq2delete = 0;
25 my $count_seq3insert = my $count_seq3delete = 0;
26 my @seq1_insert_lengths = my @seq1_delete_lengths = ();
27 my @seq2_insert_lengths = my @seq2_delete_lengths = ();
28 my @seq3_insert_lengths = my @seq3_delete_lengths = ();
29 my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = ();
30 my @seq1_insert_startOnly = my @seq1_delete_startOnly = my @seq2_insert_startOnly = my @seq2_delete_startOnly = ();
31 my @seq3_insert_startOnly = my @seq3_delete_startOnly = ();
32 my @indels = ();
33
34 # check to make sure correct files
35 my $usage = "usage: parseMAF_smallIndels.pl [MAF.in] [small_Indels_summary.out] [outgroup]\n";
36 die $usage unless @ARGV == 3;
37
38 # perform some standard subroutines
39 $fh = open_file($library);
40
41 $offset = tell($fh);
42
43 #my $ofile = $ARGV[2];
44 #unless (open(OFILE, ">$ofile")){
45 # print "Cannot open output file \"$ofile\"\n\n";
46 # exit;
47 #}
48
49 my $ofile2 = $ARGV[1];
50 unless (open(OFILE2, ">$ofile2")){
51 print "Cannot open output file \"$ofile2\"\n\n";
52 exit;
53 }
54
55
56 # header line for output files
57 #print OFILE "# small indel events, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n";
58 #print OFILE "#align\tingroup1\tingroup1_coord\tingroup1_orient\tingroup2\tingroup2_coord\tingroup2_orient\toutgroup\toutgroup_coord\toutgroup_orient\tindel_type\n";
59
60 #print OFILE2 "# small indels summary, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n";
61 print OFILE2 "#block\tindel_type\tindel_length\tingroup1\tingroup1_start\tingroup1_end\tingroup1_alignSize\tingroup1_orient\tingroup2\tingroup2_start\tingroup2_end\tingroup2_alignSize\tingroup2_orient\toutgroup\toutgroup_start\toutgroup_end\toutgroup_alignSize\toutgroup_orient\n";
62
63 # main body of program
64 while ($record = get_next_record($fh) ){
65 if ($record=~ m/\s*##maf(.*)\s*# maf/s){
66 next;
67 }
68
69 my @sequences = get_sequences_within_block($record);
70 my @seq_info = get_indels_within_block(@sequences);
71 get_indels_lengths(@seq_info);
72
73 $offset = tell($fh);
74 $count++;
75
76 }
77
78 get_starts_only(@seq1_insert);
79 get_starts_only(@seq1_delete);
80 get_starts_only(@seq2_insert);
81 get_starts_only(@seq2_delete);
82 get_starts_only(@seq3_insert);
83 get_starts_only(@seq3_delete);
84
85 # print some things to keep track of progress
86 #print "# $library\n";
87 #print "# number of records = $count\n";
88 #print "# number of sequence \"s\" lines = $count2\n";
89 if ($count3 != 0){
90 print "Skipped $count3 blocks with only 2 seqs;\n";
91 }
92 #print "# number of records with only h-m = $count4\n\n";
93
94 print "Ingroup1 = $ingroup1; Ingroup2 = $ingroup2; Outgroup = $outgroup;\n";
95 print "# of ingroup1 inserts = $count_seq1insert;\n";
96 print "# of ingroup1 deletes = $count_seq1delete;\n";
97 print "# of ingroup2 inserts = $count_seq2insert;\n";
98 print "# of ingroup2 deletes = $count_seq2delete;\n";
99 print "# of outgroup3 inserts = $count_seq3insert;\n";
100 print "# of outgroup3 deletes = $count_seq3delete\n";
101
102
103 #close OFILE;
104
105 if ($count == $count3){
106 print STDERR "Skipped all blocks since none of them contain 3-way alignments.\n";
107 exit -1;
108 }
109
110 ###################SUBROUTINES#####################################
111
112 # subroutine to open file
113 sub open_file {
114 my($filename) = @_;
115 my $fh;
116
117 unless (open($fh, $filename)){
118 print "Cannot open file $filename\n";
119 exit;
120 }
121 return $fh;
122 }
123
124 # get next record
125 sub get_next_record {
126 my ($fh) = @_;
127 my ($offset);
128 my ($record) = "";
129 my ($save_input_separator) = $/;
130
131 $/ = "a score";
132
133 $record = <$fh>;
134
135 $/ = $save_input_separator;
136 return $record;
137 }
138
139 # get header info
140 sub get_sequences_within_block{
141 my (@alignment) = @_;
142 my @lines = ();
143
144 my @sequences = ();
145
146 @lines = split ("\n", $record);
147 foreach (@lines){
148 chomp($_);
149 if (m/^\s*$/){
150 next;
151 }
152 elsif (m/^\s*=(\d+\.*\d*)/){
153
154 }elsif (m/^\s*s(.*)$/){
155 $count2++;
156
157 push (@sequences,$_);
158 }
159 }
160 return @sequences;
161 }
162
163 sub get_indels_within_block{
164 my (@sequences) = @_;
165 my $line1 = my $line2 = my $line3 = "";
166 my @line1 = my @line2 = my @line3 = ();
167 my $score = 0;
168 my $start1 = my $align_length1 = my $end1 = my $seq_length1 = 0;
169 my $start2 = my $align_length2 = my $end2 = my $seq_length2 = 0;
170 my $start3 = my $align_length3 = my $end3 = my $seq_length3 = 0;
171 my $seq1 = my $orient1 = "";
172 my $seq2 = my $orient2 = "";
173 my $seq3 = my $orient3 = "";
174 my $start1_plus = my $end1_plus = 0;
175 my $start2_plus = my $end2_plus = 0;
176 my $start3_plus = my $end3_plus = 0;
177 my @test = ();
178 my $test = "";
179 my $header = "";
180 my @header = ();
181 my $sequence1 = my $sequence2 = my $sequence3 ="";
182 my @array_return = ();
183 my $test1 = 0;
184 my $line1_stat = my $line2_stat = my $line3_stat = "";
185
186 # process 3-way blocks only
187 if (scalar(@sequences) == 3){
188 $line1 = $sequences[0];
189 chomp ($line1);
190 $line2 = $sequences[1];
191 chomp ($line2);
192 $line3 = $sequences[2];
193 chomp ($line3);
194 # check order of sequences and assign uniformly seq1= human, seq2 = chimp, seq3 = macaque
195 if ($line1 =~ m/$outgroup/){
196 $line1_stat = "out";
197 $line2=~ s/^\s*//;
198 $line2 =~ s/\s+/\t/g;
199 @line2 = split(/\t/, $line2);
200 if (($ingroup1 eq "") || ($line2[1] =~ m/$ingroup1/)){
201 $line2_stat = "in1";
202 $line3_stat = "in2";
203 }
204 else{
205 $line3_stat = "in1";
206 $line2_stat = "in2"; }
207 }
208 elsif ($line2 =~ m/$outgroup/){
209 $line2_stat = "out";
210 $line1=~ s/^\s*//;
211 $line1 =~ s/\s+/\t/g;
212 @line1 = split(/\t/, $line1);
213 if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){
214 $line1_stat = "in1";
215 $line3_stat = "in2";
216 }
217 else{
218 $line3_stat = "in1";
219 $line1_stat = "in2"; }
220 }
221 elsif ($line3 =~ m/$outgroup/){
222 $line3_stat = "out";
223 $line1=~ s/^\s*//;
224 $line1 =~ s/\s+/\t/g;
225 @line1 = split(/\t/, $line1);
226 if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){
227 $line1_stat = "in1";
228 $line2_stat = "in2";
229 }
230 else{
231 $line2_stat = "in1";
232 $line1_stat = "in2"; }
233 }
234
235 #print "# l1 = $line1_stat\n";
236 #print "# l2 = $line2_stat\n";
237 #print "# l3 = $line3_stat\n";
238 my $linei1 = my $linei2 = my $lineo = "";
239 my @linei1 = my @linei2 = my @lineo = ();
240
241 if ($line1_stat eq "out"){
242 $lineo = $line1;
243 }
244 elsif ($line1_stat eq "in1"){
245 $linei1 = $line1;
246 }
247 else{
248 $linei2 = $line1;
249 }
250
251 if ($line2_stat eq "out"){
252 $lineo = $line2;
253 }
254 elsif ($line2_stat eq "in1"){
255 $linei1 = $line2;
256 }
257 else{
258 $linei2 = $line2;
259 }
260
261 if ($line3_stat eq "out"){
262 $lineo = $line3;
263 }
264 elsif ($line3_stat eq "in1"){
265 $linei1 = $line3;
266 }
267 else{
268 $linei2 = $line3;
269 }
270
271 $linei1=~ s/^\s*//;
272 $linei1 =~ s/\s+/\t/g;
273 @linei1 = split(/\t/, $linei1);
274 $end1 =($linei1[2]+$linei1[3]-1);
275 $seq1 = $linei1[1].":".$linei1[3];
276 $ingroup1 = (split(/\./, $seq1))[0];
277 $start1 = $linei1[2];
278 $align_length1 = $linei1[3];
279 $orient1 = $linei1[4];
280 $seq_length1 = $linei1[5];
281 $sequence1 = $linei1[6];
282 $test1 = length($sequence1);
283 my $total_length1 = $test1+$start1;
284 my @array1 = ($start1,$end1,$orient1,$seq_length1);
285 ($start1_plus, $end1_plus) = convert_coords(@array1);
286
287 $linei2=~ s/^\s*//;
288 $linei2 =~ s/\s+/\t/g;
289 @linei2 = split(/\t/, $linei2);
290 $end2 =($linei2[2]+$linei2[3]-1);
291 $seq2 = $linei2[1].":".$linei2[3];
292 $ingroup2 = (split(/\./, $seq2))[0];
293 $start2 = $linei2[2];
294 $align_length2 = $linei2[3];
295 $orient2 = $linei2[4];
296 $seq_length2 = $linei2[5];
297 $sequence2 = $linei2[6];
298 my $test2 = length($sequence2);
299 my $total_length2 = $test2+$start2;
300 my @array2 = ($start2,$end2,$orient2,$seq_length2);
301 ($start2_plus, $end2_plus) = convert_coords(@array2);
302
303 $lineo=~ s/^\s*//;
304 $lineo =~ s/\s+/\t/g;
305 @lineo = split(/\t/, $lineo);
306 $end3 =($lineo[2]+$lineo[3]-1);
307 $seq3 = $lineo[1].":".$lineo[3];
308 $start3 = $lineo[2];
309 $align_length3 = $lineo[3];
310 $orient3 = $lineo[4];
311 $seq_length3 = $lineo[5];
312 $sequence3 = $lineo[6];
313 my $test3 = length($sequence3);
314 my $total_length3 = $test3+$start3;
315 my @array3 = ($start3,$end3,$orient3,$seq_length3);
316 ($start3_plus, $end3_plus) = convert_coords(@array3);
317
318 #print "# l1 = $ingroup1\n";
319 #print "# l2 = $ingroup2\n";
320 #print "# l3 = $outgroup\n";
321
322 my $ABC = "";
323 my $coord1 = my $coord2 = my $coord3 = 0;
324 $coord1 = $start1_plus;
325 $coord2 = $start2_plus;
326 $coord3 = $start3_plus;
327
328 for (my $position = 0; $position < $test1; $position++) {
329 my $indelType = "";
330 my $indel_line = "";
331 # seq1 deletes
332 if ((substr($sequence1,$position,1) eq "-")
333 && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
334 && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
335 $ABC = join("",($ABC,"X"));
336 my @s = split(/:/, $seq1);
337 $indelType = $s[0]."_delete";
338
339 #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
340 $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
341 push (@indels,$indel_line);
342 push (@seq1_delete,$indel_line);
343 $coord2++; $coord3++;
344 }
345 # seq2 deletes
346 elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
347 && (substr($sequence2,$position,1) eq "-")
348 && (substr($sequence3,$position,1) !~ m/[-*\$?^]/)){
349 $ABC = join("",($ABC,"Y"));
350 my @s = split(/:/, $seq2);
351 $indelType = $s[0]."_delete";
352 #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
353 $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
354 push (@indels,$indel_line);
355 push (@seq2_delete,$indel_line);
356 $coord1++;
357 $coord3++;
358
359 }
360 # seq1 inserts
361 elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
362 && (substr($sequence2,$position,1) eq "-")
363 && (substr($sequence3,$position,1) eq "-")){
364 $ABC = join("",($ABC,"Z"));
365 my @s = split(/:/, $seq1);
366 $indelType = $s[0]."_insert";
367 #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
368 $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
369 push (@indels,$indel_line);
370 push (@seq1_insert,$indel_line);
371 $coord1++;
372 }
373 # seq2 inserts
374 elsif ((substr($sequence1,$position,1) eq "-")
375 && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
376 && (substr($sequence3,$position,1) eq "-")){
377 $ABC = join("",($ABC,"W"));
378 my @s = split(/:/, $seq2);
379 $indelType = $s[0]."_insert";
380 #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
381 $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
382 push (@indels,$indel_line);
383 push (@seq2_insert,$indel_line);
384 $coord2++;
385 }
386 # seq3 deletes
387 elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
388 && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
389 && (substr($sequence3,$position,1) eq "-")){
390 $ABC = join("",($ABC,"S"));
391 my @s = split(/:/, $seq3);
392 $indelType = $s[0]."_delete";
393 #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
394 $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
395 push (@indels,$indel_line);
396 push (@seq3_delete,$indel_line);
397 $coord1++; $coord2++;
398 }
399 # seq3 inserts
400 elsif ((substr($sequence1,$position,1) eq "-")
401 && (substr($sequence2,$position,1) eq "-")
402 && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
403 $ABC = join("",($ABC,"T"));
404 my @s = split(/:/, $seq3);
405 $indelType = $s[0]."_insert";
406 #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
407 $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
408 push (@indels,$indel_line);
409 push (@seq3_insert,$indel_line);
410 $coord3++;
411 }else{
412 $ABC = join("",($ABC,"N"));
413 $coord1++; $coord2++; $coord3++;
414 }
415
416 }
417 @array_return=($seq1,$seq2,$seq3,$ABC);
418 return (@array_return);
419
420 }
421 # ignore pairwise cases for now, just count the number of blocks
422 elsif (scalar(@sequences) == 2){
423 my $ABC = "";
424 my $coord1 = my $coord2 = my $coord3 = 0;
425 $count3++;
426
427 $line1 = $sequences[0];
428 $line2 = $sequences[1];
429 chomp ($line1);
430 chomp ($line2);
431
432 if ($line2 !~ m/$ingroup2/){
433 $count4++;
434 }
435 }
436 }
437
438
439 sub get_indels_lengths{
440 my (@array) = @_;
441 if (scalar(@array) == 4){
442 my $seq1 = $array[0]; my $seq2 = $array[1]; my $seq3 = $array[2]; my $ABC = $array[3];
443
444 while ($ABC =~ m/(X+)/g) {
445 push (@seq1_delete_lengths,length($1));
446 $count_seq1delete++;
447 }
448
449 while ($ABC =~ m/(Y+)/g) {
450 push (@seq2_delete_lengths,length($1));
451 $count_seq2delete++;
452 }
453 while ($ABC =~ m/(S+)/g) {
454 push (@seq3_delete_lengths,length($1));
455 $count_seq3delete++;
456 }
457 while ($ABC =~ m/(Z+)/g) {
458 push (@seq1_insert_lengths,length($1));
459 $count_seq1insert++;
460 }
461 while ($ABC =~ m/(W+)/g) {
462 push(@seq2_insert_lengths,length($1));
463 $count_seq2insert++;
464 }
465 while ($ABC =~ m/(T+)/g) {
466 push (@seq3_insert_lengths,length($1));
467 $count_seq3insert++;
468 }
469 }
470 elsif (scalar(@array) == 0){
471 next;
472 }
473
474 }
475 # convert to coordinates to + strand if align orient = -
476 sub convert_coords{
477 my (@array) = @_;
478 my $s = $array[0];
479 my $e = $array[1];
480 my $o = $array[2];
481 my $l = $array[3];
482 my $start_plus = my $end_plus = 0;
483
484 if ($o eq "-"){
485 $start_plus = ($l - $e);
486 $end_plus = ($l - $s);
487 }elsif ($o eq "+"){
488 $start_plus = $s;
489 $end_plus = $e;
490 }
491
492 return ($start_plus, $end_plus);
493 }
494
495 # get first line only for each event
496 sub get_starts_only{
497 my (@test) = @_;
498 my $seq1 = my $seq2 = my $seq3 = my $indelType = my $old_seq1 = my $old_seq2 = my $old_seq3 = my $old_indelType = my $old_line = "";
499 my $coord1 = my $coord2 = my $coord3 = my $old_coord1 = my $old_coord2 = my $old_coord3 = 0;
500
501 my @matches = ();
502 my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = ();
503
504
505 foreach my $line (@test){
506 chomp($line);
507 $line =~ s/^\s*//;
508 $line =~ s/\s+/\t/g;
509 my @line1 = split(/\t/, $line);
510 $seq1 = $line1[1];
511 $coord1 = $line1[2];
512 $seq2 = $line1[4];
513 $coord2 = $line1[5];
514 $seq3 = $line1[7];
515 $coord3 = $line1[8];
516 $indelType = $line1[10];
517 if ($indelType =~ m/$ingroup1/ && $indelType =~ m/insert/){
518 if ($coord1 != $old_coord1+1 || ($coord2 != $old_coord2 || $coord3 != $old_coord3)){
519 $start1++;
520 push (@seq1_insert_startOnly,$line);
521 }
522 }
523 elsif ($indelType =~ m/$ingroup1/ && $indelType =~ m/delete/){
524 if ($coord1 != $old_coord1 || ($coord2 != $old_coord2+1 || $coord3 != $old_coord3+1)){
525 $start2++;
526 push(@seq1_delete_startOnly,$line);
527 }
528 }
529 elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/insert/){
530 if ($coord2 != $old_coord2+1 || ($coord1 != $old_coord1 || $coord3 != $old_coord3)){
531 $start3++;
532 push(@seq2_insert_startOnly,$line);
533 }
534 }
535 elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/delete/){
536 if ($coord2 != $old_coord2 || ($coord1 != $old_coord1+1 || $coord3 != $old_coord3+1)){
537 $start4++;
538 push(@seq2_delete_startOnly,$line);
539 }
540 }
541 elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/insert/){
542 if ($coord3 != $old_coord3+1 || ($coord1 != $old_coord1 || $coord2 != $old_coord2)){
543 $start5++;
544 push(@seq3_insert_startOnly,$line);
545 }
546 }
547 elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/delete/){
548 if ($coord3 != $old_coord3 || ($coord1 != $old_coord1+1 ||$coord2 != $old_coord2+1)){
549 $start6++;
550 push(@seq3_delete_startOnly,$line);
551 }
552 }
553 $old_indelType = $indelType;
554 $old_seq1 = $seq1;
555 $old_coord1 = $coord1;
556 $old_seq2 = $seq2;
557 $old_coord2 = $coord2;
558 $old_seq3 = $seq3;
559 $old_coord3 = $coord3;
560 $old_line = $line;
561 }
562 }
563 # append lengths to each event start line
564 my $counter1; my $counter2; my $counter3; my $counter4; my $counter5; my $counter6;
565 my @final1 = my @final2 = my @final3 = my @final4 = my @final5 = my @final6 = ();
566 my $final_line1 = my $final_line2 = my $final_line3 = my $final_line4 = my $final_line5 = my $final_line6 = "";
567
568
569 for ($counter1 = 0; $counter1 < @seq3_insert_startOnly; $counter1++){
570 $final_line1 = join("\t",($seq3_insert_startOnly[$counter1],$seq3_insert_lengths[$counter1]));
571 push (@final1,$final_line1);
572 }
573
574 for ($counter2 = 0; $counter2 < @seq3_delete_startOnly; $counter2++){
575 $final_line2 = join("\t",($seq3_delete_startOnly[$counter2],$seq3_delete_lengths[$counter2]));
576 push(@final2,$final_line2);
577 }
578
579 for ($counter3 = 0; $counter3 < @seq2_insert_startOnly; $counter3++){
580 $final_line3 = join("\t",($seq2_insert_startOnly[$counter3],$seq2_insert_lengths[$counter3]));
581 push(@final3,$final_line3);
582 }
583
584 for ($counter4 = 0; $counter4 < @seq2_delete_startOnly; $counter4++){
585 $final_line4 = join("\t",($seq2_delete_startOnly[$counter4],$seq2_delete_lengths[$counter4]));
586 push(@final4,$final_line4);
587 }
588
589 for ($counter5 = 0; $counter5 < @seq1_insert_startOnly; $counter5++){
590 $final_line5 = join("\t",($seq1_insert_startOnly[$counter5],$seq1_insert_lengths[$counter5]));
591 push(@final5,$final_line5);
592 }
593
594 for ($counter6 = 0; $counter6 < @seq1_delete_startOnly; $counter6++){
595 $final_line6 = join("\t",($seq1_delete_startOnly[$counter6],$seq1_delete_lengths[$counter6]));
596 push(@final6,$final_line6);
597 }
598
599 # format final output
600 # # if inserts, increase coords for the sequence inserted, other sequences give coords for 5' and 3' bases flanking the gap
601 # # for deletes, increase coords for other 2 sequences and the one deleted give coords for 5' and 3' bases flanking the gap
602
603 get_final_format(@final5);
604 get_final_format(@final6);
605 get_final_format(@final3);
606 get_final_format(@final4);
607 get_final_format(@final1);
608 get_final_format(@final2);
609
610 sub get_final_format{
611 my (@final) = @_;
612 foreach (@final){
613 my $event_line = $_;
614 my @events = split(/\t/, $event_line);
615 my $event_type = $events[10];
616 my @name_align1 = split(/:/, $events[1]);
617 my @name_align2 = split(/:/, $events[4]);
618 my @name_align3 = split(/:/, $events[7]);
619 my $seq1_event_start = my $seq1_event_end = my $seq2_event_start = my $seq2_event_end = my $seq3_event_start = my $seq3_event_end = 0;
620 my $final_event_line = "";
621 # seq1_insert
622 if ($event_type =~ m/$ingroup1/ && $event_type =~ m/insert/){
623 # only increase coord for human
624 # remember that other two sequnences, the gap spans (coord - 1) --> coord
625 $seq1_event_start = ($events[2]);
626 $seq1_event_end = ($events[2]+$events[11]-1);
627 $seq2_event_start = ($events[5]-1);
628 $seq2_event_end = ($events[5]);
629 $seq3_event_start = ($events[8]-1);
630 $seq3_event_end = ($events[8]);
631 $final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
632 }
633 # seq1_delete
634 elsif ($event_type =~ m/$ingroup1/ && $event_type =~ m/delete/){
635 # only increase coords for seq2 and seq3
636 # remember for seq1, the gap spans (coord - 1) --> coord
637 $seq1_event_start = ($events[2]-1);
638 $seq1_event_end = ($events[2]);
639 $seq2_event_start = ($events[5]);
640 $seq2_event_end = ($events[5]+$events[11]-1);
641 $seq3_event_start = ($events[8]);
642 $seq3_event_end = ($events[8]+$events[11]-1);
643 $final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
644 }
645 # seq2_insert
646 elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/insert/){
647 # only increase coords for seq2
648 # remember that other two sequnences, the gap spans (coord - 1) --> coord
649 $seq1_event_start = ($events[2]-1);
650 $seq1_event_end = ($events[2]);
651 $seq2_event_start = ($events[5]);
652 $seq2_event_end = ($events[5]+$events[11]-1);
653 $seq3_event_start = ($events[8]-1);
654 $seq3_event_end = ($events[8]);
655 $final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
656 }
657 # seq2_delete
658 elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/delete/){
659 # only increase coords for seq1 and seq3
660 # remember for seq2, the gap spans (coord - 1) --> coord
661 $seq1_event_start = ($events[2]);
662 $seq1_event_end = ($events[2]+$events[11]-1);
663 $seq2_event_start = ($events[5]-1);
664 $seq2_event_end = ($events[5]);
665 $seq3_event_start = ($events[8]);
666 $seq3_event_end = ($events[8]+$events[11]-1);
667 $final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
668 }
669 # start testing w/seq3_insert
670 elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/insert/){
671 # only increase coord for rheMac
672 # remember that other two sequnences, the gap spans (coord - 1) --> coord
673 $seq1_event_start = ($events[2]-1);
674 $seq1_event_end = ($events[2]);
675 $seq2_event_start = ($events[5]-1);
676 $seq2_event_end = ($events[5]);
677 $seq3_event_start = ($events[8]);
678 $seq3_event_end = ($events[8]+$events[11]-1);
679 $final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
680 }
681 # seq3_delete
682 elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/delete/){
683 # only increase coords for seq1 and seq2
684 # remember for seq3, the gap spans (coord - 1) --> coord
685 $seq1_event_start = ($events[2]);
686 $seq1_event_end = ($events[2]+$events[11]-1);
687 $seq2_event_start = ($events[5]);
688 $seq2_event_end = ($events[5]+$events[11]-1);
689 $seq3_event_start = ($events[8]-1);
690 $seq3_event_end = ($events[8]);
691 $final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
692
693 }
694
695 print OFILE2 "$final_event_line\n";
696 }
697 }
698 close OFILE2;