blast_parse_1.pl

1	#!/usr/local/bin/perl -w
2	
3	# Parsing BLAST reports with BioPerl's Bio::Tools::BPlite module
4	# WI Biocomputing course - Bioinformatics for Biologists - October 2003
5	
6	# See documentation at http://www.bioperl.org/Core/POD/Bio/Tools/BPlite.html
7	
8	use Bio::Tools::BPlite;
9	
10	# Prompt the user for the file name if it's not an argument
11	# NOTE: BLAST file must be in text (not html) format
12	if (! $ARGV[0])
13	{
14	   print "What is the BLAST file to parse? ";
15	
16	   # Get input and remove the newline character at the end
17	   chomp ($inFile = <STDIN>);
18	}
19	else
20	{
21	   $inFile = $ARGV[0];
22	}
23	
24	$report = new Bio::Tools::BPlite(-file=>"$inFile");
25	
26	# Count hits
27	$i = 1;
28	
29	print "HIT_NUM\tSCORE\tBITS\tPERCENT\tEXPECT\tSUBJECT_NAME\tIDENTITIES\tMATCH_LENGTH\t";
30	print "QUERY_START\tQUERY_END\tSUBJECT_START\tSUBJECT_END\n";
31	print "\n**** DATA FOR QUERY ", $report->query, "****\n\n";
32	
33	while(my $sbjct = $report->nextSbjct)
34	{
35	   while (my $hsp = $sbjct->nextHSP)
36	   {
37	      print "$i\t";
38	      print $hsp->score, "\t",
39	      $hsp->bits,   "\t",
40	      $hsp->percent, "\t",
41	      $hsp->P, "\t",
42	      $hsp->subject->seq_id, "\t",
43	      $hsp->match, "\t",
44	      $hsp->length, "\t";
45	      # $hsp->querySeq,   "\t",     
46	      # $hsp->sbjctSeq,   "\t",     
47	      # $hsp->homologySeq, "\t",     
48	      print $hsp->query->start, "\t",
49	      $hsp->query->end, "\t",
50	      $hsp->subject->start, "\t",
51	      $hsp->subject->end, "\n";
52	   }
53	   $i++;
54	}
55