blast_summary.v2.pl

#!/usr/bin/perl
use strict;
use Bio::SearchIO;

my $usage = '
perl blast_summary.pl <full_path> <blast_file>
';

die $usage unless scalar @ARGV == 2;
my ($full_path, $blast_file) = @ARGV;

# get all the viral read sequences
my %viral_reads_blastn = ();
my %viral_reads_blastx = ();
my %best_e_blastn = (); # viral_read_ID => best_e value for this read in blastn
my %best_e_blastx = (); # viral_read_ID => best_e value for this read in blastx
my @blast_files_blastn = (); # all blastn.out files
my @blast_files_blastx = (); # all blastx.out files
my @unassigned_reads = ();
my @ambiguous_reads = (); #cai added 12/2010

# category => num of sequence assigned to this category by blastn
my %blastn = (
    "Bacteria"   => 0,
    "Fungi"      => 0,
    "Homo"       => 0,
    "Mus"        => 0,
    "Phage"      => 0,
    "Viruses"    => 0,
    "other"      => 0,
    "unassigned" => 0,
    "Ambiguous"  => 0,
);

# category => num of sequence assigned to this category by blastn of Reference genome
my %blastn_RefG = ();
foreach my $key (keys %blastn) {
    $blastn_RefG{$key} = 0;
}

# category => num of sequence assigned to this category by tblastx of viral genome
my %blastx = ();
foreach my $key (keys %blastn) {
    $blastx{$key} = 0;
}

my %num_reads = ();
my %blast_readinfo = ();   # readID => information about this read
my %lineage_blastn = ();   # lineage => [read ID]
my %lineage_gi = ();
my %lineage_blastx = ();   # lineage => [read ID]

my $blast_out = $blast_file;
$blast_out =~ s/\.blastn\.parsed/\.blastn\.out/;
$blast_out = $full_path . "/" . $blast_out;

my $blast_s = $blast_file;
print $blast_s, "\n";
$blast_s =~ s/\.blastn\.parsed/\.blastn\.summary/;
$blast_s = $full_path . "/" . $blast_s;

push @blast_files_blastn, $blast_out;

my $parsed = $full_path . "/" . $blast_file;
collect_information(
    $parsed,
    \%blastn,
    \%viral_reads_blastn,
    \%best_e_blastn,
    \%lineage_blastn,
    \%lineage_gi,
    \%num_reads,
    \@unassigned_reads,
    \@ambiguous_reads
);

get_viral_read_info(\@blast_files_blastn, "blastn", \%viral_reads_blastn, \%best_e_blastn, \%blast_readinfo);

open(OUT, ">$blast_s") or die "Cannot open $blast_s\n";
foreach my $id (keys %blast_readinfo) {
    print OUT $blast_readinfo{$id};
}
print OUT "Finished summary\n";
close OUT;

#####################################################################################
# collect information from parsed file
sub collect_information {
    my ($infile, $category_hash_ref, $viral_reads_hash_ref, $best_e_hash_ref,
        $lineage_hash_ref, $lineage_hash_gi, $num_reads_hash_ref,
        $unassigned_reads_arr_ref, $ambiguous_reads_arr_ref) = @_;

    open(IN, $infile) or die "can not open file $infile!\n";
    while (<IN>) {
        if ($_ =~ /#/) { # skip comment line
            next;
        }
        chomp;
        my ($read_ID, $length, $category, $lineage, $hit_name, $e_value) = split("\t", $_);

        my $gid = 0;
        if ($hit_name =~ /gi\|(\d+)\|/) {
            $gid = $1;
            $lineage_hash_gi->{$gid} = $lineage;
        }

        # Replace Switch/case with a simple counter
        if (exists $category_hash_ref->{$category}) {
            $category_hash_ref->{$category}++;
        } else {
            $category_hash_ref->{"other"}++;
        }

        if (($category eq "Viruses" || $category eq "Bacteria") && $gid != 0) {

            $viral_reads_hash_ref->{$read_ID} = 1;
            $best_e_hash_ref->{$read_ID} = $e_value;

            if (!(defined $lineage_hash_ref->{$gid})) {
                $lineage_hash_ref->{$gid} = [$read_ID];
            } else {
                push @{$lineage_hash_ref->{$gid}}, $read_ID;
            }

            if (defined $num_reads_hash_ref->{$gid}) {
                $num_reads_hash_ref->{$gid}++;
            } else {
                $num_reads_hash_ref->{$gid} = 1;
            }

        } elsif ($category eq "Ambiguous") {
            push @{$ambiguous_reads_arr_ref}, $read_ID;

        } elsif ($category eq "unassigned") {
            push @{$unassigned_reads_arr_ref}, $read_ID;
        }
    }
    close IN;
}

#############################################################################
# get detailed information about each viral read
sub get_viral_read_info {
    my ($report_file_ref, $report_type, $viral_reads_hash_ref, $best_e_hash_ref, $blast_readinfo_hash_ref) = @_;

    foreach my $file (@{$report_file_ref}) {
        my $report = Bio::SearchIO->new(-format => 'blast', -file => $file, -report_type => $report_type);

        while (my $result = $report->next_result) {
            my $read_ID = $result->query_name;

            if (defined $viral_reads_hash_ref->{$read_ID}) {
                my $desc = "";
                my $hit_count = 0;

                while (my $hit = $result->next_hit()) {
                    if ($hit->significance() == $best_e_hash_ref->{$read_ID}) {
                        $hit_count++;
                        # only take first best hit (your original stops at 2)
                        if ($hit_count == 2) {
                            last;
                        }

                        $desc .= $result->query_name() . "\t";
                        $desc .= $result->query_length() . "\t";
                        $desc .= $hit->name() . "\t";
                        $desc .= $hit->length() . "\t";
                        $desc .= $hit->description(60) . "\t";

                        while (my $hsp = $hit->next_hsp()) {
                            $desc .= $hsp->length('hit') . "\t";
                            my $percent_id = sprintf("%4.1f", $hsp->percent_identity());
                            $desc .= $percent_id . "\%\t[";
                            $desc .= $hsp->start('hit') . "\t";
                            $desc .= $hsp->end('hit') . "]\t";
                            $desc .= $hsp->evalue() . "\n";
                            last;
                        }
                    }
                }

                $blast_readinfo_hash_ref->{$read_ID} = $desc;
            }
        }
    }
}