BLASTn_VIRUSDB_parser.pl.sqlite

################################################################################
##  Author:
##      Guoyan Zhao <gzhao@pathology.wustl.edu>
#*
#* Copyright (C) Washington University in St. Louis Developed by Guoyan Zhao,
#* Washington University in St. Louis, Department of Pathology and Immunology.
#*
###############################################################################

#!/usr/bin/perl -w

use strict;
use Bio::SearchIO;
use Bio::Taxon;
use Bio::DB::Taxonomy;
use Bio::Tree::Tree;
use DBI();
use File::Temp qw/ tempfile tempdir /;

my $Usage = '
This script accepts a BLASTn output file and parse the information

perl script <dir><blast output file> 
<dir> = directory that blast output file resides in, without last "/"
<blast output file> = name of the BLASTx output file

';

die $Usage unless scalar @ARGV == 2;
my ($dir, $blastout) = @ARGV;

################################################################################
# Modify script to use local configuration:
my $database_dir = "/scratch/dwlab/databases/taxdump_20160802";

# open a connection to mysql database
my $dbh_sqlite = DBI->connect("dbi:SQLite:dbname=/scratch/dwlab/databases/taxdump_20160802/vhunter.db","", "",{'RaiseError'=>1}) or die "Unable to connect $DBI::errstr\n";
#my $dbh_mysql = DBI->connect("DBI:mysql:database=vhunter;host=192.168.123.1","kchen", "k123",{'RaiseError'=>1}) or die "Unable to connect $DBI::errstr\n";

#################################################################################
# Don't need to change anything below
# cutoff value for having a significant hit
my $E_cutoff = 1e-5;
my $TMPDIR = tempdir();

my %assignment = ();
my %assignment_NotBestE = ();

# get file name prefix
my $file_name_prefix = $blastout;
$file_name_prefix =~ s/\.blastnv\.out//;

# create ouput file
my $outFile = $dir."/".$file_name_prefix.".blastnv.parsed";
open (OUT, ">$outFile") or die "can not open file $outFile!\n";


# get a Taxon from a Bio::DB::Taxonomy object
my $dbh = Bio::DB::Taxonomy->new(-source => 'flatfile',
		-directory=> "$TMPDIR",
		-nodesfile=> "$database_dir/nodes.dmp",
		-namesfile=> "$database_dir/names.dmp",
		);

my @unassigned = (); # query should be kept for further analysis
my @virusseq = (); 
my @phageseq = (); 
my $total_records = 0;

print "parsing blast output files...\n\n";

my $input_file = $dir."/".$blastout;
my $report = new Bio::SearchIO(-format => 'blast', -file => $input_file, -report_type => 'blastn');

# Go through BLAST reports one by one        
while(my $result = $report->next_result) {# next query output
	$total_records++;
#	print $result->query_name, " \n";
	my $haveHit = 0;
	my $keep_for_next_step = 1;
	%assignment = ();
	%assignment_NotBestE = ();

	# only take the best hits
	my $best_e = 100;
	my $hit_count = 0;
	my $determined = 0;
	while(my $hit = $result->next_hit) {
		# from hit name get hit gi number
		my $hit_name = $hit->name; # gi|num|database|accessionNum|
		my @temp_arr = split(/\|/, $hit_name);
		my $gi = $temp_arr[1];
		if ($temp_arr[2] eq "pdb") { # skip data from pdb database
			next;
		}
		$haveHit = 1;
		$hit_count++;
		if ($hit_count == 1) {
			$best_e = $hit->significance;
		}

		# check whether the hit should be kept for further analysis
		if ($best_e <= $E_cutoff) { # similar to known, need Phylotyped
#			print $result->query_name, " similar to known, output description!\n\n";
#			print "the $hit_count hit, $best_e \n"; 

			if ($hit->significance == $best_e ) { # only get best hits
				# from gi get taxonomy lineage
				my $sth = $dbh_sqlite->prepare("SELECT * FROM gi_taxid_nucl where gi = $gi");
				$sth->execute();
				my $ref = $sth->fetchrow_hashref();
#				print "gi = $ref->{'gi'}, taxid = $ref->{'tax_id'}\n";
				
				$sth->finish();
				my $taxID = $ref->{'tax_id'};
				if ($taxID) { # some gi don't have record in gi_taxid_nucl
#					print "taxID is $taxID\n";
					my $taxon_obj = $dbh->get_taxon(-taxonid => $taxID);
#					print "taxon object is \\$taxon_obj\\\n";

					if (!(defined $taxon_obj)) {
						my $description = "undefined taxon ".$hit->description."\t".$hit->name."\t".$hit->significance;
						$assignment{"other"} = $description;
					}
					else {
						my $tree_function = Bio::Tree::Tree->new();
						my @lineage = $tree_function->get_lineage_nodes($taxon_obj);
						# each lineage node is a Bio::Tree::NodeI object
#						print "hit gi is $gi\n";
#						print "id is ", $taxon_obj->id, "\n";
#						print "rank is ", $taxon_obj->rank, "\n";
#						print "divison is ", $taxon_obj->division, "\n\n";
#						print "lineage is @lineage\n";
						if (scalar @lineage) {	
							$determined = 1;
#							print "PhyloTyped, don't save for further analysis\n";
							&PhyloType(\@lineage,$hit, $dbh, \%assignment);
						}
					}
				}	
				else { # for situations that gi does not have corresponding taxid
#					print $result->query_name, " ", $hit->name, "\n";
#					print "gi = $ref->{'gi'}, taxid = $ref->{'tax_id'}\n";
#					print "hit gi is $gi\n";
					my $description = $hit->description."\t".$hit->name."\t";
#					print $result->query_name, "\t", $desc, "\n";
					$determined = 1;
					$description .= $hit->length()."\t";
					$description .= $hit->description(60)."\t";
					while (my $hsp = $hit->next_hsp()) {
					$description .= $hsp->length('hit')."\t";
					my $percent_id = sprintf("%4.1f", $hsp->percent_identity());
						$description .= $percent_id."\%\t[";
						$description .= $hsp->start('hit')."\t";
						$description .= $hsp->end('hit')."]\t";
#						$description .= $hsp->evalue()."\n";
						last;
					}
					$description .= $hit->significance;

					$assignment{"other"} = $description;
				} 
			}
			elsif ($hit->significance <= $E_cutoff) { # significant but is not the same e value as best hit 
				# from gi get taxonomy lineage
				my $sth = $dbh_sqlite->prepare("SELECT * FROM gi_taxid_nucl where gi = $gi");
				$sth->execute();
				my $ref = $sth->fetchrow_hashref();
#				print "gi = $ref->{'gi'}, taxid = $ref->{'tax_id'}\n";
				$sth->finish();
				my $taxID = $ref->{'tax_id'};
				if ($taxID) { # some gi don't have record in gi_taxid_nucl
#					print "taxID is $taxID\n";
					my $taxon_obj = $dbh->get_taxon(-taxonid => $taxID);

					if (!(defined $taxon_obj)) {
						my $description = "undefined taxon ".$hit->description."\t".$hit->name."\t";
						$description .= $hit->length()."\t";
						$description .= $hit->description(60)."\t";
						while (my $hsp = $hit->next_hsp()) {
						$description .= $hsp->length('hit')."\t";
						my $percent_id = sprintf("%4.1f", $hsp->percent_identity());
							$description .= $percent_id."\%\t[";
							$description .= $hsp->start('hit')."\t";
							$description .= $hsp->end('hit')."]\t";
#							$description .= $hsp->evalue()."\n";
							last;
						}
						$description .= $hit->significance;

						$assignment_NotBestE{"other"} = $description;
					}
					else {
						my $tree_function = Bio::Tree::Tree->new();
						my @lineage = $tree_function->get_lineage_nodes($taxon_obj);
						# each lineage node is a Bio::Tree::NodeI object
#						print "hit gi is $gi\n";
#						print "id is ", $taxon_obj->id, "\n";
#						print "rank is ", $taxon_obj->rank, "\n";
#						print "divison is ", $taxon_obj->division, "\n\n";
#						print "lineage is @lineage\n";
						if (scalar @lineage) {				
							&PhyloType(\@lineage,$hit,  $dbh, \%assignment_NotBestE);
						}

						#################################################################
						# If the sequence also hit any other species with significant e value skip all the rest hits.
						if (((defined $assignment_NotBestE{"Bacteria"}) || (defined $assignment_NotBestE{"Fungi"}) || (defined $assignment_NotBestE{"Homo"}) || (defined $assignment_NotBestE{"Mus"}) || (defined $assignment_NotBestE{"Phage"}) || (defined $assignment_NotBestE{"other"})) ) {
							last;
						}
					}
				}	
			}
		} # finish phylotype for given hit
		else { # e value is not significant 
			if($determined){ # skip the rest hits that are not significant
				last;
			}
			else {
				my $desc = "hit not significant\t".$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;
				}
				$desc .= $hit->significance;

				$assignment{"unassigned"} = $desc;
				last;
			}
		}
	}  # finish all hits

	if (!$haveHit) {
		$assignment{"unassigned"} = "no hit";
	}

#	foreach my $key (keys %assignment) {
#		print "after parsing ", $key, "\t", $assignment{$key},"\n";
#	}

	# remove duplicate assignment
	# If a query is assigned both Homo and Primates, it will be reported as Homo only
	# If a query is assigned a real taxon name and "other" for reason like"other sequences;
	# artificial sequences", or no taxon id in taxon database it will be reported only as 
	# the real taxon name
	my $num_assignment = keys %assignment;
	if ($num_assignment > 1) { # have multiple assignment, can only be virus or phage
		###############################################
		# determine phage hits
		# If a sequence hits virus and  phage, the sequence is assigned to "Phage" category. 
		if (defined $assignment{"Viruses"}) {
			if (defined $assignment{"Phage"})  {
#				print "has phage hits!!!!!!!!!!!!!!\n";
				delete $assignment{"Viruses"};
			}
		}
	}
	elsif ($num_assignment == 1) { # have exactly one assignment
		if (defined $assignment{"Viruses"})  { # it's virus assignment
			if (defined $assignment_NotBestE{"Phage"})  { # but has phage as significant (not best) hit
#				print "has phage hits!!!!!!!!!!!!!!\n";
				$assignment{"Phage"} = $assignment_NotBestE{"Phage"}; # considered to be phage seq
				delete $assignment{"Viruses"};
			}
		}
	}

=head1
	print " now print significant but not best hit assignment:\n";
	foreach my $key (keys %assignment_NotBestE) {
		print "significant but not best hit assignment \\", $key, "\\\t", $assignment_NotBestE{$key},"\n";
	}
	print "end\n\n";

	foreach my $key (keys %assignment) {
		print "after consolidateion ", $key, "\t", $assignment{$key},"\n";
	}
	print "end******************************\n\n";
=cut

	# print out assignment for this query
	my $assign_to_virus=0; 

	foreach my $assign (keys %assignment) {
		if (($assign eq "Viruses") ) {
			push @virusseq, $result->query_name; 
			$assign_to_virus=1; 
		}
		if (($assign eq "Phage") ) {
			push @phageseq, $result->query_name; 
			$assign_to_virus=1; 
		}

	    print OUT $result->query_name, "\t", $result->query_length, "\t", $assign, "\t", $assignment{$assign}, "\n"; 
	}

	if ($assign_to_virus==0) {
		push @unassigned, $result->query_name;
	}

}
print OUT "# Summary: ", scalar @unassigned, " out of $total_records ", (scalar @unassigned)*100/$total_records, "% are unassigned.\n";

close OUT;

# generate a fasta file that contains all the sequences that are saved for downstream analysis 
# read in blast input sequences
my $file = $dir."/".$file_name_prefix.".fasta";
my %seq = &read_FASTA_data($file);

$outFile = $dir."/".$file_name_prefix.".BNVIRUSDB_hit.fa";
open (OUT2, ">$outFile") or die "can not open file $outFile!\n";
foreach my $seq_name (@virusseq) {
	print OUT2 ">$seq_name\n";
	print OUT2 $seq{$seq_name}, "\n";
}
close OUT2;

$outFile = $dir."/".$file_name_prefix.".BNVIRUSDB_phage.fa";
open (OUT2, ">$outFile") or die "can not open file $outFile!\n";
foreach my $seq_name (@phageseq) {
	print OUT2 ">$seq_name\n";
	print OUT2 $seq{$seq_name}, "\n";
}
close OUT2;

$outFile = $dir."/".$file_name_prefix.".BNVIRUSDB_filtered.fa";
open (OUT2, ">$outFile") or die "can not open file $outFile!\n";
foreach my $seq_name (@unassigned) {
        print OUT2 ">$seq_name\n";
        print OUT2 $seq{$seq_name}, "\n";
}
close OUT2;

$dbh_sqlite->disconnect();
#$dbh_mysql->disconnect();

exit;


############################################################################
sub read_FASTA_data () {
	my $fastaFile = shift @_;

	#keep old read seperator and set new read seperator to ">"
	my $oldseperator = $/;
	$/ = ">";
	 
	my %fastaSeq;	 
	open (FastaFile, $fastaFile) or die "Can't Open FASTA file: $fastaFile";

	while (my $line = <FastaFile>){
		# Discard blank lines
		if ($line =~ /^\s*$/) {
			next;
		}
		# discard comment lines
		elsif ($line =~ /^\s*#/) {
			next;
		}
		# discard the first line which only has ">", keep the rest
		elsif ($line ne ">") {
			chomp $line;
			my @rows = ();
			@rows = split (/\n/, $line);	
			my $firstLine = shift @rows;
			my @temp = split (/\s/, $firstLine);
			my $contigName = shift @temp;
			my $contigSeq = join("", @rows);
			$contigSeq =~ s/\s//g; #remove white space
			$fastaSeq{$contigName} = $contigSeq;
		}
    }

    # execute right
#	 foreach my $key (keys %fastaSeq){
#	      print "Here is the key for fasta seq: $key \t $fastaSeq{$key}\n";
#	 }

    #reset the read seperator
    $/ = $oldseperator;
   	close FastaFile;
 
    return %fastaSeq;
}

		
###############################################################################
# subroutine to determine the taxonomy lineage for a given blast hit
sub PhyloType {
	my ($lineage_ref, $hit_ref,  $dbh_taxonomy, $assignment_ref) = @_;
	my $description = "";
	my $node_id; 
	my $obj;
	my $name;
	my $assigned = 0;

	my $Lineage = "";
	for (my $i = 0; $i <= $#$lineage_ref; $i++) { 
		my $temp_node_id = $lineage_ref->[$i]->id;
		my $temp_obj = $dbh_taxonomy->get_taxon(-taxonid=>$temp_node_id);
		my $temp_name = $temp_obj->scientific_name;
		$Lineage .= $temp_name.";";
	}					
   
	# check to see if it is a human sequence
	if (scalar @{$lineage_ref} >= 4) {
		$node_id = $lineage_ref->[3]->id;
		$obj = $dbh_taxonomy->get_taxon(-taxonid=>$node_id);
		$name = $obj->scientific_name;
		if ($name eq "Metazoa") {
			# make assignment
			for (my $i = 0; $i <= $#$lineage_ref; $i++) { 
				my $temp_node_id = $lineage_ref->[$i]->id;
				my $temp_obj = $dbh_taxonomy->get_taxon(-taxonid=>$temp_node_id);
				my $temp_name = $temp_obj->scientific_name;
#				print "name = $temp_name\n";
				if ($temp_name eq "Homo") {
                    if(!defined $assignment_ref->{"Homo"}) {  
#						print "assigned to Homo\n\n";
						$description .= "Homo\t".$hit_ref->name."\t";

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

						$assignment_ref->{"Homo"} = $description; 
					}
					$assigned = 1;
					last; 
				}
			}
			if (!$assigned) {
				for (my $i = 0; $i <= $#$lineage_ref; $i++) { 
					my $temp_node_id = $lineage_ref->[$i]->id;
					my $temp_obj = $dbh_taxonomy->get_taxon(-taxonid=>$temp_node_id);
					my $temp_name = $temp_obj->scientific_name;
#					print "name = $temp_name\n";
	
					if ($temp_name eq "Mus") {
                       if(!defined $assignment_ref->{"Mus"}) {  
#							print "assigned to Mus\n\n";
							$description .= "Mus\t".$hit_ref->name."\t";

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

							$assignment_ref->{"Mus"} = $description; 
						}
						$assigned = 1;
						last; 
					}
				}
			}
			if (!$assigned) {
				if(!defined $assignment_ref->{"other"})  { 
					$description .= $Lineage."\t".$hit_ref->name."\t";
#					print "assigned to other\n\n";

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

					$assignment_ref->{"other"} = $description; 
				}
				$assigned = 1; 
			}
		}
	}

	# check to see if it is bacteria sequence
	if ((scalar @{$lineage_ref} >= 2)&&(!$assigned)) {
		$node_id = $lineage_ref->[1]->id;
		$obj = $dbh_taxonomy->get_taxon(-taxonid=>$node_id);
		$name = $obj->scientific_name;
		if ($name eq "Bacteria") {
			if(!defined $assignment_ref->{"Bacteria"})  {    
				$description = $Lineage."\t".$hit_ref->name."\t";

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

				$assignment_ref->{"Bacteria"} = $description;
			}
			$assigned = 1; 
		}
	}

	# check to see if it is a phage virus sequence
	if (!$assigned) {
		$node_id = $lineage_ref->[0]->id;
		$obj = $dbh_taxonomy->get_taxon(-taxonid=>$node_id);
		$name = $obj->scientific_name;
		if ($name eq "Viruses") {
			for (my $i = 0; $i <= $#$lineage_ref; $i++) { 
				my $temp_node_id = $lineage_ref->[$i]->id;
				my $temp_obj = $dbh_taxonomy->get_taxon(-taxonid=>$temp_node_id);
				my $temp_name = $temp_obj->scientific_name;
				$description .= $temp_name.";";
				if (($temp_name =~ /Lipothrixviridae/i)||($temp_name =~ /Caudovirales/i)||($temp_name =~ /Corticoviridae/i)||($temp_name =~ /Cystoviridae/i)||($temp_name =~ /Inoviridae/i)||($temp_name =~ /Leviviridae/i)||($temp_name =~ /Microviridae/i)||($temp_name =~ /Tectiviridae/i)||($temp_name =~ /phage/i)||($temp_name =~ /Plasmaviridae/i) ) {
#					print "assigned to phage\n\n";
					if(!defined $assignment_ref->{"Phage"}) { # take the first best hit description
						$description = $Lineage."\t".$hit_ref->name."\t";

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

						$assignment_ref->{"Phage"} = $description;
					}
					$assigned = 1;
					last;
				}
			}
		}
	}

	# check to see if it is a virus sequence
	$description = "";
	if (!$assigned) {
		$node_id = $lineage_ref->[0]->id;
		$obj = $dbh_taxonomy->get_taxon(-taxonid=>$node_id);
		$name = $obj->scientific_name;
		if ($name eq "Viruses") {
			if(!defined $assignment_ref->{"Viruses"}) { 
				$description = $Lineage."\t".$hit_ref->name."\t";

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

				$assignment_ref->{"Viruses"} = $description;
			}
			$assigned = 1;
		}
	}

	# check to see if it is a fungi sequence
	if ((scalar @{$lineage_ref} >= 4)&&(!$assigned)) {
		$node_id = $lineage_ref->[3]->id;
		$obj = $dbh->get_taxon(-taxonid=>$node_id);
		$name = $obj->scientific_name;
		if ($name eq "Fungi") {
			if(!defined $assignment_ref->{"Fungi"}) { 
				$description = $Lineage."\t".$hit_ref->name."\t";

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

				$assignment_ref->{"Fungi"} = $description;
			}
			$assigned = 1;
		}
	}

	# if still not assigned, assigned to "other" category
	if (!$assigned) {
		if(!defined $assignment_ref->{"other"}) {
			$description = $Lineage."\t".$hit_ref->name."\t";

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

			$assignment_ref->{"other"} = $description;
		}
		$assigned = 1;
	}
 
	return $assigned;
}