Added "coverage" virtual track function that calculates the fraction …

.gitignore

@@ -10,3 +10,4 @@ docs
 .DS_Store
 mm10
 add_include.sh
+.aider*

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: misha
 Title: Toolkit for Analysis of Genomic Data
-Version: 4.3.2
+Version: 4.3.3
 Authors@R: c(
     person("Misha", "Hoichman", , "[email protected]", role = "aut"),
     person("Aviezer", "Lifshitz", , "[email protected]", role = c("aut", "cre")),

NEWS.md

@@ -1,3 +1,7 @@
+# misha 4.3.3
+
+* Added "coverage" virtual track function that calculates the fraction that iterator intervals are covered by given source intervals.
+
 # misha 4.3.2 
 
 * Added requirement for C++14 standard.

R/vtrack.R

@@ -162,6 +162,12 @@
 #' strand (-1 or 1) of the interval. Distance is always positive if 'strand' is
 #' '0' or if 'strand' column is missing.
 #'
+#' \emph{src = [1D intervals], func = "coverage", params = NULL} \cr
+#' For each iterator interval, calculates the fraction of its length that is covered by the
+#' source intervals. Returns a value between 0 and 1. For example, if an iterator interval is [100,200]
+#' and the source intervals cover positions 120-140 and 160-170, the coverage would be 0.3
+#' ((20 + 10) / 100 = 0.3). Overlapping source intervals are first unified.
+#'
 #' \emph{func = "pwm", params = list(pssm = matrix, bidirect = TRUE,
 #' prior = 0.01, extend = TRUE)} \cr
 #' Calculates total log-likelihood score of DNA sequence against PSSM.
@@ -230,6 +236,9 @@
 #'     seq(-500, 500, 200)
 #' )
 #'
+#' gvtrack.create("cov", "annotations", "coverage")
+#' gextract("cov", gintervals(1, 0, 1000), iterator = 100)
+#'
 #' pssm <- matrix(
 #'     c(
 #'         0.7, 0.1, 0.1, 0.1, # Example PSSM

src/TrackExpressionVars.cpp

@@ -22,7 +22,7 @@
 const char *TrackExpressionVars::Track_var::FUNC_NAMES[TrackExpressionVars::Track_var::NUM_FUNCS] = {
 		"avg", "min", "max", "nearest", "stddev", "sum", "quantile", "global.percentile", "global.percentile.min", "global.percentile.max", "weighted.sum", "area", "pwm", "pwm.max", "pwm.max.pos"
 };
-const char *TrackExpressionVars::Interv_var::FUNC_NAMES[TrackExpressionVars::Interv_var::NUM_FUNCS] = { "distance", "distance.center" };
+const char *TrackExpressionVars::Interv_var::FUNC_NAMES[TrackExpressionVars::Interv_var::NUM_FUNCS] = { "distance", "distance.center", "coverage" };
 
 using namespace rdb;
 
@@ -621,8 +621,20 @@ TrackExpressionVars::Interv_var &TrackExpressionVars::add_vtrack_var_src_interv(
 			if (iinterv->do_touch(*(iinterv - 1)))
 				verror("Virtual track %s: intervals are overlapping and hence incompatible with %s function", vtrack.c_str(), func.c_str());
 		}
-	} else
+	} else if (!strcmp(func.c_str(), Interv_var::FUNC_NAMES[Interv_var::COVERAGE])) {
+        var.val_func = Interv_var::COVERAGE;
+
+        if (!Rf_isNull(rparams))
+            verror("Virtual track %s: function %s does not accept any parameters", vtrack.c_str(), func.c_str());
+
+        var.dist_margin = 0.;
+        var.sintervs.swap(intervs1d);
+        var.sintervs.sort();
+        var.sintervs.unify_overlaps(); // Unify overlaps since we want total coverage
+
+	} else {
 		verror("Virtual track %s: invalid function %s used with intervals", vtrack.c_str(), func.c_str());
+	}
 
 	return var;
 }
@@ -1243,5 +1255,70 @@ void TrackExpressionVars::set_vars(unsigned idx)
 				ivar->var[idx] = dist;
 			}
 		}
+		else if (ivar->val_func == Interv_var::COVERAGE)
+		{
+			const GInterval &interval = ivar->imdf1d ? ivar->imdf1d->interval : m_interval1d;
+
+			if (ivar->imdf1d && ivar->imdf1d->out_of_range)
+			{
+				ivar->var[idx] = 0;
+				continue;
+			}
+
+			int64_t total_overlap = 0;
+			GIntervals::const_iterator iinterv;
+
+			// For non-sequential access or first access
+			if (ivar->imdf1d || ivar->siinterv == ivar->sintervs.end())
+			{
+				iinterv = lower_bound(ivar->sintervs.begin(), ivar->sintervs.end(), interval,
+									  GIntervals::compare_by_start_coord);
+
+				// Check previous interval too
+				if (iinterv != ivar->sintervs.begin())
+				{
+					auto prev = iinterv - 1;
+					if (prev->chromid == interval.chromid && prev->end > interval.start)
+					{
+						int64_t overlap_start = max(interval.start, prev->start);
+						int64_t overlap_end = min(interval.end, prev->end);
+						total_overlap += overlap_end - overlap_start;
+					}
+				}
+			}
+			else
+			{
+				// For sequential access, start from last position
+				iinterv = ivar->siinterv;
+				// But check if we need to back up
+				while (iinterv != ivar->sintervs.begin() &&
+					   (iinterv - 1)->chromid == interval.chromid &&
+					   (iinterv - 1)->end > interval.start)
+				{
+					--iinterv;
+				}
+			}
+
+			// Check forward intervals
+			while (iinterv != ivar->sintervs.end() &&
+				   iinterv->chromid == interval.chromid &&
+				   iinterv->start < interval.end)
+			{
+				if (iinterv->end > interval.start)
+				{
+					int64_t overlap_start = max(interval.start, iinterv->start);
+					int64_t overlap_end = min(interval.end, iinterv->end);
+					total_overlap += overlap_end - overlap_start;
+				}
+				++iinterv;
+			}
+
+			if (!ivar->imdf1d)
+			{
+				ivar->siinterv = iinterv;
+			}
+
+			ivar->var[idx] = (double)total_overlap / (interval.end - interval.start);
+		}
 	}
-	}
+}

src/TrackExpressionVars.h

@@ -106,7 +106,7 @@ class TrackExpressionVars {
     typedef vector<Track_var> Track_vars;
 
     struct Interv_var {
-        enum Val_func { DIST, DIST_CENTER, NUM_FUNCS };
+        enum Val_func { DIST, DIST_CENTER, COVERAGE, NUM_FUNCS };
 
         static const char *FUNC_NAMES[NUM_FUNCS];