| Example Program Blast Reports Parsing the output of BLAST call. File "blast_report.cpp"
1 | #include <iostream>
2 | #include <seqan/blast.h>
3 |
4 | using namespace std;
5 | using namespace seqan;
6 |
7 |
8 | template <typename TFile>
9 | void read_blast_report(TFile & strm)
10 | {
11 |
First the type of Blast report needs to 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)
12 |
13 | typedef BlastHsp<BlastN, FullInfo> TBlastHsp;
14 |
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)
15 |
16 | typedef BlastReport<TBlastHsp, StreamReport<TFile> > TBlastReport;
17 |
some more types that will be needed
18 | typedef typename Hit<TBlastReport>::Type TBlastHit;
19 | typedef typename Iterator<TBlastReport,HitIterator>::Type THitIterator;
20 | typedef typename Iterator<TBlastHit,HspIterator>::Type THspIterator;
21 |
22 |
our Blast report oject
23 | TBlastReport blast;
24 |
25 | // counters
26 | unsigned hspcount = 0;
27 | unsigned hitcount = 0;
28 | unsigned highsignif = 0;
29 |
30 | while(!atEnd(strm,blast))
31 | {
32 | /// get the current Blast report (there can be multiple reports in one file)
33 | read(strm,blast,Blast());
34 | std::cout << "Query: "<<getQueryName(blast) <<"\n";
35 | std::cout << "Database: "<<getDatabaseName(blast) <<"\n\n";
36 |
37 | /// iterate over hits
38 | THitIterator hit_it(blast);
39 | for(; !atEnd(strm,hit_it); goNext(strm,hit_it))
40 | {
41 | ++hitcount;
42 | TBlastHit hit = getValue(strm,hit_it);
43 | std::cout << " Hit: " <<name(hit) <<"\n\n";
44 |
45 | /// iterate over alignments (HSPs)
46 | THspIterator hsp_it(hit);
47 | for(; !atEnd(strm,hsp_it); goNext(strm,hsp_it))
48 | {
49 | ++hspcount;
50 | TBlastHsp hsp = getValue(strm,hsp_it);
51 |
52 | /// do something with the alignment, e.g.
53 | /// output score and length of alignment
54 | std::cout << " Score = " << getBitScore(hsp) << "\n";
55 | std::cout << " Length = " << length(hsp) << "\n\n";
56 | /// and count alignments with highly significant e-values
57 | if(getEValue(hsp)<0.01)
58 | ++highsignif;
59 |
60 | }
61 | }
62 | }
63 | std::cout <<"Total number of Hits: "<< hitcount<<std::endl;
64 | std::cout <<"Total number of HSPs: "<< hspcount<<std::endl;
65 | std::cout <<"Number of highly significant HSPs: "<< highsignif<<std::endl;
66 |
67 | }
68 |
69 |
70 | int main()
71 | {
72 |
73 | std::fstream strm;
74 | strm.open("ecoln.out",ios_base::in | ios_base::binary);
75 | read_blast_report(strm);
76 |
77 | return 0;
78 | }
See BlastHsp, BlastReport, getBitScore, getDatabaseName, getEValue, getQueryName, Hit SeqAn - Sequence Analysis Library - www.seqan.de