| Example Program Blast Reports Parsing the output of BLAST call. File "blast_report.cpp" A tutorial about parsing a blast report.
1 | #include <iostream>
2 | #include <seqan/blast.h>
3 |
4 | using namespace seqan;
5 |
6 |
7 | template <typename TFile>
8 | void read_blast_report(TFile & strm)
9 | {
10 |
First the type of Blast report needs to be specified. In this case, a BlastN report is parsed and for each alignment in the report (i.e. each HSP) we choose to parse all the possible information delivered with the alignment (e.g. scores, percentage of gaps, orientation... see BlastHsp).
11 |
12 | typedef BlastHsp<BlastN, FullInfo> TBlastHsp;
13 |
The StreamReport specialization determines that the alignments are parsed when iterating over them (i.e. only one alignment is stored at a time, as opposed to the StoreReport specialization).
14 |
15 | typedef BlastReport<TBlastHsp, StreamReport<TFile> > TBlastReport;
16 | typedef typename Hit<TBlastReport>::Type TBlastHit;
17 | typedef typename Iterator<TBlastReport,HitIterator>::Type THitIterator;
18 | typedef typename Iterator<TBlastHit,HspIterator>::Type THspIterator;
19 |
20 |
Our Blast report object
21 | TBlastReport blast;
22 |
Counters
23 | unsigned hspcount = 0;
24 | unsigned hitcount = 0;
25 | unsigned highsignif = 0;
26 |
27 | while(!atEnd(strm,blast))
28 | {
Get the current Blast report (there can be multiple reports in one file).
29 | read(strm,blast,Blast());
30 | ::std::cout << "Query: "<<getQueryName(blast) <<"\n";
31 | ::std::cout << "Database: "<<getDatabaseName(blast) <<"\n\n";
32 |
Iterate over hits
33 | THitIterator hit_it(blast);
34 | for(; !atEnd(strm,hit_it); goNext(strm,hit_it))
35 | {
36 | ++hitcount;
37 | TBlastHit hit = getValue(strm,hit_it);
38 | ::std::cout << " Hit: " <<name(hit) <<"\n\n";
39 |
40 | /// iterate over alignments (HSPs)
41 | THspIterator hsp_it(hit);
42 | for(; !atEnd(strm,hsp_it); goNext(strm,hsp_it))
43 | {
44 | ++hspcount;
45 | TBlastHsp hsp = getValue(strm,hsp_it);
46 |
Do something with the alignment, e.g. output score and length of alignment.
47 | ::std::cout << " Score = " << getBitScore(hsp) << "\n";
48 | ::std::cout << " Length = " << length(hsp) << "\n\n";
Count alignments with highly significant e-values.
49 | if(getEValue(hsp)<0.01)
50 | ++highsignif;
51 |
52 | }
53 | }
54 | }
55 | ::std::cout <<"Total number of Hits: "<< hitcount<<::std::endl;
56 | ::std::cout <<"Total number of HSPs: "<< hspcount<<::std::endl;
57 | ::std::cout <<"Number of highly significant HSPs: "<< highsignif<<::std::endl;
58 |
59 | }
60 |
61 |
62 | int main()
63 | {
64 | ::std::fstream strm;
65 | strm.open("ecoln.out", ::std::ios_base::in | ::std::ios_base::binary);
66 | read_blast_report(strm);
67 | return 0;
68 | }
See BlastHsp, BlastReport, getBitScore, getDatabaseName, getEValue, getQueryName, Hit SeqAn - Sequence Analysis Library - www.seqan.de