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 | typedef BlastHsp<BlastN, FullInfo> TBlastHsp;
12 |

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).

13 | typedef BlastReport<TBlastHsp, StreamReport<TFile> > TBlastReport;
14 | typedef typename Hit<TBlastReport>::Type TBlastHit;
15 | typedef typename Iterator<TBlastReport,HitIterator>::Type THitIterator;
16 | typedef typename Iterator<TBlastHit,HspIterator>::Type THspIterator;
17 |
18 |

Our Blast report object

19 | TBlastReport blast;
20 |

Counters

21 | unsigned hspcount = 0;
22 | unsigned hitcount = 0;
23 | unsigned highsignif = 0;
24 |
25 | while(!atEnd(strm,blast))
26 | {

Get the current Blast report (there can be multiple reports in one file).

27 | read(strm,blast,Blast());
28 | ::std::cout << "Query: "<<getQueryName(blast) <<"\n";
29 | ::std::cout << "Database: "<<getDatabaseName(blast) <<"\n\n";
30 |

Iterate over hits

31 | THitIterator hit_it(blast);
32 | for(; !atEnd(strm,hit_it); goNext(strm,hit_it))
33 | {
34 | ++hitcount;
35 | TBlastHit hit = getValue(strm,hit_it);
36 | ::std::cout << " Hit: " <<name(hit) <<"\n\n";
37 |
38 | /// iterate over alignments (HSPs)
39 | THspIterator hsp_it(hit);
40 | for(; !atEnd(strm,hsp_it); goNext(strm,hsp_it))
41 | {
42 | ++hspcount;
43 | TBlastHsp hsp = getValue(strm,hsp_it);
44 |

Do something with the alignment, e.g. output score and length of alignment.

45 | ::std::cout << " Score = " << getBitScore(hsp) << "\n";
46 | ::std::cout << " Length = " << length(hsp) << "\n\n";

Count alignments with highly significant e-values.

47 | if(getEValue(hsp)<0.01)
48 | ++highsignif;
49 |
50 | }
51 | }
52 | }
53 | ::std::cout <<"Total number of Hits: "<< hitcount<<::std::endl;
54 | ::std::cout <<"Total number of HSPs: "<< hspcount<<::std::endl;
55 | ::std::cout <<"Number of highly significant HSPs: "<< highsignif<<::std::endl;
56 |
57 | }
58 |
59 |
60 | int main()
61 | {
62 | ::std::fstream strm;
63 | strm.open("ecoln.out", ::std::ios_base::in | ::std::ios_base::binary);
64 | read_blast_report(strm);
65 | return 0;
66 | }

See
getEValue, Hit, getQueryName, BlastReport, getDatabaseName, getBitScore, BlastHsp
SeqAn - Sequence Analysis Library - www.seqan.de

Page built @2012/03/14 16:34:45