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
4using namespace seqan;
5
6
7template <typename TFile>
8void 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
60int 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}
SeqAn - Sequence Analysis Library - www.seqan.de
 

Page built @2011/02/08 21:37:01