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