C++ ContigPath类代码示例

/** Find every path that overlaps with the specified path. */
static void findOverlaps(const Graph& g,
		const Paths& paths, const SeedMap& seedMap,
		const Vertex& v, Overlaps& overlaps)
{
	ContigPath rc;
	if (v.sense) {
		rc = paths[v.id];
		reverseComplement(rc.begin(), rc.end());
	}
	const ContigPath& path = v.sense ? rc : paths[v.id];

	for (ContigPath::const_iterator it = path.begin();
			it != path.end(); ++it) {
		if (it->ambiguous())
			continue;

		pair<SeedMap::const_iterator, SeedMap::const_iterator>
			range = seedMap.equal_range(*it);
		for (SeedMap::const_iterator seed = range.first;
				seed != range.second; ++seed) {
			if (v == seed->second)
				continue;
			int distance = 0;
			unsigned overlap = findOverlap(g, paths, it, path.end(),
					   seed->second, distance);
			if (overlap > 0)
				overlaps.push_back(Overlap(v, seed->second,
					overlap, distance));

		}
	}
}

/** Read contig paths from the specified file.
 * @param g the contig adjacency graph
 * @param inPath the file of contig paths
 * @param[out] pathIDs the path IDs
 * @return the paths
 */
static Paths readPaths(Graph& g,
		const string& inPath, vector<string>& pathIDs)
{
	typedef graph_traits<Graph>::vertex_descriptor V;

	assert(pathIDs.empty());
	ifstream fin(inPath.c_str());
	if (opt::verbose > 0)
		cerr << "Reading `" << inPath << "'..." << endl;
	if (inPath != "-")
		assert_good(fin, inPath);
	istream& in = inPath == "-" ? cin : fin;

	assert_good(in, inPath);
	Paths paths;
	string id;
	ContigPath path;
	while (in >> id >> path) {
		if (path.empty()) {
			// Remove this contig from the graph.
			V u = find_vertex(id, false, g);
			clear_vertex(u, g);
			remove_vertex(u, g);
		} else {
			pathIDs.push_back(id);
			paths.push_back(path);
		}
	}
	assert(in.eof());
	return paths;
}

/** Attempt to merge the paths specified in mergeQ with path.
 * @return the number of paths merged
 */
static unsigned mergePaths(const Lengths& lengths,
		ContigPath& path,
		deque<ContigNode>& mergeQ, set<ContigNode>& seen,
		const ContigPathMap& paths)
{
	unsigned merged = 0;
	deque<ContigNode> invalid;
	for (ContigNode pivot; !mergeQ.empty(); mergeQ.pop_front()) {
		pivot = mergeQ.front();
		ContigPathMap::const_iterator path2It
			= paths.find(pivot.contigIndex());
		if (path2It == paths.end())
			continue;

		ContigPath path2 = path2It->second;
		if (pivot.sense())
			reverseComplement(path2.begin(), path2.end());
		ContigPath consensus = align(lengths, path, path2, pivot);
		if (consensus.empty()) {
			invalid.push_back(pivot);
			continue;
		}

		appendToMergeQ(mergeQ, seen, path2);
		path.swap(consensus);
		if (gDebugPrint)
#pragma omp critical(cout)
			cout << get(g_contigNames, pivot)
				<< '\t' << path2 << '\n'
				<< '\t' << path << '\n';
		merged++;
	}
	mergeQ.swap(invalid);
	return merged;
}

/** Return a map of contig IDs to their distance along this path.
 * Repeat contigs, which would have more than one position, are not
 * represented in this map.
 */
map<ContigNode, int> makeDistanceMap(const Graph& g,
		const ContigNode& origin, const ContigPath& path)
{
	map<ContigNode, int> distances;
	int distance = 0;
	for (ContigPath::const_iterator it = path.begin();
			it != path.end(); ++it) {
		vertex_descriptor u = it == path.begin() ? origin : *(it - 1);
		vertex_descriptor v = *it;
		distance += getDistance(g, u, v);

		bool inserted = distances.insert(
				make_pair(v, distance)).second;
		if (!inserted) {
			// Mark this contig as a repeat.
			distances[v] = INT_MIN;
		}

		distance += g[v].length;
	}

	// Remove the repeats.
	for (map<ContigNode, int>::iterator it = distances.begin();
			it != distances.end();)
		if (it->second == INT_MIN)
			distances.erase(it++);
		else
			++it;
	return distances;
}

/** Remove ambiguous edges that overlap by only a small amount.
 * Remove the edge (u,v) if deg+(u) > 1 and deg-(v) > 1 and the
 * overlap of (u,v) is small.
 */
static void removeSmallOverlaps(PathGraph& g,
		const ContigPathMap& paths)
{
	typedef graph_traits<PathGraph>::edge_descriptor E;
	typedef graph_traits<PathGraph>::out_edge_iterator Eit;
	typedef graph_traits<PathGraph>::vertex_descriptor V;
	typedef graph_traits<PathGraph>::vertex_iterator Vit;

	vector<E> edges;
	pair<Vit, Vit> urange = vertices(g);
	for (Vit uit = urange.first; uit != urange.second; ++uit) {
		V u = *uit;
		if (out_degree(u, g) < 2)
			continue;
		ContigPath pathu = getPath(paths, u);
		pair<Eit, Eit> uvits = out_edges(u, g);
		for (Eit uvit = uvits.first; uvit != uvits.second; ++uvit) {
			E uv = *uvit;
			V v = target(uv, g);
			assert(v != u);
			if (in_degree(v, g) < 2)
				continue;
			ContigPath pathv = getPath(paths, v);
			if (pathu.back() == pathv.front()
					&& paths.count(pathu.back().contigIndex()) > 0)
				edges.push_back(uv);
		}
	}
	remove_edges(g, edges.begin(), edges.end());
	if (opt::verbose > 0)
		cout << "Removed " << edges.size()
			<< " small overlap edges.\n";
	if (!opt::db.empty())
		addToDb(db, "Edges_removed_small_overlap", edges.size());
}

static void appendToMergeQ(deque<ContigNode>& mergeQ,
	set<ContigNode>& seen, const ContigPath& path)
{
	for (ContigPath::const_iterator it = path.begin();
			it != path.end(); ++it)
		if (!it->ambiguous() && seen.insert(*it).second)
			mergeQ.push_back(*it);
}

/** Mark every contig in path as seen. */
static void markSeen(vector<bool>& seen, const ContigPath& path,
		bool flag)
{
	for (Path::const_iterator it = path.begin();
			it != path.end(); ++it)
		if (!it->ambiguous() && it->id() < seen.size())
			seen[it->id()] = flag;
}

/** Return the specified path. */
static ContigPath getPath(const ContigPathMap& paths, ContigNode u)
{
	ContigPathMap::const_iterator it = paths.find(u.contigIndex());
	assert(it != paths.end());
	ContigPath path = it->second;
	if (u.sense())
		reverseComplement(path.begin(), path.end());
	return path;
}

/** Check whether path starts with the sequence [first, last). */
static bool startsWith(ContigPath path, bool rc,
		ContigPath::const_iterator first,
		ContigPath::const_iterator last)
{
	if (rc)
		reverseComplement(path.begin(), path.end());
	assert(*first == path.front());
	assert(first < last);
	return unsigned(last - first) > path.size() ? false
		: equal(first, last, path.begin());
}

/** Return a FASTA comment for the specified path. */
static void pathToComment(ostream& out,
		const Graph& g, const ContigPath& path)
{
	assert(path.size() > 1);
	out << get(vertex_name, g, path.front());
	if (path.size() == 3)
		out << ',' << get(vertex_name, g, path[1]);
	else if (path.size() > 3)
		out << ",...";
	out << ',' << get(vertex_name, g, path.back());
}

/** Remove ambiguous contigs from the ends of the path. */
static void removeAmbiguousContigs(ContigPath& path)
{
	if (!path.empty() && path.back().ambiguous())
		path.erase(path.end() - 1);
	if (!path.empty() && path.front().ambiguous())
		path.erase(path.begin());
}

/** Return the length of the specified path in k-mer. */
static unsigned calculatePathLength(const Graph& g,
		const ContigNode& origin,
		const ContigPath& path, size_t prefix = 0, size_t suffix = 0)
{
	if (prefix + suffix == path.size())
		return 0;
	assert(prefix + suffix < path.size());
	int length = addProp(g, path.begin() + prefix,
			path.end() - suffix).length;

	// Account for the overlap on the left.
	vertex_descriptor u = prefix == 0 ? origin : path[prefix - 1];
	length += getDistance(g, u, path[prefix]);
	assert(length > 0);
	return length;
}

/** Merge a sequence of overlapping paths. */
static ContigPath mergePaths(const Paths& paths,
		const OverlapMap& overlaps, const ContigPath& merge)
{
	assert(!merge.empty());
	ContigNode u = merge.front();
	ContigPath path(getPath(paths, u));
	for (ContigPath::const_iterator it = merge.begin() + 1;
			it != merge.end(); ++it) {
		ContigNode v = *it;
		ContigPath vpath(getPath(paths, v));
		unsigned overlap = getOverlap(overlaps, u, v);
		assert(path.size() > overlap);
		assert(vpath.size() > overlap);
		assert(equal(path.end() - overlap, path.end(),
					vpath.begin()));
		path.insert(path.end(), vpath.begin() + overlap, vpath.end());
		u = v;
	}
	return path;
}

/** Merge the specified path. */
static Contig mergePath(const Graph& g, const Contigs& contigs,
		const ContigPath& path)
{
	Sequence seq;
	unsigned coverage = 0;
	for (ContigPath::const_iterator it = path.begin();
			it != path.end(); ++it) {
		if (!it->ambiguous())
			coverage += g[*it].coverage;
		if (seq.empty()) {
			seq = sequence(contigs, *it);
		} else {
			assert(it != path.begin());
			mergeContigs(g, contigs, *(it-1), *it, seq, path);
		}
	}
	ostringstream ss;
	ss << seq.size() << ' ' << coverage << ' ';
	pathToComment(ss, g, path);
	return Contig(ss.str(), seq);
}

/** Add an edge if the two paths overlap.
 * @param pivot the pivot at which to seed the alignment
 * @return whether an overlap was found
 */
static bool addOverlapEdge(const Lengths& lengths,
		PathGraph& gout, ContigNode pivot,
		ContigNode seed1, const ContigPath& path1,
		ContigNode seed2, const ContigPath& path2)
{
	assert(seed1 != seed2);

	// Determine the orientation of the overlap edge.
	dir_type orientation = DIR_X;
	ContigPath consensus = align(lengths,
			path1, path2, pivot, orientation);
	if (consensus.empty())
		return false;
	assert(orientation != DIR_X);
	if (orientation == DIR_B) {
		// One of the paths subsumes the other. Use the order of the
		// seeds to determine the orientation of the edge.
		orientation = find(consensus.begin(), consensus.end(), seed1)
			< find(consensus.begin(), consensus.end(), seed2)
			? DIR_F : DIR_R;
	}
	assert(orientation == DIR_F || orientation == DIR_R);

	// Add the edge.
	ContigNode u = orientation == DIR_F ? seed1 : seed2;
	ContigNode v = orientation == DIR_F ? seed2 : seed1;
	bool added = false;
#pragma omp critical(gout)
	if (!edge(u, v, gout).second) {
		add_edge(u, v, gout);
		added = true;
	}
	return added;
}