New upstream version 1.7.1+dfsg (608e890a) · Commits · Debian Med / canu

src/AS_UTL/stddev.H

+6 −6

Original line number	Diff line number	Diff line
		@@ -286,7 +286,7 @@ computeExponentialMovingAverage(TT alpha, TT ema, TT value) {




		#if 0

		template<typename TT>
		class genericStatistics {
		@@ -338,12 +338,12 @@ public:

		vector<uint64> &histogram(void) { // Returns pointer to private histogram data
		finalizeData();
		return(&_histogram);
		return(_histogram);
		};

		vector<uint64> &Nstatistics(void) { // Returns pointer to private N data
		finalizeData();
		return(&_Nstatistics);
		return(_Nstatistics);
		};

		void finalizeData(void) {
		@@ -375,7 +375,7 @@ private:
		vector<uint64> _Nstatistics;
		};


		#endif



		@@ -422,12 +422,12 @@ public:
		#if 0
		vector<uint64> &histogram(void) { // Returns pointer to private histogram data
		finalizeData();
		return(&_histogram);
		return(_histogram);
		};

		vector<uint64> &Nstatistics(void) { // Returns pointer to private N data
		finalizeData();
		return(&_Nstatistics);
		return(_Nstatistics);
		};
		#endif

src/bogart/AS_BAT_BestOverlapGraph.C

+70 −0

Original line number	Diff line number	Diff line
		@@ -130,6 +130,27 @@ BestOverlapGraph::removeHighErrorBestEdges(void) {

		if (b5->readId() != 0) edgeStats.insert(erates[eratesLen++] = b5->erate());
		if (b3->readId() != 0) edgeStats.insert(erates[eratesLen++] = b3->erate());

		// If there are NO best edges, find the overlap with the most matches and use that.

		if ((b5->readId() == 0) &&
		(b3->readId() == 0)) {
		uint32 no = 0;
		BAToverlap *ovl = OC->getOverlaps(fi, no);
		uint32 bestM = 0;
		double bestE = 0.0;

		for (uint32 oo=0; oo<no; oo++) {
		double matches = (1 - ovl[oo].erate()) * RI->overlapLength(ovl[oo].a_iid, ovl[oo].b_iid, ovl[oo].a_hang, ovl[oo].b_hang);
		if (bestM < matches) {
		bestM = matches;
		bestE = ovl[oo].erate();
		}
		}

		if (no > 0)
		edgeStats.insert(erates[eratesLen++] = bestE);
		}
		}

		_mean = edgeStats.mean();
		@@ -357,6 +378,53 @@ BestOverlapGraph::removeSpurs(const char *prefix) {



		// Mark zombie masters. Any read that has only contained overlaps (it is the container) and is the
		// smallest ID of those with no hangs, is a master. These get promoted to unitigs.
		//
		void
		BestOverlapGraph::findZombies(const char *prefix) {
		uint32 fiLimit = RI->numReads();
		uint32 numThreads = omp_get_max_threads();
		uint32 blockSize = (fiLimit < 100 * numThreads) ? numThreads : fiLimit / 99;

		#pragma omp parallel for schedule(dynamic, blockSize)
		for (uint32 fi=1; fi <= fiLimit; fi++) {
		uint32 no = 0;
		BAToverlap *ovl = OC->getOverlaps(fi, no);
		uint32 nc = 0;

		if (no == 0)
		continue;

		for (uint32 ii=0; ii<no; ii++, nc++) // If any overlap makes A not
		if (ovl[ii].AisContainer() == false) // a container, it's not a zombie
		break;

		if (nc < no)
		continue;

		nc = UINT32_MAX;

		for (uint32 ii=0; ii<no; ii++) // Find the smallest ID
		if ((ovl[ii].a_hang == 0) && // with no hangs.
		(ovl[ii].b_hang == 0) &&
		(ovl[ii].b_iid < nc))
		nc = ovl[ii].b_iid;

		if (fi < nc) { // If we're smaller, we're a
		#pragma omp critical (suspInsert) // Zombie Master!
		{
		writeLog("read %u is a zombie.\n", fi);
		_zombie.insert(fi);
		}
		}
		}

		writeStatus("BestOverlapGraph()-- detected " F_SIZE_T " zombie reads.\n", _zombie.size());
		}



		void
		BestOverlapGraph::findEdges(void) {
		uint32 fiLimit = RI->numReads();
		@@ -513,6 +581,8 @@ BestOverlapGraph::BestOverlapGraph(double erateGraph,
		findEdges();
		}

		findZombies(prefix);

		reportBestEdges(prefix, logFileFlagSet(LOG_ALL_BEST_EDGES) ? "best.3.final" : "best");

		// One more pass, to find any ambiguous best edges.

src/bogart/AS_BAT_BestOverlapGraph.H

+6 −0

Original line number	Diff line number	Diff line
		@@ -183,6 +183,7 @@ private:
		void removeSuspicious(const char *prefix);
		void removeSpurs(const char *prefix);
		void removeLopsidedEdges(const char *prefix);
		void findZombies(const char *prefix);

		void findEdges(void);

		@@ -240,6 +241,10 @@ public:
		return(_suspicious.count(readid) > 0);
		};

		bool isZombie(const uint32 readid) {
		return(_zombie.count(readid) > 0);
		};

		void reportEdgeStatistics(const char prefix, const char label);
		void reportBestEdges(const char prefix, const char label);

		@@ -283,6 +288,7 @@ private:
		set<uint32> _suspicious;
		set<uint32> _singleton;
		set<uint32> _spur;
		set<uint32> _zombie;

		map<uint32, BestOverlaps> _bestM;
		map<uint32, BestScores> _scorM;

src/bogart/AS_BAT_DropDeadEnds.C

+66 −34

Original line number	Diff line number	Diff line
		@@ -133,17 +133,16 @@ findPrevRead(Unitig *tig,

		uint32
		dropDeadFirstRead(AssemblyGraph *AG,
		Unitig *tig) {
		Unitig *tig,
		bool isForward) {

		ufNode *fn = tig->firstRead();
		ufNode *sn = findNextRead(tig, fn);

		// No next read, keep fn in the tig.

		if (sn == NULL) {
		writeLog("dropDead()- read %8u no sn\n", fn->ident);
		if (sn == NULL)
		return(0);
		}

		// Over all edges from the first read, look for any edge to something else.
		//
		@@ -153,23 +152,35 @@ dropDeadFirstRead(AssemblyGraph *AG,
		// the tig. We assume that this is always the first read, which is OK, because the function name
		// says so. Any edge to anywhere means the read is good and should be kept.

		if (AG->getForward(fn->ident).size() == 0)
		writeLog("dropDead()-- (%s) 1st read %8u has no edges\n", (isForward) ? "fwd" : "rev", fn->ident);

		for (uint32 pp=0; pp<AG->getForward(fn->ident).size(); pp++) {
		BestPlacement &pf = AG->getForward(fn->ident)[pp];

		writeLog("dropDead()-- 1st read %8u %s pf %3u/%3u best5 %8u best3 %8u bestC %8u\n",

		if (pf.bestC.b_iid > 0) {
		writeLog("dropDead()-- (%s) 1st read %8u %s edge %4u/%-4u - contained in %u - confirmed\n",
		(isForward) ? "fwd" : "rev",
		fn->ident,
		fn->position.isForward() ? "->" : "<-",
		pp, AG->getForward(fn->ident).size(),
		pf.best5.b_iid, pf.best3.b_iid, pf.bestC.b_iid);

		if (pf.bestC.b_iid > 0) {
		pf.bestC.b_iid);
		return(0);
		}

		if (((fn->position.isForward() == true) && (pf.best5.b_iid != 0)) \|\|
		((fn->position.isForward() == false) && (pf.best3.b_iid != 0)))
		((fn->position.isForward() == false) && (pf.best3.b_iid != 0))) {
		writeLog("dropDead()-- (%s) 1st read %8u %s edge %4u/%-4u - 5:%u 3:%u - confirmed\n",
		(isForward) ? "fwd" : "rev",
		fn->ident,
		fn->position.isForward() ? "->" : "<-",
		pp, AG->getForward(fn->ident).size(),
		pf.best5.b_iid, pf.best3.b_iid);
		return(0);
		}
		}


		// But no edge means we need to check the second read. If it has an edge, then we infer the
		// first read is bogus and should be removed. If it also has no edge (except to the first read,
		@@ -180,26 +191,41 @@ dropDeadFirstRead(AssemblyGraph *AG,
		// first read. Well, and that if the second read has an edge we declare the first read to be
		// junk. That's also a bit of a difference from the previous loop.

		if (AG->getForward(sn->ident).size() == 0) {
		writeLog("dropDead()-- (%s) 2nd read %8u has no edges - keep first\n", (isForward) ? "fwd" : "rev", sn->ident);
		return(0);
		}

		for (uint32 pp=0; pp<AG->getForward(sn->ident).size(); pp++) {
		BestPlacement &pf = AG->getForward(sn->ident)[pp];

		writeLog("dropDead()-- 2nd read %8u %s pf %3u/%3u best5 %8u best3 %8u bestC %8u\n",
		if ((pf.bestC.b_iid > 0) && (pf.bestC.b_iid != fn->ident)) {
		writeLog("dropDead()-- (%s) 2nd read %8u %s edge %4u/%-4u - contained in %u - drop first\n",
		(isForward) ? "fwd" : "rev",
		sn->ident,
		sn->position.isForward() ? "->" : "<-",
		pp, AG->getForward(sn->ident).size(),
		pf.best5.b_iid, pf.best3.b_iid, pf.bestC.b_iid);

		if ((pf.bestC.b_iid > 0) && (pf.bestC.b_iid != fn->ident))
		pf.bestC.b_iid);
		return(fn->ident);
		}

		if (((sn->position.isForward() == true) && (pf.best5.b_iid != 0) && (pf.best5.b_iid != fn->ident)) \|\|
		((sn->position.isForward() == false) && (pf.best3.b_iid != 0) && (pf.best3.b_iid != fn->ident)))
		((sn->position.isForward() == false) && (pf.best3.b_iid != 0) && (pf.best3.b_iid != fn->ident))) {
		writeLog("dropDead()-- (%s) 2nd read %8u %s edge %4u/%-4u - 5:%u 3:8u - drop first\n",
		(isForward) ? "fwd" : "rev",
		sn->ident,
		sn->position.isForward() ? "->" : "<-",
		pp, AG->getForward(sn->ident).size(),
		pf.best5.b_iid, pf.best3.b_iid);
		return(fn->ident);
		}
		}

		// Otherwise, the second read had only edges to the first read, and we should keep the first
		// read.

		writeLog("dropDead()-- (%s) 2nd read %8u has no useful external edges - keep first\n", (isForward) ? "fwd" : "rev", sn->ident);

		return(0);
		}

		@@ -222,20 +248,32 @@ dropDeadEnds(AssemblyGraph *AG,
		(tig->_isUnassembled == true))
		continue;

		uint32 fn = dropDeadFirstRead(AG, tig); // Decide if the first read is junk.
		writeLog("\n");
		writeLog("dropDead()-- testing tig %u length %u\n", ti, tig->getLength());

		uint32 fn = dropDeadFirstRead(AG, tig, true); // Decide if the first read is junk.

		//fprintf(stderr, "drop dead for tig %u (flipped)\n", ti);

		tig->reverseComplement(); // Flip.
		uint32 ln = dropDeadFirstRead(AG, tig); // Decide if the last (now first) read is junk.
		uint32 ln = dropDeadFirstRead(AG, tig, false); // Decide if the last (now first) read is junk.
		tig->reverseComplement(); // Flip back.

		if ((fn == 0) && (ln == 0)) // Nothing to remove, just get out of here.
		if ((fn == 0) && (ln == 0)) { // Nothing to remove, nothing to do.
		writeLog("dropDead()-- both ends confirmed\n");
		continue;
		}

		if (fn == ln) { // The same thing to remove, ignore it.
		writeLog("dropDead()-- retaining spanning read %u\n", fn);
		continue;
		}

		// At least one read needs to be kicked out. Make new tigs for everything.

		char fnMsg[80] = {0}; Unitig *fnTig = NULL;
		char nnMsg[80] = {0}; Unitig *nnTig = NULL; int32 nnOff = INT32_MAX;
		char lnMsg[80] = {0}; Unitig *lnTig = NULL;
		Unitig *fnTig = NULL;
		Unitig *nnTig = NULL; int32 nnOff = INT32_MAX;
		Unitig *lnTig = NULL;

		if (fn > 0)
		fnTig = tigs.newUnitig(false);
		@@ -256,32 +294,26 @@ dropDeadEnds(AssemblyGraph *AG,

		// Move reads to their new unitig.

		strcpy(fnMsg, " ");
		strcpy(nnMsg, " ");
		strcpy(lnMsg, "");

		for (uint32 cc=0, tt=0; tt<tig->ufpath.size(); tt++) {
		ufNode &read = tig->ufpath[tt];

		if (read.ident == fn) {
		sprintf(fnMsg, "first read %9u to tig %7u --", read.ident, fnTig->id());
		writeLog("dropDead()-- tig %u gets first read %u\n", fnTig->id(), read.ident);
		fnTig->addRead(read, -read.position.min(), false);

		} else if (read.ident == ln) {
		sprintf(lnMsg, "-- last read %9u to tig %7u", read.ident, lnTig->id());
		writeLog("dropDead()-- tig %u gets last read %u\n", lnTig->id(), read.ident);
		lnTig->addRead(read, -read.position.min(), false);

		} else {
		if (nnOff == INT32_MAX) {
		sprintf(nnMsg, "other reads to tig %7u", nnTig->id());
		writeLog("dropDead()-- tig %u gets all other reads\n", nnTig->id());
		nnOff = read.position.min();
		}
		nnTig->addRead(read, -nnOff, false);
		}
		}

		writeLog("dropDeadEnds()-- tig %7u --> %s %s %s\n", tig->id(), fnMsg, nnMsg, lnMsg);

		if (fnTig) fnTig->cleanUp(); // Probably not neeeded, but cheap.
		if (lnTig) lnTig->cleanUp(); // Probably not neeeded, but cheap.
		if (nnTig) nnTig->cleanUp(); // Most likely needed.
		@@ -292,6 +324,6 @@ dropDeadEnds(AssemblyGraph *AG,
		tigs[ti] = NULL;
		}

		writeStatus("dropDeadEnds()-- Modified %u tigs. Dropped %u first and %u last reads, %u tig%s had both reads dropped.\n",
		writeStatus("dropDead()-- Modified %u tigs. Dropped %u first and %u last reads, %u tig%s had both reads dropped.\n",
		numT, numF, numL, numB, (numB == 1) ? "" : "s");
		}

src/bogart/AS_BAT_MarkRepeatReads.C

+4 −1

Original line number	Diff line number	Diff line
		@@ -708,6 +708,9 @@ findConfusedEdges(TigVector &tigs,
		continue;

		// Skip if this overlap is the best we're trying to match.
		//
		// NOTE. This doesn't care about potential duplicate overlaps between a pair of reads,
		// as we're looking for thicker overlaps off only one end of the read.
		if ((rdBid == b5->readId()) \|\|
		(rdBid == b3->readId()))
		continue;