Skip to content

Instantly share code, notes, and snippets.

@radaniba
Created November 29, 2012 16:50
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save radaniba/4170305 to your computer and use it in GitHub Desktop.
Save radaniba/4170305 to your computer and use it in GitHub Desktop.
parse BLAST output files with BioPerl's BPlite
#!/usr/local/bin/perl -w
# Parsing BLAST reports with BioPerl's Bio::Tools::BPlite module
# WI Bioinformatics course - Feb 2002 - Lecture 6
# See documentation at http://www.bioperl.org/Core/POD/Bio/Tools/BPlite.html
use Bio::Tools::BPlite;
# Prompt the user for the file name if it's not an argument
# NOTE: BLAST file must be in text (not html) format
if (! $ARGV[0])
{
print "What is the BLAST file to parse? ";
# Get input and remove the newline character at the end
chomp ($inFile = <STDIN>);
}
else
{
$inFile = $ARGV[0];
}
$report = new Bio::Tools::BPlite(-file=>"$inFile");
# Count hits
$i = 1;
print "HIT_NUM\tSCORE\tBITS\tPERCENT\tEXPECT\tSUBJECT_NAME\tIDENTITIES\tMATCH_LENGTH\t";
print "QUERY_START\tQUERY_END\tSUBJECT_START\tSUBJECT_END\n";
print "\n**** DATA FOR QUERY ", $report->query, "****\n\n";
while(my $sbjct = $report->nextSbjct)
{
while (my $hsp = $sbjct->nextHSP)
{
print "$i\t";
print $hsp->score, "\t",
$hsp->bits, "\t",
$hsp->percent, "\t",
$hsp->P, "\t",
$hsp->subject->seqname, "\t",
$hsp->match, "\t",
$hsp->length, "\t";
# $hsp->querySeq, "\t",
# $hsp->sbjctSeq, "\t",
# $hsp->homologySeq, "\t",
print $hsp->query->start, "\t",
$hsp->query->end, "\t",
$hsp->subject->start, "\t",
$hsp->subject->end, "\n";
}
$i++;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment