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

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