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 }