Mercurial > repos > bgruening > cmsearch_deoverlap

comparison cmsearch-deoverlap.pl @ 0:571e644de888 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/rna_tools/cmsearch_deoverlap commit 675c458c910bb85bd472e25209fa1df43179efb6
author bgruening
date Sat, 19 Aug 2023 05:25:46 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:571e644de888
1 #!/usr/bin/env perl
2 #
3 # cmsearch-deoverlap.pl: remove lower scoring overlaps from cmsearch
4 # --tblout files.
5 #
6 # EPN, Mon May 8 08:41:36 2017
7 #
8 # The EMG Team took a SNAPSHOT of this script on Thurs July 13 2017 from src:
9 # https://raw.githubusercontent.com/nawrockie/cmsearch_tblout_deoverlap/master/cmsearch-deoverlap.pl
10 #
11 #
12 # Without --maxkeep this script will exactly reproduce how cmscan removes overlapping hits.
13 # With --maxkeep, it won't. Here's an example to explain the difference:
14 #
15 ##target name accession query name accession mdl mdl from mdl to seq from seq to strand trunc pass gc bias score E-value inc description of target
16 ##------------------- --------- -------------------- --------- --- -------- -------- -------- -------- ------ ----- ---- ---- ----- ------ --------- --- ---------------------
17 #1. contig--151565 - LSU_rRNA_eukarya RF02543 hmm 632 2329 410 1883 + - 6 0.48 1.1 726.0 1.1e-216 ! -
18 #2. contig--151565 - LSU_rRNA_archaea RF02540 hmm 170 2084 12 1882 + - 6 0.47 3.0 490.0 2.3e-145 ! -
19 #3. contig--151565 - LSU_rRNA_bacteria RF02541 hmm 187 2005 10 1883 + - 6 0.47 5.7 331.4 1.8e-97 ! -
20 #4. contig--151565 - LSU_rRNA_eukarya RF02543 hmm 29 431 12 366 + - 6 0.41 7.8 100.0 6.3e-28 ! -
21 #
22 # hit 1: 410..1883 euk
23 # hit 2: 12..1182 arc
24 # hit 3: 10..1883 bac
25 # hit 4: 12..366 euk
26 #
27 # Without --maxkeep (default behavior) hits 2, 3, and 4 will be removed because for all 3, there is a better scoring hit
28 # that overlaps.
29 #
30 # With --maxkeep only hits 2 and 3 will be removed because only those 2 have
31 # higher scoring hits that overlap with them AND are not removed. If we remove only hits 2 and 3
32 # hits 1 and 4 no longer overlap.
33 #
34 use strict;
35 use warnings;
36 use Getopt::Long;
37 use File::Basename;
38
39 my $in_tblout = ""; # name of input tblout file
40
41 my $usage;
42 $usage = "cmsearch-deoverlap.pl [OPTIONS] <tblout file>\n\tOR\n";
43 $usage .= "cmsearch-deoverlap.pl -l [OPTIONS] <list of tblout files>\n\n";
44 $usage .= "\tOPTIONS:\n";
45 $usage .= "\t\t-l : single command line argument is a list of tblout files, not a single tblout file\n";
46 $usage .= "\t\t-s : sort hits by bit score [default: sort by E-value]\n";
47 $usage .= "\t\t-d : run in debugging mode (prints extra info)\n";
48 $usage .= "\t\t--clanin <s> : only remove overlaps within clans, read clan info from file <s> [default: remove all overlaps]\n";
49 $usage .= "\t\t--maxkeep : keep hits that only overlap with other hits that are not kept [default: remove all hits with higher scoring overlap]\n";
50 $usage .= "\t\t--dirty : keep intermediate files (sorted tblout files)\n\n";
51
52 my $do_listfile = 0; # set to '1' if -l used
53 my $rank_by_score = 0; # set to '1' if -s used, rank by score, not evalues
54 my $do_debug = 0; # set to '1' if -d used
55 my $in_clanin = undef; # defined if --clanin option used
56 my $do_maxkeep = 0; # set to '1' if --maxkeep, only remove hits that have
57 # higher scoring overlap that is not removed
58 my $do_dirty = 0; # set to '1' if --dirty used, keep intermediate files
59
60 &GetOptions( "-l" => \$do_listfile,
61 "-s" => \$rank_by_score,
62 "-d" => \$do_debug,
63 "clanin=s" => \$in_clanin,
64 "maxkeep" => \$do_maxkeep,
65 "keep" => \$do_dirty);
66
67 if(scalar(@ARGV) != 1) { die $usage; }
68 ($in_tblout) = @ARGV;
69
70 my @tblout_file_A = ();
71
72 if($do_listfile) {
73 # read the list file
74 my $list_file = $in_tblout;
75 open(IN, $list_file) || die "ERROR unable to open $list_file for reading";
76 while(my $line = <IN>) {
77 if($line =~ m/\w/ && $line !~ m/^\#/) {
78 chomp $line;
79 if(! -e $line) { die "ERROR file $line read from $list_file does not exist"; }
80 if(! -s $line) { die "ERROR file $line read from $list_file is empty"; }
81 push(@tblout_file_A, $line);
82 }
83 }
84 close(IN);
85 }
86 else {
87 $tblout_file_A[0] = $in_tblout;
88 }
89
90 my %clan_H = (); # key: model name, value: clan name
91 if(defined $in_clanin) {
92 %clan_H = ();
93 parse_claninfo($in_clanin, \%clan_H)
94 }
95
96 my $sort_cmd = undef; # command used to sort the tblout file
97 my $sorted_tblout_file = undef; # sorted tblout file to create, temporarily
98 my $output_file = undef; # name of output file we create
99 my $out_FH = undef; # output file handle
100 my $nkept = undef; # number of sequences kept from a file
101 my $nremoved = undef; # number of sequences removed from a file
102
103 foreach my $tblout_file (@tblout_file_A) {
104 if(! -e $tblout_file) { die "ERROR tblout file $tblout_file does not exist"; }
105 if(! -s $tblout_file) { die "ERROR tblout file $tblout_file is empty"; }
106
107 # sort tblout file by target sequence name
108 $sorted_tblout_file = $tblout_file . ".sort";
109 $sort_cmd = ((defined $rank_by_score) && ($rank_by_score == 1)) ?
110 "grep -v ^\# $tblout_file | sort -k 1,1 -k 15,15rn -k 16,16g > $sorted_tblout_file" :
111 "grep -v ^\# $tblout_file | sort -k 1,1 -k 16,16g -k 15,15rn > $sorted_tblout_file";
112 run_command($sort_cmd, $do_debug);
113 $output_file = basename($tblout_file) . ".deoverlapped";
114 $out_FH = undef;
115 open($out_FH, ">", $output_file) || die "ERROR unable to open $output_file for writing";
116 ($nkept, $nremoved) = parse_sorted_tblout_file($sorted_tblout_file, (defined $in_clanin) ? \%clan_H : undef, $rank_by_score, $do_maxkeep, $do_debug, $out_FH);
117 close $out_FH;
118 if(! $do_dirty) {
119 #unlink $sorted_tblout_file;
120 }
121 printf("Saved %5d hits (%5d removed) to $output_file\n", $nkept, $nremoved)
122 }
123
124 #################################################################
125 # Subroutine : parse_sorted_tblout_file()
126 # Incept: EPN, Mon May 8 08:57:54 2017
127 #
128 # Purpose: Parse a sorted tabular output file, and output
129 # all hits that do not have a higher scoring overlap.
130 #
131 # Arguments:
132 # $sorted_tbl_file: file with sorted tabular search results
133 # $clan_HR: ref to hash of clan info, key is model, value is clan
134 # $rank_by_score: '1' if rank is determined by score, '0' if
135 # determined by E-value
136 # $do_maxkeep: '1' if --maxkeep option used
137 # only remove hits with higher scoring overlaps
138 # *THAT ARE NOT THEMSELVES REMOVED*
139 # $do_debug; '1' if we should print debugging output
140 # $out_FH: file handle to output to
141 #
142 # Returns: Two values:
143 # $nkept: number of hits saved and output
144 # $nremoved: number of hits removed and not output
145 #
146 # Dies: If we see a line of input that is an an unexpected
147 # format
148 #
149 #################################################################
150 sub parse_sorted_tblout_file {
151 my $nargs_expected = 6;
152 my $sub_name = "parse_sorted_tblout_file";
153 if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); }
154
155 my ($sorted_tbl_file, $clan_HR, $rank_by_score, $do_maxkeep, $do_debug, $out_FH) = @_;
156
157 my $prv_target = undef; # target name of previous line
158 my $prv_score = undef; # bit score of previous line
159 my $prv_evalue = undef; # E-value of previous line
160 my $clan = undef; # clan of current model
161 my $nkept = 0; # number of hits kept and output
162 my $nremoved = 0; # number of hits removed and not output
163 my $do_clans = (defined $clan_HR) ? 1 : 0;
164
165 open(IN, $sorted_tbl_file) || die "ERROR unable to open sorted tabular file $sorted_tbl_file for reading";
166
167 my ($target, $model, $domain, $mdlfrom, $mdlto, $seqfrom, $seqto, $strand, $score, $evalue) =
168 (undef, undef, undef, undef, undef, undef, undef, undef, undef, undef);
169
170 my @line_A = (); # array of output lines for kept hits for current sequence
171 my @seqfrom_A = (); # array of seqfroms for kept hits for current sequence
172 my @seqto_A = (); # array of seqtos for kept hits for current sequence
173 my @strand_A = (); # array of strands for kept hits for current sequence
174 my @clan_A = (); # array of clans for kept hits for current sequence
175 my @keepme_A = (); # array of '1' and '0', '1' if we should keep this hit, '0' if it had a higher scoring overlap
176 my $nhits = 0; # number of hits for current sequence (size of all arrays)
177
178 my %already_seen_H = (); # hash, key is sequence name, value is '1' if we have output info for this sequence
179
180 $prv_evalue = 0.;
181 while(my $line = <IN>) {
182 ######################################################
183 # Parse the data on this line, this differs depending
184 # on our annotation method
185 chomp $line;
186 $line =~ s/^\s+//; # remove leading whitespace
187
188 if($line =~ m/^\#/) {
189 die "ERROR, found line that begins with #, input should have these lines removed and be sorted by the first column:$line.";
190 }
191 my @el_A = split(/\s+/, $line);
192
193 ##target name accession query name accession mdl mdl from mdl to seq from seq to strand trunc pass gc bias score E-value inc description of target
194 ##----------------------- --------- -------------------- --------- --- -------- -------- -------- -------- ------ ----- ---- ---- ----- ------ --------- --- ---------------------
195 #lcl|dna_BP444_24.8k:251 - SSU_rRNA_archaea RF01959 hmm 3 1443 2 1436 + - 6 0.53 6.0 1078.9 0 ! -
196 if(scalar(@el_A) < 18) { die "ERROR found less than 18 columns in cmsearch tabular output at line: $line"; }
197 # ($target, $model, $mdlfrom, $mdlto, $seqfrom, $seqto, $strand, $score, $evalue) =
198 ($target, $model, $seqfrom, $seqto, $strand, $score, $evalue) =
199 ($el_A[0], $el_A[2], $el_A[7], $el_A[8], $el_A[9], $el_A[14], $el_A[15]);
200
201 $clan = undef;
202 if(defined $clan_HR) {
203 if(exists $clan_HR->{$model}) {
204 $clan = $clan_HR->{$model};
205 }
206 }
207
208 # make sure we didn't output information for this sequence already
209 if(exists $already_seen_H{$target}) {
210 die "ERROR found line with target $target previously output, did you sort by sequence name?";
211 }
212
213 ##############################################################
214 # Are we now finished with the previous sequence?
215 # Yes, if target sequence we just read is different from it
216 # If yes, output info for it, and re-initialize data structures
217 # for new sequence just read
218 if((defined $prv_target) && ($prv_target ne $target)) {
219 output_one_target($out_FH, \@line_A, \@keepme_A);
220 $already_seen_H{$prv_target} = 1;
221 @line_A = ();
222 @seqfrom_A = ();
223 @seqto_A = ();
224 @strand_A = ();
225 @clan_A = ();
226 @keepme_A = ();
227 $nhits = 0;
228 }
229 else { # this is not a new sequence
230 # make sure that our current score or E-value is less than previous
231 if($rank_by_score && ($score > $prv_score)) {
232 die "ERROR found lines with same target [$target] incorrectly sorted by score, did you sort by sequence name and score?";
233 }
234 elsif((! $rank_by_score) && ($evalue < $prv_evalue)) {
235 die "ERROR found lines with same target [$target] incorrectly sorted by E-value, did you sort by sequence name and E-value?";
236 }
237 }
238 ##############################################################
239
240 # look through all other hits $i on the same strand and see if any of them overlap with this one
241 # if (! $do_maxkeep) we look at all hits
242 # if ( $do_maxkeep) we only look at hits that haven't been removed yet
243 my $keep_me = 1; # set to '0' below if we find a better scoring hit
244 my $overlap_idx = -1;
245 for(my $i = 0; $i < $nhits; $i++) {
246 if(($strand_A[$i] eq $strand) && # same strand
247 (determine_if_clans_match($do_clans, $clan, $clan_A[$i])) && # clans match, or --clanin not used
248 ((! $do_maxkeep) || ($keepme_A[$i] == 1))) { # either --maxkeep not enabled, or this hit is a keepter
249 if(($strand eq "+") && (get_overlap($seqfrom, $seqto, $seqfrom_A[$i], $seqto_A[$i]) > 0)) {
250 $keep_me = 0;
251 $overlap_idx = $i;
252 $i = $nhits; # breaks for loop
253 }
254 elsif(($strand eq "-") && (get_overlap($seqto, $seqfrom, $seqto_A[$i], $seqfrom_A[$i]) > 0)) {
255 $keep_me = 0;
256 $overlap_idx = $i;
257 $i = $nhits; # breaks for loop
258 }
259 }
260 }
261 # add hit to list of hits we have for this sequence
262 $line_A[$nhits] = $line . "\n";
263 $seqfrom_A[$nhits] = $seqfrom;
264 $seqto_A[$nhits] = $seqto;
265 $strand_A[$nhits] = $strand;
266 $clan_A[$nhits] = $clan;
267 if($keep_me) {
268 $nkept++;
269 $keepme_A[$nhits] = 1;
270 }
271 else {
272 $nremoved++;
273 $keepme_A[$nhits] = 0;
274 if($do_debug) {
275 printf("Removing $seqfrom..$seqto, it overlapped with $seqfrom_A[$overlap_idx]..$seqto_A[$overlap_idx]\n");
276 }
277 }
278 $nhits++;
279 $prv_target = $target;
280 $prv_score = $score;
281 $prv_evalue = $evalue;
282 }
283
284 # output data for final sequence
285 output_one_target($out_FH, \@line_A, \@keepme_A);
286
287 # close file handle
288 close(IN);
289
290 return ($nkept, $nremoved);
291 }
292
293 #################################################################
294 # Subroutine : output_one_target()
295 # Incept: EPN, Mon May 8 10:20:40 2017
296 #
297 # Purpose: Output all hits for a target. Overlapping hits
298 # are not included, they've been skipped.
299 #
300 # Arguments:
301 # $out_FH: file handle to output short output to (can be undef to not output short output)
302 # $line_AR: array of lines to output
303 # $keepme_AR: array of '1', '0', $keepme_AR->[$i]==1 indicates we should output $line_AR->[$i]
304 #
305 # Returns: Nothing.
306 #
307 # Dies: Never.
308 #
309 #################################################################
310 sub output_one_target {
311 my $nargs_expected = 3;
312 my $sub_name = "output_one_target";
313 if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); }
314
315 my ($out_FH, $line_AR, $keepme_AR) = @_;
316
317 for(my $i = 0; $i < scalar(@{$line_AR}); $i++) {
318 if($keepme_AR->[$i]) {
319 print $out_FH $line_AR->[$i];
320 }
321 }
322
323 return;
324 }
325
326 #################################################################
327 # Subroutine : determine_if_clans_match()
328 # Incept: EPN, Mon May 8 10:28:41 2017
329 #
330 # Purpose: Given two clan values return true if they
331 # either match or if the --clanin option is
332 # not used.
333 #
334 # Arguments:
335 # $do_clans: '1' if the --clanin option was used
336 # $clan1: clan 1, can be undef
337 # $clan2: clan 2, can be undef
338 #
339 # Returns: '1' if the clans match or if $do_clans is FALSE
340 #
341 # Dies: Never.
342 #
343 #################################################################
344 sub determine_if_clans_match {
345 my $nargs_expected = 3;
346 my $sub_name = "determine_if_clans_match";
347 if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); }
348
349 my ($do_clans, $clan1, $clan2) = @_;
350
351 if(! $do_clans) {
352 return 1;
353 }
354 elsif((! defined $clan1) || (! defined $clan2)) {
355 # 1 or both of $clan1 and $clan2 are undefined, can't be a match
356 return 0;
357 }
358 elsif($clan1 eq $clan2) {
359 # equal clans
360 return 1;
361 }
362 else {
363 return 0;
364 }
365
366 }
367
368 #################################################################
369 # Subroutine: get_overlap()
370 # Incept: EPN, Mon Mar 14 13:47:57 2016 [dnaorg_scripts:dnaorg.pm:getOverlap()]
371 #
372 # Purpose: Calculate number of nucleotides of overlap between
373 # two regions.
374 #
375 # Args:
376 # $start1: start position of hit 1 (must be <= $end1)
377 # $end1: end position of hit 1 (must be >= $end1)
378 # $start2: start position of hit 2 (must be <= $end2)
379 # $end2: end position of hit 2 (must be >= $end2)
380 #
381 # Returns: $noverlap: Number of nucleotides of overlap between hit1 and hit2,
382 # 0 if none
383 #
384 # Dies: if $end1 < $start1 or $end2 < $start2.
385 sub get_overlap {
386 my $sub_name = "get_overlap";
387 my $nargs_exp = 4;
388 if(scalar(@_) != 4) { die "ERROR $sub_name entered with wrong number of input args"; }
389
390 my ($start1, $end1, $start2, $end2) = @_;
391
392 # printf("in $sub_name $start1..$end1 $start2..$end2\n");
393
394 if($start1 > $end1) { die "ERROR in $sub_name start1 > end1 ($start1 > $end1)"; }
395 if($start2 > $end2) { die "ERROR in $sub_name start2 > end2 ($start2 > $end2)"; }
396
397 # Given: $start1 <= $end1 and $start2 <= $end2.
398
399 # Swap if nec so that $start1 <= $start2.
400 if($start1 > $start2) {
401 my $tmp;
402 $tmp = $start1; $start1 = $start2; $start2 = $tmp;
403 $tmp = $end1; $end1 = $end2; $end2 = $tmp;
404 }
405
406 # 3 possible cases:
407 # Case 1. $start1 <= $end1 < $start2 <= $end2 Overlap is 0
408 # Case 2. $start1 <= $start2 <= $end1 < $end2
409 # Case 3. $start1 <= $start2 <= $end2 <= $end1
410 if($end1 < $start2) { return (0); } # case 1
411 if($end1 < $end2) { return ($end1 - $start2 + 1); } # case 2
412 if($end2 <= $end1) { return ($end2 - $start2 + 1); } # case 3
413 die "ERROR in $sub_name, unforeseen case in $start1..$end1 and $start2..$end2";
414
415 return; # NOT REACHED
416 }
417
418 #################################################################
419 # Subroutine: run_command()
420 # Incept: EPN, Mon Dec 19 10:43:45 2016
421 #
422 # Purpose: Runs a command using system() and exits in error
423 # if the command fails. If $be_verbose, outputs
424 # the command to stdout.
425 #
426 # Arguments:
427 # $cmd: command to run, with a "system" command;
428 # $be_verbose: '1' to output command to stdout before we run it, '0' not to
429 #
430 # Returns: nothing
431 #
432 # Dies: if $cmd fails
433 #################################################################
434 sub run_command {
435 my $sub_name = "run_command()";
436 my $nargs_expected = 2;
437 if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); }
438
439 my ($cmd, $be_verbose) = @_;
440
441 if($be_verbose) {
442 print ("Running cmd: $cmd\n");
443 }
444
445 system($cmd);
446
447 if($? != 0) {
448 die "ERROR in $sub_name, the following command failed:\n$cmd\n";
449 }
450
451 return;
452 }
453
454 #################################################################
455 # Subroutine: parse_claninfo()
456 # Incept: EPN, Wed May 10 15:01:07 2017
457 #
458 # Purpose: Parse a claninfo file and fill %{$clan_HR}.
459 #
460 # Arguments:
461 # $claninfo_file: clan info file
462 # $clan_HR: ref to hash of clan info, key: model name, value: clan name
463 #
464 # Returns: nothing
465 #
466 # Dies: if $cmd fails
467 #################################################################
468 sub parse_claninfo {
469 my $sub_name = "parse_claninfo()";
470 my $nargs_expected = 2;
471 if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); }
472
473 my ($claninfo_file, $clan_HR) = @_;
474
475 open(IN, $claninfo_file) || die "ERROR unable to open clan info file";
476
477 %{$clan_HR} = ();
478
479 while(my $line = <IN>) {
480 if($line !~ m/^#/) {
481 chomp $line;
482 my @el_A = split(/\s+/, $line);
483 # first element is clan name, all other elements are model names in that clan
484 #CL00111 SSU_rRNA_bacteria SSU_rRNA_archaea SSU_rRNA_eukarya SSU_rRNA_microsporidia SSU_trypano_mito
485 for(my $i = 1; $i < scalar(@el_A); $i++) {
486 if(exists $clan_H{$el_A[$i]}) {
487 die "ERROR in $sub_name, parsing clan info file $claninfo_file, read model $el_A[$i] more than once";
488 }
489 $clan_HR->{$el_A[$i]} = $el_A[0];
490 }
491 }
492 }
493
494 return;
495 }