1 #include <iostream>
2 #include <seqan/basic.h>
3 #ifndef STDLIB_VS
4 #include <seqan/blast.h>
5
6 using namespace seqan;
7
main(int argc,char ** argv)8 int main(int argc, char ** argv)
9 {
10 if (argc != 2)
11 {
12 std::cerr << "USAGE: FILE_OUT\n";
13 return 0;
14 }
15
16 typedef String<AminoAcid> TSequence;
17 typedef std::string TId;
18 typedef Gaps<TSequence, ArrayGaps> TGaps;
19 typedef BlastMatch<TGaps, TGaps> TBlastMatch;
20 typedef BlastRecord<TBlastMatch> TBlastRecord;
21 typedef BlastIOContext<Blosum62> TContext;
22
23 StringSet<TSequence> queries;
24 StringSet<TSequence> subjects;
25 StringSet<TId> qIds;
26 StringSet<TId> sIds;
27
28 appendValue(queries, "VAYAQPRKLCYP");
29 appendValue(queries, "NAYPRUTEIN");
30 appendValue(queries, "AVITSFTQ");
31
32 appendValue(subjects,
33 "SSITEEKHIPHKEQDKDAEFLSKEALKTHMTENVLQMDRRAVQDPSTSFLQLLKAKGLLG"
34 "LPDYEVNLADVNSPGFRKVAYAQTKPRRLCFPNGGTRRGSFIMDTAVVVMVSLRYVNIGK"
35 "VIFPGATDVSEGEDEFWAGLPQAYGCLATEFLCIHIAIYSWIHVQSSRYDDMNASVIRAK"
36 "LNLAVITSWTQLIQAEKETI");
37
38 appendValue(subjects,
39 "GATRDSKGNAVITSFTQARLRVYADLLGPYWIILHVIELTGVGNTGQKCTLNHMGTYAVF"
40 "DLKQPPATNDLGLPKPCFIGFDIQNELAIGTVGHSEAVIAAFTQRDRLEERAESKQSLAR"
41 "PVISPKLIAEVSTVLESALNQMYSSLGFYRVERAEDYAQPRKLCVVKKKSFNCLNADIWL"
42 "EYRMEDQKSVPKVFKIMMDD");
43
44 appendValue(qIds, "Query_Numero_Uno with args");
45 appendValue(qIds, "Query_Numero_Dos with args");
46 appendValue(qIds, "Query_Numero_Tres with args");
47
48 appendValue(sIds, "Subject_Numero_Uno");
49 appendValue(sIds, "Subject_Numero_Dos");
50
51 BlastTabularFileOut<TContext> outfile(argv[1]);
52 String<TBlastRecord> records;
53
54 // protein vs protein search is BLASTP
55 context(outfile).blastProgram = BlastProgram::BLASTP;
56
57 // set gap parameters in blast notation
58 setScoreGapOpenBlast(context(outfile).scoringScheme, -11);
59 setScoreGapExtend(context(outfile).scoringScheme, -1);
60 SEQAN_ASSERT(isValid(context(outfile).scoringScheme));
61
62 // set the database properties in the context
63 context(outfile).dbName = "The Foo Database";
64 context(outfile).dbTotalLength = length(concat(subjects));
65 context(outfile).dbNumberOfSeqs = length(subjects);
66
67 writeHeader(outfile); // write file header
68
69 for (unsigned q = 0; q < length(queries); ++q)
70 {
71 appendValue(records, TBlastRecord(qIds[q]));
72 TBlastRecord & r = back(records);
73
74 r.qLength = length(queries[q]);
75
76 for (unsigned s = 0; s < length(subjects); ++s)
77 {
78 appendValue(r.matches, TBlastMatch(sIds[s]));
79 TBlastMatch & m = back(records[q].matches);
80
81 assignSource(m.alignRow0, queries[q]);
82 assignSource(m.alignRow1, subjects[s]);
83
84 localAlignment(m.alignRow0, m.alignRow1, seqanScheme(context(outfile).scoringScheme));
85
86 m.qStart = beginPosition(m.alignRow0);
87 m.qEnd = endPosition(m.alignRow0);
88 m.sStart = beginPosition(m.alignRow1);
89 m.sEnd = endPosition(m.alignRow1);
90
91 m.sLength = length(subjects[s]);
92
93 computeAlignmentStats(m, context(outfile));
94 computeBitScore(m, context(outfile));
95 computeEValue(m, r.qLength, context(outfile));
96
97 if (m.eValue > 1)
98 eraseBack(records[q].matches);
99 }
100
101 r.matches.sort(); // sort by bitscore
102
103 writeRecord(outfile, r);
104 }
105
106 writeFooter(outfile);
107
108 return 0;
109 }
110 #else
main()111 int main()
112 {
113 std::cerr << "USAGE: FILE_OUT\n";
114 return 0;
115 }
116 #endif
117