« back to all changes in this revision
Viewing changes to MergePaths/MergeContigs.cpp
- browse files at revision 7
- compare with another revision
- download diff
- download tarball
- view history from revision 7
- Committer: Package Import Robot
- Author(s): Shaun Jackman
- Date: 2011-12-13 17:04:24 UTC
- mfrom: (1.1.3)
- Revision ID: package-import@ubuntu.com-20111213170424-5lenznlcdridsj81
Tags: 1.3.2-1
New upstream release.
- files added:
- Common/PMF.h
- files removed:
- Common/HashMap.h
- files modified:
- .pc/abyss-pe/bin/abyss-pe
added
removed
MergePaths/MergeContigs.cpp
1
1
#include "config.h"
2
2
#include "Common/Options.h"
3
#include "ContigGraph.h"
4
3
#include "ContigNode.h"
5
4
#include "ContigPath.h"
6
5
#include "ContigProperties.h"
7
6
#include "DataLayer/Options.h"
8
7
#include "Dictionary.h"
9
#include "DirectedGraph.h"
10
8
#include "FastaReader.h"
11
#include "GraphIO.h"
12
9
#include "IOUtil.h"
10
#include "MemoryUtil.h"
13
11
#include "smith_waterman.h"
14
12
#include "StringUtil.h"
15
13
#include "Uncompress.h"
14
#include "Graph/ContigGraph.h"
15
#include "Graph/ContigGraphAlgorithms.h"
16
#include "Graph/DirectedGraph.h"
17
#include "Graph/GraphIO.h"
18
#include "Graph/GraphUtil.h"
19
#include "Graph/Options.h"
16
20
#include <algorithm>
17
21
#include <cstdlib>
18
22
#include <fstream>
43
47
"\n"
44
48
" -k, --kmer=KMER_SIZE k-mer size\n"
45
49
" -o, --out=FILE output the merged contigs to FILE [stdout]\n"
46
" -p, --path=PATH_FILE paths output by SimpleGraph\n"
50
" -g, --graph=FILE write the contig overlap graph to FILE\n"
47
51
" --merged output only merged contigs\n"
52
" --adj output the graph in adj format\n"
53
" --dot output the graph in dot format [default]\n"
54
" --dot-meancov same as above but give the mean coverage\n"
55
" --sam output the graph in SAM format\n"
48
56
" -v, --verbose display verbose output\n"
49
57
" --help display this help and exit\n"
50
58
" --version output version information and exit\n"
53
61
54
62
namespace opt {
55
63
unsigned k; // used by ContigProperties
64
65
/** Output FASTA path. */
56
66
static string out = "-";
57
static string path;
67
68
/** Output graph path. */
69
static string graphPath;
70
71
/** Output graph format. */
72
int format = DOT;
58
73
59
74
/** Output only merged contigs. */
60
75
int onlyMerged;
66
81
static float minIdentity = 0.9;
67
82
}
68
83
69
static const char shortopts[] = "k:o:p:v";
84
static const char shortopts[] = "g:k:o:v";
70
85
71
86
enum { OPT_HELP = 1, OPT_VERSION };
72
87
73
88
static const struct option longopts[] = {
89
{ "adj", no_argument, &opt::format, ADJ },
90
{ "dot", no_argument, &opt::format, DOT },
91
{ "dot-meancov", no_argument, &opt::format, DOT_MEANCOV },
92
{ "sam", no_argument, &opt::format, SAM },
93
{ "graph", required_argument, NULL, 'g' },
74
94
{ "kmer", required_argument, NULL, 'k' },
75
95
{ "merged", no_argument, &opt::onlyMerged, 1 },
76
96
{ "out", required_argument, NULL, 'o' },
82
102
};
83
103
84
104
/* A contig sequence. */
85
typedef FastaRecord Contig;
105
struct Contig {
106
Contig(const string& comment, const string& seq)
107
: comment(comment), seq(seq) { }
108
Contig(const FastaRecord& o) : comment(o.comment), seq(o.seq) { }
109
string comment;
110
string seq;
111
};
86
112
87
113
/** The contig sequences. */
88
static vector<Contig> g_contigs;
114
typedef vector<Contig> Contigs;
89
115
90
116
/** Return the sequence of the specified contig node. The sequence
91
117
* may be ambiguous or reverse complemented.
92
118
*/
93
static Sequence sequence(const ContigNode& id)
119
static Sequence sequence(const Contigs& contigs, const ContigNode& id)
94
120
{
95
121
if (id.ambiguous()) {
96
122
string s(id.ambiguousSequence());
98
124
transform(s.begin(), s.end(), s.begin(), ::tolower);
99
125
return string(opt::k - 1, 'N') + s;
100
126
} else {
101
const Sequence& seq = g_contigs[id.id()].seq;
127
const Sequence& seq = contigs[id.id()].seq;
102
128
return id.sense() ? reverseComplement(seq) : seq;
103
129
}
104
130
}
131
157
132
158
typedef ContigGraph<DirectedGraph<ContigProperties, Distance> > Graph;
133
159
typedef graph_traits<Graph>::vertex_descriptor vertex_descriptor;
134
typedef ContigPath Path;
160
161
/** Return the properties of the specified vertex, unless u is
162
* ambiguous, in which case return the length of the ambiguous
163
* sequence.
164
*/
165
static inline
166
ContigProperties get(vertex_bundle_t, const Graph& g, ContigNode u)
167
{
168
return u.ambiguous()
169
? ContigProperties(u.length() + opt::k - 1, 0)
170
: g[u];
171
}
135
172
136
173
/** Append the sequence of contig v to seq. */
137
static void mergeContigs(const Graph& g,
174
static void mergeContigs(const Graph& g, const Contigs& contigs,
138
175
vertex_descriptor u, vertex_descriptor v,
139
Sequence& seq, const Path& path)
176
Sequence& seq, const ContigPath& path)
140
177
{
141
178
int d = get(edge_bundle, g, u, v).distance;
142
179
assert(d < 0);
143
180
unsigned overlap = -d;
144
const Sequence& s = sequence(v);
181
const Sequence& s = sequence(contigs, v);
145
182
assert(s.length() > overlap);
146
183
Sequence ao;
147
184
Sequence bo(s, 0, overlap);
159
196
} while (chomp(seq, 'n'));
160
197
161
198
// Try an overlap alignment.
162
if (opt::verbose > 1)
199
if (opt::verbose > 2)
163
200
cerr << '\n';
164
201
vector<overlap_align> overlaps;
165
alignOverlap(ao, bo, 0, overlaps, false, opt::verbose > 1);
202
alignOverlap(ao, bo, 0, overlaps, false, opt::verbose > 2);
166
203
bool good = false;
167
204
if (!overlaps.empty()) {
168
205
assert(overlaps.size() == 1);
172
209
float identity = (float)matches / consensus.size();
173
210
good = matches >= opt::minOverlap
174
211
&& identity >= opt::minIdentity;
175
if (opt::verbose > 1)
212
if (opt::verbose > 2)
176
213
cerr << matches << " / " << consensus.size()
177
214
<< " = " << identity
178
215
<< (matches < opt::minOverlap ? " (too few)"
194
231
}
195
232
}
196
233
197
static Contig mergePath(const Graph& g, const Path& path)
234
/** Return a FASTA comment for the specified path. */
235
static void pathToComment(ostream& out, const ContigPath& path)
236
{
237
assert(path.size() > 1);
238
out << path.front();
239
if (path.size() == 3)
240
out << ',' << path[1];
241
else if (path.size() > 3)
242
out << ",...";
243
out << ',' << path.back();
244
}
245
246
/** Merge the specified path. */
247
static Contig mergePath(const Graph& g, const Contigs& contigs,
248
const ContigPath& path)
198
249
{
199
250
Sequence seq;
200
251
unsigned coverage = 0;
201
for (Path::const_iterator it = path.begin();
252
for (ContigPath::const_iterator it = path.begin();
202
253
it != path.end(); ++it) {
203
254
if (!it->ambiguous())
204
255
coverage += g[*it].coverage;
205
256
if (seq.empty()) {
206
seq = sequence(*it);
257
seq = sequence(contigs, *it);
207
258
} else {
208
259
assert(it != path.begin());
209
mergeContigs(g, *(it-1), *it, seq, path);
260
mergeContigs(g, contigs, *(it-1), *it, seq, path);
210
261
}
211
262
}
212
263
ostringstream ss;
213
ss << seq.size() << ' ' << coverage;
214
return Contig("", ss.str(), seq);
264
ss << seq.size() << ' ' << coverage << ' ';
265
pathToComment(ss, path);
266
return Contig(ss.str(), seq);
215
267
}
216
268
217
/** Return a FASTA comment for the specified path. */
218
static inline string pathToComment(const ContigPath& path)
219
{
220
assert(path.size() > 1);
221
ostringstream s;
222
s << path.front();
223
if (path.size() == 3)
224
s << ',' << path[1];
225
else if (path.size() > 3)
226
s << ",...";
227
s << ',' << path.back();
228
return s.str();
229
}
269
/** A container of ContigPath. */
270
typedef vector<ContigPath> ContigPaths;
230
271
231
272
/** Read contig paths from the specified file.
232
273
* @param ids [out] the string ID of the paths
233
274
*/
234
static vector<Path> readPaths(const string& inPath,
275
static ContigPaths readPaths(const string& inPath,
235
276
vector<string>* ids = NULL)
236
277
{
237
278
if (ids != NULL)
243
284
assert_good(fin, inPath);
244
285
istream& in = inPath == "-" ? cin : fin;
245
286
246
vector<Path> paths;
287
unsigned count = 0;
288
ContigPaths paths;
247
289
string id;
248
Path path;
290
ContigPath path;
249
291
while (in >> id >> path) {
250
292
paths.push_back(path);
251
293
if (ids != NULL)
252
294
ids->push_back(id);
295
296
++count;
297
if (opt::verbose > 1 && count % 1000000 == 0)
298
cerr << "Read " << count << " paths. "
299
"Using " << toSI(getMemoryUsage())
300
<< "B of memory.\n";
253
301
}
302
if (opt::verbose > 0)
303
cerr << "Read " << count << " paths. "
304
"Using " << toSI(getMemoryUsage()) << "B of memory.\n";
254
305
assert(in.eof());
255
306
return paths;
256
307
}
257
308
258
309
/** Finds all contigs used in each path in paths, and
259
310
* marks them as seen in the vector seen. */
260
static void seenContigs(vector<bool>& seen, const vector<Path>& paths)
311
static void seenContigs(vector<bool>& seen, const ContigPaths& paths)
261
312
{
262
for (vector<Path>::const_iterator it = paths.begin();
313
for (ContigPaths::const_iterator it = paths.begin();
263
314
it != paths.end(); ++it)
264
for (Path::const_iterator itc = it->begin();
315
for (ContigPath::const_iterator itc = it->begin();
265
316
itc != it->end(); ++itc)
266
317
if (itc->id() < seen.size())
267
318
seen[itc->id()] = true;
271
322
* should be removed.
272
323
*/
273
324
static void markRemovedContigs(vector<bool>& marked,
274
const vector<string>& pathIDs, const vector<Path>& paths)
325
const vector<string>& pathIDs, const ContigPaths& paths)
275
326
{
276
for (vector<Path>::const_iterator it = paths.begin();
327
for (ContigPaths::const_iterator it = paths.begin();
277
328
it != paths.end(); ++it)
278
329
if (it->empty())
279
330
marked[ContigID(pathIDs[it - paths.begin()])] = true;
280
331
}
281
332
333
/** Output the updated overlap graph. */
334
static void outputGraph(Graph& g,
335
const vector<string>& pathIDs, const ContigPaths& paths,
336
const string& commandLine)
337
{
338
// Add the path vertices.
339
for (ContigPaths::const_iterator it = paths.begin();
340
it != paths.end(); ++it) {
341
const ContigPath& path = *it;
342
const string& id = pathIDs[it - paths.begin()];
343
if (!path.empty()) {
344
ContigID::insert(id);
345
merge(g, path.begin(), path.end());
346
}
347
}
348
349
// Remove the vertices that are used in paths.
350
for (ContigPaths::const_iterator it = paths.begin();
351
it != paths.end(); ++it) {
352
const ContigPath& path = *it;
353
const string& id = pathIDs[it - paths.begin()];
354
if (path.empty()) {
355
remove_vertex(ContigNode(id, false), g);
356
} else {
357
remove_vertex_if(g, path.begin(), path.end(),
358
not1(std::mem_fun_ref(&ContigNode::ambiguous)));
359
}
360
}
361
362
// Output the graph.
363
const string& graphPath = opt::graphPath;
364
assert(!graphPath.empty());
365
if (opt::verbose > 0)
366
cerr << "Writing `" << graphPath << "'..." << endl;
367
ofstream fout(graphPath.c_str());
368
assert_good(fout, graphPath);
369
write_graph(fout, g, PROGRAM, commandLine);
370
assert_good(fout, graphPath);
371
if (opt::verbose > 0)
372
printGraphStats(cerr, g);
373
}
374
282
375
int main(int argc, char** argv)
283
376
{
284
377
opt::trimMasked = false;
285
378
379
string commandLine;
380
{
381
ostringstream ss;
382
char** last = argv + argc - 1;
383
copy(argv, last, ostream_iterator<const char *>(ss, " "));
384
ss << *last;
385
commandLine = ss.str();
386
}
387
286
388
bool die = false;
287
389
for (int c; (c = getopt_long(argc, argv,
288
390
shortopts, longopts, NULL)) != -1;) {
289
391
istringstream arg(optarg != NULL ? optarg : "");
290
392
switch (c) {
291
393
case '?': die = true; break;
394
case 'g': arg >> opt::graphPath; assert(arg.eof()); break;
292
395
case 'k': arg >> opt::k; break;
293
396
case 'o': arg >> opt::out; break;
294
case 'p': arg >> opt::path; break;
295
397
case 'v': opt::verbose++; break;
296
398
case OPT_HELP:
297
399
cout << USAGE_MESSAGE;
317
419
die = true;
318
420
}
319
421
422
if (argc - optind > 3) {
423
cerr << PROGRAM ": too many arguments\n";
424
die = true;
425
}
426
320
427
if (die) {
321
428
cerr << "Try `" << PROGRAM
322
429
<< " --help' for more information.\n";
325
432
326
433
const char* contigFile = argv[optind++];
327
434
string adjPath, mergedPathFile;
328
if (argc - optind > 1)
435
Graph g;
436
if (argc - optind > 1) {
329
437
adjPath = string(argv[optind++]);
438
439
// Read the contig adjacency graph.
440
if (opt::verbose > 0)
441
cerr << "Reading `" << adjPath << "'..." << endl;
442
ifstream fin(adjPath.c_str());
443
assert_good(fin, adjPath);
444
fin >> g;
445
assert(fin.eof());
446
if (opt::verbose > 0)
447
cerr << "Read " << num_vertices(g) << " vertices. "
448
"Using " << toSI(getMemoryUsage()) << "B of memory.\n";
449
}
330
450
mergedPathFile = string(argv[optind++]);
331
451
332
452
// Read the contig sequence.
333
vector<Contig>& contigs = g_contigs;
453
Contigs contigs;
334
454
{
455
if (opt::verbose > 0)
456
cerr << "Reading `" << contigFile << "'..." << endl;
457
unsigned count = 0;
335
458
FastaReader in(contigFile, FastaReader::NO_FOLD_CASE);
336
459
for (FastaRecord rec; in >> rec;) {
337
ContigID::insert(rec.id);
460
if (!adjPath.empty()
461
&& ContigID::count(rec.id) == 0)
462
continue;
463
ContigID id = adjPath.empty()
464
? ContigID::insert(rec.id)
465
: ContigID(rec.id);
466
assert(id == contigs.size());
338
467
contigs.push_back(rec);
468
469
++count;
470
if (opt::verbose > 1 && count % 1000000 == 0)
471
cerr << "Read " << count << " sequences. "
472
"Using " << toSI(getMemoryUsage())
473
<< "B of memory.\n";
339
474
}
475
if (opt::verbose > 0)
476
cerr << "Read " << count << " sequences. "
477
"Using " << toSI(getMemoryUsage())
478
<< "B of memory.\n";
340
479
assert(in.eof());
341
480
assert(!contigs.empty());
342
481
opt::colourSpace = isdigit(contigs[0].seq[0]);
343
if (optind == argc)
344
ContigID::lock();
482
ContigID::lock();
345
483
}
346
484
347
485
vector<string> pathIDs;
348
vector<Path> paths = readPaths(mergedPathFile, &pathIDs);
349
if (opt::verbose > 0)
350
cerr << "Number of paths: " << paths.size() << '\n';
486
ContigPaths paths = readPaths(mergedPathFile, &pathIDs);
351
487
352
488
// Record all the contigs that are in a path.
353
489
vector<bool> seen(contigs.size());
354
490
seenContigs(seen, paths);
355
491
markRemovedContigs(seen, pathIDs, paths);
356
492
357
// Record all the contigs that were in a previous path.
358
if (argc - optind > 0) {
359
unsigned count = 0;
360
for (; optind < argc; optind++) {
361
vector<Path> prevPaths = readPaths(argv[optind]);
362
seenContigs(seen, prevPaths);
363
count += prevPaths.size();
364
}
365
if (opt::verbose > 0)
366
cerr << "Number of previous paths: " << count << '\n';
367
}
368
369
// Record all the contigs that are seeds.
370
if (!opt::path.empty()) {
371
vector<bool> seenPivots(contigs.size());
372
ifstream fin(opt::path.c_str());
373
assert_good(fin, opt::path);
374
for (ContigID id; fin >> id;) {
375
fin.ignore(numeric_limits<streamsize>::max(), '\n');
376
assert(id < contigs.size());
377
// Only count a pivot as seen if it was in a final path.
378
if (seen[id])
379
seenPivots[id] = true;
380
}
381
assert(fin.eof());
382
seen = seenPivots;
383
}
384
385
493
// Output those contigs that were not seen in a path.
386
494
ofstream fout;
387
495
ostream& out = opt::out == "-" ? cout
388
496
: (fout.open(opt::out.c_str()), fout);
389
497
assert_good(out, opt::out);
390
498
if (!opt::onlyMerged) {
391
for (vector<Contig>::const_iterator it = contigs.begin();
392
it != contigs.end(); ++it)
393
if (!seen[it - contigs.begin()])
394
out << *it;
499
for (Contigs::const_iterator it = contigs.begin();
500
it != contigs.end(); ++it) {
501
ContigID id(it - contigs.begin());
502
if (!seen[id]) {
503
const Contig& contig = *it;
504
out << '>' << id;
505
if (!contig.comment.empty())
506
out << ' ' << contig.comment;
507
out << '\n' << contig.seq << '\n';
508
}
509
}
395
510
}
396
511
397
512
if (adjPath.empty())
398
513
return 0;
399
514
400
// Read the contig adjacency graph.
401
ifstream fin(adjPath.c_str());
402
assert_good(fin, adjPath);
403
Graph g;
404
fin >> g;
405
assert(fin.eof());
406
407
515
unsigned npaths = 0;
408
for (vector<Path>::const_iterator it = paths.begin();
516
for (ContigPaths::const_iterator it = paths.begin();
409
517
it != paths.end(); ++it) {
410
const Path& path = *it;
518
const ContigPath& path = *it;
411
519
if (path.empty())
412
520
continue;
413
Contig contig = mergePath(g, path);
414
contig.id = pathIDs[it - paths.begin()];
415
FastaRecord rec(contig);
416
rec.comment += ' ' + pathToComment(path);
417
out << rec;
521
Contig contig = mergePath(g, contigs, path);
522
out << '>' << pathIDs[it - paths.begin()]
523
<< ' ' << contig.comment << '\n'
524
<< contig.seq << '\n';
418
525
assert_good(out, opt::out);
419
526
npaths++;
420
527
}
435
542
minCovUsed = min(minCovUsed, cov);
436
543
}
437
544
545
if (!opt::graphPath.empty())
546
outputGraph(g, pathIDs, paths, commandLine);
547
438
548
cerr << "The minimum coverage of single-end contigs is "
439
549
<< minCov << ".\n"
440
550
<< "The minimum coverage of merged contigs is "
Loggerhead is a web-based interface for Breezy
Version: 2.0.1