PageRenderTime 118ms CodeModel.GetById 21ms app.highlight 90ms RepoModel.GetById 2ms app.codeStats 0ms

/tools/regVariation/parseMAF_smallIndels.pl

https://bitbucket.org/cistrome/cistrome-harvard/
Perl | 698 lines | 552 code | 73 blank | 73 comment | 101 complexity | 52a8c0d056cf5c5467079a7707586383 MD5 | raw file
  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 
  7use strict;
  8use warnings;
  9
 10# declare and initialize variables
 11my $fh; # variable to store filehandle
 12my $record;
 13my $offset;
 14my $library = $ARGV[0]; 
 15my $count = 0;
 16my $count2 = 0;
 17my $count3 = 0;
 18my $count4 = 0;
 19my $start1 = my $start2 = my $start3 = my $start4 = my $start5 = my $start6 = 0;
 20my $orient = "";
 21my $outgroup = $ARGV[2];
 22my $ingroup1 = my $ingroup2 = "";
 23my $count_seq1insert = my $count_seq1delete = 0;
 24my $count_seq2insert = my $count_seq2delete = 0;
 25my $count_seq3insert = my $count_seq3delete = 0;
 26my @seq1_insert_lengths = my @seq1_delete_lengths = ();
 27my @seq2_insert_lengths = my @seq2_delete_lengths = ();
 28my @seq3_insert_lengths = my @seq3_delete_lengths = ();
 29my @seq1_insert =  my @seq1_delete =  my @seq2_insert =  my @seq2_delete =  my @seq3_insert =  my @seq3_delete = ();
 30my @seq1_insert_startOnly = my @seq1_delete_startOnly = my @seq2_insert_startOnly = my @seq2_delete_startOnly = ();
 31my @seq3_insert_startOnly = my @seq3_delete_startOnly = ();
 32my @indels = (); 
 33
 34# check to make sure correct files
 35my $usage = "usage: parseMAF_smallIndels.pl [MAF.in] [small_Indels_summary.out] [outgroup]\n";
 36die $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
 49my $ofile2 = $ARGV[1];
 50unless (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";
 61print 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
 64while ($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
 78get_starts_only(@seq1_insert);
 79get_starts_only(@seq1_delete);
 80get_starts_only(@seq2_insert);
 81get_starts_only(@seq2_delete);
 82get_starts_only(@seq3_insert);
 83get_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";
 89if ($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
 94print "Ingroup1 = $ingroup1; Ingroup2 = $ingroup2; Outgroup = $outgroup;\n";
 95print "# of ingroup1 inserts = $count_seq1insert;\n";
 96print "# of ingroup1 deletes = $count_seq1delete;\n";
 97print "# of ingroup2 inserts = $count_seq2insert;\n";
 98print "# of ingroup2 deletes = $count_seq2delete;\n";
 99print "# of outgroup3 inserts = $count_seq3insert;\n";
100print "# of outgroup3 deletes = $count_seq3delete\n";
101
102
103#close OFILE;
104
105if ($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
113sub 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
125sub 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
140sub 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	
163sub 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		
439sub 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 = -
476sub 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
496sub 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
564my $counter1; my $counter2; my $counter3; my $counter4; my $counter5; my $counter6; 
565my @final1 = my @final2 = my @final3 = my @final4 = my @final5 = my @final6 = ();
566my $final_line1 = my $final_line2 = my $final_line3 = my $final_line4 = my $final_line5 = my $final_line6 = "";
567
568
569for ($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
574for ($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
579for ($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
584for ($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		
589for ($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
594for ($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
603get_final_format(@final5);
604get_final_format(@final6);
605get_final_format(@final3);
606get_final_format(@final4);
607get_final_format(@final1);
608get_final_format(@final2);
609
610sub 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}
698close OFILE2;