PageRenderTime 71ms CodeModel.GetById 27ms RepoModel.GetById 0ms app.codeStats 0ms

/tools/regVariation/multispecies_MicrosatDataGenerator_interrupted_GALAXY.pl

https://bitbucket.org/chapmanb/galaxy-central/
Perl | 5606 lines | 5065 code | 245 blank | 296 comment | 163 complexity | e64ba560e907a6b8d7818123ba23613a MD5 | raw file

Large files files are truncated, but you can click here to view the full file

    

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

  2. use strict;
    3. 

  3. use warnings;
    4. 

  4. use Term::ANSIColor;
    5. 

  5. use File::Basename;
    6. 

  6. use IO::Handle;
    7. 

  7. use Cwd;
    8. 

  8. use File::Path;
    9. 

  9. use vars qw($distance @thresholds @tags $species_set @allspecies $printer $treeSpeciesNum $focalspec $mergestarts $mergeends $mergemicros $interrtypecord $microscanned $interrcord $interr_poscord $no_of_interruptionscord $infocord $typecord $startcord $strandcord $endcord $microsatcord $motifcord $sequencepos $no_of_species $gapcord $prinkter);
    10. 

  10. use File::Path qw(make_path remove_tree);
    11. 

  11. use File::Temp qw/ tempfile tempdir /;
    12. 

  12. my $tdir = tempdir( CLEANUP => 1 );
    13. 

  13. chdir $tdir;
    14. 

  14. my $dir = getcwd;  
    15. 

  15. #print "dir = $dir\n";
    16. 

  16. 

  17. #$ENV{'PATH'} .= ':' . dirname($0);
    18. 

  18. my $date = `date`;
    19. 

  19. 

  20. my ($mafile, $mafile_sputt, $orthfile, $threshold_array,  $allspeciesin, $tree_definition_all, $separation) = @ARGV;
    21. 

  21. if (!$mafile or !$mafile_sputt or !$orthfile or !$threshold_array or !$separation or !$tree_definition_all or !$allspeciesin) { die "missing arguments\n"; }
    22. 

  22. 

  23. $tree_definition_all =~ s/\s+//g;
    24. 

  24. $threshold_array =~ s/\s+//g;
    25. 

  25. $allspeciesin =~ s/\s+//g;
    26. 

  26. #-------------------------------------------------------------------------------
    27. 

  27. # WHICH SPUTNIK USED?
    28. 

  28. my $sputnikpath = ();
    29. 

  29. $sputnikpath = "sputnik_lowthresh_MATCH_MIN_SCORE3" ;
    30. 

  30. #$sputnikpath = "/Users/ydk/work/rhesus_microsat/codes/./sputnik_Mac-PowerPC";
    31. 

  31. #print "sputnik_Mac-PowerPC non-existant\n" if !-e $sputnikpath;
    32. 

  32. #exit if !-e $sputnikpath;
    33. 

  33. #$sputnikpath = "bx-sputnik" ;
    34. 

  34. #print "ARGV input = @ARGV\n";
    35. 

  35. #print "ARGV input :\n mafile=$mafile\n orthfile=$orthfile\n threshold_array=$threshold_array\n  species_set=$species_set\n tree_definition=$tree_definition\n separation=$separation\n";
    36. 

  36. #-------------------------------------------------------------------------------
    37. 

  37. # RUNFILE
    38. 

  38. #-------------------------------------------------------------------------------
    39. 

  39. $distance = 1; #bp
    40. 

  40. $distance++;
    41. 

  41. my @tree_definitions=MakeTrees($tree_definition_all);
    42. 

  42. my $allspeciesset = $tree_definition_all;
    43. 

  43. $allspeciesset =~ s/[\(\) ]+//g;
    44. 

  44. @allspecies = split(/,/,$allspeciesset);
    45. 

  45. 

  46. my @outputfiles = ();
    47. 

  47. my $round = 0;
    48. 

  48. #my $tdir = tempdir( CLEANUP => 0 );
    49. 

  49. #chdir $tdir;
    50. 

  50. 

  51. foreach my $tree_definition (@tree_definitions){
    52. 

  52. 	my @commas = ($tree_definition =~ /,/g) ;
    53. 

  53. 	#print "commas = @commas\n"; <STDIN>;
    54. 

  54. 	next if scalar(@commas) <= 1;
    55. 

  55. 	#print "species_set = $species_set\n";
    56. 

  56. 	$treeSpeciesNum = scalar(@commas) + 1;
    57. 

  57. 	$species_set = $tree_definition;
    58. 

  58. 	$species_set =~ s/[\)\( ;]+//g;
    59. 

  59. 	#print "species_set = $species_set\n"; <STDIN>;
    60. 

  60. 

  61. 	$round++;
    62. 

  62. 	#-------------------------------------------------------------------------------
    63. 

  63. 	# MICROSATELLITE THRESHOLD SETTINGS (LENGTH, BP)
    64. 

  64. 	$threshold_array=~ s/,/_/g;
    65. 

  65. 	my @thresharr = split("_",$threshold_array);
    66. 

  66. 	@thresholds=@thresharr;
    67. 

  67. 	#my $threshold_array = join("_",($mono_threshold, $di_threshold, $tri_threshold, $tetra_threshold));
    68. 

  68. 	#print "current dit=$dir\n";
    69. 

  69. 	#-------------------------------------------------------------------------------
    70. 

  70. 	# CREATE AXT FILES IN FORWARD AND REVERSE ORDERS IF NECESSARY
    71. 

  71. 	my @chrfiles=();
    72. 

  72. 	
    73. 

  73. 	#my $mafile =  "/Users/ydk/work/rhesus_microsat/results/galay/align.txt"; #$ARGV[0];
    74. 

  74. 	my $chromt=int(rand(10000));
    75. 

  75. 	my $p_chr=$chromt;
    76. 

  76. 	
    77. 

  77. 	
    78. 

  78. 	my @exactspeciesset_unarranged = split(/,/,$species_set);
    79. 

  79. 	$tree_definition=~s/[\)\(, ]/\t/g;
    80. 

  80. 	my @treespecies=split(/\t+/,$tree_definition);
    81. 

  81. 	my @exactspecies=();
    82. 

  82. 	
    83. 

  83. 	foreach my $spec (@treespecies){
    84. 

  84. 		foreach my $espec (@exactspeciesset_unarranged){
    85. 

  85. 			push @exactspecies, $spec if $spec eq $espec;
    86. 

  86. 		}
    87. 

  87. 	}
    88. 

  88. 	#print "exactspecies=@exactspecies\n";
    89. 

  89. 	$focalspec = $exactspecies[0];
    90. 

  90. 	my $arranged_species_set=join(".",@exactspecies);
    91. 

  91. 	my $chr_name = join(".",("chr".$p_chr),$arranged_species_set, "net", "axt");
    92. 

  92. 	my $chr_name_sputt = join(".",("chr".$p_chr),$arranged_species_set, "net", "axt_sputt");
    93. 

  93. 	#print "sending to maftoAxt_multispecies: $mafile, $tree_definition, $chr_name, $species_set .. focalspec=$focalspec \n"; 
    94. 

  94. 	maftoAxt_multispecies($mafile, $tree_definition, $chr_name, $species_set);
    95. 

  95. 	maftoAxt_multispecies($mafile_sputt, $tree_definition, $chr_name_sputt, $species_set);
    96. 

  96. 	#print "done maf to axt conversion\n";
    97. 

  97. 	my $reverse_chr_name = join(".",("chr".$p_chr."r"),$arranged_species_set, "net", "axt");
    98. 

  98. 	artificial_axdata_inverter ($chr_name, $reverse_chr_name);
    99. 

  99. 	#print "reverse_chr_name=$reverse_chr_name\n"; 
    100. 

  100. 	#-------------------------------------------------------------------------------
    101. 

  101. 	# FIND THE CORRESPONDING CHIMP CHROMOSOME FROM FILE ORTp_chrS.TXT
    102. 

  102. 	foreach my $direct ("reverse_direction","forward_direction"){
    103. 

  103. 		$p_chr=$chromt;
    104. 

  104. 		#print "direction = $direct\n";
    105. 

  105. 		$p_chr = $p_chr."r" if $direct eq "reverse_direction";
    106. 

  106. 		$p_chr = $p_chr if $direct eq "forward_direction";
    107. 

  107. 		my $config = $species_set;
    108. 

  108. 		$config=~s/,/./g;
    109. 

  109. 		my @orgs = split(/\./,$arranged_species_set);
    110. 

  110. 		#print "ORGS= @orgs\n";
    111. 

  111. 		my @tag=@orgs;
    112. 

  112. 			
    113. 

  113. 		
    114. 

  114. 		my $tags = join(",", @tag);
    115. 

  115. 		my @tags=@tag;
    116. 

  116. 		chomp $p_chr;
    117. 

  117. 		$tags = join("_", split(/,/, $tags));
    118. 

  118. 		my $pchr = "chr".$p_chr;
    119. 

  119. 		
    120. 

  120. 		my $ptag = $orgs[0]."-".$pchr.".".join(".",@orgs[1 ... scalar(@orgs)-1])."-".$threshold_array;
    121. 

  121. 		my @sp_tags = ();
    122. 

  122. 		
    123. 

  123. #		print "$ptag _ orthfile\n"; <STDIN>;
    124. 

  124. 		#print "orgs=@orgs, pchr=$pchr, hence, ptag = $ptag\n";
    125. 

  125. 		foreach my $sp (@tag){
    126. 

  126. 			push(@sp_tags, ($sp.".".$ptag));
    127. 

  127. 		}
    128. 

  128. 	
    129. 

  129. 		my $preptag = $orgs[0]."-".$pchr.".".join(".",@orgs[1 ... scalar(@orgs)-1]);
    130. 

  130. 		my @presp_tags = ();
    131. 

  131. 		
    132. 

  132. 		foreach my $sp (@tag){
    133. 

  133. 			push(@presp_tags, ($sp.".".$preptag));
    134. 

  134. 		}
    135. 

  135. 	
    136. 

  136. 		my $resultdir = "";
    137. 

  137. 		my $orthdir = "";
    138. 

  138. 		my $filtereddir = "";
    139. 

  139. 		my $pipedir = "";
    140. 

  140. 		
    141. 

  141. 		my @title_queries = ();
    142. 

  142. 		push(@title_queries, "^[0-9]+");
    143. 

  143. 		my $sep="\\s";
    144. 

  144. 		for my $or (0 ... $#orgs){
    145. 

  145. 			my $title =  join($sep, ($orgs[$or],  "[A-Za-z_]+[0-9a-zA-Z]+", "[0-9]+", "[0-9]+", "[\\-\\+]"));
    146. 

  146. 			#$title =~ s/chr\\+\\s+\+/chr/g;
    147. 

  147. 			push(@title_queries, $title);
    148. 

  148. 		}
    149. 

  149. 		my $title_query = join($sep, @title_queries);
    150. 

  150. 		#print "title_queries=@title_queries\n";
    151. 

  151. 		#print "query = >$title_query<\n"; 
    152. 

  152. 		#print "orgs = @orgs\n"; 
    153. 

  153. 		#-------------------------------------------------------------------------------
    154. 

  154. 		# GET AXTNET FILES, EDIT THEM AND SPLIT THEM INTO HUMAN AND CHIMP INPUT FILES
    155. 

  155. 		my $t1input = $pchr.".".$arranged_species_set.".net.axt";
    156. 

  156. 		
    157. 

  157. 		my @t1outputs = ();
    158. 

  158. 		
    159. 

  159. 		foreach my $sp (@presp_tags){
    160. 

  160. 			push(@t1outputs, $sp."_gap_op");
    161. 

  161. 		}
    162. 

  162. 		
    163. 

  163. 		
    164. 

  164. 		
    165. 

  165. 		multi_species_t1($t1input,$tags,(join(",", @t1outputs)), $title_query); 
    166. 

  166. 		#print "t1outputs=@t1outputs\n";
    167. 

  167. 		#print "done t1\n"; <STDIN>;
    168. 

  168. 		#-------------------------------------------------------------------------------
    169. 

  169. 		#START T2.PL
    170. 

  170. 		
    171. 

  171. 		my $stag  = (); my $tag1 = (); my $tag2 = ();  my $schrs = ();
    172. 

  172. 		
    173. 

  173. 		for my $t (0 ... scalar(@tags)-1){
    174. 

  174. 			multi_species_t2($t1outputs[$t], $tag[$t]);		
    175. 

  175. 		}
    176. 

  176. 		#-------------------------------------------------------------------------------
    177. 

  177. 		#START T2.2.PL
    178. 

  178. 		
    179. 

  179. 		my @temp_tags = @tag;
    180. 

  180. 		
    181. 

  181. 		foreach my $sp (@presp_tags){
    182. 

  182. 			my $t2input =  $sp."_nogap_op_unrand";
    183. 

  183. 			multi_species_t2_2($t2input, shift(@temp_tags));
    184. 

  184. 		}
    185. 

  185. 		undef (@temp_tags);
    186. 

  186. 		
    187. 

  187. 		#-------------------------------------------------------------------------------
    188. 

  188. 		#START SPUTNIK
    189. 

  189. 		
    190. 

  190. 		my @jobIDs = ();
    191. 

  191. 		@temp_tags = @tag;
    192. 

  192. 		my @sput_filelist = ();
    193. 

  193. 		
    194. 

  194. 		foreach my $sp (@presp_tags){
    195. 

  195. 			#print "sp = $sp\n";
    196. 

  196. 			my $sputnikoutput = $pipedir.$sp."_sput_op0";
    197. 

  197. 			my $sputnikinput = $pipedir.$sp."_nogap_op_unrand";
    198. 

  198. 			push(@sput_filelist, $sputnikinput);
    199. 

  199. 			my $sputnikcommand = $sputnikpath." ".$sputnikinput." > ".$sputnikoutput;
    200. 

  200. 		#	print "$sputnikcommand\n";
    201. 

  201. 			my @sputnikcommand_system = $sputnikcommand;
    202. 

  202. 			system(@sputnikcommand_system);
    203. 

  203. 		}
    204. 

  204. 	
    205. 

  205. 		#-------------------------------------------------------------------------------
    206. 

  206. 		#START SPUTNIK OUTPUT CORRECTOR
    207. 

  207. 		
    208. 

  208. 		foreach my $sp (@presp_tags){
    209. 

  209. 			my $corroutput = $pipedir.$sp."_sput_op1";
    210. 

  210. 			my $corrinput = $pipedir.$sp."_sput_op0";
    211. 

  211. 			sputnikoutput_corrector($corrinput,$corroutput);
    212. 

  212. 			
    213. 

  213. 			my $t4output = $pipedir.$sp."_sput_op2";
    214. 

  214. 			multi_species_t4($corroutput,$t4output);
    215. 

  215. 		
    216. 

  216. 			my $t5output = $pipedir.$sp."_sput_op3";
    217. 

  217. 			multi_species_t5($t4output,$t5output);
    218. 

  218. 			#print "done t5.pl for $sp\n";
    219. 

  219. 		
    220. 

  220. 			my $t6output = $pipedir.$sp."_sput_op4";
    221. 

  221. 			multi_species_t6($t5output,$t6output,scalar(@orgs));
    222. 

  222. 		}
    223. 

  223. 		#-------------------------------------------------------------------------------
    224. 

  224. 		#START T9.PL FOR T10.PL AND FOR INTERRUPTED HUNTING
    225. 

  225. 		
    226. 

  226. 		foreach my $sp (@presp_tags){
    227. 

  227. 			my $t9output = $pipedir.$sp."_gap_op_unrand_match";
    228. 

  228. 			my $t9sequence = $pipedir.$sp."_gap_op_unrand2";
    229. 

  229. 			my $t9micro = $pipedir.$sp."_sput_op4";
    230. 

  230. 			t9($t9micro,$t9sequence,$t9output);
    231. 

  231. 			
    232. 

  232. 			my $t9output2 = $pipedir.$sp."_nogap_op_unrand2_match";
    233. 

  233. 			my $t9sequence2 = $pipedir.$sp."_nogap_op_unrand2";
    234. 

  234. 			t9($t9micro,$t9sequence2,$t9output2);
    235. 

  235. 		}
    236. 

  236. 		#print "done both t9.pl for all orgs\n";
    237. 

  237. 	
    238. 

  238. 		#-------------------------------------------------------------------------------
    239. 

  239. 		# FIND COMPOUND MICROSATELLITES
    240. 

  240. 		
    241. 

  241. 		@jobIDs = ();
    242. 

  242. 		my $species_counter = 0;
    243. 

  243. 		
    244. 

  244. 		foreach my $sp (@presp_tags){
    245. 

  245. 			my $simple_microsats=$pipedir.$sp."_sput_op4_simple";
    246. 

  246. 			my $compound_microsats=$pipedir.$sp."_sput_op4_compound";
    247. 

  247. 			my $input_micro = $pipedir.$sp."_sput_op4";
    248. 

  248. 			my $input_seq = $pipedir.$sp."_nogap_op_unrand2_match";
    249. 

  249. 			multiSpecies_compound_microsat_hunter3($input_micro,$input_seq,$simple_microsats,$compound_microsats,$orgs[$species_counter], scalar(@sp_tags), $threshold_array );
    250. 

  250. 			$species_counter++;
    251. 

  251. 		}
    252. 

  252. 		
    253. 

  253. 		#-------------------------------------------------------------------------------
    254. 

  254. 		# READING  AND FILTERING SIMPLE MICROSATELLITES
    255. 

  255. 		my $spcounter2=0;
    256. 

  256. 		foreach my $sp (@sp_tags){
    257. 

  257. 			my $presp = $presp_tags[$spcounter2];
    258. 

  258. 			$spcounter2++;
    259. 

  259. 			my $simple_microsats=$pipedir.$presp."_sput_op4_simple";
    260. 

  260. 			my $simple_filterout = $pipedir.$sp."_sput_op4_simple_filtered";
    261. 

  261. 			my $simple_residue = $pipedir.$sp."_sput_op4_simple_residue";
    262. 

  262. 			multiSpecies_filtering_interrupted_microsats($simple_microsats, $simple_filterout, $simple_residue,$threshold_array,$threshold_array,scalar(@sp_tags));
    263. 

  263. 		}
    264. 

  264. 		
    265. 

  265. 		#-------------------------------------------------------------------------------
    266. 

  266. 		# ANALYZE  COMPOUND MICROSATELLITES FOR BEING INTERRUPTED MICROSATS
    267. 

  267. 		
    268. 

  268. 		$species_counter = 0;
    269. 

  269. 		foreach my $sp (@sp_tags){
    270. 

  270. 			my $presp = $presp_tags[$species_counter];
    271. 

  271. 			my $compound_microsats = $pipedir.$presp."_sput_op4_compound";
    272. 

  272. 			my $analyzed_simple_microsats=$pipedir.$presp."_sput_op4_compound_interrupted";
    273. 

  273. 			my $analyzed_compound_microsats=$pipedir.$presp."_sput_op4_compound_pure";
    274. 

  274. 			my $seq_file = $pipedir.$presp."_nogap_op_unrand2_match";
    275. 

  275. 			multiSpecies_compound_microsat_analyzer($compound_microsats,$seq_file,$analyzed_simple_microsats,$analyzed_compound_microsats,$orgs[$species_counter], scalar(@sp_tags));
    276. 

  276. 			$species_counter++;
    277. 

  277. 		}
    278. 

  278. 		#-------------------------------------------------------------------------------
    279. 

  279. 		# REANALYZE COMPOUND MICROSATELLITES FOR PRESENCE OF SIMPLE ONES WITHIN THEM.. 
    280. 

  280. 		$species_counter = 0;
    281. 

  281. 		
    282. 

  282. 		foreach my $sp (@sp_tags){
    283. 

  283. 			my $presp = $presp_tags[$species_counter];
    284. 

  284. 			my $compound_microsats = $pipedir.$presp."_sput_op4_compound_pure";
    285. 

  285. 			my $compound_interrupted = $pipedir.$presp."_sput_op4_compound_clarifiedInterrupted";
    286. 

  286. 			my $compound_compound = $pipedir.$presp."_sput_op4_compound_compound";
    287. 

  287. 			my $seq_file = $pipedir.$presp."_nogap_op_unrand2_match";
    288. 

  288. 			multiSpecies_compoundClarifyer($compound_microsats,$seq_file,$compound_interrupted,$compound_compound,$orgs[$species_counter], scalar(@sp_tags), "2_4_6_8", "3_4_6_8", "2_4_6_8");
    289. 

  289. 			$species_counter++;
    290. 

  290. 		}
    291. 

  291. 		#-------------------------------------------------------------------------------
    292. 

  292. 		# READING  AND FILTERING SIMPLE AND COMPOUND MICROSATELLITES
    293. 

  293. 		$species_counter = 0;
    294. 

  294. 		
    295. 

  295. 		foreach my $sp (@sp_tags){
    296. 

  296. 			my $presp = $presp_tags[$species_counter];
    297. 

  297. 		
    298. 

  298. 			my $simple_microsats=$pipedir.$presp."_sput_op4_compound_clarifiedInterrupted";
    299. 

  299. 			my $simple_filterout = $pipedir.$sp."_sput_op4_compound_clarifiedInterrupted_filtered";
    300. 

  300. 			my $simple_residue = $pipedir.$sp."_sput_op4_compound_clarifiedInterrupted_residue";
    301. 

  301. 			multiSpecies_filtering_interrupted_microsats($simple_microsats, $simple_filterout, $simple_residue,$threshold_array,$threshold_array,scalar(@sp_tags));
    302. 

  302. 		
    303. 

  303. 			my $simple_microsats2 = $pipedir.$presp."_sput_op4_compound_interrupted";
    304. 

  304. 			my $simple_filterout2 = $pipedir.$sp."_sput_op4_compound_interrupted_filtered";
    305. 

  305. 			my $simple_residue2 = $pipedir.$sp."_sput_op4_compound_interrupted_residue";
    306. 

  306. 			multiSpecies_filtering_interrupted_microsats($simple_microsats2, $simple_filterout2, $simple_residue2,$threshold_array,$threshold_array,scalar(@sp_tags));
    307. 

  307. 		
    308. 

  308. 			my $compound_microsats=$pipedir.$presp."_sput_op4_compound_compound";
    309. 

  309. 			my $compound_filterout = $pipedir.$sp."_sput_op4_compound_compound_filtered";
    310. 

  310. 			my $compound_residue = $pipedir.$sp."_sput_op4_compound_compound_residue";
    311. 

  311. 			multispecies_filtering_compound_microsats($compound_microsats, $compound_filterout, $compound_residue,$threshold_array,$threshold_array,scalar(@sp_tags));
    312. 

  312. 			$species_counter++;
    313. 

  313. 		}
    314. 

  314. 		#print "done filtering both simple and compound microsatellites \n";
    315. 

  315. 		
    316. 

  316. 		#-------------------------------------------------------------------------------
    317. 

  317. 		
    318. 

  318. 		my @combinedarray = ();
    319. 

  319. 		my @combinedarray_indicators = ("mononucleotide", "dinucleotide", "trinucleotide", "tetranucleotide");
    320. 

  320. 		my @combinedarray_tags = ("mono", "di", "tri", "tetra");
    321. 

  321. 		$species_counter = 0;
    322. 

  322. 		
    323. 

  323. 		foreach my $sp (@sp_tags){
    324. 

  324. 			my $simple_interrupted = $pipedir.$sp."_simple_analyzed_simple";
    325. 

  325. 			push @{$combinedarray[$species_counter]}, $pipedir.$sp."_simple_analyzed_simple_mono", $pipedir.$sp."_simple_analyzed_simple_di", $pipedir.$sp."_simple_analyzed_simple_tri", $pipedir.$sp."_simple_analyzed_simple_tetra";
    326. 

  326. 			$species_counter++;
    327. 

  327. 		}
    328. 

  328. 		
    329. 

  329. 		#-------------------------------------------------------------------------------
    330. 

  330. 		# PUT TOGETHER THE INTERRUPTED AND SIMPLE MICROSATELLITES BASED ON THEIR MOTIF SIZE FOR FURTHER EXTENTION	
    331. 

  331. 		my $sp_counter = 0;
    332. 

  332. 		foreach my $sp (@sp_tags){	
    333. 

  333. 			my $analyzed_simple = $pipedir.$sp."_sput_op4_compound_interrupted_filtered";
    334. 

  334. 			my $clarifyed_simple = $pipedir.$sp."_sput_op4_compound_clarifiedInterrupted_filtered";
    335. 

  335. 			my $simple = $pipedir.$sp."_sput_op4_simple_filtered";
    336. 

  336. 			my $simple_analyzed_simple = $pipedir.$sp."_simple_analyzed_simple";
    337. 

  337. 			`cat $analyzed_simple $clarifyed_simple $simple > $simple_analyzed_simple`;
    338. 

  338. 			for my $i (0 ... 3){
    339. 

  339. 				`grep "$combinedarray_indicators[$i]" $simple_analyzed_simple > $combinedarray[$sp_counter][$i]`;	
    340. 

  340. 			}
    341. 

  341. 			$sp_counter++;		
    342. 

  342. 		}
    343. 

  343. 		#print "\ndone grouping interrupted & simple microsats based on their motif size for further extention\n";
    344. 

  344. 		
    345. 

  345. 		#-------------------------------------------------------------------------------
    346. 

  346. 		# BREAK CHROMOSOME INTO PARTS OF CERTAIN NO. CONTIGS EACH, FOR FUTURE SEARCHING OF INTERRUPTED MICROSATELLITES
    347. 

  347. 		# ESPECIALLY DI, TRI AND TETRANUCLEOTIDE MICROSATELLITES
    348. 

  348. 		@temp_tags = @sp_tags;
    349. 

  349. 		my $increment = 1000000;
    350. 

  350. 		my @splist = ();
    351. 

  351. 		my $targetdir = $pipedir;
    352. 

  352. 		$species_counter=0;
    353. 

  353. 		
    354. 

  354. 		foreach my $sp (@sp_tags){
    355. 

  355. 			my $presp = $presp_tags[$species_counter];
    356. 

  356. 			$species_counter++;
    357. 

  357. 			my $localtag = shift @temp_tags;
    358. 

  358. 			my $locallist = $targetdir.$localtag."_".$p_chr."_list";
    359. 

  359. 			push(@splist, $locallist);
    360. 

  360. 			my $input = $pipedir.$presp."_nogap_op_unrand2_match";
    361. 

  361. 			chromosome_unrand_breaker($input,$targetdir,$locallist,$increment, $localtag, $pchr);
    362. 

  362. 		}
    363. 

  363. 		
    364. 

  364. 	
    365. 

  365. 		my @unionarray = ();
    366. 

  366. 		#print "splist=@splist\n";
    367. 

  367. 		#-------------------------------------------------------------------------------
    368. 

  368. 		# FIND INTERRUPTED MICROSATELLITES
    369. 

  369. 		
    370. 

  370. 		$species_counter = 0;		
    371. 

  371. 		
    372. 

  372. 		for my $i (0 .. $#combinedarray){
    373. 

  373. 			
    374. 

  374. 			@jobIDs = ();
    375. 

  375. 			open (JLIST1, "$splist[$i]") or die "Cannot open file $splist[$i]: $!";
    376. 

  376. 			
    377. 

  377. 			while (my $sp1  = <JLIST1>){
    378. 

  378. 				#print "$splist[$i]: sp1=$sp1\n";
    379. 

  379. 				chomp $sp1;
    380. 

  380. 				
    381. 

  381. 				for my $j (0 ... $#combinedarray_tags){
    382. 

  382. 					my $interr  = $sp1."_interr_".$combinedarray_tags[$j];
    383. 

  383. 					my $simple  = $sp1."_simple_".$combinedarray_tags[$j];
    384. 

  384. 					push @{$unionarray[$i]}, $interr, $simple;	
    385. 

  385. 					multiSpecies_interruptedMicrosatHunter($combinedarray[$i][$j],$sp1,$interr ,$simple, $orgs[$species_counter], scalar(@sp_tags), "3_4_6_8"); 
    386. 

  386. 				}		
    387. 

  387. 			}
    388. 

  388. 			$species_counter++;
    389. 

  389. 		}
    390. 

  390. 		close JLIST1;	
    391. 

  391. 		#-------------------------------------------------------------------------------
    392. 

  392. 		#	REUNION AND ZIPPING BEFORE T10.PL
    393. 

  393. 		
    394. 

  394. 		my @allarray = ();
    395. 

  395. 		
    396. 

  396. 		for my $i (0 ... $#sp_tags){	
    397. 

  397. 			my $localfile = $pipedir.$sp_tags[$i]."_allmicrosats";
    398. 

  398. 			unlink $localfile if -e $localfile;
    399. 

  399. 			push(@allarray, $localfile);
    400. 

  400. 		
    401. 

  401. 			my $unfiltered_localfile= $localfile."_unfiltered";
    402. 

  402. 			my $residue_localfile= $localfile."_residue";
    403. 

  403. 			
    404. 

  404. 			unlink $unfiltered_localfile;
    405. 

  405. 			#unlink $unfiltered_localfile;
    406. 

  406. 			for my $j (0 ... $#{$unionarray[$i]}){
    407. 

  407. 				#print "listing files for species $i  and list number $j= \n$unionarray[$i][$j] \n";
    408. 

  408. 				`cat $unionarray[$i][$j] >> $unfiltered_localfile`;
    409. 

  409. 				unlink $unionarray[$i][$j];
    410. 

  410. 			}
    411. 

  411. 	
    412. 

  412. 			multiSpecies_filtering_interrupted_microsats($unfiltered_localfile, $localfile, $residue_localfile,$threshold_array,$threshold_array,scalar(@sp_tags) );
    413. 

  413. 			my $analyzed_compound = $pipedir.$sp_tags[$i]."_sput_op4_compound_compound_filtered";
    414. 

  414. 			my $simple_residue = $pipedir.$sp_tags[$i]."_sput_op4_simple_residue";
    415. 

  415. 			my $compound_residue = $pipedir.$sp_tags[$i]."_sput_op4_compound_residue";
    416. 

  416. 			
    417. 

  417. 			`cat $analyzed_compound >> $localfile`;
    418. 

  418. 		}
    419. 

  419. 		#-------------------------------------------------------------------------------
    420. 

  420. 		# MERGING MICROSATELLITES THAT ARE VERY CLOSE TO EACH OTHER, INCLUDING THOSE FOUND BY SEARCHING IN 2 OPPOSIT DIRECTIONS
    421. 

  421. 		
    422. 

  422. 		my $toescape=0;
    423. 

  423. 		
    424. 

  424. 		
    425. 

  425. 		for my $i (0 ... $#sp_tags){	
    426. 

  426. 			my $localfile = $pipedir.$sp_tags[$i]."_allmicrosats";
    427. 

  427. 			$localfile =~ /$focalspec\-(chr[0-9a-zA-Z]+)\./;
    428. 

  428. 			my $direction = $1;
    429. 

  429. 			#print "localfile = $localfile , direction = $direction\n"; 
    430. 

  430. 	#		`gzip $reverse_chr_name` if $direction =~ /chr[0-9a-zA-Z]+r/ && $switchboard{"deleting_processFiles"} != 1;
    431. 

  431. 			$toescape =1 if $direction =~ /chr[0-9a-zA-Z]+r/;
    432. 

  432. 			last if $direction =~ /chr[0-9a-zA-Z]+r/;
    433. 

  433. 			my $nogap_sequence = $pipedir.$presp_tags[$i]."_nogap_op_unrand2_match";
    434. 

  434. 			my $gap_sequence = $pipedir.$presp_tags[$i]."_gap_op_unrand_match";
    435. 

  435. 			my $reverselocal = $localfile;
    436. 

  436. 			$reverselocal =~ s/\-chr([0-9a-zA-Z]+)\./-chr$1r./g;
    437. 

  437. 			merge_interruptedMicrosats($nogap_sequence,$localfile, $reverselocal ,scalar(@sp_tags)); 
    438. 

  438. 			#-------------------------------------------------------------------------------
    439. 

  439. 			my $forward_separate = $localfile."_separate";
    440. 

  440. 			my $reverse_separate = $reverselocal."_separate";
    441. 

  441. 			my $diff = $forward_separate."_diff";
    442. 

  442. 			my $miss = $forward_separate."_miss";
    443. 

  443. 			my $common = $forward_separate."_common";
    444. 

  444. 			forward_reverse_sputoutput_comparer($nogap_sequence,$forward_separate, $reverse_separate, $diff, $miss, $common ,scalar(@sp_tags)); 
    445. 

  445. 			#-------------------------------------------------------------------------------
    446. 

  446. 			my $symmetrical_file = $localfile."_symmetrical";
    447. 

  447. 			my $merged_file = $localfile."_merged";
    448. 

  448. 			#print "cating: $merged_file $common  into -> $symmetrical_file \n"; 
    449. 

  449. 			`cat $merged_file $common > $symmetrical_file`;
    450. 

  450. 			#-------------------------------------------------------------------------------
    451. 

  451. 			my $t10output = $symmetrical_file."_fin_hit_all_2";
    452. 

  452. 			new_multispecies_t10($gap_sequence, $symmetrical_file, $t10output, join(".", @orgs));
    453. 

  453. 			#-------------------------------------------------------------------------------		
    454. 

  454. 		}
    455. 

  455. 		next if $toescape == 1;
    456. 

  456. 		#------------------------------------------------------------------------------------------------
    457. 

  457. 		# BRINGING IT ALL TOGETHER: FINDING ORTHOLOGOUS MICROSATELLITES AMONG THE SPECIES
    458. 

  458. 		
    459. 

  459. 		
    460. 

  460. 		my @micros_array = ();
    461. 

  461. 		my $sampletag = ();
    462. 

  462. 		for my $i (0 ... $#sp_tags){
    463. 

  463. 			my $finhitFile = $pipedir.$sp_tags[$i]."_allmicrosats_symmetrical_fin_hit_all_2";
    464. 

  464. 			push(@micros_array, $finhitFile);
    465. 

  465. 			$sampletag = $sp_tags[$i];
    466. 

  466. 		}
    467. 

  467. 		#$sampletag =~ s/^([A-Z]+\.)/ORTH_/;
    468. 

  468. 		#$sampletag = $sampletag."_monoThresh-".$mono_threshold."bp";
    469. 

  469. 		my $orthfiletemp = $ptag."_orthfile";
    470. 

  470. 		my $orthanswer = multiSpecies_orthFinder4($t1input, join(":",@micros_array), $orthfiletemp, join(":", @orgs), $separation);
    471. 

  471. 

  472. 		my $maskedorthfiletemp = $ptag."_orthfile_masked";
    473. 

  473. 		qualityFilter ($orthfiletemp, $chr_name_sputt, $maskedorthfiletemp);
    474. 

  474. 

  475. 		push @outputfiles , $maskedorthfiletemp;
    476. 

  476. 	}
    477. 

  477. 	$date = `date`;
    478. 

  478. }
    479. 

  479. 

  480. `cat @outputfiles > $orthfile`;
    481. 

  481. 

  482. my $rootdir = $dir;
    483. 

  483. $rootdir =~ s/\/[A-Za-z0-9\-_]+$//;
    484. 

  484. chdir $rootdir;
    485. 

  485. remove_tree($dir);
    486. 

  486. 

  487. #print "date = $date\n";
    488. 

  488. #remove_tree($tdir);
    489. 

  489. #------------------------------------------------------------------------------------------------
    490. 

  490. #------------------------------------------------------------------------------------------------
    491. 

  491. #------------------------------------------------------------------------------------------------
    492. 

  492. #------------------------------------------------------------------------------------------------
    493. 

  493. 

  494. #xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx 
    495. 

  495. 

  496. sub maftoAxt_multispecies {
    497. 

  497. 	#print "in maftoAxt_multispecies : got @_\n";
    498. 

  498. 	my $fname=$_[0];
    499. 

  499. 	open(IN,"<$_[0]") or die "Cannot open $_[0]: $! \n";
    500. 

  500. 	my $treedefinition = $_[1];
    501. 

  501. 	open(OUT,">$_[2]") or die "Cannot open $_[2]: $! \n";
    502. 

  502. 	my $counter = 0;
    503. 

  503. 	my $exactspeciesset = $_[3];
    504. 

  504. 	my @exactspeciesset_unarranged = split(/,/,$exactspeciesset);
    505. 

  505. 	
    506. 

  506. 	$treedefinition=~s/[\)\(, ]/\t/g;
    507. 

  507. 	my @species=split(/\t+/,$treedefinition);
    508. 

  508. 	my @exactspecies=();
    509. 

  509. 	
    510. 

  510. 	foreach my $spec (@species){
    511. 

  511. 		foreach my $espec (@exactspeciesset_unarranged){
    512. 

  512. 			push @exactspecies, $spec if $spec eq $espec;
    513. 

  513. 		}
    514. 

  514. 	}
    515. 

  515. 	#print "exactspecies=@exactspecies\n";
    516. 

  516. 	
    517. 

  517. 	###########
    518. 

  518. 	my $select = 2;  
    519. 

  519. 	#select = 1 if all species need sequences to be present for each block otherwise, it is 0
    520. 

  520. 	#select = 2 only the allowed set make up the alignment. use the removeset
    521. 

  521. 	# information to detect alignmenets that have other important genomes aligned.
    522. 

  522. 	###########
    523. 

  523. 	my @allowedset = ();
    524. 

  524. 	@allowedset = split(/;/,allowedSetOfSpecies(join("_",@species))) if $select == 0;
    525. 

  525. 	@allowedset = join("_",0,@species) if $select == 1;
    526. 

  526. 	#print "species = @species , allowedset =",join("\n", @allowedset) ," \n"; 
    527. 

  527. 	@allowedset = join("_",0,@exactspecies) if $select == 2;
    528. 

  528. 	#print "allowedset = @allowedset and exactspecies = @exactspecies\n";
    529. 

  529. 	
    530. 

  530. 	my $start = 0;
    531. 

  531. 	my @sequences = ();
    532. 

  532. 	my @titles = ();
    533. 

  533. 	my $species_counter = "0";
    534. 

  534. 	my $countermatch = 0;
    535. 

  535. 	my $outsideSpecies=0;
    536. 

  536. 	
    537. 

  537. 	while(my $line = <IN>){
    538. 

  538. #		print $line;
    539. 

  539. 		next if $line =~ /^#/;
    540. 

  540. 		next if $line =~ /^i/;
    541. 

  541. 		chomp $line;
    542. 

  542. 		my @fields = split(/\s+/,$line);
    543. 

  543. 		chomp $line;
    544. 

  544. 		if ($line =~ /^a /){
    545. 

  545. 			$start = 1;
    546. 

  546. 		}
    547. 

  547. 		
    548. 

  548. 		if ($line =~ /^s /){
    549. 

  549. 		
    550. 

  550. 			foreach my $sp (@allspecies){
    551. 

  551. #				print "checking species $sp\n";
    552. 

  552. 				if ($fields[1] =~ /$sp/){
    553. 

  553. 					$species_counter = $species_counter."_".$sp;
    554. 

  554. 					push(@sequences, $fields[6]);
    555. 

  555. 					my @sp_info = split(/\./,$fields[1]);
    556. 

  556. 					my $title = join(" ",@sp_info, $fields[2], ($fields[2]+$fields[3]), $fields[4]);
    557. 

  557. 					push(@titles, $title);				
    558. 

  558. #					print "species_counter = $species_counter\n";
    559. 

  559. 				}
    560. 

  560. 			}
    561. 

  561. 		}
    562. 

  562. 		
    563. 

  563. 		if (($line !~ /^a/) && ($line !~ /^s/) && ($line !~ /^#/) && ($line !~ /^i/) && ($start = 1)){
    564. 

  564. #			print "species_counter = $species_counter\n"; 
    565. 

  565. 			my $arranged = reorderSpecies($species_counter, @allspecies);
    566. 

  566. 			my $stopper = 1;
    567. 

  567. 			my $arrno = 0;
    568. 

  568. 			
    569. 

  569. #			print "checking if ", scalar(@sequences), " match @exactspecies allowedset=@allowedset\n";
    570. 

  570. 			if (scalar(@sequences) == scalar(@exactspecies)){
    571. 

  571. 				foreach my $set (@allowedset){
    572. 

  572. #					print "testing $arranged against $set\n";
    573. 

  573. 					if ($arranged eq $set){
    574. 

  574. 						$stopper = 0; last;
    575. 

  575. 					}
    576. 

  576. 					$arrno++;
    577. 

  577. 				}
    578. 

  578. 			}
    579. 

  579. 			else{
    580. 

  580. 				$stopper = 1;
    581. 

  581. 			}
    582. 

  582. 			
    583. 

  583. 			
    584. 

  584. 			if ($stopper == 0) {
    585. 

  585. 				@titles = split ";", orderInfo(join(";", @titles), $species_counter, $arranged) if $species_counter ne $arranged;				
    586. 

  586. 				@sequences = split ";", orderInfo(join(";", @sequences), $species_counter, $arranged) if $species_counter ne $arranged;				
    587. 

  587. 				my $filteredseq = filter_gaps(@sequences);
    588. 

  588. 				
    589. 

  589. 				if ($filteredseq ne "SHORT"){
    590. 

  590. 					#print "printing"; <STDIN>;
    591. 

  591. 					$counter++;
    592. 

  592. 					print OUT join (" ",$counter, @titles), "\n";
    593. 

  593. 					print OUT $filteredseq, "\n";
    594. 

  594. 					print OUT "\n"; 
    595. 

  595. 					$countermatch++;
    596. 

  596. 				}
    597. 

  597. 			}
    598. 

  598. 			else{ #print "nexting\n";<STDIN>;
    599. 

  599. 			}
    600. 

  600. 	
    601. 

  601. 			@sequences = (); @titles = (); $start = 0;$species_counter = "0";
    602. 

  602. 			next;		
    603. 

  603. 			
    604. 

  604. 		}
    605. 

  605. 	}
    606. 

  606. #	print "countermatch = $countermatch\n";
    607. 

  607. }
    608. 

  608. 

  609. sub reorderSpecies{
    610. 

  610. 	my @inarr=@_;
    611. 

  611. 	my $currSpecies = shift (@inarr);
    612. 

  612. 	my $ordered_species = 0;
    613. 

  613. 	my @species=@inarr;
    614. 

  614. 	#print "species = @species\n";
    615. 

  615. 	foreach my $order (@species){
    616. 

  616. 		$ordered_species = $ordered_species."_".$order	if	$currSpecies=~ /$order/;
    617. 

  617. 	}
    618. 

  618. 	return $ordered_species;
    619. 

  619. 

  620. }
    621. 

  621. 

  622. sub filter_gaps{
    623. 

  623. 	my @sequences = @_;
    624. 

  624. #	print "sequences sent are @sequences\n";
    625. 

  625. 	my $seq_length = length($sequences[0]);
    626. 

  626. 	my $seq_no = scalar(@sequences);
    627. 

  627. 	my $allgaps = ();
    628. 

  628. 	for (1 ... $seq_no){
    629. 

  629. 		$allgaps = $allgaps."-";
    630. 

  630. 	}
    631. 

  631. 	
    632. 

  632. 	my @seq_array = ();
    633. 

  633. 	my $seq_counter = 0;
    634. 

  634. 	foreach my $seq (@sequences){
    635. 

  635. #		my @sequence = split(/\s*/,$seq);
    636. 

  636. 		$seq_array[$seq_counter] = [split(/\s*/,$seq)];
    637. 

  637. #		push @seq_array, [@sequence];
    638. 

  638. 		$seq_counter++;
    639. 

  639. 	}
    640. 

  640. 	my $g = 0;
    641. 

  641. 	while ( $g < $seq_length){
    642. 

  642. 		last if (!exists $seq_array[0][$g]);
    643. 

  643. 		my $bases = ();
    644. 

  644. 		for my $u (0 ... $#seq_array){
    645. 

  645. 			$bases = $bases.$seq_array[$u][$g];
    646. 

  646. 		}	
    647. 

  647. #		print $bases, "\n";
    648. 

  648. 		if ($bases eq $allgaps){
    649. 

  649. #			print "bases are $bases, position is $g \n";
    650. 

  650. 			for my $seq (@seq_array){
    651. 

  651. 				splice(@$seq , $g, 1);
    652. 

  652. 			}
    653. 

  653. 		}
    654. 

  654. 		else {
    655. 

  655. 			$g++;
    656. 

  656. 		}
    657. 

  657. 	}
    658. 

  658. 	
    659. 

  659. 	my @outs = ();
    660. 

  660. 	
    661. 

  661. 	foreach my $seq (@seq_array){
    662. 

  662. 		push(@outs, join("",@$seq));
    663. 

  663. 	}
    664. 

  664. 	return "SHORT" if length($outs[0]) <=100;
    665. 

  665. 	return (join("\n", @outs));	
    666. 

  666. }
    667. 

  667. 

  668. 

  669. sub allowedSetOfSpecies{
    670. 

  670. 	my @allowed_species = split(/_/,$_[0]);
    671. 

  671. 	unshift @allowed_species, 0;
    672. 

  672. #	print "allowed set = @allowed_species \n";
    673. 

  673. 	my @output = ();
    674. 

  674. 	for (0 ... scalar(@allowed_species) - 4){
    675. 

  675. 		push(@output, join("_",@allowed_species));
    676. 

  676. 		pop @allowed_species;
    677. 

  677. 	}
    678. 

  678. 	return join(";",reverse(@output));
    679. 

  679. 

  680. }
    681. 

  681. 

  682. 

  683. sub orderInfo{
    684. 

  684. 	my @info = split(/;/,$_[0]);
    685. 

  685. #	print "info = @info";
    686. 

  686. 	my @old = split(/_/,$_[1]);
    687. 

  687. 	my @new = split(/_/,$_[2]);
    688. 

  688. 	shift @old; shift @new;
    689. 

  689. 	my @outinfo = ();
    690. 

  690. 	foreach my $spe (@new){
    691. 

  691. 		for my $no (0 ... $#old){
    692. 

  692. 			if ($spe eq $old[$no]){
    693. 

  693. 				push(@outinfo, $info[$no]);
    694. 

  694. 			}
    695. 

  695. 		}
    696. 

  696. 	}
    697. 

  697. #	print "outinfo = @outinfo \n"; 
    698. 

  698. 	return join(";", @outinfo);
    699. 

  699. }
    700. 

  700. 

  701. #xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx xxxxxxx maftoAxt_multispecies xxxxxxx 
    702. 

  702. 

  703. #xxxxxxx artificial_axdata_inverter xxxxxxx xxxxxxx artificial_axdata_inverter xxxxxxx xxxxxxx artificial_axdata_inverter xxxxxxx 
    704. 

  704. sub artificial_axdata_inverter{
    705. 

  705. 	open(IN,"<$_[0]") or die "Cannot open file $_[0]: $!";
    706. 

  706. 	open(OUT,">$_[1]") or die "Cannot open file $_[1]: $!";
    707. 

  707. 	my $linecounter=0;
    708. 

  708. 	while (my $line = <IN>){
    709. 

  709. 		$linecounter++;
    710. 

  710. 		#print "$linecounter\n";
    711. 

  711. 		chomp $line;	
    712. 

  712. 		my $final_line = $line;	
    713. 

  713. 		my $trycounter = 0;
    714. 

  714. 		if ($line =~ /^[a-zA-Z\-]/){
    715. 

  715. 		#	while ($final_line eq $line){
    716. 

  716. 				my @fields = split(/\s*/,$line);
    717. 

  717. 				
    718. 

  718. 				$final_line = join("",reverse(@fields));
    719. 

  719. 		#		print colored ['red'], "$line\n$final_line\n" if $final_line eq $line && $line !~ /chr/ && $line =~ /[a-zA-Z]/;
    720. 

  720. 		#		$trycounter++;
    721. 

  721. 		#		print "trying again....$trycounter : $final_line\n" if $final_line eq $line;
    722. 

  722. 		#	}
    723. 

  723. 		}
    724. 

  724. 		
    725. 

  725. 	#	print colored ['yellow'], "$line\n$final_line\n" if $final_line eq $line && $line !~ /chr/ && $line =~ /[a-zA-Z]/;
    726. 

  726. 		if ($line =~ /^[0-9]/){
    727. 

  727. 			$line =~ s/chr([A-Z0-9a-b]+)/chr$1r/g;
    728. 

  728. 			$final_line = $line;
    729. 

  729. 		}
    730. 

  730. 		print OUT $final_line,"\n";
    731. 

  731. 		#print "$line\n$final_line\n" if $final_line eq $line && $line !~ /chr/ && $line =~ /[a-zA-Z]/;
    732. 

  732. 	}
    733. 

  733. 	close OUT;
    734. 

  734. }
    735. 

  735. #xxxxxxx artificial_axdata_inverter xxxxxxx xxxxxxx artificial_axdata_inverter xxxxxxx xxxxxxx artificial_axdata_inverter xxxxxxx 
    736. 

  736. 

  737. 

  738. #xxxxxxx multi_species_t1 xxxxxxx xxxxxxx multi_species_t1 xxxxxxx xxxxxxx multi_species_t1 xxxxxxx 
    739. 

  739. 

  740. sub multi_species_t1 {
    741. 

  741. 

  742. 	my $input1 = $_[0];
    743. 

  743. 	#print "@_\n"; <STDIN>;
    744. 

  744. 	my @tags = split(/_/, $_[1]);
    745. 

  745. 	my @outputs = split(/,/, $_[2]);
    746. 

  746. 	my $title_query = $_[3];
    747. 

  747. 	my @handles = ();
    748. 

  748. 	
    749. 

  749. 	open(FILEB,"<$input1")or die "Cannot open file: $input1 $!";
    750. 

  750. 	my $i = 0;
    751. 

  751. 	foreach my $path (@outputs){
    752. 

  752. 		$handles[$i] = IO::Handle->new();
    753. 

  753. 		open ($handles[$i], ">$path") or die "Can't open $path : $!";
    754. 

  754. 		$i++;
    755. 

  755. 	}
    756. 

  756. 	
    757. 

  757. 	my $curdef;
    758. 

  758. 	my $start = 0;
    759. 

  759. 	
    760. 

  760. 	while (my $line = <FILEB> ) {
    761. 

  761. 		if ($line =~ /^\d/){
    762. 

  762. 			$line =~ s/ +/\t/g;
    763. 

  763. 			my @fields = split(/\s+/, $line);
    764. 

  764. 			if (($line =~ /$title_query/)){
    765. 

  765. 				my $title = $line;
    766. 

  766. 				my $counter = 0;
    767. 

  767. 				foreach my $tag (@tags){
    768. 

  768. 					$line = <FILEB>;
    769. 

  769. 					print {$handles[$counter]} ">",$tag,"\t",$title, " ",$line;  	
    770. 

  770. 					$counter++;
    771. 

  771. 				}
    772. 

  772. 			}
    773. 

  773. 			else{
    774. 

  774. 					foreach my $tag (@tags){
    775. 

  775. 					my $tine = <FILEB>;
    776. 

  776. 				}		
    777. 

  777. 			}
    778. 

  778. 		
    779. 

  779. 		}
    780. 

  780. 	}
    781. 

  781. 	
    782. 

  782. 	foreach my $hand (@handles){
    783. 

  783. 		$hand->close(); 
    784. 

  784. 	}
    785. 

  785. 	
    786. 

  786. 	close FILEB;
    787. 

  787. }
    788. 

  788. 

  789. #xxxxxxx multi_species_t1 xxxxxxx xxxxxxx multi_species_t1 xxxxxxx xxxxxxx multi_species_t1 xxxxxxx 
    790. 

  790. 

  791. #xxxxxxx multi_species_t2 xxxxxxx xxxxxxx multi_species_t2 xxxxxxx xxxxxxx multi_species_t2 xxxxxxx 
    792. 

  792. 

  793. sub multi_species_t2{
    794. 

  794. 	
    795. 

  795. 	my $input = $_[0];
    796. 

  796. 	my $species = $_[1];
    797. 

  797. 	my $output1 = $input."_unr";
    798. 

  798. 	
    799. 

  799. 	#------------------------------------------------------------------------------------------
    800. 

  800. 	open (FILEF1, "<$input") or die "Cannot open file $input :$!";
    801. 

  801. 	open (FILEF2, ">$output1") or die "Cannot open file $output1 :$!";
    802. 

  802. 	
    803. 

  803. 	my $line1 = <FILEF1>;
    804. 

  804. 	
    805. 

  805. 	while($line1){
    806. 

  806. 	{
    807. 

  807. 	#    chomp($line);
    808. 

  808. 		if ($line1 =~ (m/^\>$species/)){
    809. 

  809. 		chomp($line1);
    810. 

  810. 		print FILEF2 $line1;
    811. 

  811. 		$line1 = <FILEF1>;
    812. 

  812. 		chomp($line1);
    813. 

  813. 		print FILEF2 "\t", $line1,"\n";
    814. 

  814. 	   }   
    815. 

  815. 	}
    816. 

  816. 	$line1 = <FILEF1>;
    817. 

  817. 	}
    818. 

  818. 	
    819. 

  819. 	close FILEF1;
    820. 

  820. 	close FILEF2;
    821. 

  821. 	#------------------------------------------------------------------------------------------
    822. 

  822. 	
    823. 

  823. 	my $output2 = $output1."and";
    824. 

  824. 	my $output3 = $output1."and2";
    825. 

  825. 	open(IN,"<$output1");
    826. 

  826. 	open (FILEF3, ">$output2");
    827. 

  827. 	open (FILEF4, ">$output3");
    828. 

  828. 	
    829. 

  829. 	
    830. 

  830. 	while (<IN>){
    831. 

  831. 		my $line = $_;
    832. 

  832. 		chomp($line);
    833. 

  833. 		my @fields=split (/\t/, $line);
    834. 

  834. 	  #  print $line,"\n"; <STDIN>;
    835. 

  835. 		if($line !~ /random/){
    836. 

  836. 			print FILEF3 join ("\t",@fields[0 ... scalar(@fields)-2]), "\n", $fields[scalar(@fields)-1], "\n";
    837. 

  837. 			print FILEF4 join ("\t",@fields[0 ... scalar(@fields)-2]), "\t", $fields[scalar(@fields)-1], "\n";	
    838. 

  838. 		}
    839. 

  839. 	}
    840. 

  840. 	
    841. 

  841. 	
    842. 

  842. 	close IN;
    843. 

  843. 	close FILEF3;
    844. 

  844. 	close FILEF4;
    845. 

  845. 	unlink $output1;
    846. 

  846. 	
    847. 

  847. 	#------------------------------------------------------------------------------------------
    848. 

  848. 	# OLD T3.PL RUDIMENT
    849. 

  849. 	
    850. 

  850. 	my $t3output = $output2;
    851. 

  851. 	$t3output =~ s/gap_op_unrand/nogap_op_unrand/g;
    852. 

  852. 	
    853. 

  853. 	open(IN,"<$output2");
    854. 

  854. 	open(OUTA,">$t3output");
    855. 

  855. 	
    856. 

  856. 	
    857. 

  857. 	while (<IN>){
    858. 

  858. 		s/-//g unless /^>/;
    859. 

  859. 		print OUTA;
    860. 

  860. 	}
    861. 

  861. 	
    862. 

  862. 	close IN;
    863. 

  863. 	close OUTA;
    864. 

  864. 	#------------------------------------------------------------------------------------------
    865. 

  865. }
    866. 

  866. #xxxxxxx multi_species_t2 xxxxxxx xxxxxxx multi_species_t2 xxxxxxx xxxxxxx multi_species_t2 xxxxxxx 
    867. 

  867. 

  868. 

  869. #xxxxxxx multi_species_t2_2 xxxxxxx xxxxxxx multi_species_t2_2 xxxxxxx xxxxxxxmulti_species_t2_2 xxxxxxx 
    870. 

  870. sub multi_species_t2_2{
    871. 

  871. 	#print "IN multi_species_t2_2 : @_\n";
    872. 

  872. 	my $input = $_[0];
    873. 

  873. 	my $species = $_[1];
    874. 

  874. 	my $output1 = $input."2";
    875. 

  875. 	
    876. 

  876. 	
    877. 

  877. 	open (FILEF1, "<$input");
    878. 

  878. 	open (FILEF2, ">$output1");
    879. 

  879. 	
    880. 

  880. 	my $line1 = <FILEF1>;
    881. 

  881. 	
    882. 

  882. 	while($line1){
    883. 

  883. 	{
    884. 

  884. 	#    chomp($line);
    885. 

  885. 		if ($line1 =~ (m/^\>$species/)){
    886. 

  886. 		chomp($line1);
    887. 

  887. 		print FILEF2 $line1;
    888. 

  888. 		$line1 = <FILEF1>;
    889. 

  889. 		chomp($line1);
    890. 

  890. 		print FILEF2 "\t", $line1,"\n";
    891. 

  891. 	   }   
    892. 

  892. 	}
    893. 

  893. 	$line1 = <FILEF1>;
    894. 

  894. 	}
    895. 

  895. 	
    896. 

  896. 	close FILEF1;
    897. 

  897. 	close FILEF2;
    898. 

  898. }
    899. 

  899. 

  900. #xxxxxxx multi_species_t2_2 xxxxxxx xxxxxxx multi_species_t2_2 xxxxxxx xxxxxxx multi_species_t2_2 xxxxxxx 
    901. 

  901. 

  902. 

  903. #xxxxxxx sputnikoutput_corrector xxxxxxx xxxxxxx sputnikoutput_corrector xxxxxxx xxxxxxx sputnikoutput_corrector xxxxxxx 
    904. 

  904. sub sputnikoutput_corrector{
    905. 

  905. 	my $input = $_[0];
    906. 

  906. 	my $output = $_[1];
    907. 

  907. 	open(IN,"<$input") or die "Cannot open file $input :$!";
    908. 

  908. 	open(OUT,">$output") or die "Cannot open file $output :$!";
    909. 

  909. 	my $tine;
    910. 

  910. 	while (my $line=<IN>){
    911. 

  911. 		if($line =~/length /){
    912. 

  912. 			$tine = $line;
    913. 

  913. 			$tine =~ s/\s+/\t/g;
    914. 

  914. 			my @fields = split(/\t/,$tine);	
    915. 

  915. 			if ($fields[6] > 60){
    916. 

  916. 				print OUT $line;
    917. 

  917. 				$line = <IN>;
    918. 

  918. 			
    919. 

  919. 				while (($line !~ /nucleotide/) && ($line !~ /^>/)){
    920. 

  920. 					chomp $line;
    921. 

  921. 					print OUT $line;
    922. 

  922. 					$line = <IN>;
    923. 

  923. 				}
    924. 

  924. 				print OUT "\n";
    925. 

  925. 				print OUT $line;
    926. 

  926. 			}
    927. 

  927. 			else{
    928. 

  928. 				print OUT $line;
    929. 

  929. 			}
    930. 

  930. 		}
    931. 

  931. 		else{
    932. 

  932. 			print OUT $line;
    933. 

  933. 		}
    934. 

  934. 	}
    935. 

  935. 	close IN;
    936. 

  936. 	close OUT;
    937. 

  937. }
    938. 

  938. #xxxxxxx sputnikoutput_corrector xxxxxxx xxxxxxx sputnikoutput_corrector xxxxxxx xxxxxxx sputnikoutput_corrector xxxxxxx 
    939. 

  939. 

  940. 

  941. #xxxxxxx multi_species_t4 xxxxxxx xxxxxxx multi_species_t4 xxxxxxx xxxxxxx multi_species_t4 xxxxxxx 
    942. 

  942. sub multi_species_t4{
    943. 

  943. #	print "multi_species_t4 : @_\n";
    944. 

  944. 	my $input = $_[0];
    945. 

  945. 	my $output = $_[1];
    946. 

  946. 	open (FILEA, "<$input");
    947. 

  947. 	open (FILEB, ">$output");
    948. 

  948. 	
    949. 

  949. 	my $line = <FILEA>;
    950. 

  950. 	
    951. 

  951. 	while ($line) {
    952. 

  952. 	   # chomp $line;
    953. 

  953. 		if ($line =~ />/) {
    954. 

  954. 			chomp $line;
    955. 

  955. 			print FILEB $line, "\n"; 
    956. 

  956. 		}
    957. 

  957. 	
    958. 

  958. 		
    959. 

  959. 		if ($line =~ /^m/ | $line =~ /^d/ | $line =~ /^t/ | $line =~ /^p/){
    960. 

  960. 		chomp $line;
    961. 

  961. 		print FILEB $line, " " ;
    962. 

  962. 		$line = <FILEA>;
    963. 

  963. 		chomp $line;
    964. 

  964. 		print FILEB $line,"\n";
    965. 

  965. 		}
    966. 

  966. 	
    967. 

  967. 		$line = <FILEA>;
    968. 

  968. 	}
    969. 

  969. 	
    970. 

  970. 	
    971. 

  971. 	close FILEA;
    972. 

  972. 	close FILEB;
    973. 

  973. 

  974. }
    975. 

  975. 

  976. #xxxxxxx multi_species_t4 xxxxxxx xxxxxxx multi_species_t4 xxxxxxx xxxxxxx multi_species_t4 xxxxxxx 
    977. 

  977. 

  978. 

  979. #xxxxxxx multi_species_t5 xxxxxxx xxxxxxx multi_species_t5 xxxxxxx xxxxxxx multi_species_t5 xxxxxxx 
    980. 

  980. sub multi_species_t5{
    981. 

  981. 	
    982. 

  982. 	my $input = $_[0];
    983. 

  983. 	my $output = $_[1];
    984. 

  984. 	
    985. 

  985. 	open(FILEB,"<$input");
    986. 

  986. 	open(FILEC,">$output");
    987. 

  987. 	
    988. 

  988. 	my $curdef;
    989. 

  989. 	
    990. 

  990. 	while (my $line = <FILEB> ) {
    991. 

  991. 	
    992. 

  992. 		if ($line =~ /^>/){
    993. 

  993. 		chomp $line;
    994. 

  994. 		$curdef = $line;
    995. 

  995. 		next;
    996. 

  996. 	}
    997. 

  997. 	
    998. 

  998. 	if ($line =~ /^m/ | $line =~ /^d/ | $line =~ /^t/ | $line =~ /^p/){
    999. 

  999. 		print  FILEC $curdef," ",$line;
    1000. 

  1000. 	}
    1001. 

  1001. 	
    1002. 

  1002. 	}
    1003. 

  1003. 	
    1004. 

  1004. 	
    1005. 

  1005. 	close FILEB;
    1006. 

  1006. 	close FILEC;
    1007. 

  1007. 

  1008. }
    1009. 

  1009. #xxxxxxx multi_species_t5 xxxxxxx xxxxxxx multi_species_t5 xxxxxxx xxxxxxx multi_species_t5 xxxxxxx 
    1010. 

  1010. 

  1011. 

  1012. #xxxxxxx multi_species_t6 xxxxxxx xxxxxxx multi_species_t6 xxxxxxx xxxxxxx multi_species_t6 xxxxxxx 
    1013. 

  1013. sub multi_species_t6{
    1014. 

  1014. 	my $input = $_[0];
    1015. 

  1015. 	my $output = $_[1];	
    1016. 

  1016. 	my $focalstrand=$_[3];
    1017. 

  1017. #	print "inpput = @_\n"; 
    1018. 

  1018. 	open (FILE, "<$input");
    1019. 

  1019. 	open (FILE_MICRO, ">$output");
    1020. 

  1020. 	my $linecounter=0;	
    1021. 

  1021. 	while (my $line = <FILE>){
    1022. 

  1022. 		$linecounter++;
    1023. 

  1023. 		chomp $line;
    1024. 

  1024. 		#print "line = $line\n";
    1025. 

  1025. 		#MONO#
    1026. 

  1026. 		$line =~ /$focalspec\s[a-zA-Z]+[0-9a-zA-Z]+\s[0-9]+\s[0-9]+\s([+\-])/;
    1027. 

  1027. 		my $strand=$1;
    1028. 

  1028. 		my $no_of_species = ($line =~ s/\s+[+\-]\s+/ /g);
    1029. 

  1029. 		#print "line = $line\n";
    1030. 

  1030. 		my $specfieldsend = 2 + ($no_of_species*4) - 1;
    1031. 

  1031. 		my @fields = split(/\s+/, $line);
    1032. 

  1032. 		my @speciesdata = @fields[0 ... $specfieldsend];
    1033. 

  1033. 		$line =~ /([a-z]+nucleotide)\s([0-9]+)\s:\s([0-9]+)/;
    1034. 

  1034. 		my ($tide, $start, $end) = ($1, $2, $3);
    1035. 

  1035. 		#print "no_of_species=$no_of_species.. speciesdata = @speciesdata and ($tide, $start, $end)\n";
    1036. 

  1036. 		if($line  =~ /mononucleotide/){
    1037. 

  1037. 			print FILE_MICRO join("\t",@speciesdata, $tide, $start, $strand,$end, $fields[$#fields],  mono($fields[$#fields]),),"\n";
    1038. 

  1038. 		}
    1039. 

  1039. 		#DI#	
    1040. 

  1040. 		elsif($line =~ /dinucleotide/){
    1041. 

  1041. 			print FILE_MICRO join("\t",@speciesdata, $tide, $start, $strand,$end, $fields[$#fields],  di($fields[$#fields]),),"\n";
    1042. 

  1042. 		}
    1043. 

  1043. 		#TRI#	
    1044. 

  1044. 		elsif($line =~ /trinucleotide/ ){
    1045. 

  1045. 			print FILE_MICRO join("\t",@speciesdata, $tide, $start, $strand,$end, $fields[$#fields],  tri($fields[$#fields]),),"\n";
    1046. 

  1046. 		}
    1047. 

  1047. 		#TETRA#
    1048. 

  1048. 		elsif($line =~ /tetranucleotide/){
    1049. 

  1049. 			print FILE_MICRO join("\t",@speciesdata, $tide, $start, $strand,$end, $fields[$#fields],  tetra($fields[$#fields]),),"\n";
    1050. 

  1050. 		}
    1051. 

  1051. 		#PENTA#
    1052. 

  1052. 		elsif($line =~ /pentanucleotide/){
    1053. 

  1053. 			#print FILE_MICRO join("\t",@speciesdata, $tide, $start, $strand,$end, $fields[$#fields],  penta($fields[$#fields]),),"\n";
    1054. 

  1054. 		}
    1055. 

  1055. 		else{
    1056. 

  1056. 		#	print "not: @fields\n";
    1057. 

  1057. 		}
    1058. 

  1058. 	}
    1059. 

  1059. #	print "linecounter=$linecounter\n"; 
    1060. 

  1060. 	close FILE;
    1061. 

  1061. 	close FILE_MICRO;
    1062. 

  1062. }
    1063. 

  1063. 

  1064. sub mono {
    1065. 

  1065. 	my $st = $_[0];
    1066. 

  1066. 	my $tp = unpack "A1"x(length($st)/1),$st;
    1067. 

  1067. 	my $var1 = substr($tp, 0, 1);
    1068. 

  1068. 	return join ("\t", $var1);
    1069. 

  1069. }
    1070. 

  1070. sub di {
    1071. 

  1071. 	my $st = $_[0];
    1072. 

  1072. 	my $tp = unpack "A2"x(length($st)/2),$st;
    1073. 

  1073. 	my $var1 = substr($tp, 0, 2);
    1074. 

  1074. 	return join ("\t", $var1);
    1075. 

  1075. }
    1076. 

  1076. sub tri {
    1077. 

  1077. 	my $st = $_[0];
    1078. 

  1078. 	my $tp = unpack "A3"x(length($st)/3),$st;
    1079. 

  1079. 	my $var1 = substr($tp, 0, 3);
    1080. 

  1080. 	return join ("\t", $var1);
    1081. 

  1081. }
    1082. 

  1082. sub tetra {
    1083. 

  1083. 	my $st = $_[0];
    1084. 

  1084. 	my $tp = unpack "A4"x(length($st)/4),$st;
    1085. 

  1085. 	my $var1 = substr($tp, 0, 4);
    1086. 

  1086. 	return join ("\t", $var1);
    1087. 

  1087. }
    1088. 

  1088. sub penta {
    1089. 

  1089. 	my $st = $_[0];
    1090. 

  1090. 	my $tp = unpack "A5"x(length($st)/5),$st;
    1091. 

  1091. 	my $var1 = substr($tp, 0, 5);
    1092. 

  1092. 	return join ("\t", $var1);
    1093. 

  1093. }
    1094. 

  1094. 

  1095. #xxxxxxx multi_species_t6 xxxxxxx xxxxxxx multi_species_t6 xxxxxxx xxxxxxx multi_species_t6 xxxxxxx 
    1096. 

  1096. 

  1097. 

  1098. #xxxxxxxxxxxxxx t9 xxxxxxxxxxxxxx xxxxxxxxxxxxxx t9 xxxxxxxxxxxxxx xxxxxxxxxxxxxx t9 xxxxxxxxxxxxxx 
    1099. 

  1099. sub t9{
    1100. 

  1100. 	my $input1 = $_[0];
    1101. 

  1101. 	my $input2 = $_[1];
    1102. 

  1102. 	my $output = $_[2];
    1103. 

  1103. 	
    1104. 

  1104. 	
    1105. 

  1105. 	open(IN1,"<$input1") if -e $input1; 
    1106. 

  1106. 	open(IN2,"<$input2") or die "cannot open file $_[1] : $!";
    1107. 

  1107. 	open(OUT,">$output") or die "cannot open file $_[2] : $!"; 
    1108. 

  1108. 	
    1109. 

  1109. 	
    1110. 

  1110. 	my %seen = ();
    1111. 

  1111. 	my $prevkey = 0;
    1112. 

  1112. 	
    1113. 

  1113. 	if (-e $input1){
    1114. 

  1114. 		while (my $line = <IN1>){
    1115. 

  1115. 			chomp($line);
    1116. 

  1116. 			my @fields = split(/\t/,$line);
    1117. 

  1117. 			my $key1 = join ("_K10K1_",@fields[0,1,3,4,5]);
    1118. 

  1118. 		#	print "key in t9 = $key1\n";
    1119. 

  1119. 			$seen{$key1}++	if ($prevkey ne $key1) ;
    1120. 

  1120. 			$prevkey = $key1;
    1121. 

  1121. 		}
    1122. 

  1122. #		print "done first hash\n";
    1123. 

  1123. 		close IN1;
    1124. 

  1124. 	}
    1125. 

  1125. 	
    1126. 

  1126. 	while (my $line = <IN2>){
    1127. 

  1127. 	#	print $line, "**\n";
    1128. 

  1128. 		if (-e $input1){
    1129. 

  1129. 			chomp($line);
    1130. 

  1130. 			my @fields = split(/\t/,$line);
    1131. 

  1131. 			my $key2 = join ("_K10K1_",@fields[0,1,3,4,5]);
    1132. 

  1132. 			if (exists $seen{$key2}){
    1133. 

  1133. 				print OUT "$line\n"	;
    1134. 

  1134. 				delete $seen{$key2};
    1135. 

  1135. 			}
    1136. 

  1136. 		}
    1137. 

  1137. 		else {
    1138. 

  1138. 			print OUT "$line\n"	;		
    1139. 

  1139. #			print  "$line\n"	;		
    1140. 

  1140. 		}
    1141. 

  1141. 	}
    1142. 

  1142. 	
    1143. 

  1143. 	close IN2;
    1144. 

  1144. 	close OUT;
    1145. 

  1145. }
    1146. 

  1146. #xxxxxxxxxxxxxx t9 xxxxxxxxxxxxxx xxxxxxxxxxxxxx t9 xxxxxxxxxxxxxx xxxxxxxxxxxxxx t9 xxxxxxxxxxxxxx 
    1147. 

  1147. 

  1148. 

  1149. #xxxxxxxxxxxxxx multiSpecies_compound_microsat_hunter3 xxxxxxxxxxxxxx  multiSpecies_compound_microsat_hunter3 xxxxxxxxxxxxxx  multiSpecies_compound_microsat_hunter3 xxxxxxxxxxxxxx 
    1150. 

  1150. 

  1151. 

  1152. sub multiSpecies_compound_microsat_hunter3{	
    1153. 

  1153. 	
    1154. 

  1154. 	my $input1 = $_[0];  ###### the *_sput_op4_ii file
    1155. 

  1155. 	my $input2 = $_[1];  ###### looks like this: my $t8humanoutput = $pipedir.$ptag."_nogap_op_unrand2"
    1156. 

  1156. 	my $output1 = $_[2]; ###### plain microsatellite file
    1157. 

  1157. 	my $output2 = $_[3]; ###### compound microsatellite file
    1158. 

  1158. 	my $org = $_[4]; ###### 1 or 2
    1159. 

  1159. 	$no_of_species = $_[5];
    1160. 

  1160. 	#print "IN multiSpecies_compound_microsat_hunter3: @_\n"; 
    1161. 

  1161. 	#my @tags = split(/\t/,$info);
    1162. 

  1162. 	sub compoundify;
    1163. 

  1163. 	open(IN,"<$input1") or die "Cannot open file $input1 $!";
    1164. 

  1164. 	open(SEQ,"<$input2") or die "Cannot open file $input2 $!";
    1165. 

  1165. 	open(OUT,">$output1") or die "Cannot open file $output1 $!";
    1166. 

  1166. 	open(OUT2,">$output2") or die "Cannot open file $output2 $!";
    1167. 

  1167. 	$infocord = 2 + (4*$no_of_species) - 1;
    1168. 

  1168. 	$startcord = 2 + (4*$no_of_species) + 2 - 1;
    1169. 

  1169. 	$strandcord = 2 + (4*$no_of_species) + 3 - 1;
    1170. 

  1170. 	$endcord = 2 + (4*$no_of_species) + 4 - 1;
    1171. 

  1171. 	$microsatcord = 2 + (4*$no_of_species) + 5 - 1;
    1172. 

  1172. 	$motifcord = 2 + (4*$no_of_species) + 6 - 1;
    1173. 

  1173. 	my $sequencepos = 2 + (5*$no_of_species) + 1 -1 ; 
    1174. 

  1174. 	
    1175. 

  1175. 	my @thresholds = ("0");
    1176. 

  1176. 	push(@thresholds, split(/_/,$_[6]));
    1177. 

  1177. 	sub thresholdCheck;
    1178. 

  1178. 	my %micros = ();
    1179. 

  1179. 	while (my $line = <IN>){
    1180. 

  1180. 	#	print "$org\t(chr[0-9]+)\t([0-9]+)\t([0-9])+\t \n";
    1181. 

  1181. 		next if $line =~ /\t\t/;
    1182. 

  1182. 		if ($line =~ /^>[A-Za-z0-9_]+\s+([0-9]+)\s+([a-zA-Z0-9]+)\s([a-zA-Z]+[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)\s/ ) {
    1183. 

  1183. 			my $key = join("\t",$1, $2, $3, $4, $5);
    1184. 

  1184. 		#	print $key, "#-#-#-#-#-#-#-#\n";
    1185. 

  1185. 			push (@{$micros{$key}},$line);		
    1186. 

  1186. 		}
    1187. 

  1187. 		else{
    1188. 

  1188. 		}
    1189. 

  1189. 	}
    1190. 

  1190. 	close IN;
    1191. 

  1191. 	my @deletedlines = ();
    1192. 

  1192. 	
    1193. 

  1193. 	my $linecount = 0;
    1194. 

  1194. 	
    1195. 

  1195. 	while(my $sine = <SEQ>){
    1196. 

  1196. 		my %microstart=();
    1197. 

  1197. 		my %microend=();
    1198. 

  1198. 	
    1199. 

  1199. 		my @sields = split(/\t/,$sine);
    1200. 

  1200. 	
    1201. 

  1201. 		my $key = ();
    1202. 

  1202. 	
    1203. 

  1203. 		if ($sine =~ /^>[A-Za-z0-9]+\s+([0-9]+)\s+([a-zA-Z0-9]+)\s([a-zA-Z]+[0-9a-zA-Z]+)\s([0-9]+)\s([0-9]+)\s/ ) {
    1204. 

  1204. 			$key = join("\t",$1, $2, $3, $4, $5);
    1205. 

  1205. 		#	print $key, "<-<-<-<-<-<-<-<\n";		
    1206. 

  1206. 		}
    1207. 

  1207. 		else{
    1208. 

  1208. 		}
    1209. 

  1209. 	
    1210. 

  1210. 		if (exists $micros{$key}){
    1211. 

  1211. 			$linecount++;	
    1212. 

  1212. 			my @microstring = @{$micros{$key}};
    1213. 

  1213. 			my @tempmicrostring = @{$micros{$key}};
    1214. 

  1214. 			
    1215. 

  1215. 			foreach my $line (@tempmicrostring){
    1216. 

  1216. 				my @fields = split(/\t/,$line);
    1217. 

  1217. 				my $start = $fields[$startcord];
    1218. 

  1218. 				my $end = $fields[$endcord];
    1219. 

  1219. 				push (@{$microstart{$start}},$line);
    1220. 

  1220. 				push (@{$microend{$end}},$line);	
    1221. 

  1221. 			}
    1222. 

  1222. 			my $firstflag = 'down';
    1223. 

  1223. 			while( my $line =shift(@microstring)){
    1224. 

  1224. 	#			print "-----------\nline = $line ";
    1225. 

  1225. 				chomp $line;
    1226. 

  1226. 				my @fields = split(/\t/,$line);
    1227. 

  1227. 				my $start = $fields[$startcord];
    1228. 

  1228. 				my $end = $fields[$endcord];
    1229. 

  1229. 				my $startmicro = $line;
    1230. 

  1230. 				my $endmicro = $line;
    1231. 

  1231. 	
    1232. 

  1232. 			#	print "fields=@fields, start = $start end=$end, startcord=$startcord, endcord=$endcord\n";
    1233. 

  1233. 	
    1234. 

  1234. 				delete ($microstart{$start});
    1235. 

  1235. 				delete ($microend{$end});
    1236. 

  1236. 				my $flag = 'down';	
    1237. 

  1237. 				my $startflag = 'down';
    1238. 

  1238. 				my $endflag = 'down';
    1239. 

  1239. 				my $prestart = $start - $distance;
    1240. 

  1240. 				my $postend = $end + $distance;
    1241. 

  1241. 				my @compoundlines = ();
    1242. 

  1242. 				my %compoundhash = ();
    1243. 

  1243. 				push (@compoundlines, $line);
    1244. 

  1244. 				push (@{$compoundhash{$line}},$line);
    1245. 

  1245. 				my $startrank = 1;
    1246. 

  1246. 				my $endrank = 1;
    1247. 

  1247. 				
    1248. 

  1248. 				while( ($startflag eq "down") || ($endflag eq "down") ){
    1249. 

  1249. 				if ((($prestart < 0) && $firstflag eq "up") || (($postend > length($sields[$sequencepos])) && $firstflag eq "up") ) {
    1250. 

  1250. #					print "coming to the end of sequence,prestart = $prestart &  post end = $postend and sequence length =", length($sields[$sequencepos])," so exiting\n";
    1251. 

  1251. 					last;
    1252. 

  1252. 				}
    1253. 

  1253. 					
    1254. 

  1254. 				$firstflag = "up";
    1255. 

  1255. 				if ($startflag eq "down"){		
    1256. 

  1256. 					for my $i ($prestart ... $start){
    1257. 

  1257. 					
    1258. 

  1258. 						if(exists $microend{$i}){	
    1259. 

  1259. 							chomp $microend{$i}[0];
    1260. 

  1260. 							if(exists $compoundhash{$microend{$i}[0]}) {next;}
    1261. 

  1261. 	#						print "sending from microend $startmicro, $microend{$i}[0] |||\n";
    1262. 

  1262. 							if (identityMatch_thresholdCheck($startmicro, $microend{$i}[0], $startrank) eq "proceed"){
    1263. 

  1263. 								push(@compoundlines, $microend{$i}[0]);
    1264. 

  1264. 	#							print "accepted\n";
    1265. 

  1265. 								my @tields = split(/\t/,$microend{$i}[0]);
    1266. 

  1266. 								$startmicro = $microend{$i}[0];
    1267. 

  1267. 								chomp $startmicro;
    1268. 

  1268. 								$start = $tields[$startcord];
    1269. 

  1269. 								$flag = 'down';
    1270. 

  1270. 								$startrank++;
    1271. 

  1271. 	#							print "startcompund = $microend{$i}[0]\n";
    1272. 

  1272. 								delete $microend{$i};
    1273. 

  1273. 								delete $microstart{$start};
    1274. 

  1274. 								$startflag = 'down';
    1275. 

  1275. 								$prestart = $start - $distance;
    1276. 

  1276. 								last;
    1277. 

  1277. 							}	
    1278. 

  1278. 							else{
    1279. 

  1279. 								$flag = 'up';
    1280. 

  1280. 								$startflag = 'up';							
    1281. 

  1281. 							}
    1282. 

  1282. 						}
    1283. 

  1283. 						else{
    1284. 

  1284. 							$flag = 'up';
    1285. 

  1285. 							$startflag = 'up';
    1286. 

  1286. 						}
    1287. 

  1287. 					}
    1288. 

  1288. 				}
    1289. 

  1289. 					
    1290. 

  1290. 				$endrank = $startrank;
    1291. 

  1291. 				
    1292. 

  1292. 				if ($endflag eq "down"){				
    1293. 

  1293. 					for my $i ($end ... $postend){
    1294. 

  1294. 					
    1295. 

  1295. 						if(exists $microstart{$i} ){
    1296. 

  1296. 							chomp $microstart{$i}[0];
    1297. 

  1297. 							if(exists $compoundhash{$microstart{$i}[0]}) {next;}	
    1298. 

  1298. 	#						print "sending from microstart $endmicro, $microstart{$i}[0] |||\n";
    1299. 

  1299. 	
    1300. 

  1300. 							if(identityMatch_thresholdCheck($endmicro,$microstart{$i}[0], $endrank) eq "proceed"){
    1301. 

  1301. 								push(@compoundlines, $microstart{$i}[0]);
    1302. 

  1302. 	#								print "accepted\n";
    1303. 

  1303. 								my @tields = split(/\t/,$microstart{$i}[0]);
    1304. 

  1304. 								$end = $tields[$endcord]-0;
    1305. 

  1305. 								$endmicro = $microstart{$i}[0];
    1306. 

  1306. 								$endrank++;
    1307. 

  1307. 								chomp $endmicro;
    1308. 

  1308. 								$flag = 'down';
    1309. 

  1309. 	#							print "endcompund = $microstart{$i}[0]\n";
    1310. 

  1310. 								delete $microstart{$i};
    1311. 

  1311. 								delete $microend{$end};
    1312. 

  1312. 								shift @microstring;
    1313. 

  1313. 								$postend = $end + $distance;
    1314. 

  1314. 								$endflag = 'down';
    1315. 

  1315. 								last;								
    1316. 

  1316. 							}
    1317. 

  1317. 							else{
    1318. 

  1318. 								$flag = 'up';
    1319. 

  1319. 								$endflag = 'up';							
    1320. 

  1320. 							}
    1321. 

  1321. 						}
    1322. 

  1322. 						else{
    1323. 

  1323. 							$flag = 'up';
    1324. 

  1324. 							$endflag = 'up';
    1325. 

  1325. 						}
    1326. 

  1326. 					}
    1327. 

  1327. 				}
    1328. 

  1328. 	#			print "for next turn, flag status: startflag = $startflag and endflag = $endflag \n";
    1329. 

  1329. 			} 														#end while( $flag eq "down")
    1330. 

  1330. 	#			print "compoundlines = @compoundlines \n";
    1331. 

  1331. 			if (scalar (@compoundlines) == 1){
    1332. 

  1332. 				print OUT $line,"\n";			
    1333. 

  1333. 			}
    1334. 

  1334. 			if (scalar (@compoundlines) > 1){
    1335. 

  1335. 				my $compoundline = compoundify(\@compoundlines, $sields[$sequencepos]);
    1336. 

  1336. 	#				print $compoundline,"\n";
    1337. 

  1337. 				print OUT2 $compoundline,"\n";
    1338. 

  1338. 			}
    1339. 

  1339. 			} #end foreach my $line (@microstring){
    1340. 

  1340. 		}	#if (exists $micros{$key}){
    1341. 

  1341. 	
    1342. 

  1342. 	
    1343. 

  1343. 	}
    1344. 

  1344. 	
    1345. 

  1345. 	close OUT;
    1346. 

  1346. 	close OUT2;
    1347. 

  1347. }
    1348. 

  1348. 

  1349. 

  1350. #------------------------------------------------------------------------------------------------
    1351. 

  1351. sub compoundify{
    1352. 

  1352. 	my ($compoundlines, $sequence)  = @_;
    1353. 

  1353. #	print "\nfound to compound : @$compoundlines and$sequence \n";
    1354. 

  1354. 	my $noOfComps = @$compoundlines;
    1355. 

  1355. #	print "Number of elements in hash is $noOfComps\n";
    1356. 

  1356. 	my @starts;
    1357. 

  1357. 	my @ends;
    1358. 

  1358. 	foreach my $line (@$compoundlines){
    1359. 

  1359. #		print "compoundify.. line = $line \n";
    1360. 

  1360. 		chomp $line;
    1361. 

  1361. 		my @fields = split(/\t/,$line);
    1362. 

  1362. 		my $start = $fields[$startcord];
    1363. 

  1363. 		my $end = $fields[$endcord];
    1364. 

  1364. 	#	print "start = $start, end = $end \n";
    1365. 

  1365. 		push(@starts, $start);
    1366. 

  1366. 		push(@ends,$end);		
    1367. 

  1367. 	}
    1368. 

  1368. 	my @temp = @$compoundlines;
    1369. 

  1369. 	my $startline=$temp[0];
    1370. 

  1370. 	my @mields  = split(/\t/,$startline);
    1371. 

  1371. 	my $startcoord = $mields[$startcord];
    1372. 

  1372. 	my $startgapsign=$mields[$endcord];
    1373. 

  1373. 	my @startsorted = sort { $a <=> $b } @starts;
    1374. 

  1374. 	my @endsorted = sort { $a <=> $b } @ends;
    1375. 

  1375. 	my @intervals;
    1376. 

  1376. 	for my $end (0 ... (scalar(@endsorted)-2)){
    1377. 

  1377. 		my $interval = substr($sequence,($endsorted[$end]+1),(($startsorted[$end+1])-($endsorted[$end])-1));
    1378. 

  1378. 		push(@intervals,$interval);
    1379. 

  1379. 	#	print "interval = $interval =\n";
    1380. 

  1380. 	#	print "substr(sequence,($endsorted[$end]+1),(($startsorted[$end+1])-($endsorted[$end])-1))\n";
    1381. 

  1381. 	}
    1382. 

  1382. 	push(@intervals,"");
    1383. 

  1383. 	my $compoundmicrosat=();
    1384. 

  1384. 	my $multiunit="";
    1385. 

  1385. 	foreach my $line (@$compoundlines){
    1386. 

  1386. 		my @fields = split(/\t/,$line);
    1387. 

  1387. 		my $component="[".$fields[$microsatcord]."]".shift(@intervals);
    1388. 

  1388. 		$compoundmicrosat=$compoundmicrosat.$component;
    1389. 

  1389. 		$multiunit=$multiunit."[".$fields[$motifcord]."]";
    1390. 

  1390. #		print "multiunit = $multiunit\n";
    1391. 

  1391. 	}
    1392. 

  1392. 	my $compoundcopy = $compoundmicrosat;
    1393. 

  1393. 	$compoundcopy =~ s/\[|\]//g;
    1394. 

  1394. 	my $compoundlength = $mields[$startcord] + length($compoundcopy) - 1;
    1395. 

  1395. 	
    1396. 

  1396. 	
    1397. 

  1397. 	my $compoundline = join("\t",(@mields[0 ... $infocord], "compound",@mields[$startcord ... $startcord+1],$compoundlength,$compoundmicrosat, $multiunit));
    1398. 

  1398. 	return $compoundline;
    1399. 

  1399. }	
    1400. 

  1400. 

  1401. #------------------------------------------------------------------------------------------------
    1402. 

  1402. 

  1403. sub identityMatch_thresholdCheck{
    1404. 

  1404. 	my $line1 = $_[0];
    1405. 

  1405. 	my $line2 = $_[1];
    1406. 

  1406. 	my $rank = $_[2];
    1407. 

  1407. 	my @lields1 = split(/\t/,$line1);
    1408. 

  1408. 	my @lields2 = split(/\t/,$line2);
    1409. 

  1409. #	print "recieved $line1 && $line2\n motif comparison: ", length($lields1[$motifcord])," : ",length($lields2[$motifcord]),"\n";
    1410. 

  1410. 	
    1411. 

  1411. 	if (length($lields1[$motifcord]) == length($lields2[$motifcord])){
    1412. 

  1412. 		my $probe = $lields1[$motifcord].$lields1[$motifcord];
    1413. 

  1413. 			#print "$probe :: $lields2[$motifcord]\n";
    1414. 

  1414. 		return "proceed" if $probe =~ /$lields2[$motifcord]/;
    1415. 

  1415. 			#print "line recieved\n";
    1416. 

  1416. 		if ($rank ==1){
    1417. 

  1417. 			return "proceed" if thresholdCheck($line1) eq "proceed" && thresholdCheck($line2) eq "proceed";
    1418. 

  1418. 		}
    1419. 

  1419. 		else {
    1420. 

  1420. 			return "proceed" if thresholdCheck($line2) eq "proceed";
    1421. 

  1421. 			return "stop";
    1422. 

  1422. 		}
    1423. 

  1423. 	}
    1424. 

  1424. 	else{
    1425. 

  1425. 		if ($rank ==1){
    1426. 

  1426. 			return "proceed" if thresholdCheck($line1) eq "proceed" && thresholdCheck($line2) eq "proceed";
    1427. 

  1427. 		}
    1428. 

  1428. 		else {
    1429. 

  1429. 			return "proceed" if thresholdCheck($line2) eq "proceed";
    1430. 

  1430. 			return "stop";
    1431. 

  1431. 		}
    1432. 

  1432. 	}
    1433. 

  1433. 	return "stop";
    1434. 

  1434. }
    1435. 

  1435. #------------------------------------------------------------------------------------------------
    1436. 

  1436. 

  1437. sub thresholdCheck{
    1438. 

  1438. 	my @checkthresholds=(0,@thresholds);
    1439. 

  1439. 	#print "IN thresholdCheck: @_\n";
    1440. 

  1440. 	my $line = $_[0];
    1441. 

  1441. 	my @lields = split(/\t/,$line);
    1442. 

  1442. 	return "proceed" if length($lields[$microsatcord]) >= $checkthresholds[length($lields[$motifcord])];
    1443. 

  1443. 	return "stop";
    1444. 

  1444. }
    1445. 

  1445. #xxxxxxxxxxxxxx multiSpecies_compound_microsat_hunter3 xxxxxxxxxxxxxx  multiSpecies_compound_microsat_hunter3 xxxxxxxxxxxxxx  multiSpecies_compound_microsat_hunter3 xxxxxxxxxxxxxx 
    1446. 

  1446. 

  1447. 

  1448. #xxxxxxxxxxxxxx multiSpecies_filtering_interrupted_microsats xxxxxxxxxxxxxx  multiSpecies_filtering_interrupted_microsats xxxxxxxxxxxxxx  multiSpecies_filtering_interrupted_microsats xxxxxxxxxxxxxx 
    1449. 

  1449. 

  1450. sub multiSpecies_filtering_interrupted_microsats{
    1451. 

  1451. #	print "IN multiSpecies_filtering_interrupted_microsats: @_\n";
    1452. 

  1452. 	my $unfiltered = $_[0];
    1453. 

  1453. 	my $filtered = $_[1];
    1454. 

  1454. 	my $residue = $_[2];
    1455. 

  1455. 	my $no_of_species = $_[5];
    1456. 

  1456. 	open(UNF,"<$unfiltered") or die "Cannot open file $unfiltered: $!";
    1457. 

  1457. 	open(FIL,">$filtered") or die "Cannot open file $filtered: $!";
    1458. 

  1458. 	open(RES,">$residue") or die "Cannot open file $residue: $!";
    1459. 

  1459. 		
    1460. 

  1460. 	$infocord = 2 + (4*$no_of_species) - 1;
    1461. 

  1461. 	$startcord = 2 + (4*$no_of_species) + 2 - 1;
    1462. 

  1462. 	$strandcord = 2 + (4*$no_of_species) + 3 - 1;
    1463. 

  1463. 	$endcord = 2 + (4*$no_of_species) + 4 - 1;
    1464. 

  1464. 	$microsatcord = 2 + (4*$no_of_species) + 5 - 1;
    1465. 

  1465. 	$motifcord = 2 + (4*$no_of_species) + 6 - 1;
    1466. 

  1466. 	
    1467. 

  1467. 	
    1468. 

  1468. 	my @sub_thresholds = (0);
    1469. 

  1469. 	
    1470. 

  1470. 	push(@sub_thresholds, split(/_/,$_[3]));
    1471. 

  1471. 	my @thresholds = (0);
    1472. 

  1472. 	
    1473. 

  1473. 	push(@thresholds, split(/_/,$_[4]));
    1474. 

  1474. 	
    1475. 

  1475. 	while (my $line = <UNF>) {
    1476. 

  1476. 		next if $line !~ /[a-z]/;
    1477. 

  1477. 		#print $line;
    1478. 

  1478. 		chomp $line;
    1479. 

  1479. 		my @fields = split(/\t/,$line);
    1480. 

  1480. 		my $motif = $fields[$motifcord];
    1481. 

  1481. 		my $realmotif = $motif;
    1482. 

  1482. 		#print "motif = $motif\n";
    1483. 

  1483. 		if ($motif =~ /^\[/){
    1484. 

  1484. 			$motif =~ s/^\[//g;
    1485. 

  1485. 			my @motifs = split(/\]/,$motif);
    1486. 

  1486. 			$realmotif = $motifs[0];
    1487. 

  1487. 		}
    1488. 

  1488. #		print "realmotif = $realmotif";
    1489. 

  1489. 		my $motif_size = length($realmotif);
    1490. 

  1490. 		
    1491. 

  1491. 		my $microsat = $fields[$microsatcord];
    1492. 

  1492. #		print "microsat = $microsat\n";
    1493. 

  1493. 		$microsat =~ s/^\[|\]$//sg;
    1494. 

  1494. 		my @microsats = split(/\][a-zA-Z|-]*\[/,$microsat);
    1495. 

  1495. 		
    1496. 

  1496. 		$microsat = join("",@microsats);
    1497. 

  1497. 		if (length($microsat) < $thresholds[$motif_size]) {
    1498. 

  1498. 		#	print length($microsat)," < ",$thresholds[$motif_size],"\n"; 
    1499. 

  1499. 			print RES $line,"\n"; next;
    1500. 

  1500. 		}
    1501. 

  1501. 		my @lengths = ();
    1502. 

  1502. 		foreach my $mic (@microsats){
    1503. 

  1503. 			push(@lengths, length($mic));
    1504. 

  1504. 		}
    1505. 

  1505. 		if (largest_microsat(@lengths) < $sub_thresholds[$motif_size]) {
    1506. 

  1506. 	#		print largest_microsat(@lengths)," < ",$sub_thresholds[$motif_size],"\n"; 
    1507. 

  1507. 			print RES $line,"\n"; next;}
    1508. 

  1508. 			else {print FIL $line,"\n"; next;
    1509. 

  1509. 		}
    1510. 

  1510. 	}
    1511. 

  1511. 	close FIL;
    1512. 

  1512. 	close RES;
    1513. 

  1513. 

  1514. }
    1515. 

  1515. 

  1516. sub largest_microsat{
    1517. 

  1517. 	my $counter = 0;
    1518. 

  1518. 	my($max) = shift(@_);
    1519. 

  1519.     foreach my $temp (@_) {
    1520. 

  1520.     	#print "finding largest array: $maxcounter \n";
    1521. 

  1521.     	if($temp > $max){
    1522. 

  1522.         	$max = $temp;
    1523. 

  1523.         }
    1524. 

  1524.     }
    1525. 

  1525.     return($max);
    1526. 

  1526. }
    1527. 

  1527. 

  1528. #xxxxxxxxxxxxxx multiSpecies_filtering_interrupted_microsats xxxxxxxxxxxxxx  multiSpecies_filtering_interrupted_microsats xxxxxxxxxxxxxx  multiSpecies_filtering_interrupted_microsats xxxxxxxxxxxxxx 
    1529. 

  1529. 

  1530. 

  1531. #xxxxxxxxxxxxxx multiSpecies_compound_microsat_analyzer xxxxxxxxxxxxxx  multiSpecies_compound_microsat_analyzer xxxxxxxxxxxxxx  multiSpecies_compound_microsat_analyzer xxxxxxxxxxxxxx…
    1532.

Large files files are truncated, but you can click here to view the full file