fixed build, enabled constraining flux and k changes at once

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: dycone
 Title: dycone - analyze the dynamic landscape of metabolome data
-Version: 1.3.2
+Version: 1.4.0
 Authors@R: person("Christian", "Diener", , "[email protected]", role = c("aut", "cre"))
 Description: Dycone ("dynamic cone") allows you infer enzymatic regulation from
     metabolome mesurements. It employs formalisms based on flux and k-cone analysis

NAMESPACE

@@ -25,6 +25,7 @@ export(kcone_null)
 export(ma_terms)
 export(make_irreversible)
 export(mass_action)
+export(minimal_perturbation)
 export(minimal_perturbation_lp)
 export(occupation)
 export(patch)

R/kcone.R

@@ -497,11 +497,7 @@ eigendynamics <- function(basis, n = 1) {
 #'  and log2-fold changes of k2 relative to k1. Log-fold changes where either of the
 #'  constants is 0 are evaluated to 0.
 #' @examples
-#' data(eryth)
-#' n <- length(make_irreversible(eryth))
-#' k1 <- runif(n)
-#' k2 <- runif(n)
-#' single_hyp(k1, k2, eryth)
+#'  NULL
 single_hyp <- function(k1, k2, reacts) {
     reacts <- make_irreversible(reacts)
     logfold <- k2 - k1

R/minimal_perturbation.R

@@ -2,6 +2,9 @@
 #
 # MIT license. See LICENSE for more information.
 
+# to make R CMD CHECK happy :/
+utils::globalVariables(c("disease", "normal"))
+
 #' Assemble the LP problem for minimal perturbation.
 #'
 #' @param ma_terms A matrix or data frame containing the metabolic terms in the
@@ -10,14 +13,19 @@
 #'  "disease" entries.
 #' @param S The stoichiometric matrix or a list containing one matrix for each
 #'  sample.
-#' @param bounds The bounds for the fluxes. Either a single value for the max
-#'  bound or a list each containing a vector specifying the bounds for the
-#'  specific sample.
-#' @param tradeoff Double in [0, 1] specifiying the equlibrium between
-#'  minimizing differences in k changes versus flux changes.
-#' @param growth Constraints fro growth rate. If not NA must be a named vector
-#'  with entries "idx", "min", "max" denoting the index of the biomass reaction
-#'  and its minimum and maximum flux respectively.
+#' @param v_min The lower bounds for the fluxes. Either a single value for the
+#'  min bounds for all reactions or a matrix/data.frame with two columns and as
+#'  many rows as irreversible reactions where the columns contain the lower
+#'  bounds for each factor level in `samples`.
+#' @param v_max The upper bounds for the fluxes. Either a single value for the
+#'  max bounds for all reactions or a matrix/data.frame with two columns and as
+#'  many rows as irreversible reactions where the columns contain the upper
+#'  bounds for each factor level in `samples`.
+#' @param min_obj Constraints for growth rate. If a vector with two elements
+#'  the first will be interpreted as the index of the objective reaction and
+#'  the second value as its lower bound. If only a single value will constrain
+#'  the absolute sum of fluxes (sum |v_i|) and take the value as the lower
+#'  bound.
 #' @return A list with the following components:
 #'  \describe{
 #'  \item{coefficients}{The coefficients for the equalities/inequalities in
@@ -45,9 +53,9 @@ minimal_perturbation_lp <- function(ma_terms, samples, S, v_min = 0,
         }
         S <- new_S
     }
-    n_var <- ncol(S[[1]])
+    n_vars <- ncol(S[[1]])
     coefficients <- bdiag(S)
-    coefficients <- cbind(coefficients, Matrix(0, ncol = n_var,
+    coefficients <- cbind(coefficients, Matrix(0, ncol = 2 * n_vars,
                           nrow = nrow(coefficients)))
     type <- rep("equal", nrow(coefficients))
     row_bounds <- rep(0, nrow(coefficients))
@@ -67,9 +75,9 @@ minimal_perturbation_lp <- function(ma_terms, samples, S, v_min = 0,
         v_min[idx] <- min_obj[2]
         v_max[idx] <- pmax(v_max[idx], min_obj[2])
     } else {
-        coefs <- matrix(0, nrow = 2, ncol = 3 * n_var)
+        coefs <- matrix(0, nrow = 2, ncol = 4 * n_vars)
         for (i in 1:2) {
-            coefs[i, ((i - 1) * n_var + 1):(i * n_var)] <- 1
+            coefs[i, ((i - 1) * n_vars + 1):(i * n_vars)] <- 1
         }
         coefficients <- rbind(coefficients, coefs)
         type <- append(type, rep("larger", 2))
@@ -91,14 +99,24 @@ minimal_perturbation_lp <- function(ma_terms, samples, S, v_min = 0,
 #' @param samples A factor or character string with either "normal" and
 #'  "disease" entries.
 #' @param reacts The reactions for the respective model.
-#' @param bounds The bounds for the fluxes. Either a single value for the max
-#'  bound or a list each containing a vector specifying the bounds for the
-#'  specific sample.
+#' @param permutations The maximum number of permutations. If this number is
+#'  smaller than all possible permutation will check all permutations,
+#'  otherwise will sample that many permutations *with replacements*.
+#' @param v_min The lower bounds for the fluxes. Either a single value for the
+#'  min bounds for all reactions or a matrix/data.frame with two columns and as
+#'  many rows as irreversible reactions where the columns contain the lower
+#'  bounds for each factor level in `samples`.
+#' @param v_max The upper bounds for the fluxes. Either a single value for the
+#'  max bounds for all reactions or a matrix/data.frame with two columns and as
+#'  many rows as irreversible reactions where the columns contain the upper
+#'  bounds for each factor level in `samples`.
 #' @param tradeoff Double in [0, 1] specifiying the equlibrium between
 #'  minimizing differences in k changes versus flux changes.
-#' @param growth Constraints fro growth rate. If not NA must be a named vector
-#'  with entries "idx", "min", "max" denoting the index of the biomass reaction
-#'  and its minimum and maximum flux respectively.
+#' @param min_obj Constraints for growth rate. If a vector with two elements
+#'  the first will be interpreted as the index of the objective reaction and
+#'  the second value as its lower bound. If only a single value will constrain
+#'  the absolute sum of fluxes (sum |v_i|) and take the value as the lower
+#'  bound.
 #' @return A list with the following components:
 #'  \describe{
 #'  \item{fluxes}{The alterations in fluxes.}

src/min_pertub.cpp

@@ -10,28 +10,44 @@ using namespace Rcpp;
 
 const double ztol = 1.0e-9;
 
-
+/*
+ * Update the terms for |k^i_d - k^i_n| using the current metabolic terms.
+ */
 void update_problem(glp_prob* lp, NumericVector cond_terms,
                     NumericVector norm_terms, int neq, double tradeoff) {
     int size = cond_terms.size();
     for(int i=1; i <= size; i++) {
         NumericVector vals = NumericVector::create(cond_terms[i-1],
                                                    norm_terms[i-1]);
-        double obj_coef = tradeoff;
-        if(is_true(any(is_na(vals)))) {
-            vals[0] = 1.0;
-            vals[1] = 1.0;
-            obj_coef = 1.0 - tradeoff;
+
+        if (is_false(any(is_na(vals)))) {
+            int idx[4] = {0, i, i + size, i + 3*size};
+            double coefs[4] = {0.0, 1.0/vals[0], -1.0/vals[1], 1.0};
+            glp_set_mat_row(lp, 2*i - 1 + neq, 3, idx, coefs);
+            glp_set_row_bnds(lp, 2*i - 1 + neq, GLP_LO, 0.0, 0.0);
+            coefs[3] = -1.0;
+            glp_set_mat_row(lp, 2*i + neq, 3, idx, coefs);
+            glp_set_row_bnds(lp, 2*i + neq, GLP_UP, 0.0, 0.0);
+            glp_set_obj_coef(lp, 3*size + i, tradeoff);
         }
+    }
+}
+
+/*
+ * Add the terms for |v^i_d - v^i_n|.
+ */
+void add_flux_diffs(glp_prob* lp, int neq, double weight) {
+    int size = glp_get_num_cols(lp) / 4;
 
+    for(int i=1; i <= size; i++) {
         int idx[4] = {0, i, i + size, i + 2*size};
-        double coefs[4] = {0.0, 1.0/vals[0], -1.0/vals[1], 1.0};
+        double coefs[4] = {0.0, 1.0, -1.0, 1.0};
         glp_set_mat_row(lp, 2*i - 1 + neq, 3, idx, coefs);
         glp_set_row_bnds(lp, 2*i - 1 + neq, GLP_LO, 0.0, 0.0);
         coefs[3] = -1.0;
         glp_set_mat_row(lp, 2*i + neq, 3, idx, coefs);
         glp_set_row_bnds(lp, 2*i + neq, GLP_UP, 0.0, 0.0);
-        glp_set_obj_coef(lp, 2*size + i, obj_coef);
+        glp_set_obj_coef(lp, 2*size + i, weight);
     }
 }
 
@@ -47,7 +63,7 @@ SEXP glpk_min_perturb(IntegerVector ridx, IntegerVector cidx,
     glp_term_out(GLP_OFF);
     glp_prob* lp = glp_create_prob();
     glp_set_prob_name(lp, "minimal perturbation");
-    glp_add_rows(lp, nrows + 2*n_vars);
+    glp_add_rows(lp, nrows + 4*n_vars);
     glp_add_cols(lp, ncols);
 
     double bnd;
@@ -87,6 +103,8 @@ SEXP glpk_min_perturb(IntegerVector ridx, IntegerVector cidx,
     NumericMatrix cond_fluxes(n_vars, perms.nrow());
     NumericMatrix norm_fluxes(n_vars, perms.nrow());
     NumericVector obj_vals(perms.nrow());
+    add_flux_diffs(lp, nrows, 1.0 - tradeoff);
+    nrows += 2 * n_vars;
     Rcout<<"Running permutations";
     for(int p = 0; p < perms.nrow(); p++) {
         Rcout<<".";

tests/testthat/test_min_perturb.R

@@ -0,0 +1,14 @@
+# Copyright 2015 Christian Diener <[email protected]>
+# MIT license. See LICENSE for more information.
+
+context("minimum perturbation")
+
+test_that("random data gives no significant results", {
+    mats <- matrix(rnorm(6 * 68, 10, 1), ncol = 6)
+    samples <- factor(rep(c("disease", "normal"), each = 3))
+    perturb <- minimal_perturbation(mats, samples, eryth)
+    expect_true(all(perturb$k$P.value < 0.05))
+    expect_true(all(perturb$k$adj.P.Val > 0.99))
+    expect_true(all(perturb$fluxes$P.value < 0.05))
+    expect_true(all(perturb$fluxes$adj.P.Val > 0.99))
+})

vignettes/eryth.Rmd

@@ -103,10 +103,10 @@ set.seed(42)
 exp_concs = cbind(noise(d_concs[1,]), noise)
 mats <- cbind(replicate(4, ma_terms(S, noise(d_concs[1,]))),
               replicate(4, ma_terms(S, noise(d_concs[3,]))))
-samples <- rep(c("normal", "disease"), each = 4)
+samples <- factor(rep(c("normal", "disease"), each = 4))
 
-h = hyp(reacts, samples, mats)
-print(h[h$pval<0.05,])
+p = minimal_perturbation(mats, samples, reacts)
+print(p$k[p$k$adj.P.Val < 0.05,])
 ```
 
 We see that the most prominent change is an increased reaction rate for two reactions,