| Example Program Blast Reports Parsing the output of BLAST call. 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 |
|
11 |
| 12 | typedef BlastHsp<BlastN, FullInfo> TBlastHsp;
| 13 |
|
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 |
|
21 | TBlastReport blast;
| 22 |
|
23 | unsigned hspcount = 0;
| 24 | unsigned hitcount = 0;
| 25 | unsigned highsignif = 0;
| 26 |
| 27 | while(!atEnd(strm,blast))
| 28 | {
|
29 | read(strm,blast,Blast());
| 30 | ::std::cout << "Query: "<<getQueryName(blast) <<"\n";
| 31 | ::std::cout << "Database: "<<getDatabaseName(blast) <<"\n\n";
| 32 |
|
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 |
|
47 | ::std::cout << " Score = " << getBitScore(hsp) << "\n";
| 48 | ::std::cout << " Length = " << length(hsp) << "\n\n";
|
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 | }
|
|