← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to scripts/taxa/bp_classify_hits_kingdom.pl

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/taxa/bp_classify_hits_kingdom.pl
 
1
#!/usr/bin/perl
 
2
 
 
3
=head1 NAME
 
4
 
 
5
bp_classify_hits_kingdom - classify BLAST hits by taxonomic kingdom
 
6
 
 
7
=head2 USAGE
 
8
 
 
9
bp_classify_hits_kingdom [-i tab_file] [-i second_BLAST_file] [-e evalue_cutoff]
 
10
                      [-t dir_where_TAXONOMY_files_are] [-g gi2taxid] 
 
11
                      [-z PATH_TO_zcat] [-v]
 
12
 
 
13
=head2 DESCRIPTION
 
14
 
 
15
Will print out the taxonomic distribution (at the kingdom level) for a
 
16
set of hits against the NR database.  By default, this script assumes you
 
17
did a search against the protein database (gi_taxid_nuc.dump file).
 
18
 
 
19
This expects BLAST files in tabbed -m9 or -m8 format.  Output with -m
 
20
8 or use blast2table.pl to convert (or fastam9_to_table.PLS if using
 
21
FASTA).
 
22
 
 
23
  Input values:
 
24
    -t/--taxonomy Directory where the taxonomy .dmp files are (from NCBI)
 
25
    -g/--gi       File path of the gi2taxid file (gi_taxid_prot.dmp for proteins
 
26
                  or gi_taxid_nucl.dmp if the search was against a nucleid database)
 
27
    -i/--in       The name of the tab delimited -m8/-m9 output files to process
 
28
    -e/--evalue   Provide an E-value cutoff for hits to be considered
 
29
    -z/--zcat     Path to the 'zcat' executable, can also be 'gunzip -c'
 
30
                  if no zcat on your system.
 
31
  Flags:
 
32
    -v/--verbose  To turn on verbose messages
 
33
    -h/--help     Display this helpful information
 
34
 
 
35
This is intended to be useful starting script, but users may want to
 
36
customize the output and parameters. Note that I am summarizing the
 
37
kingdoms here and Eukaryota not falling into Metazoa, Viridiplantae, or
 
38
Fungi gets grouped into the general superkingdom Eukaryota for simplicity.
 
39
There are comments in the code directing you to where changes can be made
 
40
if you wanted to display hits by phylum for example.  Note that you must
 
41
wipe out the cache file 'gi2class' that is created in your directory after
 
42
making these changes.
 
43
 
 
44
=head2 AUTHOR
 
45
 
 
46
Jason Stajich jason_at_bioperl_dot_org
 
47
 
 
48
=cut
 
49
 
 
50
use strict;
 
51
use warnings;
 
52
use Bio::DB::Taxonomy;
 
53
use DBI;
 
54
use Env;
 
55
use File::Spec;
 
56
use vars qw($SEP);
 
57
my $DEBUG = 0;
 
58
use Getopt::Long;
 
59
$SEP = '_';
 
60
 
 
61
my $evalue_filter = 1e-3;
 
62
my @files;
 
63
my $zcat = 'zcat'; # or gunzip -c 
 
64
my $prefix = File::Spec->catfile($HOME,'taxonomy');
 
65
my $gi2taxidfile = "$prefix/gi_taxid_prot.dmp.gz";
 
66
my $force = 0; # don't use the cached gi2taxid file
 
67
GetOptions(
 
68
       'v|verbose|debug' => \$DEBUG,
 
69
       'force!'          => \$force,
 
70
       'z|zcat:s'        => \$zcat,
 
71
       'i|in:s'          => \@files,
 
72
       'e|evalue:f'      => \$evalue_filter,
 
73
       't|taxonomy:s'    => \$prefix,
 
74
       'g|gi|gi2taxid:s' => \$gi2taxidfile,
 
75
       'h|help'          => sub { system('perldoc', $0); exit },
 
76
       );
 
77
 
 
78
# insure idx location is created
 
79
mkdir(File::Spec->catfile($prefix,'idx')) 
 
80
    unless -d File::Spec->catfile($prefix,'idx');
 
81
 
 
82
# these files came from ftp://ftp.ncbi.nih.gov/pub/taxonomy
 
83
my $taxdb = Bio::DB::Taxonomy->new
 
84
    (-source => 'flatfile',
 
85
     -directory => File::Spec->catfile
 
86
     ($prefix, 'idx'), 
 
87
     -nodesfile => File::Spec->catfile($prefix,'nodes.dmp'),
 
88
     -namesfile => File::Spec->catfile($prefix,'names.dmp')
 
89
     );
 
90
my %query;
 
91
 
 
92
my (%taxid4gi,%gi2node);
 
93
my $dbh = tie(%gi2node, 'DB_File', 'gi2class');
 
94
my $giidxfile = File::Spec->catfile($prefix,'idx','gi2taxid');
 
95
my $done = -f $giidxfile;
 
96
$done = 0 if $force;
 
97
my $dbh2 = $dbh = DBI->connect("dbi:SQLite:dbname=$giidxfile","","");
 
98
if( ! $done ) {
 
99
    $dbh2->do("CREATE TABLE gi2taxid ( gi integer PRIMARY KEY,
 
100
                      taxid integer NOT NULL)");
 
101
    $dbh2->{AutoCommit} = 0;
 
102
    my $fh;
 
103
    # this file came from ftp://ftp.ncbi.nih.gov/pub/taxonomy
 
104
    # I'm interested in protein hits therefor _prot file.
 
105
    if (not -f $gi2taxidfile) {
 
106
        die "Error: File $gi2taxidfile does not exist\n";
 
107
    }
 
108
    if( $gi2taxidfile =~ /\.gz$/ ) {
 
109
        open($fh, "$zcat $gi2taxidfile |" ) || die "$zcat $gi2taxidfile: $!";
 
110
    } else {
 
111
        open($fh, $gi2taxidfile ) || die "Error: could not read file $gi2taxidfile: $!";
 
112
    }
 
113
    my $i = 0;
 
114
    my $sth = $dbh2->prepare("INSERT INTO gi2taxid (gi,taxid) VALUES (?,?)");
 
115
 
 
116
    while(<$fh>) {
 
117
        my ($gi,$taxid) = split;
 
118
        $sth->execute($gi,$taxid);
 
119
        $i++;
 
120
        if( $i % 500000 == 0 ) {
 
121
            $dbh->commit;
 
122
            warn("$i\n") if $DEBUG;
 
123
        } 
 
124
    }
 
125
    $dbh->commit;
 
126
    $sth->finish;
 
127
}
 
128
 
 
129
for my $file ( @files ) {
 
130
    warn("$file\n");
 
131
    my $gz;
 
132
    if( $file =~ /\.gz$/) {
 
133
        $gz = 1;
 
134
    }
 
135
    my ($spname) = split(/\./,$file); 
 
136
    my ($fh,$i);
 
137
    if( $gz ) {
 
138
        open($fh, "$zcat $file |")  || die "$zcat $file: $!";
 
139
    } else {
 
140
        open($fh, $file) || die "$file: $!";
 
141
    }
 
142
    my $sth = $dbh->prepare("SELECT taxid from gi2taxid WHERE gi=?");
 
143
    while(<$fh>) {
 
144
        next if /^\#/;
 
145
        my ($qname,$hname,$pid,$qaln,$mismatch,$gaps,
 
146
            $qstart,$qend,$hstart,$hend,
 
147
            $evalue,$bits,$score) = split(/\t/,$_);        
 
148
        next if( $evalue > $evalue_filter );
 
149
        if( ! exists $query{$spname}->{$qname} ) {
 
150
            $query{$spname}->{$qname} = {};
 
151
        }
 
152
        
 
153
        if( $hname =~ /gi\|(\d+)/) {
 
154
            my $gi = $1;
 
155
            if( ! $gi2node{$gi} ){ # see if we cached the results from before
 
156
                $sth->execute($gi);
 
157
                my $taxid;
 
158
                $sth->bind_columns(\$taxid);
 
159
                if( ! $sth->fetch ) {
 
160
                    warn("no taxid for $gi\n");
 
161
                    next;
 
162
                }
 
163
                my $node = $taxdb->get_Taxonomy_Node($taxid);
 
164
                if( ! $node ) {
 
165
                    warn("cannot find node for gi=$gi ($hname) (taxid=$taxid)\n");
 
166
                    next;
 
167
                }
 
168
                my $parent = $taxdb->get_Taxonomy_Node($node->parent_id);
 
169
 
 
170
                # THIS IS WHERE THE KINGDOM DECISION IS MADE
 
171
                # DON'T FORGET TO WIPE OUT YOUR CACHE FILE
 
172
                # gi2class after you make changes here
 
173
                while( defined $parent && $parent->node_name ne 'root' ) { 
 
174
                    # this is walking up the taxonomy hierarchy
 
175
                    # can be a little slow, but works...
 
176
                    #warn( "\t",$parent->rank, " ", $parent->node_name, "\n");
 
177
                    # deal with Eubacteria, Archea separate from 
 
178
                    # Metazoa, Fungi, Viriplantae differently
 
179
                    # (everything else Eukaryotic goes in Eukaryota)
 
180
                    if( $parent->rank eq 'kingdom') {
 
181
                        # caching in ... 
 
182
                        ($gi2node{$gi}) = $parent->node_name;
 
183
                        last;
 
184
                    } elsif( $parent->rank eq 'superkingdom' ) {
 
185
                        # caching in ... 
 
186
                        ($gi2node{$gi}) = $parent->node_name;
 
187
                        $gi2node{$gi} =~ s/ \<(bacteria|archaea)\>//g;
 
188
                        last;
 
189
                    }
 
190
                    $parent = $taxdb->get_Taxonomy_Node($parent->parent_id);
 
191
                }
 
192
            } 
 
193
            my ($kingdom) = $gi2node{$gi};
 
194
            #warn("$gi2node{$gi}\n");
 
195
            unless( defined $kingdom && length($kingdom) ) {
 
196
                #warn("no kingdom for $hname\n");
 
197
            } else {
 
198
                $query{$spname}->{$qname}->{$kingdom}++;
 
199
            }
 
200
        } else {
 
201
            warn("no GI in $hname\n");
 
202
        }
 
203
    }
 
204
    last if ( $DEBUG && $i++ > 10000);
 
205
    $sth->finish;
 
206
}
 
207
 
 
208
# print out the taxonomic distribution
 
209
while( my ($sp,$d) = each %query ) {
 
210
    my $total = scalar keys %$d;
 
211
    print "$sp total=$total\n";
 
212
    my %seen;
 
213
    for my $v ( values %$d ) {
 
214
        my $tag = join(",",sort keys %$v );
 
215
        $seen{$tag}++;
 
216
    }
 
217
    for my $t ( sort { $seen{$a} <=> $seen{$b} } keys %seen ) {
 
218
        printf " %-20s\t%d\t%.2f%%\n",
 
219
        $t,$seen{$t}, 100 * $seen{$t} / $total;
 
220
    }
 
221
    print "\n\n";
 
222
}