-
Notifications
You must be signed in to change notification settings - Fork 47
/
taxonomy_lookup.pl
executable file
·254 lines (209 loc) · 7.03 KB
/
taxonomy_lookup.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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
#!/usr/bin/perl -w
#
# taxonomy_lookup.pl
#
# This script will parse a SAM/blast file and generate an annotated sam /blast file.
# It will append the family/genus/species name and the taxonomic lineage tree to the end of each hit/line.
# Chiu Laboratory
# University of California, San Francisco
#
# Copyright (C) 2014 Scot Federman - All Rights Reserved
# SURPI has been released under a modified BSD license.
# Please see license file for details.
use strict;
#use diagnostics;
use Time::HiRes qw[gettimeofday tv_interval];
use DBI;
if ( @ARGV != 5 ) {
print "USAGE: taxonomy_lookup.pl <blast_file/sam_file> <file_type:blast/sam> <nucl/prot> <cores> <taxonomy_reference_directory>\n";
exit;
}
my ($inputfile, $filetype, $seq_type, $cores, $database_directory) = @ARGV;
my $sql_taxdb_loc_nucl = "$database_directory/gi_taxid_nucl.db";
my $sql_taxdb_loc_prot = "$database_directory/gi_taxid_prot.db";
my $names_nodes = "$database_directory/names_nodes_scientific.db";
my $begintime = [gettimeofday()];
my $gi_table;
my $sql_taxdb_loc;
my $basef_inputfile = $inputfile;
$basef_inputfile =~ s{\.[^.]+$}{};
# First, extract gi from SAM/BLAST inputfile
if ($filetype eq "sam") {
system ("awk {'print \$3'} $inputfile | awk -F \"|\" {'print \$2'} > $basef_inputfile.gi");
} elsif ($filetype eq "blast") {
system ("awk {'print \$2'} $inputfile | awk -F \"|\" {'print \$2'} > $basef_inputfile.gi");
}
if ($seq_type eq "nucl") {
$gi_table = "GI_Taxa_nucl";
$sql_taxdb_loc = $sql_taxdb_loc_nucl;
}
elsif ($seq_type eq "prot") {
$gi_table = "GI_Taxa_prot";
$sql_taxdb_loc = $sql_taxdb_loc_prot;
}
else {
print "\nImproper database specified.\n\n";
exit;
}
my $extracttime = tv_interval($begintime);
print localtime()."\t$0\tTime to extract gi: $extracttime seconds\n";
#create a unique list of gi, make a hash from that list, use the hash to populate the original nonunique list with tax info
#starting point -> file containing all gi to look up
# 1. uniq the list
# 2. foreach unique gi
# look up tax info
# put tax info into hash of hashes like below
# %taxonomy = (
# 149408158 => {
# family => "Hominidae",
# genus => "Homo",
# species => "sapiens"
# },
# 06292007 => {
# family => "Muridae",
# genus => "Mus",
# species => "Mus musculus",
# },
# );
# 3. lookup taxonomy for original gi within hash & output to .gi.taxonomy file
# sort/uniq gi file
my $startsort = [gettimeofday()];
system ("sort --parallel=$cores -u $basef_inputfile.gi > $basef_inputfile.gi.uniq");
my $sorttime = tv_interval($startsort);
print localtime()."\t$0\tTime to sort -u gi file: $sorttime seconds\n";
# Parallelization can occur at this point in the code. Since the file in now sorted, it can be split into n chunks
# with no overlap.
my %taxonomy;
# prepare database connections and statements
my $db = DBI->connect("dbi:SQLite:dbname=$sql_taxdb_loc", "", "", {RaiseError => 1, AutoCommit => 1}) or die $DBI::errstr;
my $names_nodes_db = DBI->connect("dbi:SQLite:dbname=$names_nodes", "", "", {RaiseError => 1, AutoCommit => 1}) or die $DBI::errstr;
# prepare gi to taxid lookup query
my $taxid_return = $db->prepare("SELECT taxid FROM gi_taxid WHERE gi = ?");
my $taxidrow;
# prepare scientific name query
my $name_return = $names_nodes_db->prepare("SELECT name FROM names WHERE taxid = ?");
my $namerow;
# prepare taxonomy query
my $sth = $names_nodes_db->prepare("SELECT * FROM nodes WHERE taxid = ?");
my @noderow;
open (UNIQGI, "$basef_inputfile.gi.uniq") or die $!;
my $starthash = [gettimeofday()];
while (<UNIQGI>) {
chomp;
my ($family, $genus, $species, $lineage) = ("") x 4;
my $result = taxonomy_fgsl($_, $seq_type);
if ($result =~ /family--(.*?)\t/) {
$family = $1;
}
if ($result =~ /genus--(.*?)\t/) {
$genus = $1;
}
if ($result =~ /species--(.*?)\t/) {
$species = $1;
}
if ($result =~ /lineage--(.*)$/) {
$lineage = $1;
}
$taxonomy{$_}{"family"} = $family;
$taxonomy{$_}{"genus"} = $genus;
$taxonomy{$_}{"species"} = $species;
$taxonomy{$_}{"lineage"} = $lineage;
#provide feedback on hash construction
#my $count = scalar(keys %taxonomy);
#if ($count % 500 == 0) {print "$count\n";}
}
close (UNIQGI);
my $endhash = tv_interval($starthash);
print localtime()."\t$0\tTime to create hash: $endhash seconds\n";
my $starttaxwrite = [gettimeofday()];
open (ALLGI, "$basef_inputfile.gi") or die $!;
open (FINALTAXOUTPUT, ">$basef_inputfile.gi.taxonomy") or die $!;
while (my $gi = <ALLGI>) {
chomp $gi;
print FINALTAXOUTPUT "$gi\t";
print FINALTAXOUTPUT "family--$taxonomy{$gi}{\"family\"}\t";
print FINALTAXOUTPUT "genus--$taxonomy{$gi}{\"genus\"}\t";
print FINALTAXOUTPUT "species--$taxonomy{$gi}{\"species\"}\t";
print FINALTAXOUTPUT "lineage--$taxonomy{$gi}{\"lineage\"}\n";
}
close (ALLGI);
close (FINALTAXOUTPUT);
my $endtaxwrite = tv_interval($starttaxwrite);
print localtime()."\t$0\tTime to write taxonomy file: $endtaxwrite seconds\n";
my $start_annotated_write = [gettimeofday()];
open (TAXONOMY, "$basef_inputfile.gi.taxonomy");
open (SAMFILE, "$inputfile");
open (OUTALL, ">$basef_inputfile.all.annotated");
while (my $sam_line = <SAMFILE>) {
chomp($sam_line);
my $tax_line = <TAXONOMY>;
chomp($tax_line);
print OUTALL "$sam_line\t$tax_line\n";
}
close(OUTALL);
close(SAMFILE);
close (TAXONOMY);
my $end_annotated_write = tv_interval($start_annotated_write);
print localtime()."\t$0\tTime to write annotated file: $end_annotated_write seconds\n";
my $elapsedtime = tv_interval($begintime);
print localtime()."\t$0\tTotal time: $elapsedtime seconds\n";
#remove intermediate files
unlink ("$basef_inputfile".".gi");
unlink ("$basef_inputfile".".gi.uniq");
unlink ("$basef_inputfile".".gi.taxonomy");
exit;
sub trim ($){
my $str = shift;
$str =~ s/^\s+//;
$str =~ s/\s+$//;
return $str;
}
sub taxonomy_fgsl {
my ($gi, $seq_type) = @_;
my $lineage = "";
my $name;
my $gi_count = 0;
my %rank_to_print;
my $taxonomy;
$rank_to_print{family} = "1";
$rank_to_print{genus} = "1";
$rank_to_print{species} = "1";
my $begintime = [gettimeofday()];
my $numeric_lineage ="";
# convert gi -> taxid
my $taxid;
$taxid_return->execute(($gi));
if ($taxidrow = $taxid_return->fetchrow_array) {
$taxid = $taxidrow;
}
$taxonomy = "$gi\t";
if ($taxid) {
while ($taxid > 1) {
# Obtain the scientific name corresponding to a taxid
$name_return->execute(($taxid));
# Obtain the parent taxa taxid
# nodes table: taxid - parent_tax id - rank
$sth->execute(($taxid));
if (@noderow = $sth->fetchrow_array) {
(my $tid, my $parent, my $rank) = @noderow;
if ($namerow = $name_return->fetchrow_array) {
$name = trim($namerow);
}
# print the rank if specified on the command line
if (exists($rank_to_print{$rank})) {
$taxonomy = $taxonomy."$rank--$name\t";
}
# Build the taxonomy path
$lineage = "$name;$lineage";
$numeric_lineage = "$tid "."$numeric_lineage ";
$taxid=$parent;
} else {
last;
}
}
}
$taxonomy = $taxonomy."lineage--$lineage";
# print "$numeric_lineage" if $test;
my $elapsedtime = tv_interval($begintime);
return $taxonomy;
}