PageRenderTime 70ms CodeModel.GetById 30ms RepoModel.GetById 0ms app.codeStats 0ms

/tools/regVariation/parseMAF_smallIndels.pl

https://bitbucket.org/chapmanb/galaxy-central/
Perl | 698 lines | 552 code | 73 blank | 73 comment | 101 complexity | 52a8c0d056cf5c5467079a7707586383 MD5 | raw file
    

  1. #!/usr/bin/perl -w
    2. 

  2. # a program to get indels 
    3. 

  3. # input is a MAF format 3-way alignment file
    4. 

  4. # from 3-way blocks only at this time
    5. 

  5. # translate seq2, seq3, etc coordinates to + if align orient is reverse complement
    6. 

  6.  
    7. 

  7. use strict;
    8. 

  8. use warnings;
    9. 

  9. 

  10. # declare and initialize variables
    11. 

  11. my $fh; # variable to store filehandle
    12. 

  12. my $record;
    13. 

  13. my $offset;
    14. 

  14. my $library = $ARGV[0]; 
    15. 

  15. my $count = 0;
    16. 

  16. my $count2 = 0;
    17. 

  17. my $count3 = 0;
    18. 

  18. my $count4 = 0;
    19. 

  19. my $start1 = my $start2 = my $start3 = my $start4 = my $start5 = my $start6 = 0;
    20. 

  20. my $orient = "";
    21. 

  21. my $outgroup = $ARGV[2];
    22. 

  22. my $ingroup1 = my $ingroup2 = "";
    23. 

  23. my $count_seq1insert = my $count_seq1delete = 0;
    24. 

  24. my $count_seq2insert = my $count_seq2delete = 0;
    25. 

  25. my $count_seq3insert = my $count_seq3delete = 0;
    26. 

  26. my @seq1_insert_lengths = my @seq1_delete_lengths = ();
    27. 

  27. my @seq2_insert_lengths = my @seq2_delete_lengths = ();
    28. 

  28. my @seq3_insert_lengths = my @seq3_delete_lengths = ();
    29. 

  29. my @seq1_insert =  my @seq1_delete =  my @seq2_insert =  my @seq2_delete =  my @seq3_insert =  my @seq3_delete = ();
    30. 

  30. my @seq1_insert_startOnly = my @seq1_delete_startOnly = my @seq2_insert_startOnly = my @seq2_delete_startOnly = ();
    31. 

  31. my @seq3_insert_startOnly = my @seq3_delete_startOnly = ();
    32. 

  32. my @indels = (); 
    33. 

  33. 

  34. # check to make sure correct files
    35. 

  35. my $usage = "usage: parseMAF_smallIndels.pl [MAF.in] [small_Indels_summary.out] [outgroup]\n";
    36. 

  36. die $usage unless @ARGV == 3;
    37. 

  37. 

  38. # perform some standard subroutines 
    39. 

  39. $fh = open_file($library);
    40. 

  40. 

  41. $offset = tell($fh);
    42. 

  42. 

  43. #my $ofile = $ARGV[2];
    44. 

  44. #unless (open(OFILE, ">$ofile")){
    45. 

  45. #	 print "Cannot open output file \"$ofile\"\n\n";
    46. 

  46. #	 exit;
    47. 

  47. #}
    48. 

  48. 

  49. my $ofile2 = $ARGV[1];
    50. 

  50. unless (open(OFILE2, ">$ofile2")){
    51. 

  51.          print "Cannot open output file \"$ofile2\"\n\n";
    52. 

  52.          exit;
    53. 

  53. }
    54. 

  54. 

  55. 

  56. # header line for output files
    57. 

  57. #print OFILE "# small indel events, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n";
    58. 

  58. #print OFILE "#align\tingroup1\tingroup1_coord\tingroup1_orient\tingroup2\tingroup2_coord\tingroup2_orient\toutgroup\toutgroup_coord\toutgroup_orient\tindel_type\n";
    59. 

  59. 

  60. #print OFILE2 "# small indels summary, parsed from MAF 3-way alignment file, coords are translated from (-) to (+) if necessary\n";
    61. 

  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";
    62. 

  62. 

  63. # main body of program
    64. 

  64. while ($record = get_next_record($fh) ){
    65. 

  65. 	if ($record=~ m/\s*##maf(.*)\s*# maf/s){
    66. 

  66. 		next;
    67. 

  67. 	}
    68. 

  68. 

  69. 	my @sequences = get_sequences_within_block($record);
    70. 

  70. 	my @seq_info = get_indels_within_block(@sequences);
    71. 

  71. 	get_indels_lengths(@seq_info);
    72. 

  72. 	
    73. 

  73. 	$offset = tell($fh);
    74. 

  74.         $count++;
    75. 

  75.         
    76. 

  76. }
    77. 

  77. 

  78. get_starts_only(@seq1_insert);
    79. 

  79. get_starts_only(@seq1_delete);
    80. 

  80. get_starts_only(@seq2_insert);
    81. 

  81. get_starts_only(@seq2_delete);
    82. 

  82. get_starts_only(@seq3_insert);
    83. 

  83. get_starts_only(@seq3_delete);
    84. 

  84. 

  85. # print some things to keep track of progress
    86. 

  86. #print "# $library\n";
    87. 

  87. #print "# number of records = $count\n";
    88. 

  88. #print "# number of sequence \"s\" lines = $count2\n";
    89. 

  89. if ($count3 != 0){
    90. 

  90. 	print "Skipped $count3 blocks with only 2 seqs;\n";
    91. 

  91. }
    92. 

  92. #print "# number of records with only h-m = $count4\n\n";
    93. 

  93. 

  94. print "Ingroup1 = $ingroup1; Ingroup2 = $ingroup2; Outgroup = $outgroup;\n";
    95. 

  95. print "# of ingroup1 inserts = $count_seq1insert;\n";
    96. 

  96. print "# of ingroup1 deletes = $count_seq1delete;\n";
    97. 

  97. print "# of ingroup2 inserts = $count_seq2insert;\n";
    98. 

  98. print "# of ingroup2 deletes = $count_seq2delete;\n";
    99. 

  99. print "# of outgroup3 inserts = $count_seq3insert;\n";
    100. 

  100. print "# of outgroup3 deletes = $count_seq3delete\n";
    101. 

  101. 

  102. 

  103. #close OFILE;
    104. 

  104. 

  105. if ($count == $count3){
    106. 

  106. 	print STDERR "Skipped all blocks since none of them contain 3-way alignments.\n";
    107. 

  107.   	exit -1;
    108. 

  108. }
    109. 

  109. 

  110. ###################SUBROUTINES#####################################
    111. 

  111. 

  112. # subroutine to open file
    113. 

  113. sub open_file {
    114. 

  114.         my($filename) = @_;
    115. 

  115.         my $fh;
    116. 

  116. 

  117.         unless (open($fh, $filename)){
    118. 

  118.                 print "Cannot open file $filename\n";
    119. 

  119.                 exit;
    120. 

  120.         }
    121. 

  121.         return $fh;
    122. 

  122. }
    123. 

  123. 

  124. # get next record
    125. 

  125. sub get_next_record {
    126. 

  126.         my ($fh) = @_;
    127. 

  127.         my ($offset);
    128. 

  128.         my ($record) = "";
    129. 

  129.         my ($save_input_separator) = $/;
    130. 

  130. 	
    131. 

  131. 	$/ = "a score";
    132. 

  132. 

  133.         $record = <$fh>;
    134. 

  134. 

  135.         $/ = $save_input_separator;
    136. 

  136.         return $record;
    137. 

  137. }
    138. 

  138. 

  139. # get header info
    140. 

  140. sub get_sequences_within_block{
    141. 

  141. 	my (@alignment) = @_;
    142. 

  142. 	my @lines = ();
    143. 

  143. 	
    144. 

  144. 	my @sequences = ();
    145. 

  145. 		
    146. 

  146. 	@lines = split ("\n", $record);
    147. 

  147. 	foreach (@lines){
    148. 

  148. 		chomp($_);
    149. 

  149. 		if (m/^\s*$/){
    150. 

  150. 			next;
    151. 

  151. 		}
    152. 

  152. 		elsif (m/^\s*=(\d+\.*\d*)/){
    153. 

  153. 

  154. 		}elsif (m/^\s*s(.*)$/){
    155. 

  155. 			$count2++;
    156. 

  156. 			
    157. 

  157. 			push (@sequences,$_);
    158. 

  158. 		}
    159. 

  159. 	}
    160. 

  160. 	return @sequences;
    161. 

  161. }
    162. 

  162. 	
    163. 

  163. sub get_indels_within_block{
    164. 

  164. 	my (@sequences) = @_;	
    165. 

  165. 	my $line1 = my $line2 = my $line3 = "";
    166. 

  166.     my @line1 = my @line2 = my @line3 = ();
    167. 

  167. 	my $score = 0;
    168. 

  168.         my $start1 = my $align_length1 = my $end1 = my $seq_length1 = 0;
    169. 

  169.         my $start2 = my $align_length2 = my $end2 = my $seq_length2 = 0;
    170. 

  170.         my $start3 = my $align_length3 = my $end3 = my $seq_length3 = 0;
    171. 

  171.         my $seq1 = my $orient1 = "";
    172. 

  172.         my $seq2 = my $orient2 = "";
    173. 

  173.         my $seq3 = my $orient3 = "";
    174. 

  174.         my $start1_plus = my $end1_plus = 0;
    175. 

  175.         my $start2_plus = my $end2_plus = 0;
    176. 

  176.         my $start3_plus = my $end3_plus = 0;
    177. 

  177.    	my @test = ();
    178. 

  178.         my $test = "";
    179. 

  179.         my $header = "";
    180. 

  180.         my @header = ();
    181. 

  181.         my $sequence1 = my $sequence2 = my $sequence3 ="";
    182. 

  182. 	my @array_return = ();	
    183. 

  183. 	my $test1 = 0;
    184. 

  184. 	my $line1_stat = my $line2_stat = my $line3_stat = "";
    185. 

  185.     
    186. 

  186. 	# process 3-way blocks only
    187. 

  187. 	if (scalar(@sequences) == 3){
    188. 

  188. 		$line1 = $sequences[0]; 
    189. 

  189. 		chomp ($line1);
    190. 

  190. 		$line2 = $sequences[1]; 
    191. 

  191. 		chomp ($line2);
    192. 

  192. 		$line3 = $sequences[2];	
    193. 

  193. 		chomp ($line3);
    194. 

  194. 		# check order of sequences and assign uniformly seq1= human, seq2 = chimp, seq3 = macaque
    195. 

  195. 		if ($line1 =~ m/$outgroup/){
    196. 

  196. 			$line1_stat = "out";
    197. 

  197. 			$line2=~ s/^\s*//;
    198. 

  198.             $line2 =~ s/\s+/\t/g;
    199. 

  199.             @line2 = split(/\t/, $line2);
    200. 

  200. 			if (($ingroup1 eq "") || ($line2[1] =~ m/$ingroup1/)){
    201. 

  201. 			 $line2_stat = "in1";
    202. 

  202. 			 $line3_stat = "in2";
    203. 

  203. 			 }
    204. 

  204. 			 else{
    205. 

  205.              $line3_stat = "in1";
    206. 

  206.              $line2_stat = "in2";              }
    207. 

  207. 			}
    208. 

  208. 		elsif ($line2 =~ m/$outgroup/){
    209. 

  209. 			$line2_stat = "out";
    210. 

  210. 			$line1=~ s/^\s*//;
    211. 

  211.             $line1 =~ s/\s+/\t/g;
    212. 

  212.             @line1 = split(/\t/, $line1);
    213. 

  213.             if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){
    214. 

  214.              $line1_stat = "in1";
    215. 

  215.              $line3_stat = "in2";
    216. 

  216.              }
    217. 

  217.              else{
    218. 

  218.              $line3_stat = "in1";
    219. 

  219.              $line1_stat = "in2";              }
    220. 

  220.             }
    221. 

  221. 		elsif ($line3 =~ m/$outgroup/){
    222. 

  222. 			$line3_stat = "out";
    223. 

  223. 			$line1=~ s/^\s*//;
    224. 

  224.             $line1 =~ s/\s+/\t/g;
    225. 

  225.             @line1 = split(/\t/, $line1);
    226. 

  226.             if (($ingroup1 eq "") || ($line1[1] =~ m/$ingroup1/)){
    227. 

  227.              $line1_stat = "in1";
    228. 

  228.              $line2_stat = "in2";
    229. 

  229.              }
    230. 

  230.              else{
    231. 

  231.              $line2_stat = "in1";
    232. 

  232.              $line1_stat = "in2";              }
    233. 

  233. 			}
    234. 

  234. 

  235. 		#print "# l1 = $line1_stat\n";
    236. 

  236. 		#print "# l2 = $line2_stat\n";
    237. 

  237. 		#print "# l3 = $line3_stat\n";
    238. 

  238. 		my $linei1 = my $linei2 = my $lineo = "";
    239. 

  239. 		my @linei1 = my @linei2 = my @lineo = ();
    240. 

  240. 		
    241. 

  241.         if ($line1_stat eq "out"){
    242. 

  242.             $lineo = $line1;
    243. 

  243.         }
    244. 

  244.         elsif ($line1_stat eq "in1"){
    245. 

  245.             $linei1 = $line1;
    246. 

  246.         }
    247. 

  247.         else{
    248. 

  248.             $linei2 = $line1;
    249. 

  249.         }
    250. 

  250.         
    251. 

  251.         if ($line2_stat eq "out"){
    252. 

  252.             $lineo = $line2;
    253. 

  253.         }
    254. 

  254.         elsif ($line2_stat eq "in1"){
    255. 

  255.             $linei1 = $line2;
    256. 

  256.         }
    257. 

  257.         else{
    258. 

  258.             $linei2 = $line2;
    259. 

  259.         }
    260. 

  260.         
    261. 

  261.         if ($line3_stat eq "out"){
    262. 

  262.             $lineo = $line3;
    263. 

  263.         }
    264. 

  264.         elsif ($line3_stat eq "in1"){
    265. 

  265.             $linei1 = $line3;
    266. 

  266.         }
    267. 

  267.         else{
    268. 

  268.             $linei2 = $line3;
    269. 

  269.         }
    270. 

  270.         
    271. 

  271.         $linei1=~ s/^\s*//;
    272. 

  272.         $linei1 =~ s/\s+/\t/g;
    273. 

  273.         @linei1 = split(/\t/, $linei1);
    274. 

  274.         $end1 =($linei1[2]+$linei1[3]-1);
    275. 

  275.         $seq1 = $linei1[1].":".$linei1[3];
    276. 

  276.         $ingroup1 = (split(/\./, $seq1))[0];
    277. 

  277.         $start1 = $linei1[2];
    278. 

  278.         $align_length1 = $linei1[3];
    279. 

  279.         $orient1 = $linei1[4];
    280. 

  280.         $seq_length1 = $linei1[5];
    281. 

  281.         $sequence1 = $linei1[6];
    282. 

  282.         $test1 = length($sequence1);
    283. 

  283.         my $total_length1 = $test1+$start1;
    284. 

  284.         my @array1 = ($start1,$end1,$orient1,$seq_length1);
    285. 

  285.         ($start1_plus, $end1_plus) =  convert_coords(@array1);
    286. 

  286.            
    287. 

  287.         $linei2=~ s/^\s*//;
    288. 

  288.         $linei2 =~ s/\s+/\t/g;
    289. 

  289.         @linei2 = split(/\t/, $linei2);
    290. 

  290.         $end2 =($linei2[2]+$linei2[3]-1);               
    291. 

  291.         $seq2 = $linei2[1].":".$linei2[3];
    292. 

  292.         $ingroup2 = (split(/\./, $seq2))[0];
    293. 

  293.         $start2 = $linei2[2];
    294. 

  294.         $align_length2 = $linei2[3];
    295. 

  295.         $orient2 = $linei2[4];
    296. 

  296.         $seq_length2 = $linei2[5];
    297. 

  297.         $sequence2 = $linei2[6];
    298. 

  298.         my $test2 = length($sequence2);
    299. 

  299.         my $total_length2 = $test2+$start2;
    300. 

  300.         my @array2 = ($start2,$end2,$orient2,$seq_length2);
    301. 

  301.         ($start2_plus, $end2_plus) = convert_coords(@array2); 
    302. 

  302.             
    303. 

  303.         $lineo=~ s/^\s*//;
    304. 

  304.         $lineo =~ s/\s+/\t/g;
    305. 

  305.         @lineo = split(/\t/, $lineo);
    306. 

  306.         $end3 =($lineo[2]+$lineo[3]-1);
    307. 

  307.         $seq3 = $lineo[1].":".$lineo[3];
    308. 

  308.         $start3 = $lineo[2];
    309. 

  309.         $align_length3 = $lineo[3];
    310. 

  310.         $orient3 = $lineo[4];
    311. 

  311.         $seq_length3 = $lineo[5];
    312. 

  312.         $sequence3 = $lineo[6];
    313. 

  313.         my $test3 = length($sequence3);
    314. 

  314.         my $total_length3 = $test3+$start3;
    315. 

  315.         my @array3 = ($start3,$end3,$orient3,$seq_length3);
    316. 

  316.         ($start3_plus, $end3_plus) = convert_coords(@array3);
    317. 

  317.         
    318. 

  318.         #print "# l1 = $ingroup1\n";
    319. 

  319. 		#print "# l2 = $ingroup2\n";
    320. 

  320. 		#print "# l3 = $outgroup\n";
    321. 

  321. 	
    322. 

  322. 		my $ABC = "";
    323. 

  323. 		my $coord1 = my $coord2 = my $coord3 = 0;
    324. 

  324.                 $coord1 = $start1_plus;
    325. 

  325.                 $coord2 = $start2_plus;
    326. 

  326.                 $coord3 = $start3_plus;
    327. 

  327. 		
    328. 

  328. 		for (my $position = 0; $position < $test1; $position++) {
    329. 

  329. 			my $indelType = "";
    330. 

  330. 			my $indel_line = "";	
    331. 

  331. 			# seq1 deletes
    332. 

  332. 			 if ((substr($sequence1,$position,1) eq "-") 
    333. 

  333. 		    		&& (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
    334. 

  334. 	       			&& (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
    335. 

  335. 					$ABC = join("",($ABC,"X"));
    336. 

  336. 					my @s = split(/:/, $seq1);
    337. 

  337. 					$indelType = $s[0]."_delete";
    338. 

  338. 

  339. 					#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";	
    340. 

  340. 					$indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
    341. 

  341. 					push (@indels,$indel_line);
    342. 

  342. 					push (@seq1_delete,$indel_line);
    343. 

  343. 		 			$coord2++; $coord3++;
    344. 

  344. 	       		}	
    345. 

  345. 			# seq2 deletes
    346. 

  346. 			elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
    347. 

  347. 				&& (substr($sequence2,$position,1) eq "-")
    348. 

  348. 				&& (substr($sequence3,$position,1) !~ m/[-*\$?^]/)){
    349. 

  349. 					$ABC = join("",($ABC,"Y"));
    350. 

  350. 					my @s = split(/:/, $seq2);
    351. 

  351. 					$indelType = $s[0]."_delete";
    352. 

  352. 					#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
    353. 

  353. 					$indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
    354. 

  354.                                         push (@indels,$indel_line);
    355. 

  355. 					push (@seq2_delete,$indel_line);
    356. 

  356. 					$coord1++;
    357. 

  357. 					$coord3++;
    358. 

  358. 

  359. 			}
    360. 

  360. 			# seq1 inserts
    361. 

  361. 			elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
    362. 

  362. 				&& (substr($sequence2,$position,1) eq "-")
    363. 

  363. 				&& (substr($sequence3,$position,1) eq "-")){
    364. 

  364. 					$ABC = join("",($ABC,"Z"));
    365. 

  365. 					my @s = split(/:/, $seq1);
    366. 

  366. 					$indelType = $s[0]."_insert";
    367. 

  367. 					#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
    368. 

  368. 					$indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
    369. 

  369. 					push (@indels,$indel_line);
    370. 

  370. 					push (@seq1_insert,$indel_line);
    371. 

  371. 					$coord1++;
    372. 

  372. 			}
    373. 

  373. 			# seq2 inserts	
    374. 

  374. 			elsif ((substr($sequence1,$position,1) eq "-")
    375. 

  375. 				&& (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
    376. 

  376. 				&& (substr($sequence3,$position,1) eq "-")){
    377. 

  377. 					$ABC = join("",($ABC,"W"));
    378. 

  378. 					my @s = split(/:/, $seq2);
    379. 

  379. 					$indelType = $s[0]."_insert";
    380. 

  380. 					#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
    381. 

  381. 					$indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
    382. 

  382. 					push (@indels,$indel_line);
    383. 

  383. 					push (@seq2_insert,$indel_line);
    384. 

  384. 					$coord2++;
    385. 

  385. 			}
    386. 

  386. 			# seq3 deletes
    387. 

  387. 			elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
    388. 

  388. 				&& (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
    389. 

  389. 				&& (substr($sequence3,$position,1) eq "-")){
    390. 

  390. 					$ABC = join("",($ABC,"S"));
    391. 

  391. 					my @s = split(/:/, $seq3);
    392. 

  392. 					$indelType = $s[0]."_delete";
    393. 

  393. 					#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
    394. 

  394. 					$indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
    395. 

  395. 					push (@indels,$indel_line);
    396. 

  396. 					push (@seq3_delete,$indel_line);
    397. 

  397. 					$coord1++; $coord2++;
    398. 

  398. 			}
    399. 

  399. 			# seq3 inserts
    400. 

  400. 			elsif ((substr($sequence1,$position,1) eq "-")
    401. 

  401. 				&& (substr($sequence2,$position,1) eq "-")
    402. 

  402. 				&& (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
    403. 

  403. 					$ABC = join("",($ABC,"T"));
    404. 

  404. 					my @s = split(/:/, $seq3);
    405. 

  405. 					$indelType = $s[0]."_insert";
    406. 

  406. 					#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
    407. 

  407. 					$indel_line = join("\t",($count,$seq1,$coord1,$orient1,$seq2,$coord2,$orient2,$seq3,$coord3,$orient3,$indelType));
    408. 

  408. 					push (@indels,$indel_line);
    409. 

  409. 					push (@seq3_insert,$indel_line);
    410. 

  410. 					$coord3++;
    411. 

  411. 			}else{
    412. 

  412. 				$ABC = join("",($ABC,"N"));
    413. 

  413. 				$coord1++; $coord2++; $coord3++;
    414. 

  414. 			}
    415. 

  415. 

  416. 		}
    417. 

  417. 		@array_return=($seq1,$seq2,$seq3,$ABC);
    418. 

  418. 		return (@array_return); 
    419. 

  419. 

  420. 	}
    421. 

  421. 	# ignore pairwise cases for now, just count the number of blocks
    422. 

  422. 	elsif (scalar(@sequences) == 2){
    423. 

  423. 		my $ABC = "";
    424. 

  424. 		my $coord1 = my $coord2 = my $coord3 = 0;
    425. 

  425. 		$count3++;
    426. 

  426. 

  427. 		$line1 = $sequences[0];
    428. 

  428. 		$line2 = $sequences[1];
    429. 

  429. 		chomp ($line1);
    430. 

  430. 		chomp ($line2);
    431. 

  431. 

  432. 		if ($line2 !~ m/$ingroup2/){
    433. 

  433. 		       $count4++;
    434. 

  434. 		}
    435. 

  435. 	}
    436. 

  436. }
    437. 

  437. 

  438. 		
    439. 

  439. sub get_indels_lengths{
    440. 

  440. 	my (@array) = @_;
    441. 

  441. 	if (scalar(@array) == 4){
    442. 

  442. 		my $seq1 = $array[0]; my $seq2 = $array[1]; my $seq3 = $array[2]; my $ABC = $array[3];
    443. 

  443. 

  444. 		while ($ABC =~ m/(X+)/g) {
    445. 

  445. 			push (@seq1_delete_lengths,length($1));
    446. 

  446. 			$count_seq1delete++;
    447. 

  447. 		}
    448. 

  448. 		
    449. 

  449. 		while ($ABC =~ m/(Y+)/g) {
    450. 

  450. 			push (@seq2_delete_lengths,length($1));
    451. 

  451. 			$count_seq2delete++;
    452. 

  452. 		}
    453. 

  453. 		while ($ABC =~ m/(S+)/g) {
    454. 

  454. 			push (@seq3_delete_lengths,length($1));
    455. 

  455. 			$count_seq3delete++;
    456. 

  456. 		}
    457. 

  457. 		while ($ABC =~ m/(Z+)/g) {	
    458. 

  458. 			push (@seq1_insert_lengths,length($1));
    459. 

  459. 			$count_seq1insert++;
    460. 

  460. 		}
    461. 

  461. 		while ($ABC =~ m/(W+)/g) {
    462. 

  462. 			push(@seq2_insert_lengths,length($1));
    463. 

  463. 			$count_seq2insert++;
    464. 

  464. 		}
    465. 

  465. 		while ($ABC =~ m/(T+)/g) {
    466. 

  466. 			push (@seq3_insert_lengths,length($1));
    467. 

  467. 			$count_seq3insert++;
    468. 

  468. 		}
    469. 

  469. 	}
    470. 

  470. 	elsif (scalar(@array) == 0){
    471. 

  471. 		next;
    472. 

  472. 	}
    473. 

  473. 		
    474. 

  474. }
    475. 

  475. # convert to coordinates to + strand if align orient = -
    476. 

  476. sub convert_coords{
    477. 

  477. 	my (@array) = @_;
    478. 

  478. 	my $s = $array[0];
    479. 

  479. 	my $e = $array[1];
    480. 

  480. 	my $o = $array[2];
    481. 

  481. 	my $l = $array[3];
    482. 

  482. 	my $start_plus = my $end_plus = 0;
    483. 

  483. 

  484. 	if ($o eq "-"){
    485. 

  485. 		$start_plus = ($l - $e);
    486. 

  486. 		$end_plus = ($l - $s);
    487. 

  487. 	}elsif ($o eq "+"){
    488. 

  488. 		$start_plus = $s;
    489. 

  489. 		$end_plus = $e;
    490. 

  490. 	}
    491. 

  491. 

  492. 	return ($start_plus, $end_plus);
    493. 

  493. }
    494. 

  494. 

  495. # get first line only for each event
    496. 

  496. sub get_starts_only{
    497. 

  497. 	my (@test) = @_;
    498. 

  498.         my $seq1 = my $seq2 = my $seq3 = my $indelType = my $old_seq1 = my $old_seq2 = my $old_seq3 = my $old_indelType = my $old_line = "";
    499. 

  499.         my $coord1 = my $coord2 = my $coord3 = my $old_coord1 = my $old_coord2 = my $old_coord3 = 0;
    500. 

  500. 

  501.         my @matches = ();
    502. 

  502.         my @seq1_insert = my @seq1_delete = my @seq2_insert = my @seq2_delete = my @seq3_insert = my @seq3_delete = ();
    503. 

  503. 				
    504. 

  504. 	
    505. 

  505.        	foreach my $line (@test){
    506. 

  506.                 chomp($line);
    507. 

  507.                 $line =~ s/^\s*//;
    508. 

  508.                 $line =~ s/\s+/\t/g;
    509. 

  509. 				my @line1 = split(/\t/, $line);
    510. 

  510.                  $seq1 = $line1[1];
    511. 

  511.                  $coord1 = $line1[2];
    512. 

  512.                 $seq2 = $line1[4];
    513. 

  513.                  $coord2 = $line1[5];
    514. 

  514.                  $seq3 = $line1[7];
    515. 

  515.                  $coord3 = $line1[8];
    516. 

  516.                  $indelType = $line1[10];
    517. 

  517.                 if ($indelType =~ m/$ingroup1/ && $indelType =~ m/insert/){
    518. 

  518.            		if ($coord1 != $old_coord1+1 || ($coord2 != $old_coord2 || $coord3 != $old_coord3)){
    519. 

  519. 	       			$start1++;
    520. 

  520.                               	push (@seq1_insert_startOnly,$line);
    521. 

  521.                      	}
    522. 

  522. 	     	}
    523. 

  523.                 elsif ($indelType =~ m/$ingroup1/ && $indelType =~ m/delete/){
    524. 

  524. 		        if ($coord1 != $old_coord1 || ($coord2 != $old_coord2+1 || $coord3 != $old_coord3+1)){
    525. 

  525. 		                $start2++;
    526. 

  526. 		                push(@seq1_delete_startOnly,$line);
    527. 

  527. 		        }
    528. 

  528. 		}
    529. 

  529.                 elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/insert/){
    530. 

  530. 	                if ($coord2 != $old_coord2+1 || ($coord1 != $old_coord1 || $coord3 != $old_coord3)){
    531. 

  531. 		                $start3++;
    532. 

  532. 		                push(@seq2_insert_startOnly,$line);
    533. 

  533. 		        }
    534. 

  534. 		}
    535. 

  535.                 elsif ($indelType =~ m/$ingroup2/ && $indelType =~ m/delete/){
    536. 

  536.                         if ($coord2 != $old_coord2 || ($coord1 != $old_coord1+1 || $coord3 != $old_coord3+1)){
    537. 

  537.                                 $start4++;
    538. 

  538.                                 push(@seq2_delete_startOnly,$line);
    539. 

  539.                         }
    540. 

  540.                 }
    541. 

  541.                 elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/insert/){
    542. 

  542.                         if ($coord3 != $old_coord3+1 || ($coord1 != $old_coord1 || $coord2 != $old_coord2)){
    543. 

  543.                                 $start5++;
    544. 

  544.                                 push(@seq3_insert_startOnly,$line);
    545. 

  545.                         }
    546. 

  546.                 }
    547. 

  547.                 elsif ($indelType =~ m/$outgroup/ && $indelType =~ m/delete/){
    548. 

  548.                         if ($coord3 != $old_coord3 || ($coord1 != $old_coord1+1 ||$coord2 != $old_coord2+1)){
    549. 

  549.                                 $start6++;
    550. 

  550.                                 push(@seq3_delete_startOnly,$line);
    551. 

  551.                         }
    552. 

  552.                 }
    553. 

  553.                  $old_indelType = $indelType;
    554. 

  554.                  $old_seq1 = $seq1;
    555. 

  555.                  $old_coord1 = $coord1;
    556. 

  556.                  $old_seq2 = $seq2;
    557. 

  557.                  $old_coord2 = $coord2;
    558. 

  558.                  $old_seq3 = $seq3;
    559. 

  559.                  $old_coord3 = $coord3;
    560. 

  560.                  $old_line = $line;
    561. 

  561.         }
    562. 

  562. }
    563. 

  563. # append lengths to each event start line
    564. 

  564. my $counter1; my $counter2; my $counter3; my $counter4; my $counter5; my $counter6; 
    565. 

  565. my @final1 = my @final2 = my @final3 = my @final4 = my @final5 = my @final6 = ();
    566. 

  566. my $final_line1 = my $final_line2 = my $final_line3 = my $final_line4 = my $final_line5 = my $final_line6 = "";
    567. 

  567. 

  568. 

  569. for ($counter1 = 0; $counter1 < @seq3_insert_startOnly; $counter1++){
    570. 

  570. 	$final_line1 = join("\t",($seq3_insert_startOnly[$counter1],$seq3_insert_lengths[$counter1]));
    571. 

  571. 	push (@final1,$final_line1);
    572. 

  572. }
    573. 

  573. 

  574. for ($counter2 = 0; $counter2 < @seq3_delete_startOnly; $counter2++){
    575. 

  575.         $final_line2 =  join("\t",($seq3_delete_startOnly[$counter2],$seq3_delete_lengths[$counter2]));
    576. 

  576.         push(@final2,$final_line2);
    577. 

  577. }
    578. 

  578. 

  579. for ($counter3 = 0; $counter3 < @seq2_insert_startOnly; $counter3++){
    580. 

  580.     $final_line3 =  join("\t",($seq2_insert_startOnly[$counter3],$seq2_insert_lengths[$counter3]));
    581. 

  581. 	push(@final3,$final_line3);
    582. 

  582. }
    583. 

  583. 

  584. for ($counter4 = 0; $counter4 < @seq2_delete_startOnly; $counter4++){
    585. 

  585.         $final_line4 =  join("\t",($seq2_delete_startOnly[$counter4],$seq2_delete_lengths[$counter4]));
    586. 

  586.         push(@final4,$final_line4);
    587. 

  587. }
    588. 

  588. 		
    589. 

  589. for ($counter5 = 0; $counter5 < @seq1_insert_startOnly; $counter5++){
    590. 

  590.         $final_line5 =  join("\t",($seq1_insert_startOnly[$counter5],$seq1_insert_lengths[$counter5]));
    591. 

  591.         push(@final5,$final_line5);
    592. 

  592. }
    593. 

  593. 

  594. for ($counter6 = 0; $counter6 < @seq1_delete_startOnly; $counter6++){
    595. 

  595.         $final_line6 =  join("\t",($seq1_delete_startOnly[$counter6],$seq1_delete_lengths[$counter6]));
    596. 

  596.         push(@final6,$final_line6);
    597. 

  597. }       
    598. 

  598. 

  599. # format final output
    600. 

  600. # # if inserts, increase coords for the sequence inserted, other sequences give coords for 5' and 3' bases flanking the gap
    601. 

  601. # # for deletes, increase coords for other 2 sequences and the one deleted give coords for 5' and 3' bases flanking the gap
    602. 

  602. 

  603. get_final_format(@final5);
    604. 

  604. get_final_format(@final6);
    605. 

  605. get_final_format(@final3);
    606. 

  606. get_final_format(@final4);
    607. 

  607. get_final_format(@final1);
    608. 

  608. get_final_format(@final2);
    609. 

  609. 

  610. sub get_final_format{
    611. 

  611. 	my (@final) = @_;
    612. 

  612. 	foreach (@final){
    613. 

  613. 		my $event_line = $_;
    614. 

  614. 		my @events = split(/\t/, $event_line);
    615. 

  615. 		my $event_type = $events[10];
    616. 

  616. 		my @name_align1 = split(/:/, $events[1]);
    617. 

  617. 		my @name_align2 = split(/:/, $events[4]);
    618. 

  618. 		my @name_align3 = split(/:/, $events[7]);
    619. 

  619. 		my $seq1_event_start = my $seq1_event_end = my $seq2_event_start = my $seq2_event_end = my $seq3_event_start = my $seq3_event_end = 0;
    620. 

  620. 		my $final_event_line = "";	
    621. 

  621. 		# seq1_insert
    622. 

  622. 		if ($event_type =~ m/$ingroup1/ && $event_type =~ m/insert/){
    623. 

  623. 			# only increase coord for human
    624. 

  624. 			# remember that other two sequnences, the gap spans (coord - 1) --> coord
    625. 

  625. 			$seq1_event_start = ($events[2]);
    626. 

  626. 			$seq1_event_end = ($events[2]+$events[11]-1);
    627. 

  627. 			$seq2_event_start = ($events[5]-1);
    628. 

  628. 			$seq2_event_end = ($events[5]);
    629. 

  629. 			$seq3_event_start = ($events[8]-1);
    630. 

  630. 			$seq3_event_end = ($events[8]);
    631. 

  631. 			$final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
    632. 

  632. 		}
    633. 

  633. 		# seq1_delete
    634. 

  634. 		elsif ($event_type =~ m/$ingroup1/ && $event_type =~ m/delete/){
    635. 

  635. 			# only increase coords for seq2 and seq3
    636. 

  636. 			# remember for seq1, the gap spans (coord - 1) --> coord
    637. 

  637. 			$seq1_event_start = ($events[2]-1);
    638. 

  638. 			$seq1_event_end = ($events[2]);
    639. 

  639.                         $seq2_event_start = ($events[5]);
    640. 

  640.                         $seq2_event_end = ($events[5]+$events[11]-1);
    641. 

  641.                         $seq3_event_start = ($events[8]);
    642. 

  642.                         $seq3_event_end = ($events[8]+$events[11]-1);
    643. 

  643. 			$final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
    644. 

  644. 		}
    645. 

  645. 		# seq2_insert
    646. 

  646. 		elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/insert/){	
    647. 

  647. 			# only increase coords for seq2 
    648. 

  648. 			# remember that other two sequnences, the gap spans (coord - 1) --> coord
    649. 

  649.                         $seq1_event_start = ($events[2]-1);
    650. 

  650.                         $seq1_event_end = ($events[2]);
    651. 

  651. 			$seq2_event_start = ($events[5]);
    652. 

  652.                         $seq2_event_end = ($events[5]+$events[11]-1);
    653. 

  653.                         $seq3_event_start = ($events[8]-1);
    654. 

  654. 			$seq3_event_end = ($events[8]);			
    655. 

  655. 			$final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
    656. 

  656. 		}
    657. 

  657. 		# seq2_delete
    658. 

  658. 		elsif ($event_type =~ m/$ingroup2/ && $event_type =~ m/delete/){
    659. 

  659. 			# only increase coords for seq1 and seq3
    660. 

  660. 			# remember for seq2, the gap spans (coord - 1) --> coord
    661. 

  661.                         $seq1_event_start = ($events[2]);
    662. 

  662. 			$seq1_event_end = ($events[2]+$events[11]-1);	
    663. 

  663.                         $seq2_event_start = ($events[5]-1);
    664. 

  664. 	                $seq2_event_end = ($events[5]);
    665. 

  665.                         $seq3_event_start = ($events[8]);
    666. 

  666.                         $seq3_event_end = ($events[8]+$events[11]-1);
    667. 

  667. 			$final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
    668. 

  668. 		}	
    669. 

  669. 		# start testing w/seq3_insert
    670. 

  670. 		elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/insert/){
    671. 

  671. 			# only increase coord for rheMac
    672. 

  672. 			# remember that other two sequnences, the gap spans (coord - 1) --> coord
    673. 

  673. 			$seq1_event_start = ($events[2]-1);
    674. 

  674. 			$seq1_event_end = ($events[2]);
    675. 

  675. 			$seq2_event_start = ($events[5]-1);
    676. 

  676. 			$seq2_event_end = ($events[5]);
    677. 

  677. 			$seq3_event_start = ($events[8]);
    678. 

  678. 			$seq3_event_end = ($events[8]+$events[11]-1);
    679. 

  679. 			$final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
    680. 

  680. 		}
    681. 

  681. 		# seq3_delete
    682. 

  682. 		elsif ($event_type =~ m/$outgroup/ && $event_type =~ m/delete/){
    683. 

  683. 			# only increase coords for seq1 and seq2
    684. 

  684. 			# remember for seq3, the gap spans (coord - 1) --> coord
    685. 

  685. 			$seq1_event_start = ($events[2]);
    686. 

  686. 			$seq1_event_end = ($events[2]+$events[11]-1);
    687. 

  687. 			$seq2_event_start = ($events[5]);
    688. 

  688. 			$seq2_event_end = ($events[5]+$events[11]-1);
    689. 

  689. 			$seq3_event_start = ($events[8]-1);
    690. 

  690. 			$seq3_event_end = ($events[8]);
    691. 

  691. 			$final_event_line = join("\t",($events[0],$event_type,$events[11],$name_align1[0],$seq1_event_start,$seq1_event_end,$name_align1[1],$events[3],$name_align2[0],$seq2_event_start,$seq2_event_end,$name_align2[1],$events[6],$name_align3[0],$seq3_event_start,$seq3_event_end,$name_align3[1],$events[9]));
    692. 

  692. 

  693. 		}
    694. 

  694. 		
    695. 

  695. 		print OFILE2 "$final_event_line\n";
    696. 

  696. 	}
    697. 

  697. }
    698. 

  698. close OFILE2;
    699. 

  699.