Add function get_ratio_variance

Author: wulingyun

NAMESPACE

@@ -7,6 +7,7 @@ export(get_adjusted_deg_diff)
 export(get_diff_ratio_net)
 export(get_ratio_distribution)
 export(get_ratio_distribution2)
+export(get_ratio_variance)
 export(get_shortest_distances)
 export(get_subnets)
 export(kappa_score)

R/deg_normalization.R

@@ -0,0 +1,9 @@
+
+
+
+#' @export
+get_ratio_variance <- function(expr.matrix, log.expr = FALSE)
+{
+  if (!log.expr) expr.matrix <- log(expr.matrix)
+  .Call(ND_RatioVariance, expr.matrix)
+}

src/Corbi.cpp

@@ -11,6 +11,7 @@ const R_CallMethodDef callMethods[] = {
   {"ND_RatioDistribution", (DL_FUNC) &ND_RatioDistribution, 2},
   {"ND_RatioDistribution2", (DL_FUNC) &ND_RatioDistribution2, 3},
   {"ND_DiffRatioNet", (DL_FUNC) &ND_DiffRatioNet, 2},
+  {"ND_RatioVariance", (DL_FUNC) &ND_RatioVariance, 1},
   {NULL, NULL, 0}
 };
 

src/Corbi.h

@@ -26,4 +26,5 @@ extern "C" {
   SEXP ND_RatioDistribution(SEXP _LogExprMatrix, SEXP _pEdge);
   SEXP ND_RatioDistribution2(SEXP _LogExprMatrix, SEXP _pEdge, SEXP _pTrim);
   SEXP ND_DiffRatioNet(SEXP _RatioLB, SEXP _LogExprVal);
+  SEXP ND_RatioVariance(SEXP _LogExprMatrix);
 }

src/NetDEG.cpp

@@ -203,3 +203,52 @@ SEXP ND_DiffRatioNet(SEXP _RatioLB, SEXP _LogExprVal)
   UNPROTECT(5);
   return(_M);
 }
+
+SEXP ND_RatioVariance(SEXP _LogExprMatrix)
+{
+  PROTECT(_LogExprMatrix = AS_NUMERIC(_LogExprMatrix));
+  double *LogExprMatrix = NUMERIC_POINTER(_LogExprMatrix);
+  int *dim = INTEGER_POINTER(AS_INTEGER(GET_DIM(_LogExprMatrix)));
+  int nGenes = dim[0];
+  int nSamples = dim[1];
+
+  SEXP _Var;
+  PROTECT(_Var = NEW_NUMERIC(nGenes * nGenes));
+  SetDim2(_Var, nGenes, nGenes);
+  double *Var = NUMERIC_POINTER(_Var);
+  SetValues(_Var, Var, 0.0);
+  
+  double *r = (double *) R_alloc(nSamples, sizeof(double));
+  double *e, m, v;
+  int n;
+  for (int i = 0; i < nGenes-1; i++)
+  {
+    for (int j = i+1; j < nGenes; j++)
+    {
+      e = LogExprMatrix;
+      n = 0;
+      m = v = 0;
+      for (int k = 0; k < nSamples; k++)
+      {
+        if (R_finite(e[i]) && R_finite(e[j]))
+        {
+          r[n] = e[i] - e[j];
+          m += r[n];
+          v += r[n] * r[n];
+          n++;
+        }
+        e += nGenes;
+      }
+      
+      if (n > 0)
+      {
+        v = (v - m * m / n)  / (n - 1);
+        Var[i+nGenes*j] = v;
+        Var[j+nGenes*i] = v;
+      }
+    }
+  }
+
+  UNPROTECT(2);
+  return(_Var);
+}