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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
#include <cassert>
#include <cmath>
#include <iostream>
#include <string>
#include <alignment/datastructures/alignment/Path.h>
#include <pbdata/defs.h>
#include <alignment/algorithms/alignment/AlignmentUtils.hpp>
#include <alignment/algorithms/alignment/SDPAlign.hpp>
#include <alignment/algorithms/alignment/SWAlign.hpp>
#include <alignment/datastructures/alignment/Alignment.hpp>
#include <alignment/format/StickAlignmentPrinter.hpp>
#include <alignment/tuples/DNATuple.hpp>
#include <alignment/tuples/TupleMetrics.hpp>
#include <pbdata/FASTAReader.hpp>
#include <pbdata/FASTASequence.hpp>
/*
* Performs sparse dynamic programming (SDP) between pairs of sequences as they
* are given in two FASTA files, one called for convenience query, the other
* target. k is the size of the k-mer used for the SDP algorithm.
*/
void PrintUsage()
{
std::cout << "usage: sdpMatcher query target k [-indelRate delta] "
"[-showalign] [-printsw] [-noRefine] [-indel i] [ -local ] "
"[-match m] [-sdpIndel i]"
<< std::endl;
}
int main(int argc, char* argv[])
{
if (argc < 4) {
PrintUsage();
std::exit(EXIT_FAILURE);
}
std::string queryName, targetName;
queryName = argv[1];
targetName = argv[2];
TupleMetrics tm;
tm.Initialize(atoi(argv[3]));
int argi = 4;
float indelRate = 0.25;
int indel = 3;
int match = 0;
int printSW = 0;
int printSimilarity = 0;
int refineAlignments = 1;
int showalign = 0;
int fixedTarget = 0;
int sdpIndel = indel;
int sdpIns = 5;
int sdpDel = 5;
(void)(sdpIndel);
(void)(sdpIns);
(void)(sdpDel); // not yet used.
AlignmentType alignType = Global;
while (argi < argc) {
if (strcmp(argv[argi], "-indelRate") == 0) {
++argi;
indelRate = atof(argv[argi]);
} else if (strcmp(argv[argi], "-printsw") == 0) {
printSW = 1;
} else if (strcmp(argv[argi], "-noRefine") == 0) {
refineAlignments = 0;
} else if (strcmp(argv[argi], "-indel") == 0) {
indel = atoi(argv[++argi]);
} else if (strcmp(argv[argi], "-sdpIndel") == 0) {
sdpIndel = atoi(argv[++argi]);
} else if (strcmp(argv[argi], "-sdpIns") == 0) {
sdpIns = atoi(argv[++argi]);
} else if (strcmp(argv[argi], "-sdpDel") == 0) {
sdpDel = atoi(argv[++argi]);
} else if (strcmp(argv[argi], "-showalign") == 0) {
showalign = 1;
} else if (strcmp(argv[argi], "-local") == 0) {
alignType = Local;
} else if (strcmp(argv[argi], "-match") == 0) {
match = atoi(argv[++argi]);
} else if (strcmp(argv[argi], "-fixedtarget") == 0) {
fixedTarget = 1;
} else if (strcmp(argv[argi], "-printSimilarity") == 0) {
printSimilarity = 1;
} else {
PrintUsage();
std::cout << "Bad option: " << argv[argi] << std::endl;
std::exit(EXIT_FAILURE);
}
++argi;
}
FASTASequence query, target;
FASTAReader queryReader, targetReader;
queryReader.Init(queryName);
targetReader.Init(targetName);
if (match != 0) {
int i;
for (i = 0; i < 4; i++) {
LocalAlignLowMutationMatrix[i][i] = match;
}
}
int seqIndex = 0;
Alignment alignment;
std::vector<int> scoreMat;
std::vector<Arrow> pathMat;
DistanceMatrixScoreFunction<DNASequence, DNASequence> distScoreFn;
distScoreFn.del = indel;
distScoreFn.ins = indel;
distScoreFn.InitializeScoreMatrix(SMRTDistanceMatrix);
if (fixedTarget) {
targetReader.GetNext(target);
}
std::cout << "qid,tid,qstart,qend,qlen,tstart,tend,tlen,score";
if (printSimilarity) std::cout << ",pctSimilarity";
std::cout << std::endl;
while (queryReader.GetNext(query) and (fixedTarget or targetReader.GetNext(target))) {
if (query.length == 0 or target.length == 0) continue;
alignment.blocks.clear();
int alignScore;
alignScore = SDPAlign(query, target, distScoreFn, tm.tupleSize, sdpIndel, sdpIndel,
indelRate, alignment, alignType, refineAlignments, false, 0);
ComputeAlignmentStats(alignment, query.seq, target.seq, distScoreFn);
if (alignScore > 0) { // in rare cases the SDP returns positive.
alignScore = 0; // this makes it more like a true local alignment
}
if (showalign) {
StickPrintAlignment(alignment, query, target, std::cout);
}
if (printSW) {
MatchedAlignment swAlignment;
std::vector<int> scoreMat;
std::vector<Arrow> pathMat;
SWAlign(query, target, scoreMat, pathMat, swAlignment, distScoreFn);
StickPrintAlignment(swAlignment, query, target, std::cout);
}
std::cout << query.GetName() << "," << target.GetName() << "," << alignment.qPos << ","
<< alignment.QEnd() << "," << query.length << "," << alignment.tPos << ","
<< alignment.TEnd() << "," << target.length << "," << alignScore;
if (printSimilarity) std::cout << "," << alignment.pctSimilarity;
std::cout << std::endl;
++seqIndex;
}
return 0;
}
|
