-
Notifications
You must be signed in to change notification settings - Fork 0
/
filter_LCA_matches.pl
165 lines (144 loc) · 5.32 KB
/
filter_LCA_matches.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
#!/usr/bin/env perl
# Filters LCA matches by rank (at least species, genus, or family), mean percent identity,
# and/or mean percent query coverage. Uses output of
# retrieve_top_blast_hits_LCA_for_each_sequence.pl as input.
# Input table is output of retrieve_top_blast_hits_LCA_for_each_sequence.pl, with
# column titles (tab-separated):
# - sequence_name
# - LCA_taxon_id
# - LCA_taxon_rank
# - LCA_taxon_species
# - LCA_taxon_genus
# - LCA_taxon_family
# - evalue_of_top_hits
# - lowest_pident_of_top_hits
# - mean_pident_of_top_hits
# - highest_pident_of_top_hits
# - lowest_qcovs_of_top_hits
# - mean_qcovs_of_top_hits
# - highest_qcovs_of_top_hits
# - number_top_hits
# Usage:
# perl filter_LCA_matches.pl
# [output of retrieve_top_blast_hits_LCA_for_each_sequence.pl for one blast search]
# [1 to requires output matches to be classified to at least species level]
# [1 to requires output matches to be classified to at least genus level]
# [1 to requires output matches to be classified to at least family level]
# [minimum mean percent identity] [maximum mean percent identity]
# [minimum mean percent query coverage] [maximum mean percent query coverage]
# [1 to print non-matches instead of matches]
# Prints to console. To print to file, use
# perl filter_LCA_matches.pl
# [output of retrieve_top_blast_hits_LCA_for_each_sequence.pl for one blast search]
# [1 to requires output matches to be classified to at least species level]
# [1 to requires output matches to be classified to at least genus level]
# [1 to requires output matches to be classified to at least family level]
# [minimum mean percent identity] [maximum mean percent identity]
# [minimum mean percent query coverage] [maximum mean percent query coverage]
# [1 to print non-matches instead of matches]
# > [output filtered LCA matches table]
use strict;
use warnings;
my $LCA_matches = $ARGV[0]; # output of retrieve_top_blast_hits_LCA_for_each_sequence.pl
my $classified_to_at_least_species = $ARGV[1]; # if 1, requires output matches to be classified to at least species level
my $classified_to_at_least_genus = $ARGV[2]; # if 1, requires output matches to be classified to at least genus level
my $classified_to_at_least_family = $ARGV[3]; # if 1, requires output matches to be classified to at least family level
my $min_mean_percent_identity = $ARGV[4]; # requires output matches to have at least this mean % identity
my $max_mean_percent_identity = $ARGV[5]; # requires output matches to have at most this mean % identity
my $min_mean_query_coverage = $ARGV[6]; # requires output matches to have at least this mean % query coverage
my $max_mean_query_coverage = $ARGV[7]; # requires output matches to have at most this mean % query coverage
my $reverse = $ARGV[8]; # if 1, prints non-matches instead of matches
my $NO_DATA = "NA";
my $NEWLINE = "\n";
my $DELIMITER = "\t";
# blast LCA table
my $sequence_name_column = 0;
my $LCA_taxon_id_column = 1;
my $LCA_taxon_rank_column = 2;
my $LCA_taxon_species_column = 3;
my $LCA_taxon_genus_column = 4;
my $LCA_taxon_family_column = 5;
my $LCA_taxon_superkingdom_column = 6;
my $evalue_of_top_hits_column = 7;
my $lowest_pident_of_top_hits_column = 8;
my $mean_pident_of_top_hits_column = 9;
my $highest_pident_of_top_hits_column = 10;
my $lowest_qcovs_of_top_hits_column = 11;
my $mean_qcovs_of_top_hits_column = 12;
my $highest_qcovs_of_top_hits_column = 13;
my $number_top_hits_column = 14;
# verifies that all input files exist and are non-empty
if(!$LCA_matches or !-e $LCA_matches or -z $LCA_matches)
{
print STDERR "Error: LCA matches file not provided, does not exist, or empty:\n\t"
.$LCA_matches."\nExiting.\n";
die;
}
# reads in and filters LCA matches
my $first_row = 1;
open LCA_MATCHES, "<$LCA_matches" || die "Could not open $LCA_matches to read\n";
while(<LCA_MATCHES>)
{
chomp;
my $line = $_;
if($line =~ /\S/)
{
if($first_row)
{
# prints header line as is
print $line.$NEWLINE;
$first_row = 0;
}
else
{
# reads in relevant lines in row
my @items = split($DELIMITER, $line);
my $species = $items[$LCA_taxon_species_column];
my $genus = $items[$LCA_taxon_genus_column];
my $family = $items[$LCA_taxon_family_column];
my $mean_percent_identity = $items[$mean_pident_of_top_hits_column];
my $mean_query_coverage = $items[$mean_qcovs_of_top_hits_column];
# determines whether or not this row should be printed
my $row_passes_thresholds = 1;
if($classified_to_at_least_species
and (!$species or $species eq $NO_DATA))
{
$row_passes_thresholds = 0;
}
if($classified_to_at_least_genus
and (!$genus or $genus eq $NO_DATA))
{
$row_passes_thresholds = 0;
}
if($classified_to_at_least_family
and (!$family or $family eq $NO_DATA))
{
$row_passes_thresholds = 0;
}
if($mean_percent_identity < $min_mean_percent_identity)
{
$row_passes_thresholds = 0;
}
if($mean_percent_identity > $max_mean_percent_identity)
{
$row_passes_thresholds = 0;
}
if($mean_query_coverage < $min_mean_query_coverage)
{
$row_passes_thresholds = 0;
}
if($mean_query_coverage > $max_mean_query_coverage)
{
$row_passes_thresholds = 0;
}
# prints line if it passes thresholds
if(!$reverse and $row_passes_thresholds
or $reverse and !$row_passes_thresholds)
{
print $line.$NEWLINE;
}
}
}
}
close LCA_MATCHES;
# December 28, 2022