ultra-deep coverage

Author: tobiasrausch

Makefile

@@ -33,7 +33,9 @@ endif
 ifeq (${EBROOTHTSLIB}, ${PWD}/src/htslib/)
 	SUBMODULES += .htslib
 endif
-
+ifdef COVTYPE
+	CXXFLAGS += -DCOVTYPE=${COVTYPE}
+endif
 
 # External sources
 HTSLIBSOURCES = $(wildcard src/htslib/*.c) $(wildcard src/htslib/*.h)

src/bamstats.h

@@ -211,7 +211,10 @@ namespace bamstats
 		if (itRg->second.bc.cov[i] == 1) ++itRg->second.bc.n1;
 		if (itRg->second.bc.cov[i] == 2) ++itRg->second.bc.n2;
 	      }
-	      if (!nrun[i]) ++itRg->second.bc.bpWithCoverage[itRg->second.bc.cov[i]];
+	      if (!nrun[i]) {
+		if (itRg->second.bc.cov[i] >= itRg->second.bc.bpWithCoverage.size()) itRg->second.bc.bpWithCoverage.resize(itRg->second.bc.cov[i] + 1, 0); 
+		++itRg->second.bc.bpWithCoverage[itRg->second.bc.cov[i]];
+	      }
 	    }
 	    itRg->second.bc.cov.clear();
 	  }
@@ -567,7 +570,10 @@ namespace bamstats
 	    if (itRg->second.bc.cov[i] == 1) ++itRg->second.bc.n1;
 	    if (itRg->second.bc.cov[i] == 2) ++itRg->second.bc.n2;
 	  }
-	  if (!nrun[i]) ++itRg->second.bc.bpWithCoverage[itRg->second.bc.cov[i]];
+	  if (!nrun[i]) {
+	    if (itRg->second.bc.cov[i] >= itRg->second.bc.bpWithCoverage.size()) itRg->second.bc.bpWithCoverage.resize(itRg->second.bc.cov[i] + 1, 0); 
+	    ++itRg->second.bc.bpWithCoverage[itRg->second.bc.cov[i]];
+	  }
 	}
 	itRg->second.bc.cov.clear();
       }

src/json.h

@@ -540,13 +540,14 @@ namespace bamstats
 
       // Coverage Histogram
       {
-	uint32_t lastValidCO = _lastNonZeroIdx(itRg->second.bc.bpWithCoverage, itRg->second.bc.maxCoverage);
-	float lastFrac = _lastPercentage(itRg->second.bc.bpWithCoverage, itRg->second.bc.maxCoverage);
+	// Ignore last bucket because it counts >= the max. coverage value
+	uint32_t lastValidCO = _lastNonZeroIdx(itRg->second.bc.bpWithCoverage, itRg->second.bc.bpWithCoverage.size() - 1);
+	float lastFrac = _lastPercentage(itRg->second.bc.bpWithCoverage);
 	nlohmann::json j;
 	j["id"] = "coverageHistogram";
 	j["title"] = "Coverage histogram";
 	if (lastFrac > 0) {
-	  std::string rlstr = boost::lexical_cast<std::string>(lastFrac) + "% of all bases with >= " + boost::lexical_cast<std::string>(itRg->second.bc.maxCoverage) + "x coverage";
+	  std::string rlstr = boost::lexical_cast<std::string>(lastFrac) + "% of all bases with >= " + boost::lexical_cast<std::string>(itRg->second.bc.bpWithCoverage.size() - 1) + "x coverage";
 	  j["subtitle"] = rlstr;
 	}
 	j["x"]["data"] = nlohmann::json::array();

src/qcstruct.h

@@ -52,10 +52,10 @@ namespace bamstats
     typedef uint32_t TCountType;
     typedef std::vector<TCountType> TCoverageBp;
 
-    typedef uint16_t TMaxCoverage;
+    typedef COVTYPE TMaxCoverage;
     typedef std::vector<TMaxCoverage> TBpCoverage;
 
-    uint32_t maxCoverage;
+    TMaxCoverage maxCoverage;
     uint32_t maxIndelSize;
     uint64_t n1;
     uint64_t n2;
@@ -76,7 +76,7 @@ namespace bamstats
     BaseCounts() : maxCoverage(std::numeric_limits<TMaxCoverage>::max()), maxIndelSize(50), n1(0), n2(0), nd(0), matchCount(0), mismatchCount(0), delCount(0), insCount(0), softClipCount(0), hardClipCount(0) {
       delHomACGTN.resize(6, 0);
       insHomACGTN.resize(6, 0);
-      bpWithCoverage.resize(maxCoverage + 1, 0);
+      bpWithCoverage.resize(std::numeric_limits<uint16_t>::max() + 1, 0);
       delSize.resize(maxIndelSize + 1, 0);
       insSize.resize(maxIndelSize + 1, 0);
       cov.clear();

src/tsv.h

@@ -57,14 +57,13 @@ namespace bamstats
       if (vec[i] > 0) lastNonZeroIdx = i;
     return lastNonZeroIdx;
   }
-  
+
   template<typename TVector>
   inline uint32_t
   _lastNonZeroIdx(TVector const& vec) {
     return _lastNonZeroIdx(vec, vec.size());
   }
 
-
   template<typename TVector>
   inline float
   _lastPercentage(TVector const& vec, uint32_t const lastIdx) {
@@ -73,6 +72,12 @@ namespace bamstats
     return ((float) (vec[lastIdx]) * 100.0) / cumsum;
   }
 
+  template<typename TVector>
+  inline float
+  _lastPercentage(TVector const& vec) {
+    return _lastPercentage(vec, vec.size() - 1);
+  }
+
   template<typename TVector>
   inline uint32_t
   _lastNonZeroIdxACGTN(TVector const& vec, uint32_t const read12) {

src/util.h

@@ -17,6 +17,10 @@
 namespace bamstats
 {
 
+  #ifndef COVTYPE
+  #define COVTYPE uint16_t
+  #endif
+
   struct Interval {
     int32_t start;
     int32_t end;