root/galaxy-central/tools/regVariation/parseMAF_smallIndels.pl

リビジョン 2, 46.1 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
1#!/usr/bin/perl -w
2# a program to get indels
3# the input is a MAF format 3-way alignment file from 3-way blocks only
4# the output is a TABULAR format file containing the coordinates of indels determined based on the outgroup
5# the program translates seq2, seq3, etc coordinates to (5'-3') strand orientation if alignment orientation
6# is reverse complement
7# the program gets only indels that satisfy the following conditions:
8# 1. two gaps can not overlap unless they are identical, which means they start and end at the same coordinates
9# 2. the indel must start after at least three bases from the beginning of the block and finish before at least three bases from the end of the block
10# 3. a low score base (masked) can not exist within the indel interval in each of the three sequences
11# 4. at lease three high score bases (ummasked) must exist on each side of the indel in each of the three sequences
12
13
14use strict;
15use warnings;
16
17# declare and initialize variables
18my $fh; # variable to store filehandle
19my $record;
20my $offset;
21my $library = $ARGV[0];
22my $count = 0;
23my $count2 = 0;
24my $count3 = 0;
25my $count4 = 0;
26my $start1 = my $start2 = my $start3 = my $start4 = my $start5 = my $start6 = 0;
27my $orient = "";
28my $outgroup = $ARGV[2];
29my $ingroup1 = my $ingroup2 = "";
30my $count_seq1insert = my $count_seq1delete = 0;
31my $count_seq2insert = my $count_seq2delete = 0;
32my $count_seq3insert = my $count_seq3delete = 0;
33my @seq1_insert_lengths = my @seq1_delete_lengths = ();
34my @seq2_insert_lengths = my @seq2_delete_lengths = ();
35my @seq3_insert_lengths = my @seq3_delete_lengths = ();
36my @seq1_insert =  my @seq1_delete =  my @seq2_insert =  my @seq2_delete =  my @seq3_insert =  my @seq3_delete = ();
37my @seq1_insert_startOnly = my @seq1_delete_startOnly = my @seq2_insert_startOnly = my @seq2_delete_startOnly = ();
38my @seq3_insert_startOnly = my @seq3_delete_startOnly = ();
39my @indels = ();
40
41# check to make sure correct files
42my $usage = "usage: parseMAF_smallIndels.pl [MAF.in] [small_Indels_summary_TABULAR.out] [outgroup]\n";
43die $usage unless @ARGV == 3;
44
45# perform some standard subroutines
46$fh = open_file($library);
47
48$offset = tell($fh);
49
50my $ofile = $ARGV[1];
51unless (open(OFILE, ">$ofile")){
52        print "Cannot open output file \"$ofile\"\n\n";
53    exit;
54}
55
56
57# header line for output files
58#print OFILE "# small indels summary, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n";
59print OFILE "#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";
60
61# main body of program
62while ($record = get_next_record($fh) ){
63        if ($record=~ m/\s*##maf(.*)\s*# maf/s){
64                next;
65        }
66   
67        my @sequences = get_sequences_within_block($record);
68        my @seq_info = get_indels_within_block(@sequences);
69       
70        get_indels_lengths(@seq_info); 
71        $offset = tell($fh);
72    $count++;
73       
74}
75
76get_starts_only(@seq1_insert);
77get_starts_only(@seq1_delete);
78get_starts_only(@seq2_insert);
79get_starts_only(@seq2_delete);
80get_starts_only(@seq3_insert);
81get_starts_only(@seq3_delete);
82
83# print some things to keep track of progress
84#print "# $library\n";
85#print "# number of records = $count\n";
86#print "# number of sequence \"s\" lines = $count2\n";
87
88if ($count3 != 0){
89        print "Skipped $count3 blocks with only 2 seqs;\n";
90}
91#print "# number of records with only h-m = $count4\n\n";
92
93print "Ingroup1 = $ingroup1; Ingroup2 = $ingroup2; Outgroup = $outgroup;\n";
94print "# of ingroup1 inserts = $count_seq1insert;\n";
95print "# of ingroup1 deletes = $count_seq1delete;\n";
96print "# of ingroup2 inserts = $count_seq2insert;\n";
97print "# of ingroup2 deletes = $count_seq2delete;\n";
98print "# of outgroup3 inserts = $count_seq3insert;\n";
99print "# of outgroup3 deletes = $count_seq3delete\n";
100
101
102#close OFILE;
103
104if ($count == $count3){
105        print STDERR "Skipped all blocks since none of them contain 3-way alignments.\n";
106        exit -1;
107}
108
109###################SUBROUTINES#####################################
110
111# subroutine to open file
112sub open_file {
113        my($filename) = @_;
114        my $fh;
115
116        unless (open($fh, $filename)){
117                print "Cannot open file $filename\n";
118                exit;
119        }
120        return $fh;
121}
122
123# get next record
124sub get_next_record {
125        my ($fh) = @_;
126        my ($offset);
127        my ($record) = "";
128        my ($save_input_separator) = $/;
129       
130            $/ = "a score";
131
132        $record = <$fh>;
133
134        $/ = $save_input_separator;
135       
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       
161        return @sequences;
162}
163
164#get the indels within the blocks       
165sub get_indels_within_block{
166        my (@sequences) = @_;   
167        my $line1 = my $line2 = my $line3 = "";
168    my @line1 = my @line2 = my @line3 = ();
169        my $score = 0;
170        my $start1 = my $align_length1 = my $end1 = my $seq_length1 = 0;
171        my $start2 = my $align_length2 = my $end2 = my $seq_length2 = 0;
172        my $start3 = my $align_length3 = my $end3 = my $seq_length3 = 0;
173        my $seq1 = my $orient1 = "";
174        my $seq2 = my $orient2 = "";
175        my $seq3 = my $orient3 = "";
176        my $start1_plus = my $end1_plus = 0;
177        my $start2_plus = my $end2_plus = 0;
178        my $start3_plus = my $end3_plus = 0;
179        my @test = ();
180        my $test = "";
181        my $header = "";
182        my @header = ();
183        my $sequence1 = my $sequence2 = my $sequence3 ="";
184        my @array_return = (); 
185        my $test1 = 0;
186        my $line1_stat = my $line2_stat = my $line3_stat = "";
187   
188        # process 3-way blocks only
189        if (scalar(@sequences) == 3){
190                $line1 = $sequences[0];
191                chomp ($line1);
192                $line2 = $sequences[1];
193                chomp ($line2);
194                $line3 = $sequences[2];
195                chomp ($line3);
196               
197                # check order of sequences and assign uniformly seq1= human, seq2 = chimp, seq3 = macaque
198                if ($line1 =~ m/$outgroup/){
199                        $line1_stat = "out";
200                        $line2=~ s/^\s*//;
201            $line2 =~ s/\s+/\t/g;
202            @line2 = split(/\t/, $line2);
203                        if (($ingroup1 eq "") || ($line2[1] =~ m/$ingroup1/)){
204                                $line2_stat = "in1";
205                                $line3_stat = "in2";
206                         }
207                         else{
208                $line3_stat = "in1";
209                $line2_stat = "in2";             
210                         }
211                        }
212                elsif ($line2 =~ m/$outgroup/){
213                        $line2_stat = "out";
214                        $line1=~ s/^\s*//;
215            $line1 =~ s/\s+/\t/g;
216            @line1 = split(/\t/, $line1);
217            if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){
218                $line1_stat = "in1";
219                $line3_stat = "in2";
220             }
221             else{
222                $line3_stat = "in1";
223                $line1_stat = "in2";             
224             }
225            }
226                elsif ($line3 =~ m/$outgroup/){
227                        $line3_stat = "out";
228                        $line1=~ s/^\s*//;
229            $line1 =~ s/\s+/\t/g;
230            @line1 = split(/\t/, $line1);
231            if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){
232                $line1_stat = "in1";
233                $line2_stat = "in2";
234             }
235             else{
236                $line2_stat = "in1";
237                $line1_stat = "in2";             
238             }
239                }
240
241                #print "# l1 = $line1_stat\n";
242                #print "# l2 = $line2_stat\n";
243                #print "# l3 = $line3_stat\n";
244               
245                my $linei1 = my $linei2 = my $lineo = "";
246                my @linei1 = my @linei2 = my @lineo = ();
247               
248        if ($line1_stat eq "out"){
249            $lineo = $line1;
250        }
251        elsif ($line1_stat eq "in1"){
252            $linei1 = $line1;
253        }
254        else{
255            $linei2 = $line1;
256        }
257       
258        if ($line2_stat eq "out"){
259            $lineo = $line2;
260        }
261        elsif ($line2_stat eq "in1"){
262            $linei1 = $line2;
263        }
264        else{
265            $linei2 = $line2;
266        }
267       
268        if ($line3_stat eq "out"){
269            $lineo = $line3;
270        }
271        elsif ($line3_stat eq "in1"){
272            $linei1 = $line3;
273        }
274        else{
275            $linei2 = $line3;
276        }
277       
278        $linei1=~ s/^\s*//;
279        $linei1 =~ s/\s+/\t/g;
280        @linei1 = split(/\t/, $linei1);
281        $end1 =($linei1[2]+$linei1[3]-1);
282        $seq1 = $linei1[1].":".$linei1[3];
283        $ingroup1 = (split(/\./, $seq1))[0];
284        $start1 = $linei1[2];
285        $align_length1 = $linei1[3];
286        $orient1 = $linei1[4];
287        $seq_length1 = $linei1[5];
288        $sequence1 = $linei1[6];
289        $test1 = length($sequence1);
290        my $total_length1 = $test1+$start1;
291        my @array1 = ($start1,$end1,$orient1,$seq_length1);
292        ($start1_plus, $end1_plus) =  convert_coords(@array1);
293           
294        $linei2=~ s/^\s*//;
295        $linei2 =~ s/\s+/\t/g;
296        @linei2 = split(/\t/, $linei2);
297        $end2 =($linei2[2]+$linei2[3]-1);               
298        $seq2 = $linei2[1].":".$linei2[3];
299        $ingroup2 = (split(/\./, $seq2))[0];
300        $start2 = $linei2[2];
301        $align_length2 = $linei2[3];
302        $orient2 = $linei2[4];
303        $seq_length2 = $linei2[5];
304        $sequence2 = $linei2[6];
305        my $test2 = length($sequence2);
306        my $total_length2 = $test2+$start2;
307        my @array2 = ($start2,$end2,$orient2,$seq_length2);
308        ($start2_plus, $end2_plus) = convert_coords(@array2);
309           
310        $lineo=~ s/^\s*//;
311        $lineo =~ s/\s+/\t/g;
312        @lineo = split(/\t/, $lineo);
313        $end3 =($lineo[2]+$lineo[3]-1);
314        $seq3 = $lineo[1].":".$lineo[3];
315        $start3 = $lineo[2];
316        $align_length3 = $lineo[3];
317        $orient3 = $lineo[4];
318        $seq_length3 = $lineo[5];
319        $sequence3 = $lineo[6];
320        my $test3 = length($sequence3);
321        my $total_length3 = $test3+$start3;
322        my @array3 = ($start3,$end3,$orient3,$seq_length3);
323        ($start3_plus, $end3_plus) = convert_coords(@array3);
324       
325        #print "# l1 = $ingroup1\n";
326                #print "# l2 = $ingroup2\n";
327                #print "# l3 = $outgroup\n";
328       
329                my $ABC = "";
330                my $coord1 = my $coord2 = my $coord3 = 0;
331                $coord1 = $start1_plus;
332                $coord2 = $start2_plus;
333                $coord3 = $start3_plus;
334               
335                #variables to store the lengths of gaps in each sequence
336                my $gapLenth = my $gapLenth1 = my $gapLenth2 = my $gapLenth3 = 0;
337                               
338                #variables to store the index of the start of a gap in each sequence
339                my $gapStartPos = my $gapStartPos1 = my $gapStartPos2 = my $gapStartPos3 = 0;   
340                my $gapEndPos = 0;
341               
342                #variable to determine if there is a valid gap
343                my $getValidGap = 0;
344               
345                #scan the block from left to right looking for gaps
346                for (my $position = 0; $position < $test1; $position++) {
347                        my $indelType = "";
348                        my $indel_line = "";
349       
350                        #print($position ." position\n");
351                        #print($sequence1 ." sequence1\n");
352                        #print($sequence2 ." sequence2\n");
353                        #print($sequence3 ." sequence3\n");
354                        #if ($position > 2 && $position < ($test1 - 3)){
355                                #print(substr($sequence1,$position - 3, 3) . " before " . substr($sequence1, $position + 1, 3) . " after\n");
356                                #print(substr($sequence2,$position - 3, 3) . " before " . substr($sequence2, $position + 1, 3) . " after\n");
357                                #print(substr($sequence3,$position - 3, 3) . " before " . substr($sequence3, $position + 1, 3) . " after\n");
358                        #}
359                                       
360                        #if the current position is the start of a gap, look for any overlapping gaps in the other two sequences
361                        if (substr($sequence1,$position,1) eq "-" || substr($sequence2,$position,1) eq "-" || substr($sequence3,$position,1) eq "-"){
362                               
363                                #print($gapStartPos . " gapStartPos\n");
364                                #print("starting a gap\n");
365                               
366                                #reset all gap related variables to 0
367                                $gapLenth = $gapLenth1 = $gapLenth2 = $gapLenth3 = 0;
368                                $gapStartPos = $gapStartPos1 = $gapStartPos2 = $gapStartPos3 = 0;
369                                $gapEndPos = 0;
370                                $getValidGap = 0;
371
372                                #if($position > 2 && $position < $test1 - 3){
373                                        #print(substr($sequence1,$position - 3, 3) . " before " . substr($sequence1, $position + 1, 3) . " after\n");
374                                        #print(substr($sequence2,$position - 3, 3) . " before " . substr($sequence2, $position + 1, 3) . " after\n");
375                                        #print(substr($sequence3,$position - 3, 3) . " before " . substr($sequence3, $position + 1, 3) . " after\n");
376                                #}
377                               
378                                #check if this position is the start of a gap
379                                if ($gapStartPos == 0){
380                                       
381                                        #the variables $gapStartPos and $gapEndPos are initialized with the current position value $position
382                                        $gapStartPos = $gapEndPos = $position;
383                                       
384                                        #iterate till the end of the gap and/or the end of any overlapping gap, whichever is greater
385                                        while (substr($sequence1, $gapEndPos, 1) eq "-" || substr($sequence2, $gapEndPos, 1) eq "-" || substr($sequence3,$gapEndPos,1) eq "-"){
386                                               
387                                                #check if the current single gap is in sequence1
388                                                if(substr($sequence1, $gapEndPos, 1) eq "-" && $gapEndPos < $test1){
389                                                        $gapLenth1++;
390                                                        if ($gapLenth1 == 1){
391                                                                $gapStartPos1 = $gapEndPos;
392                                                        }       
393                                                }
394                                               
395                                                #check if the current single gap is in sequence2
396                                                if(substr($sequence2, $gapEndPos, 1) eq "-" && $gapEndPos < $test1){
397                                                        $gapLenth2++;
398                                                        if ($gapLenth2 == 1){
399                                                                $gapStartPos2 = $gapEndPos;
400                                                        }       
401                                                }
402                                               
403                                                #check if the current single gap is in sequence3
404                                                if(substr($sequence3, $gapEndPos, 1) eq "-" && $gapEndPos < $test1){
405                                                        $gapLenth3++;
406                                                        if ($gapLenth3 == 1){
407                                                                $gapStartPos3 = $gapEndPos;
408                                                        }       
409                                                }
410                                                $gapEndPos++;
411                                        }
412                                       
413                                        #readjust the value of $gapStartPos to represent the index of the end of the gap
414                                        #and/or the end of any overlapping gaps, whichever is greater
415                                        $gapEndPos--;
416                                       
417                                        #print($gapStartPos . " - " . $gapEndPos . " gapStartPos - gapEndPos\n");
418                                        #print($gapStartPos1 . " - " . $gapLenth1 . " gapStartPos1 - gapLenth1\n");
419                                        #print($gapStartPos2 . " - " . $gapLenth2 . " gapStartPos2 - gapLenth2\n");
420                                        #print($gapStartPos3 . " - " . $gapLenth3 . " gapStartPos3 - gapLenth3\n");
421                               
422                                #variable $gapLength stores the length of the max gap among the three gaps if any
423                                if ($gapLenth < $gapLenth1){
424                                        $gapLenth = $gapLenth1;
425                                }
426                                if ($gapLenth < $gapLenth2){
427                                        $gapLenth = $gapLenth2;
428                                }
429                                if ($gapLenth < $gapLenth3){
430                                        $gapLenth = $gapLenth3;
431                                }
432                                       
433                                #print("The length of the longest gap is: " . $gapLenth . "\n");
434                               
435                                #check if there are two overlapping gaps from sequence1 and sequence2 but not from sequence3
436                                if ($gapLenth1 > 0 && $gapLenth2 > 0 && $gapLenth3 == 0){
437                                        if ($gapLenth1 != $gapLenth2 || $gapStartPos1 != $gapStartPos2){
438                                               
439                                                #print("Invalid case: gaps are overlapping starting at position: " . $position . "\n");
440                                               
441                                                # readjust the value of the variable $gapLenth to include all unequal-lenght gaps
442                                                $gapLenth = $gapEndPos - $position + 1;
443                                               
444                                                # readjust the value of the variable $position to the next non-gap character
445                                                $position = $position + $gapLenth - 1;
446                                               
447                                                #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");     
448                                               
449                                                #add as many characters "N" as the length of $gapLenth to the sequence $ABC as if there is no gap
450                                                while($gapLenth > 0){
451                                                                $ABC = join("",($ABC,"N"));
452                                                                $coord1++; $coord2++; $coord3++;
453                                                                $gapLenth--;
454                                                        }
455                                        }
456                                }
457                                else{   #check if there are two overlapping gaps from sequence1 and sequence3 but not from sequence2
458                                        if ($gapLenth1 > 0 && $gapLenth3 > 0 && $gapLenth2 == 0){
459                                               
460                                                if ($gapLenth1 != $gapLenth3 || $gapStartPos1 != $gapStartPos3){
461                                                       
462                                                        #print("Invalid case: gaps are overlapping starting at position: " . $position . "\n");
463                                                       
464                                                        # readjust the value of the variable $gapLenth to include all unequal-lenght gaps
465                                                        $gapLenth = $gapEndPos - $position + 1;
466                                       
467                                                        # readjust the value of the variable $position to the next non-gap character
468                                                        $position = $position + $gapLenth - 1;
469                                                       
470                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
471                                                       
472                                                        #add as many characters "N" as the length of $gapLenth to the sequence $ABC as if there is no gap
473                                                        while($gapLenth > 0){
474                                                                        $ABC = join("",($ABC,"N"));
475                                                                        $coord1++; $coord2++; $coord3++;
476                                                                        $gapLenth--;
477                                                                }
478                                                }
479                                        }
480                                        else{    #check if there are two overlapping gaps from sequence2 and sequence3 but not from sequence1
481                                                if ($gapLenth2 > 0 && $gapLenth3 > 0 && $gapLenth1 == 0){
482                                                        if ($gapLenth2 != $gapLenth3 || $gapStartPos2 != $gapStartPos3){
483                                                               
484                                                                #print("Invalid case: gaps are overlapping starting at position: " . $position . "\n");
485                                                               
486                                                                # readjust the value of the variable $gapLenth to include all unequal-lenght gaps
487                                                                $gapLenth = $gapEndPos - $position + 1;
488                                                               
489                                                                # readjust the value of the variable $position to the next non-gap character
490                                                                $position = $position + $gapLenth - 1;
491                                                               
492                                                                #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
493                                                       
494                                                                #add as many characters "N" as the length of $gapLenth to the sequence $ABC as if there is no gap
495                                                                while($gapLenth > 0){
496                                                                                $ABC = join("",($ABC,"N"));
497                                                                                $coord1++; $coord2++; $coord3++;
498                                                                                $gapLenth--;
499                                                                        }
500                                                        }
501                                                }
502                                                else{ #check if there are three overlapping gaps from the three sequences
503                                                        if ($gapLenth1 > 0 && $gapLenth2 > 0 && $gapLenth3 > 0){
504                                                                if ($gapLenth1 != $gapLenth2 || $gapLenth1 != $gapLenth3 ||  $gapLenth2 != $gapLenth3
505                                                                    || $gapStartPos1 != $gapStartPos2 || $gapStartPos1 != $gapStartPos3
506                                                                    || $gapStartPos2 != $gapStartPos3){
507                                                                   
508                                                                    #print("Invalid case: gaps are overlapping starting at position: " . $position . "\n");
509                                                                   
510                                                                    # readjust the value of the variable $gapLenth to include all unequal-lenght gaps
511                                                                        $gapLenth = $gapEndPos - $position + 1;
512                                       
513                                                                    # readjust the value of the variable $position to the next non-gap character
514                                                                    $position = $position + $gapLenth - 1;
515                                                                   
516                                                                    #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");                                                 
517                                                       
518                                                                    #add as many characters "N" as the length of $gapLenth to the sequence $ABC as if there is no gap
519                                                                        while($gapLenth > 0){
520                                                                                        $ABC = join("",($ABC,"N"));
521                                                                                        $coord1++; $coord2++; $coord3++;
522                                                                                        $gapLenth--;
523                                                                                }
524                                                                }
525                                                        }
526                                                }
527                                        }
528                                }
529                               
530                                #check if the gaps are overlapping
531                                if ($gapLenth == 0){
532                                               
533                                        #check if the overlapping gaps are also within less than three base from either end of the block
534                                                if ($gapStartPos <= 2 || $gapEndPos >= ($test1 - 3)){
535                                                        if ($gapStartPos <= 2){
536                                                                #print("Invalid case: at least one of the overlapping gaps is within less than 3 bases from the left end of the block starting at position: " . $gapStartPos . "\n");
537                                                        }
538                                                        elsif ($gapEndPos >= $test1 - 3){
539                                                                        #print("Invalid case: at least one of the overlapping gaps is within less than 3 bases from the right end of the block starting at position: " . $gapStartPos . "\n");
540                                                        }
541                                                        elsif ($gapStartPos <= 2 && $gapEndPos >= ($test1 - 3)){
542                                                                #print("Invalid case: at least one of the overlapping gaps is within less than 3 bases from each end of the block starting at position: " . $gapStartPos . "\n");
543                                                        }
544                                                }
545                                               
546                                                #check if the gap starts at position located within three bases from the start of the block
547                                                my $preSubSeqLength = 0;
548                                                if ($gapStartPos - 3 < 0){
549                                                        $preSubSeqLength = $gapStartPos;
550                                                }
551                                                else{
552                                                        $preSubSeqLength = 3;
553                                                }
554                                               
555                                                #check if the there is a masked character among the three bases on at least one side of the gap
556                                        if ((substr($sequence1, $gapStartPos - $preSubSeqLength, $preSubSeqLength)     =~ m/[-\#$?^@]/)
557                                                 || (substr($sequence2, $gapStartPos - $preSubSeqLength, $preSubSeqLength) =~ m/[-\#$?^@]/)
558                                                 || (substr($sequence3, $gapStartPos - $preSubSeqLength, $preSubSeqLength) =~ m/[-\#$?^@]/)
559                                                 || (substr($sequence1, $gapEndPos + 1, 3) =~ m/[-\#$?^@]/)
560                                                 || (substr($sequence2, $gapEndPos + 1, 3) =~ m/[-\#$?^@]/)
561                                                 || (substr($sequence3, $gapEndPos + 1, 3) =~ m/[-\#$?^@]/)){
562                                                       
563                                                        #print("Invalid case: at least one masked character is within 3 bases from at least one of overlapping gaps starting at position: " . $gapStartPos . "\n");
564                                                }
565                                               
566                                                #check if there is at least one base with low score (maskd) within the gap interval     
567                                                if((substr($sequence1, $gapStartPos, $gapEndPos - $gapStartPos + 1)    =~ m/[\#$?^@]/)
568                                                    || (substr($sequence2, $gapStartPos, $gapEndPos - $gapStartPos + 1) =~ m/[\#$?^@]/)
569                                                    || (substr($sequence3, $gapStartPos, $gapEndPos - $gapStartPos + 1) =~ m/[\#$?^@]/)){
570   
571                                                    #print("Invalid case: at least one masked character is within the gap interval starting at position: " . $gapStartPos . "\n");
572                                                }       
573                                               
574                                                #print($sequence1 ." sequence1\n");
575                                                #print($sequence2 ." sequence2\n");
576                                                #print($sequence3 ." sequence3\n");             
577                                    #print("________________________________________________________________\n");
578                                }           
579                               
580                                       
581                                        #if $gapLenth > 0, then there is no overlapping gaps, which means that either there is only one gap
582                                        # or two gaps that have equal lengths and start at the same position
583                                        if ($gapLenth > 0){ 
584       
585                                                #variable indicating if the gap is within less than three bases from the end of the blocks
586                                                my $blockEndGap = 0;
587                                                #variable indicating if the gap is within less than three bases from the end of the blocks
588                                                my $maskedCharWithinGap = 0;                   
589                                                #variable indicating if there is any masked character among the three bases flanking the gap from each side
590                                                my $maskedCharNearGap = 0;
591                                               
592                                                #check if there are at least there bases on each side of the gap                               
593                                                if ($position > 2 && ($position + $gapLenth - 1) < ($test1 - 3)){
594                                                        $blockEndGap = 0;
595                                                }
596                                                else{   #the gap is within less than 3 bases from the end of the block
597       
598                                                        $blockEndGap = 1;
599                                                }
600                                               
601                                                #variable $preSubSeqLength stores the number of bases available on the left side of the gap
602                                                my $preSubSeqLength = 0;
603                                                #variable $leftStartCheckingPosition stores the index of the base at which
604                                                #we start checking the number of bases available on the left side of the gap
605                                                my $leftStartCheckingPosition = 0;
606                                               
607                                                if ($gapStartPos - 3 < 0){
608                                                        $preSubSeqLength = $gapStartPos;
609                                                }
610                                                else{
611                                                        $preSubSeqLength = 3;
612                                                        $leftStartCheckingPosition = $position - 3;
613                                                }
614                                                #check if all the three bases on each side of the gap have high scores (unmasked)
615                                                if ((substr($sequence1, $leftStartCheckingPosition, $preSubSeqLength) !~ m/[-\#$?^@]/)
616                                                     && (substr($sequence2, $leftStartCheckingPosition, $preSubSeqLength) !~ m/[-\#$?^@]/)
617                                                     && (substr($sequence3, $leftStartCheckingPosition, $preSubSeqLength) !~ m/[-\#$?^@]/)
618                                                     && (substr($sequence1, $position + $gapLenth, 3) !~ m/[-\#$?^@]/)
619                                                     && (substr($sequence2, $position + $gapLenth, 3) !~ m/[-\#$?^@]/)
620                                                     && (substr($sequence3, $position + $gapLenth, 3) !~ m/[-\#$?^@]/)){
621               
622                                                        #there are three bases with high scores (unmasked) on each side of the gap
623                                                        $maskedCharNearGap = 0;
624                                                }
625                                                else{   #three is a masked character among the three bases on at least one side of the gap
626                                                        $maskedCharNearGap = 1;
627                                                }
628                                                   
629                                                #check if there is at least one base with low score (maskd) within the gap interval     
630                                                if((substr($sequence1, $position, $gapLenth)     !~ m/[\#$?^@]/)
631                                                    && (substr($sequence2, $position, $gapLenth) !~ m/[\#$?^@]/)
632                                                    && (substr($sequence3, $position, $gapLenth) !~ m/[\#$?^@]/)){
633   
634                                                    #there is no base with low score (unmasked) within the gap interval         
635                                                        $maskedCharWithinGap = 0;
636                                                }
637                                                else{   #there is at least one base with low score (maskd) within the gap interval     
638                                                        $maskedCharWithinGap = 1;
639                                                }
640                                       
641                                        #check if the gap is, valid which means:
642                                        #check if the gap is within at least 3 bases from each end of the block,
643                                                #and that all the three bases on each side of the block have high scores (unmasked),
644                                                #and that there is no base with low score (masked) within the gap interval at the same time
645                                        if ($blockEndGap == 0 && $maskedCharNearGap == 0 && $maskedCharWithinGap == 0){
646                                               
647                                                #print("The gap is valid and will be treated accordingly\n");
648                                                $getValidGap = 1
649                                        }
650                                       
651                                        #the gap is not valid, which means:
652                                        #the gap is not valid, so check if the gap is within less than 3 bases from the end of the block,
653                                                #and that there is a masked character among the three bases on at least one side of the gap,
654                                                #and that there is a masked character within the gap interval at the same time
655                                                elsif ($blockEndGap == 1 && $maskedCharNearGap == 1 && $maskedCharWithinGap == 1){
656                                                        #print("Invalid case: the gap is within less than 3 bases from the end of the block starting at position: " . $position . "\n");
657                                                        #print("Invalid case: at least one masked character is within 3 bases from the gap starting at position: " . $position . "\n");
658                                                        #print("Invalid case: at least one masked character is within the gap interval starting at position: " . $position . "\n");
659                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
660                                                }
661                                                #check if the gap is within less than 3 bases from the end of the block, and
662                                                #that there is a masked character among the three bases on at least one side of the gap at the same time
663                                                elsif ($blockEndGap == 1 && $maskedCharNearGap == 1){
664                                                        #print("Invalid case: the gap is within less than 3 bases from the end of the block starting at position: " . $position . "\n");
665                                                        #print("Invalid case: at least one masked character is within 3 bases from the gap starting at position: " . $position . "\n");
666                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
667                                                }
668                                                #check if the gap is within less than 3 bases from the end of the block and
669                                                #and that there is a masked character within the gap interval at the same time
670                                                elsif ($blockEndGap == 1 && $maskedCharWithinGap == 1){
671                                                        #print("Invalid case: the gap is within less than 3 bases from the end of the block starting at position: " . $position . "\n");
672                                                        #print("Invalid case: at least one masked character is within the gap interval starting at position: " . $position . "\n");
673                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
674                                                }
675                                                #check if there is a masked character among the three bases on at least one side of the gap
676                                                #and that there is a masked character within the gap interval at the same time
677                                                elsif ($maskedCharNearGap == 1 && $maskedCharWithinGap == 1){
678                                                        #print("Invalid case: at least one masked character is within 3 bases from the gap starting at position: " . $position . "\n");
679                                                        #print("Invalid case: at least one masked character is within the gap interval starting at position: " . $position . "\n");
680                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
681                                                }
682                                                #check if the gap is within less than 3 bases from the end of the block
683                                                elsif ($blockEndGap == 1){
684                                                                #print("Invalid case: the gap is within less than 3 bases from the end of the block starting at position: " . $position . "\n");
685                                                                #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
686                                                }
687                                                #check if there is a masked character within 3 bases on each side of the gap
688                                                elsif ($maskedCharNearGap == 1){
689                                                        #print("Invalid case: at least one masked character is within 3 bases from the gap starting at position: " . $position . "\n");
690                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
691                                                }
692                                                #check if there is a masked character within the gap interval
693                                                elsif ($maskedCharWithinGap == 1){
694                                                        #print("Invalid case: at least one masked character is within the gap interval starting at position: " . $position . "\n");
695                                                        #print("The total length of the overlapping gaps is: " . $gapLenth . "\n");
696                                                }
697                                               
698                                                #check if the gap is not valid in order to readjust the $position value and increment coordinates varaibles
699                                                if ($getValidGap == 0){
700                                                        #print($sequence1 . " sequence1\n");
701                                                        #print($sequence2 . " sequence2\n");
702                                                        #print($sequence3 . " sequence3\n");
703                                                        #print("________________________________________________________________\n");
704                                               
705                                                        # readjust the value of the variable $position to the next non-gap character
706                                                        $position = $position + $gapLenth - 1;
707                               
708                                                        #add as many characters "N" as the length of $gapLenth to the sequence $ABC as if there is no gap
709                                                        while($gapLenth > 0){
710                                                                $ABC = join("",($ABC,"N"));
711                                                                $coord1++; $coord2++; $coord3++;
712                                                                $gapLenth--;
713                                                        }
714                                                }                                       
715                                        }
716                                }                       
717                        }
718                        else{   #there is no gap at this position
719                                #if ($position <= 2){
720                                        #print("invalid case\n");
721                                        #print(substr($sequence1, 0, $position) ." before\n");
722                                        #print(substr($sequence2, 0, $position) ." before\n");
723                                        #print(substr($sequence3, 0, $position) ." before\n");
724                                #}
725                                #if ($position >= ($test1 - 3)){
726                                        #print("invalid case\n");
727                                        #print(substr($sequence1, $position - $test1) . " after\n");
728                                        #print(substr($sequence2, $position - $test1) . " after\n");
729                                        #print(substr($sequence3, $position - $test1) . " after\n");
730                                #}
731                               
732                                #add as many characters "N" as the length of $gapLenth to the sequence $ABC becasue there is no gap
733                                $ABC = join("",($ABC,"N"));
734                                $coord1++; $coord2++; $coord3++;
735                               
736                                # reset the value of the variable $gapStartPos to 0
737                                if ($gapStartPos > 0){
738                                        $gapStartPos = 0;
739                                }
740                               
741                                # reset the value of the variable $getValidGap to 0
742                                if ($getValidGap > 0){
743                                        $getValidGap = 0;
744                                }
745                        } 
746       
747                        #The gap is valid and will be treated accordingly
748                        if ($getValidGap > 0){
749               
750                                my $localGap = $gapLenth;
751                               
752                                # seq1 deletes
753                                 if ((substr($sequence1,$position,1) eq "-")
754                                        && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
755                                        && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
756                                       
757                                        while($gapLenth > 0){
758                                                        $ABC = join("",($ABC,"X"));
759                                                        my @s = split(/:/, $seq1);
760                                                        $indelType = $s[0]."_delete";
761               
762                                                        #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";     
763                                                        $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
764                                                        push (@indels,$indel_line);
765                                                        push (@seq1_delete,$indel_line);
766                                                        $coord2++;
767                                                        $coord3++;
768                                                       
769                                                        $gapLenth--;
770                                        }
771                                }       
772                                # seq2 deletes
773                                elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
774                                        && (substr($sequence2,$position,1) eq "-")
775                                        && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
776                                               
777                                                while($gapLenth > 0){
778                                                        $ABC = join("",($ABC,"Y"));
779                                                        my @s = split(/:/, $seq2);
780                                                        $indelType = $s[0]."_delete";
781                                                       
782                                                        #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
783                                                        $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
784                                    push (@indels,$indel_line);
785                                                        push (@seq2_delete,$indel_line);
786                                                        $coord1++;
787                                                        $coord3++;
788                                                       
789                                                        $gapLenth--;
790                                                }
791                                }
792                                # seq1 inserts
793                                elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
794                                        && (substr($sequence2,$position,1) eq "-")
795                                        && (substr($sequence3,$position,1) eq "-")){
796                                               
797                                                while($gapLenth > 0){
798                                                        $ABC = join("",($ABC,"Z"));
799                                                        my @s = split(/:/, $seq1);
800                                                        $indelType = $s[0]."_insert";
801                                                       
802                                                        #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
803                                                        $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
804                                                        push (@indels,$indel_line);
805                                                        push (@seq1_insert,$indel_line);
806                                                        $coord1++;
807                                                       
808                                                        $gapLenth--;
809                                                }
810                                }
811                                # seq2 inserts 
812                                elsif ((substr($sequence1,$position,1) eq "-")
813                                        && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
814                                        && (substr($sequence3,$position,1) eq "-")){
815                                               
816                                                while($gapLenth > 0){
817                                                        $ABC = join("",($ABC,"W"));
818                                                        my @s = split(/:/, $seq2);
819                                                        $indelType = $s[0]."_insert";
820                                                       
821                                                        #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
822                                                        $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
823                                                        push (@indels,$indel_line);
824                                                        push (@seq2_insert,$indel_line);
825                                                        $coord2++;
826                                                       
827                                                        $gapLenth--;
828                                                }
829                                }
830                                # seq3 deletes
831                                elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
832                                        && (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
833                                        && (substr($sequence3,$position,1) eq "-")){
834                                               
835                                                while($gapLenth > 0){
836                                                        $ABC = join("",($ABC,"S"));
837                                                        my @s = split(/:/, $seq3);
838                                                        $indelType = $s[0]."_delete";
839                                                       
840                                                        #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
841                                                        $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
842                                                        push (@indels,$indel_line);
843                                                        push (@seq3_delete,$indel_line);
844                                                        $coord1++;
845                                                        $coord2++;
846                                                       
847                                                        $gapLenth--;
848                                                }
849                                }
850                                # seq3 inserts
851                                elsif ((substr($sequence1,$position,1) eq "-")
852                                        && (substr($sequence2,$position,1) eq "-")
853                                        && (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
854                                               
855                                                while($gapLenth > 0){
856                                                        $ABC = join("",($ABC,"T"));
857                                                        my @s = split(/:/, $seq3);
858                                                        $indelType = $s[0]."_insert";
859                                                       
860                                                        #print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
861                                                        $indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
862                                                        push (@indels,$indel_line);
863                                                        push (@seq3_insert,$indel_line);
864                                                        $coord3++;
865                                                       
866                                                        $gapLenth--;
867                                                }
868                                }else{
869                                        while($gapLenth > 0){
870                                                $ABC = join("",($ABC,"N"));
871                                                $coord1++; $coord2++; $coord3++;
872                                               
873                                                $gapLenth--;
874                                        }
875                                }
876                               
877                                # readjust the value of the variable $position to indicate the last character of the gap
878                                $position = $position + $localGap - 1;
879                        }
880                       
881                }
882                @array_return=($seq1,$seq2,$seq3,$ABC);
883                return (@array_return);
884
885        }
886        # ignore pairwise cases for now, just count the number of blocks
887        elsif (scalar(@sequences) == 2){
888                my $ABC = "";
889                my $coord1 = my $coord2 = my $coord3 = 0;
890                $count3++;
891
892                $line1 = $sequences[0];
893                $line2 = $sequences[1];
894                chomp ($line1);
895                chomp ($line2);
896
897                if ($line2 !~ m/$ingroup2/){
898                        $count4++;
899                }
900        }
901}
902
903               
904sub get_indels_lengths{
905        my (@array) = @_;
906        if (scalar(@array) == 4){
907                my $seq1 = $array[0]; my $seq2 = $array[1]; my $seq3 = $array[2]; my $ABC = $array[3];
908
909                while ($ABC =~ m/(X+)/g) {
910                        push (@seq1_delete_lengths,length($1));
911                        $count_seq1delete++;
912                }
913               
914                while ($ABC =~ m/(Y+)/g) {
915                        push (@seq2_delete_lengths,length($1));
916                        $count_seq2delete++;
917                }
918                while ($ABC =~ m/(S+)/g) {
919                        push (@seq3_delete_lengths,length($1));
920                        $count_seq3delete++;
921                }
922                while ($ABC =~ m/(Z+)/g) {     
923                        push (@seq1_insert_lengths,length($1));
924                        $count_seq1insert++;
925                }
926                while ($ABC =~ m/(W+)/g) {
927                        push(@seq2_insert_lengths,length($1));
928                        $count_seq2insert++;
929                }
930                while ($ABC =~ m/(T+)/g) {
931                        push (@seq3_insert_lengths,length($1));
932                        $count_seq3insert++;
933                }
934        }
935        elsif (scalar(@array) == 0){
936                next;
937        }
938               
939}
940# convert to coordinates to + strand if align orient = -
941sub convert_coords{
942        my (@array) = @_;
943        my $s = $array[0];
944        my $e = $array[1];
945        my $o = $array[2];
946        my $l = $array[3];
947        my $start_plus = my $end_plus = 0;
948
949        if ($o eq "-"){
950                $start_plus = ($l - $e);
951                $end_plus = ($l - $s);
952        }elsif ($o eq "+"){
953                $start_plus = $s;
954                $end_plus = $e;
955        }
956
957        return ($start_plus, $end_plus);
958}
959
960# get first line only for each event
961sub get_starts_only{
962        my (@test) = @_;
963    my $seq1 = my $seq2 = my $seq3 = my $indelType = my $old_seq1 = my $old_seq2 = my $old_seq3 = my $old_indelType = my $old_line = "";
964    my $coord1 = my $coord2 = my $coord3 = my $old_coord1 = my $old_coord2 = my $old_coord3 = 0;
965
966    my @matches = ();
967    my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = ();
968                               
969       
970    foreach my $line (@test){
971        chomp($line);
972        $line =~ s/^\s*//;
973        $line =~ s/\s+/\t/g;
974                my @line1 = split(/\t/, $line);
975        $seq1 = $line1[1];
976        $coord1 = $line1[2];
977        $seq2 = $line1[4];
978        $coord2 = $line1[5];
979        $seq3 = $line1[7];
980        $coord3 = $line1[8];
981        $indelType = $line1[10];
982               
983        if ($indelType =~ m/$ingroup1/ && $indelType =~ m/insert/){
984                if ($coord1 != $old_coord1+1 || ($coord2 != $old_coord2 || $coord3 != $old_coord3)){
985                        $start1++;
986                push (@seq1_insert_startOnly,$line);
987            }
988            }
989        elsif ($indelType =~ m/$ingroup1/ && $indelType =~ m/delete/){
990                        if ($coord1 != $old_coord1 || ($coord2 != $old_coord2+1 || $coord3 != $old_coord3+1)){
991                        $start2++;
992                        push(@seq1_delete_startOnly,$line);
993                    }
994                }
995        elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/insert/){
996                if ($coord2 != $old_coord2+1 || ($coord1 != $old_coord1 || $coord3 != $old_coord3)){
997                        $start3++;
998                        push(@seq2_insert_startOnly,$line);
999                    }
1000                }
1001        elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/delete/){
1002                if ($coord2 != $old_coord2 || ($coord1 != $old_coord1+1 || $coord3 != $old_coord3+1)){
1003                $start4++;
1004                push(@seq2_delete_startOnly,$line);
1005            }
1006        }
1007        elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/insert/){
1008            if ($coord3 != $old_coord3+1 || ($coord1 != $old_coord1 || $coord2 != $old_coord2)){
1009                $start5++;
1010                 push(@seq3_insert_startOnly,$line);
1011            }
1012        }
1013        elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/delete/){
1014                if ($coord3 != $old_coord3 || ($coord1 != $old_coord1+1 ||$coord2 != $old_coord2+1)){
1015                $start6++;
1016                push(@seq3_delete_startOnly,$line);
1017            }
1018        }
1019       
1020        $old_indelType = $indelType;
1021        $old_seq1 = $seq1;
1022        $old_coord1 = $coord1;
1023        $old_seq2 = $seq2;
1024        $old_coord2 = $coord2;
1025        $old_seq3 = $seq3;
1026        $old_coord3 = $coord3;
1027        $old_line = $line;
1028        }
1029}
1030# append lengths to each event start line
1031my $counter1; my $counter2; my $counter3; my $counter4; my $counter5; my $counter6;
1032my @final1 = my @final2 = my @final3 = my @final4 = my @final5 = my @final6 = ();
1033my $final_line1 = my $final_line2 = my $final_line3 = my $final_line4 = my $final_line5 = my $final_line6 = "";
1034
1035
1036for ($counter1 = 0; $counter1 < @seq3_insert_startOnly; $counter1++){
1037        $final_line1 = join("\t",($seq3_insert_startOnly[$counter1],$seq3_insert_lengths[$counter1]));
1038        push (@final1,$final_line1);
1039}
1040
1041for ($counter2 = 0; $counter2 < @seq3_delete_startOnly; $counter2++){
1042    $final_line2 =  join("\t",($seq3_delete_startOnly[$counter2],$seq3_delete_lengths[$counter2]));
1043    push(@final2,$final_line2);
1044}
1045
1046for ($counter3 = 0; $counter3 < @seq2_insert_startOnly; $counter3++){
1047    $final_line3 =  join("\t",($seq2_insert_startOnly[$counter3],$seq2_insert_lengths[$counter3]));
1048        push(@final3,$final_line3);
1049}
1050
1051for ($counter4 = 0; $counter4 < @seq2_delete_startOnly; $counter4++){
1052    $final_line4 =  join("\t",($seq2_delete_startOnly[$counter4],$seq2_delete_lengths[$counter4]));
1053    push(@final4,$final_line4);
1054}
1055               
1056for ($counter5 = 0; $counter5 < @seq1_insert_startOnly; $counter5++){
1057    $final_line5 =  join("\t",($seq1_insert_startOnly[$counter5],$seq1_insert_lengths[$counter5]));
1058    push(@final5,$final_line5);
1059}
1060
1061for ($counter6 = 0; $counter6 < @seq1_delete_startOnly; $counter6++){
1062    $final_line6 =  join("\t",($seq1_delete_startOnly[$counter6],$seq1_delete_lengths[$counter6]));
1063    push(@final6,$final_line6);
1064}       
1065
1066# format final output
1067# # if inserts, increase coords for the sequence inserted, other sequences give coords for 5' and 3' bases flanking the gap
1068# # for deletes, increase coords for other 2 sequences and the one deleted give coords for 5' and 3' bases flanking the gap
1069
1070get_final_format(@final5);
1071get_final_format(@final6);
1072get_final_format(@final3);
1073get_final_format(@final4);
1074get_final_format(@final1);
1075get_final_format(@final2);
1076
1077sub get_final_format{
1078        my (@final) = @_;
1079        foreach (@final){
1080                my $event_line = $_;
1081                my @events = split(/\t/, $event_line);
1082                my $event_type = $events[10];
1083                my @name_align1 = split(/:/, $events[1]);
1084                my @name_align2 = split(/:/, $events[4]);
1085                my @name_align3 = split(/:/, $events[7]);
1086                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;
1087                my $final_event_line = "";     
1088                # seq1_insert
1089                if ($event_type =~ m/$ingroup1/ && $event_type =~ m/insert/){
1090                        # only increase coord for human
1091                        # remember that other two sequnences, the gap spans (coord - 1) --> coord
1092                        $seq1_event_start = ($events[2]);
1093                        $seq1_event_end = ($events[2]+$events[11]-1);
1094                        $seq2_event_start = ($events[5]-1);
1095                        $seq2_event_end = ($events[5]);
1096                        $seq3_event_start = ($events[8]-1);
1097                        $seq3_event_end = ($events[8]);
1098                        $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]));
1099                }
1100                # seq1_delete
1101                elsif ($event_type =~ m/$ingroup1/ && $event_type =~ m/delete/){
1102                        # only increase coords for seq2 and seq3
1103                        # remember for seq1, the gap spans (coord - 1) --> coord
1104                        $seq1_event_start = ($events[2]-1);
1105                        $seq1_event_end = ($events[2]);
1106            $seq2_event_start = ($events[5]);
1107            $seq2_event_end = ($events[5]+$events[11]-1);
1108            $seq3_event_start = ($events[8]);
1109            $seq3_event_end = ($events[8]+$events[11]-1);
1110                        $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]));
1111                }
1112                # seq2_insert
1113                elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/insert/){       
1114                        # only increase coords for seq2
1115                        # remember that other two sequnences, the gap spans (coord - 1) --> coord
1116            $seq1_event_start = ($events[2]-1);
1117            $seq1_event_end = ($events[2]);
1118                        $seq2_event_start = ($events[5]);
1119            $seq2_event_end = ($events[5]+$events[11]-1);
1120            $seq3_event_start = ($events[8]-1);
1121                        $seq3_event_end = ($events[8]);                 
1122                        $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]));
1123                }
1124                # seq2_delete
1125                elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/delete/){
1126                        # only increase coords for seq1 and seq3
1127                        # remember for seq2, the gap spans (coord - 1) --> coord
1128            $seq1_event_start = ($events[2]);
1129                        $seq1_event_end = ($events[2]+$events[11]-1);   
1130            $seq2_event_start = ($events[5]-1);
1131                $seq2_event_end = ($events[5]);
1132            $seq3_event_start = ($events[8]);
1133            $seq3_event_end = ($events[8]+$events[11]-1);
1134                        $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]));
1135                }       
1136                # start testing w/seq3_insert
1137                elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/insert/){
1138                        # only increase coord for rheMac
1139                        # remember that other two sequnences, the gap spans (coord - 1) --> coord
1140                        $seq1_event_start = ($events[2]-1);
1141                        $seq1_event_end = ($events[2]);
1142                        $seq2_event_start = ($events[5]-1);
1143                        $seq2_event_end = ($events[5]);
1144                        $seq3_event_start = ($events[8]);
1145                        $seq3_event_end = ($events[8]+$events[11]-1);
1146                        $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]));
1147                }
1148                # seq3_delete
1149                elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/delete/){
1150                        # only increase coords for seq1 and seq2
1151                        # remember for seq3, the gap spans (coord - 1) --> coord
1152                        $seq1_event_start = ($events[2]);
1153                        $seq1_event_end = ($events[2]+$events[11]-1);
1154                        $seq2_event_start = ($events[5]);
1155                        $seq2_event_end = ($events[5]+$events[11]-1);
1156                        $seq3_event_start = ($events[8]-1);
1157                        $seq3_event_end = ($events[8]);
1158                        $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]));
1159
1160                }
1161               
1162                print OFILE "$final_event_line\n";
1163        }
1164}
1165close OFILE;
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。