Context Navigation

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
14	use strict;
15	use warnings;
16
17	# declare and initialize variables
18	my $fh; # variable to store filehandle
19	my $record;
20	my $offset;
21	my $library = $ARGV[0];
22	my $count = 0;
23	my $count2 = 0;
24	my $count3 = 0;
25	my $count4 = 0;
26	my $start1 = my $start2 = my $start3 = my $start4 = my $start5 = my $start6 = 0;
27	my $orient = "";
28	my $outgroup = $ARGV[2];
29	my $ingroup1 = my $ingroup2 = "";
30	my $count_seq1insert = my $count_seq1delete = 0;
31	my $count_seq2insert = my $count_seq2delete = 0;
32	my $count_seq3insert = my $count_seq3delete = 0;
33	my @seq1_insert_lengths = my @seq1_delete_lengths = ();
34	my @seq2_insert_lengths = my @seq2_delete_lengths = ();
35	my @seq3_insert_lengths = my @seq3_delete_lengths = ();
36	my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = ();
37	my @seq1_insert_startOnly = my @seq1_delete_startOnly = my @seq2_insert_startOnly = my @seq2_delete_startOnly = ();
38	my @seq3_insert_startOnly = my @seq3_delete_startOnly = ();
39	my @indels = ();
40
41	# check to make sure correct files
42	my $usage = "usage: parseMAF_smallIndels.pl [MAF.in] [small_Indels_summary_TABULAR.out] [outgroup]\n";
43	die $usage unless @ARGV == 3;
44
45	# perform some standard subroutines
46	$fh = open_file($library);
47
48	$offset = tell($fh);
49
50	my $ofile = $ARGV[1];
51	unless (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";
59	print 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
62	while ($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
76	get_starts_only(@seq1_insert);
77	get_starts_only(@seq1_delete);
78	get_starts_only(@seq2_insert);
79	get_starts_only(@seq2_delete);
80	get_starts_only(@seq3_insert);
81	get_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
88	if ($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
93	print "Ingroup1 = $ingroup1; Ingroup2 = $ingroup2; Outgroup = $outgroup;\n";
94	print "# of ingroup1 inserts = $count_seq1insert;\n";
95	print "# of ingroup1 deletes = $count_seq1delete;\n";
96	print "# of ingroup2 inserts = $count_seq2insert;\n";
97	print "# of ingroup2 deletes = $count_seq2delete;\n";
98	print "# of outgroup3 inserts = $count_seq3insert;\n";
99	print "# of outgroup3 deletes = $count_seq3delete\n";
100
101
102	#close OFILE;
103
104	if ($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
112	sub 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
124	sub 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
140	sub get_sequences_within_block{
141	my (@alignment) = @_;
142	my @lines = ();
143
144	my @sequences = ();
145
146	@lines = split ("\n", $record);
147	foreach (@lines){
148	chomp($_);
149	if (m/^\s*$/){
150	next;
151	}
152	elsif (m/^\s=(\d+\.\d*)/){
153
154	}elsif (m/^\ss(.)$/){
155	$count2++;
156
157	push (@sequences,$_);
158	}
159	}
160
161	return @sequences;
162	}
163
164	#get the indels within the blocks
165	sub 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
904	sub 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 = -
941	sub 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
961	sub 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
1031	my $counter1; my $counter2; my $counter3; my $counter4; my $counter5; my $counter6;
1032	my @final1 = my @final2 = my @final3 = my @final4 = my @final5 = my @final6 = ();
1033	my $final_line1 = my $final_line2 = my $final_line3 = my $final_line4 = my $final_line5 = my $final_line6 = "";
1034
1035
1036	for ($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
1041	for ($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
1046	for ($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
1051	for ($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
1056	for ($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
1061	for ($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
1070	get_final_format(@final5);
1071	get_final_format(@final6);
1072	get_final_format(@final3);
1073	get_final_format(@final4);
1074	get_final_format(@final1);
1075	get_final_format(@final2);
1076
1077	sub 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	}
1165	close OFILE;

Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。

Context Navigation

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

異なるフォーマットでダウンロード: