Mercurial > repos > xuebing > sharplabtool
view tools/regVariation/parseMAF_smallIndels.pl @ 0:9071e359b9a3
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:37:19 -0500 |
parents | |
children |
line wrap: on
line source
#!/usr/bin/perl -w # a program to get indels # input is a MAF format 3-way alignment file # from 3-way blocks only at this time # translate seq2, seq3, etc coordinates to + if align orient is reverse complement use strict; use warnings; # declare and initialize variables my $fh; # variable to store filehandle my $record; my $offset; my $library = $ARGV[0]; my $count = 0; my $count2 = 0; my $count3 = 0; my $count4 = 0; my $start1 = my $start2 = my $start3 = my $start4 = my $start5 = my $start6 = 0; my $orient = ""; my $outgroup = $ARGV[2]; my $ingroup1 = my $ingroup2 = ""; my $count_seq1insert = my $count_seq1delete = 0; my $count_seq2insert = my $count_seq2delete = 0; my $count_seq3insert = my $count_seq3delete = 0; my @seq1_insert_lengths = my @seq1_delete_lengths = (); my @seq2_insert_lengths = my @seq2_delete_lengths = (); my @seq3_insert_lengths = my @seq3_delete_lengths = (); my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = (); my @seq1_insert_startOnly = my @seq1_delete_startOnly = my @seq2_insert_startOnly = my @seq2_delete_startOnly = (); my @seq3_insert_startOnly = my @seq3_delete_startOnly = (); my @indels = (); # check to make sure correct files my $usage = "usage: parseMAF_smallIndels.pl [MAF.in] [small_Indels_summary.out] [outgroup]\n"; die $usage unless @ARGV == 3; # perform some standard subroutines $fh = open_file($library); $offset = tell($fh); #my $ofile = $ARGV[2]; #unless (open(OFILE, ">$ofile")){ # print "Cannot open output file \"$ofile\"\n\n"; # exit; #} my $ofile2 = $ARGV[1]; unless (open(OFILE2, ">$ofile2")){ print "Cannot open output file \"$ofile2\"\n\n"; exit; } # header line for output files #print OFILE "# small indel events, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n"; #print OFILE "#align\tingroup1\tingroup1_coord\tingroup1_orient\tingroup2\tingroup2_coord\tingroup2_orient\toutgroup\toutgroup_coord\toutgroup_orient\tindel_type\n"; #print OFILE2 "# small indels summary, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n"; 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"; # main body of program while ($record = get_next_record($fh) ){ if ($record=~ m/\s*##maf(.*)\s*# maf/s){ next; } my @sequences = get_sequences_within_block($record); my @seq_info = get_indels_within_block(@sequences); get_indels_lengths(@seq_info); $offset = tell($fh); $count++; } get_starts_only(@seq1_insert); get_starts_only(@seq1_delete); get_starts_only(@seq2_insert); get_starts_only(@seq2_delete); get_starts_only(@seq3_insert); get_starts_only(@seq3_delete); # print some things to keep track of progress #print "# $library\n"; #print "# number of records = $count\n"; #print "# number of sequence \"s\" lines = $count2\n"; if ($count3 != 0){ print "Skipped $count3 blocks with only 2 seqs;\n"; } #print "# number of records with only h-m = $count4\n\n"; print "Ingroup1 = $ingroup1; Ingroup2 = $ingroup2; Outgroup = $outgroup;\n"; print "# of ingroup1 inserts = $count_seq1insert;\n"; print "# of ingroup1 deletes = $count_seq1delete;\n"; print "# of ingroup2 inserts = $count_seq2insert;\n"; print "# of ingroup2 deletes = $count_seq2delete;\n"; print "# of outgroup3 inserts = $count_seq3insert;\n"; print "# of outgroup3 deletes = $count_seq3delete\n"; #close OFILE; if ($count == $count3){ print STDERR "Skipped all blocks since none of them contain 3-way alignments.\n"; exit -1; } ###################SUBROUTINES##################################### # subroutine to open file sub open_file { my($filename) = @_; my $fh; unless (open($fh, $filename)){ print "Cannot open file $filename\n"; exit; } return $fh; } # get next record sub get_next_record { my ($fh) = @_; my ($offset); my ($record) = ""; my ($save_input_separator) = $/; $/ = "a score"; $record = <$fh>; $/ = $save_input_separator; return $record; } # get header info sub get_sequences_within_block{ my (@alignment) = @_; my @lines = (); my @sequences = (); @lines = split ("\n", $record); foreach (@lines){ chomp($_); if (m/^\s*$/){ next; } elsif (m/^\s*=(\d+\.*\d*)/){ }elsif (m/^\s*s(.*)$/){ $count2++; push (@sequences,$_); } } return @sequences; } sub get_indels_within_block{ my (@sequences) = @_; my $line1 = my $line2 = my $line3 = ""; my @line1 = my @line2 = my @line3 = (); my $score = 0; my $start1 = my $align_length1 = my $end1 = my $seq_length1 = 0; my $start2 = my $align_length2 = my $end2 = my $seq_length2 = 0; my $start3 = my $align_length3 = my $end3 = my $seq_length3 = 0; my $seq1 = my $orient1 = ""; my $seq2 = my $orient2 = ""; my $seq3 = my $orient3 = ""; my $start1_plus = my $end1_plus = 0; my $start2_plus = my $end2_plus = 0; my $start3_plus = my $end3_plus = 0; my @test = (); my $test = ""; my $header = ""; my @header = (); my $sequence1 = my $sequence2 = my $sequence3 =""; my @array_return = (); my $test1 = 0; my $line1_stat = my $line2_stat = my $line3_stat = ""; # process 3-way blocks only if (scalar(@sequences) == 3){ $line1 = $sequences[0]; chomp ($line1); $line2 = $sequences[1]; chomp ($line2); $line3 = $sequences[2]; chomp ($line3); # check order of sequences and assign uniformly seq1= human, seq2 = chimp, seq3 = macaque if ($line1 =~ m/$outgroup/){ $line1_stat = "out"; $line2=~ s/^\s*//; $line2 =~ s/\s+/\t/g; @line2 = split(/\t/, $line2); if (($ingroup1 eq "") || ($line2[1] =~ m/$ingroup1/)){ $line2_stat = "in1"; $line3_stat = "in2"; } else{ $line3_stat = "in1"; $line2_stat = "in2"; } } elsif ($line2 =~ m/$outgroup/){ $line2_stat = "out"; $line1=~ s/^\s*//; $line1 =~ s/\s+/\t/g; @line1 = split(/\t/, $line1); if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){ $line1_stat = "in1"; $line3_stat = "in2"; } else{ $line3_stat = "in1"; $line1_stat = "in2"; } } elsif ($line3 =~ m/$outgroup/){ $line3_stat = "out"; $line1=~ s/^\s*//; $line1 =~ s/\s+/\t/g; @line1 = split(/\t/, $line1); if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){ $line1_stat = "in1"; $line2_stat = "in2"; } else{ $line2_stat = "in1"; $line1_stat = "in2"; } } #print "# l1 = $line1_stat\n"; #print "# l2 = $line2_stat\n"; #print "# l3 = $line3_stat\n"; my $linei1 = my $linei2 = my $lineo = ""; my @linei1 = my @linei2 = my @lineo = (); if ($line1_stat eq "out"){ $lineo = $line1; } elsif ($line1_stat eq "in1"){ $linei1 = $line1; } else{ $linei2 = $line1; } if ($line2_stat eq "out"){ $lineo = $line2; } elsif ($line2_stat eq "in1"){ $linei1 = $line2; } else{ $linei2 = $line2; } if ($line3_stat eq "out"){ $lineo = $line3; } elsif ($line3_stat eq "in1"){ $linei1 = $line3; } else{ $linei2 = $line3; } $linei1=~ s/^\s*//; $linei1 =~ s/\s+/\t/g; @linei1 = split(/\t/, $linei1); $end1 =($linei1[2]+$linei1[3]-1); $seq1 = $linei1[1].":".$linei1[3]; $ingroup1 = (split(/\./, $seq1))[0]; $start1 = $linei1[2]; $align_length1 = $linei1[3]; $orient1 = $linei1[4]; $seq_length1 = $linei1[5]; $sequence1 = $linei1[6]; $test1 = length($sequence1); my $total_length1 = $test1+$start1; my @array1 = ($start1,$end1,$orient1,$seq_length1); ($start1_plus, $end1_plus) = convert_coords(@array1); $linei2=~ s/^\s*//; $linei2 =~ s/\s+/\t/g; @linei2 = split(/\t/, $linei2); $end2 =($linei2[2]+$linei2[3]-1); $seq2 = $linei2[1].":".$linei2[3]; $ingroup2 = (split(/\./, $seq2))[0]; $start2 = $linei2[2]; $align_length2 = $linei2[3]; $orient2 = $linei2[4]; $seq_length2 = $linei2[5]; $sequence2 = $linei2[6]; my $test2 = length($sequence2); my $total_length2 = $test2+$start2; my @array2 = ($start2,$end2,$orient2,$seq_length2); ($start2_plus, $end2_plus) = convert_coords(@array2); $lineo=~ s/^\s*//; $lineo =~ s/\s+/\t/g; @lineo = split(/\t/, $lineo); $end3 =($lineo[2]+$lineo[3]-1); $seq3 = $lineo[1].":".$lineo[3]; $start3 = $lineo[2]; $align_length3 = $lineo[3]; $orient3 = $lineo[4]; $seq_length3 = $lineo[5]; $sequence3 = $lineo[6]; my $test3 = length($sequence3); my $total_length3 = $test3+$start3; my @array3 = ($start3,$end3,$orient3,$seq_length3); ($start3_plus, $end3_plus) = convert_coords(@array3); #print "# l1 = $ingroup1\n"; #print "# l2 = $ingroup2\n"; #print "# l3 = $outgroup\n"; my $ABC = ""; my $coord1 = my $coord2 = my $coord3 = 0; $coord1 = $start1_plus; $coord2 = $start2_plus; $coord3 = $start3_plus; for (my $position = 0; $position < $test1; $position++) { my $indelType = ""; my $indel_line = ""; # seq1 deletes if ((substr($sequence1,$position,1) eq "-") && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/) && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){ $ABC = join("",($ABC,"X")); my @s = split(/:/, $seq1); $indelType = $s[0]."_delete"; #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n"; $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType)); push (@indels,$indel_line); push (@seq1_delete,$indel_line); $coord2++; $coord3++; } # seq2 deletes elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/) && (substr($sequence2,$position,1) eq "-") && (substr($sequence3,$position,1) !~ m/[-*\$?^]/)){ $ABC = join("",($ABC,"Y")); my @s = split(/:/, $seq2); $indelType = $s[0]."_delete"; #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n"; $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType)); push (@indels,$indel_line); push (@seq2_delete,$indel_line); $coord1++; $coord3++; } # seq1 inserts elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/) && (substr($sequence2,$position,1) eq "-") && (substr($sequence3,$position,1) eq "-")){ $ABC = join("",($ABC,"Z")); my @s = split(/:/, $seq1); $indelType = $s[0]."_insert"; #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n"; $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType)); push (@indels,$indel_line); push (@seq1_insert,$indel_line); $coord1++; } # seq2 inserts elsif ((substr($sequence1,$position,1) eq "-") && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/) && (substr($sequence3,$position,1) eq "-")){ $ABC = join("",($ABC,"W")); my @s = split(/:/, $seq2); $indelType = $s[0]."_insert"; #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n"; $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType)); push (@indels,$indel_line); push (@seq2_insert,$indel_line); $coord2++; } # seq3 deletes elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/) && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/) && (substr($sequence3,$position,1) eq "-")){ $ABC = join("",($ABC,"S")); my @s = split(/:/, $seq3); $indelType = $s[0]."_delete"; #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n"; $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType)); push (@indels,$indel_line); push (@seq3_delete,$indel_line); $coord1++; $coord2++; } # seq3 inserts elsif ((substr($sequence1,$position,1) eq "-") && (substr($sequence2,$position,1) eq "-") && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){ $ABC = join("",($ABC,"T")); my @s = split(/:/, $seq3); $indelType = $s[0]."_insert"; #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n"; $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType)); push (@indels,$indel_line); push (@seq3_insert,$indel_line); $coord3++; }else{ $ABC = join("",($ABC,"N")); $coord1++; $coord2++; $coord3++; } } @array_return=($seq1,$seq2,$seq3,$ABC); return (@array_return); } # ignore pairwise cases for now, just count the number of blocks elsif (scalar(@sequences) == 2){ my $ABC = ""; my $coord1 = my $coord2 = my $coord3 = 0; $count3++; $line1 = $sequences[0]; $line2 = $sequences[1]; chomp ($line1); chomp ($line2); if ($line2 !~ m/$ingroup2/){ $count4++; } } } sub get_indels_lengths{ my (@array) = @_; if (scalar(@array) == 4){ my $seq1 = $array[0]; my $seq2 = $array[1]; my $seq3 = $array[2]; my $ABC = $array[3]; while ($ABC =~ m/(X+)/g) { push (@seq1_delete_lengths,length($1)); $count_seq1delete++; } while ($ABC =~ m/(Y+)/g) { push (@seq2_delete_lengths,length($1)); $count_seq2delete++; } while ($ABC =~ m/(S+)/g) { push (@seq3_delete_lengths,length($1)); $count_seq3delete++; } while ($ABC =~ m/(Z+)/g) { push (@seq1_insert_lengths,length($1)); $count_seq1insert++; } while ($ABC =~ m/(W+)/g) { push(@seq2_insert_lengths,length($1)); $count_seq2insert++; } while ($ABC =~ m/(T+)/g) { push (@seq3_insert_lengths,length($1)); $count_seq3insert++; } } elsif (scalar(@array) == 0){ next; } } # convert to coordinates to + strand if align orient = - sub convert_coords{ my (@array) = @_; my $s = $array[0]; my $e = $array[1]; my $o = $array[2]; my $l = $array[3]; my $start_plus = my $end_plus = 0; if ($o eq "-"){ $start_plus = ($l - $e); $end_plus = ($l - $s); }elsif ($o eq "+"){ $start_plus = $s; $end_plus = $e; } return ($start_plus, $end_plus); } # get first line only for each event sub get_starts_only{ my (@test) = @_; my $seq1 = my $seq2 = my $seq3 = my $indelType = my $old_seq1 = my $old_seq2 = my $old_seq3 = my $old_indelType = my $old_line = ""; my $coord1 = my $coord2 = my $coord3 = my $old_coord1 = my $old_coord2 = my $old_coord3 = 0; my @matches = (); my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = (); foreach my $line (@test){ chomp($line); $line =~ s/^\s*//; $line =~ s/\s+/\t/g; my @line1 = split(/\t/, $line); $seq1 = $line1[1]; $coord1 = $line1[2]; $seq2 = $line1[4]; $coord2 = $line1[5]; $seq3 = $line1[7]; $coord3 = $line1[8]; $indelType = $line1[10]; if ($indelType =~ m/$ingroup1/ && $indelType =~ m/insert/){ if ($coord1 != $old_coord1+1 || ($coord2 != $old_coord2 || $coord3 != $old_coord3)){ $start1++; push (@seq1_insert_startOnly,$line); } } elsif ($indelType =~ m/$ingroup1/ && $indelType =~ m/delete/){ if ($coord1 != $old_coord1 || ($coord2 != $old_coord2+1 || $coord3 != $old_coord3+1)){ $start2++; push(@seq1_delete_startOnly,$line); } } elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/insert/){ if ($coord2 != $old_coord2+1 || ($coord1 != $old_coord1 || $coord3 != $old_coord3)){ $start3++; push(@seq2_insert_startOnly,$line); } } elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/delete/){ if ($coord2 != $old_coord2 || ($coord1 != $old_coord1+1 || $coord3 != $old_coord3+1)){ $start4++; push(@seq2_delete_startOnly,$line); } } elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/insert/){ if ($coord3 != $old_coord3+1 || ($coord1 != $old_coord1 || $coord2 != $old_coord2)){ $start5++; push(@seq3_insert_startOnly,$line); } } elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/delete/){ if ($coord3 != $old_coord3 || ($coord1 != $old_coord1+1 ||$coord2 != $old_coord2+1)){ $start6++; push(@seq3_delete_startOnly,$line); } } $old_indelType = $indelType; $old_seq1 = $seq1; $old_coord1 = $coord1; $old_seq2 = $seq2; $old_coord2 = $coord2; $old_seq3 = $seq3; $old_coord3 = $coord3; $old_line = $line; } } # append lengths to each event start line my $counter1; my $counter2; my $counter3; my $counter4; my $counter5; my $counter6; my @final1 = my @final2 = my @final3 = my @final4 = my @final5 = my @final6 = (); my $final_line1 = my $final_line2 = my $final_line3 = my $final_line4 = my $final_line5 = my $final_line6 = ""; for ($counter1 = 0; $counter1 < @seq3_insert_startOnly; $counter1++){ $final_line1 = join("\t",($seq3_insert_startOnly[$counter1],$seq3_insert_lengths[$counter1])); push (@final1,$final_line1); } for ($counter2 = 0; $counter2 < @seq3_delete_startOnly; $counter2++){ $final_line2 = join("\t",($seq3_delete_startOnly[$counter2],$seq3_delete_lengths[$counter2])); push(@final2,$final_line2); } for ($counter3 = 0; $counter3 < @seq2_insert_startOnly; $counter3++){ $final_line3 = join("\t",($seq2_insert_startOnly[$counter3],$seq2_insert_lengths[$counter3])); push(@final3,$final_line3); } for ($counter4 = 0; $counter4 < @seq2_delete_startOnly; $counter4++){ $final_line4 = join("\t",($seq2_delete_startOnly[$counter4],$seq2_delete_lengths[$counter4])); push(@final4,$final_line4); } for ($counter5 = 0; $counter5 < @seq1_insert_startOnly; $counter5++){ $final_line5 = join("\t",($seq1_insert_startOnly[$counter5],$seq1_insert_lengths[$counter5])); push(@final5,$final_line5); } for ($counter6 = 0; $counter6 < @seq1_delete_startOnly; $counter6++){ $final_line6 = join("\t",($seq1_delete_startOnly[$counter6],$seq1_delete_lengths[$counter6])); push(@final6,$final_line6); } # format final output # # if inserts, increase coords for the sequence inserted, other sequences give coords for 5' and 3' bases flanking the gap # # for deletes, increase coords for other 2 sequences and the one deleted give coords for 5' and 3' bases flanking the gap get_final_format(@final5); get_final_format(@final6); get_final_format(@final3); get_final_format(@final4); get_final_format(@final1); get_final_format(@final2); sub get_final_format{ my (@final) = @_; foreach (@final){ my $event_line = $_; my @events = split(/\t/, $event_line); my $event_type = $events[10]; my @name_align1 = split(/:/, $events[1]); my @name_align2 = split(/:/, $events[4]); my @name_align3 = split(/:/, $events[7]); 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; my $final_event_line = ""; # seq1_insert if ($event_type =~ m/$ingroup1/ && $event_type =~ m/insert/){ # only increase coord for human # remember that other two sequnences, the gap spans (coord - 1) --> coord $seq1_event_start = ($events[2]); $seq1_event_end = ($events[2]+$events[11]-1); $seq2_event_start = ($events[5]-1); $seq2_event_end = ($events[5]); $seq3_event_start = ($events[8]-1); $seq3_event_end = ($events[8]); $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])); } # seq1_delete elsif ($event_type =~ m/$ingroup1/ && $event_type =~ m/delete/){ # only increase coords for seq2 and seq3 # remember for seq1, the gap spans (coord - 1) --> coord $seq1_event_start = ($events[2]-1); $seq1_event_end = ($events[2]); $seq2_event_start = ($events[5]); $seq2_event_end = ($events[5]+$events[11]-1); $seq3_event_start = ($events[8]); $seq3_event_end = ($events[8]+$events[11]-1); $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])); } # seq2_insert elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/insert/){ # only increase coords for seq2 # remember that other two sequnences, the gap spans (coord - 1) --> coord $seq1_event_start = ($events[2]-1); $seq1_event_end = ($events[2]); $seq2_event_start = ($events[5]); $seq2_event_end = ($events[5]+$events[11]-1); $seq3_event_start = ($events[8]-1); $seq3_event_end = ($events[8]); $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])); } # seq2_delete elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/delete/){ # only increase coords for seq1 and seq3 # remember for seq2, the gap spans (coord - 1) --> coord $seq1_event_start = ($events[2]); $seq1_event_end = ($events[2]+$events[11]-1); $seq2_event_start = ($events[5]-1); $seq2_event_end = ($events[5]); $seq3_event_start = ($events[8]); $seq3_event_end = ($events[8]+$events[11]-1); $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])); } # start testing w/seq3_insert elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/insert/){ # only increase coord for rheMac # remember that other two sequnences, the gap spans (coord - 1) --> coord $seq1_event_start = ($events[2]-1); $seq1_event_end = ($events[2]); $seq2_event_start = ($events[5]-1); $seq2_event_end = ($events[5]); $seq3_event_start = ($events[8]); $seq3_event_end = ($events[8]+$events[11]-1); $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])); } # seq3_delete elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/delete/){ # only increase coords for seq1 and seq2 # remember for seq3, the gap spans (coord - 1) --> coord $seq1_event_start = ($events[2]); $seq1_event_end = ($events[2]+$events[11]-1); $seq2_event_start = ($events[5]); $seq2_event_end = ($events[5]+$events[11]-1); $seq3_event_start = ($events[8]-1); $seq3_event_end = ($events[8]); $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])); } print OFILE2 "$final_event_line\n"; } } close OFILE2;