From b85112e250eca9259b663fd6b0331d858eab3ea8 Mon Sep 17 00:00:00 2001 From: Thom Quinn Date: Sun, 5 Apr 2020 11:43:12 +1000 Subject: [PATCH] ivar = NA will now skip log-ratio transform --- DESCRIPTION | 4 +- NEWS.md | 4 ++ R/0-params.R | 3 +- R/1-propr.R | 81 +++++++++++++++++++-------------- R/6-propdViz.R | 1 - man/all.Rd | 3 +- man/calculateTheta.Rd | 11 ++++- man/getNetwork.Rd | 18 ++++++-- man/getRatios.Rd | 3 +- man/ivar2index.Rd | 3 +- man/pra.Rd | 9 +++- man/propd.Rd | 3 +- man/propr.Rd | 14 ++++-- man/proprPerb.Rd | 3 +- man/visualize.Rd | 20 +++++--- tests/testthat/test-multiomic.R | 32 +++++++++++++ 16 files changed, 151 insertions(+), 61 deletions(-) create mode 100644 tests/testthat/test-multiomic.R diff --git a/DESCRIPTION b/DESCRIPTION index 86b1816..2ed9f1c 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: propr Title: Calculating Proportionality Between Vectors of Compositional Data -Version: 4.2.6 +Version: 4.2.7 URL: http://github.com/tpq/propr BugReports: http://github.com/tpq/propr/issues Authors@R: c( @@ -22,7 +22,7 @@ Description: The bioinformatic evaluation of gene co-expression often begins wit License: GPL-2 LazyData: true VignetteBuilder: knitr -RoxygenNote: 6.1.1 +RoxygenNote: 7.0.2 Encoding: UTF-8 Depends: methods, diff --git a/NEWS.md b/NEWS.md index 85215f4..7001614 100755 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +## propr 4.2.7 +--------------------- +* Allow user to skip log-ratio transformation with `ivar = NA` + ## propr 4.2.6 --------------------- * Fix new R version `class` bug where matrix is also an array diff --git a/R/0-params.R b/R/0-params.R index 75ac091..8d368d5 100644 --- a/R/0-params.R +++ b/R/0-params.R @@ -12,7 +12,8 @@ #' for additive log-ratio transformation. The argument will also #' accept feature name(s) instead of the index position(s). #' Set to "iqlr" to use inter-quartile log-ratio transformation. -#' Ignore to use centered log-ratio transformation. +#' Ignore to use centered log-ratio transformation. If the data +#' are already transformed, skip it with \code{ivar = NA}. #' @param select Optional. Use this to subset the final #' proportionality matrix without altering the result. #' Use this argument to rearrange feature order. diff --git a/R/1-propr.R b/R/1-propr.R index 8d5912d..19e9258 100755 --- a/R/1-propr.R +++ b/R/1-propr.R @@ -128,49 +128,59 @@ propr <- function(counts, metric = c("rho", "phi", "phs", "cor", "vlr"), ivar = select, symmetrize = FALSE, alpha, p = 100){ # Clean "count matrix" - # if(any(apply(counts, 2, function(x) all(x == 0)))){ - # stop("Remove components with all zeros before proceeding.")} - if(any(counts < 0)) stop("Data may not contain negative measurements.") - if(any(is.na(counts))) stop("Remove NAs from 'counts' before proceeding.") if("data.frame" %in% class(counts)) counts <- as.matrix(counts) if(is.null(colnames(counts))) colnames(counts) <- as.character(1:ncol(counts)) if(is.null(rownames(counts))) rownames(counts) <- as.character(1:nrow(counts)) - if(missing(alpha)){ alpha <- NA - }else{ if(!is.na(alpha)) if(alpha == 0) alpha <- NA } - ct <- counts - - # Replace zeros unless alpha is provided - if(any(as.matrix(counts) == 0) & is.na(alpha)){ - message("Alert: Replacing 0s with next smallest value.") - zeros <- ct == 0 - ct[zeros] <- min(ct[!zeros]) - } + if(any(is.na(counts))) stop("Remove NAs from 'counts' before proceeding.") - # Establish reference - use <- ivar2index(ct, ivar) + if(!is.na(ivar)){ # if a log-ratio transform is performed... - # Transform counts based on reference - if(is.na(alpha)){ + if(any(counts < 0)) stop("Data may not contain negative measurements.") - # Use g(x) = Mean[log(x)] to log-ratio transform data - message("Alert: Saving log-ratio transformed counts to @logratio.") - logX <- log(ct) - logSet <- logX[, use, drop = FALSE] - ref <- rowMeans(logSet) - lr <- sweep(logX, 1, ref, "-") + if(missing(alpha)){ alpha <- NA + }else{ if(!is.na(alpha)) if(alpha == 0) alpha <- NA } + ct <- counts - }else{ + # Replace zeros unless alpha is provided + if(any(as.matrix(counts) == 0) & is.na(alpha)){ + message("Alert: Replacing 0s with next smallest value.") + zeros <- ct == 0 + ct[zeros] <- min(ct[!zeros]) + } + + # Establish reference + use <- ivar2index(ct, ivar) + + # Transform counts based on reference + if(is.na(alpha)){ - # Since y = log(x) = [x^a-1]/a, ref = Mean[log(x)] = Mean[y] - # Calculate log(x/ref) = log(x) - log(ref) = [x^a-1]/a - [ref^a-1]/a - message("Alert: Saving alpha transformed counts to @logratio.") - aX <- (ct^alpha - 1) / alpha - aSet <- aX[, use, drop = FALSE] - ref <- rowMeans(aSet) - lr <- sweep(aX, 1, ref, "-") + # Use g(x) = Mean[log(x)] to log-ratio transform data + message("Alert: Saving log-ratio transformed counts to @logratio.") + logX <- log(ct) + logSet <- logX[, use, drop = FALSE] + ref <- rowMeans(logSet) + lr <- sweep(logX, 1, ref, "-") + + }else{ + + # Since y = log(x) = [x^a-1]/a, ref = Mean[log(x)] = Mean[y] + # Calculate log(x/ref) = log(x) - log(ref) = [x^a-1]/a - [ref^a-1]/a + message("Alert: Saving alpha transformed counts to @logratio.") + aX <- (ct^alpha - 1) / alpha + aSet <- aX[, use, drop = FALSE] + ref <- rowMeans(aSet) + lr <- sweep(aX, 1, ref, "-") + } + + }else{ # if skipping the log-ratio transform... + + message("Alert: Skipping the log-ratio transformation.") + alpha <- NA + ct <- counts + lr <- ct } - # Optionally apply select for performance gain + # Optionally reduce the number of features used in matrix algebra if(!missing(select)){ message("Alert: Using 'select' disables permutation testing.") @@ -222,7 +232,8 @@ propr <- function(counts, metric = c("rho", "phi", "phs", "cor", "vlr"), ivar = # Set up @permutes result@permutes <- list(NULL) if(p > 0){ - # Sample compositions so that the clr does not change + + # Shuffle row-wise so that per-sample CLR does not change message("Alert: Fixing permutations to active random seed.") permutes <- vector("list", p) for(ins in 1:length(permutes)) permutes[[ins]] <- t(apply(ct, 1, sample)) @@ -242,7 +253,7 @@ propr <- function(counts, metric = c("rho", "phi", "phs", "cor", "vlr"), ivar = ) # Initialize @results -- Tally frequency of 0 counts - if(any(as.matrix(counts) == 0)){ + if(any(as.matrix(counts) == 0) & !is.na(ivar)){ message("Alert: Tabulating the presence of 0 counts.") result@results$Zeros <- ctzRcpp(as.matrix(counts)) # count 0s } diff --git a/R/6-propdViz.R b/R/6-propdViz.R index ebaf492..824c6fd 100755 --- a/R/6-propdViz.R +++ b/R/6-propdViz.R @@ -49,7 +49,6 @@ pals <- function(object, k){ #' where the \code{@@matrix} slot contains \eqn{1 - \theta}. #' Allows the user to interrogate theta using any #' visualization built for \code{propr} objects. -#' @rdname propd #' @export propd2propr <- function(object, ivar){ diff --git a/man/all.Rd b/man/all.Rd index 6171269..5cd00a1 100644 --- a/man/all.Rd +++ b/man/all.Rd @@ -15,7 +15,8 @@ to calculate. Choose from "rho", "phi", or "phs".} for additive log-ratio transformation. The argument will also accept feature name(s) instead of the index position(s). Set to "iqlr" to use inter-quartile log-ratio transformation. -Ignore to use centered log-ratio transformation.} +Ignore to use centered log-ratio transformation. If the data +are already transformed, skip it with \code{ivar = NA}.} \item{select}{Optional. Use this to subset the final proportionality matrix without altering the result. diff --git a/man/calculateTheta.Rd b/man/calculateTheta.Rd index 28d7432..5ded246 100755 --- a/man/calculateTheta.Rd +++ b/man/calculateTheta.Rd @@ -4,8 +4,15 @@ \alias{calculateTheta} \title{Calculate Theta} \usage{ -calculateTheta(counts, group, alpha = NA, lrv = NA, only = "all", - weighted = FALSE, weights = as.matrix(NA)) +calculateTheta( + counts, + group, + alpha = NA, + lrv = NA, + only = "all", + weighted = FALSE, + weights = as.matrix(NA) +) } \arguments{ \item{counts}{A data.frame or matrix. A "count matrix" with diff --git a/man/getNetwork.Rd b/man/getNetwork.Rd index 8dfb8f3..ed5b44d 100644 --- a/man/getNetwork.Rd +++ b/man/getNetwork.Rd @@ -4,9 +4,21 @@ \alias{getNetwork} \title{Get Network from Object} \usage{ -getNetwork(object, cutoff = NA, propr.object, propr.cutoff = NA, - thetad.object, thetad.cutoff = NA, thetae.object, thetae.cutoff = NA, - col1, col2, include = NA, or = TRUE, d3 = FALSE) +getNetwork( + object, + cutoff = NA, + propr.object, + propr.cutoff = NA, + thetad.object, + thetad.cutoff = NA, + thetae.object, + thetae.cutoff = NA, + col1, + col2, + include = NA, + or = TRUE, + d3 = FALSE +) } \arguments{ \item{object}{Any \code{propr} or \code{propd} object.} diff --git a/man/getRatios.Rd b/man/getRatios.Rd index 21c1091..099c026 100644 --- a/man/getRatios.Rd +++ b/man/getRatios.Rd @@ -4,8 +4,7 @@ \alias{getRatios} \title{Get (Log-)ratios from Object} \usage{ -getRatios(object, cutoff = NA, include = NA, or = TRUE, - melt = FALSE) +getRatios(object, cutoff = NA, include = NA, or = TRUE, melt = FALSE) } \arguments{ \item{object}{A \code{propr} or \code{propd} object.} diff --git a/man/ivar2index.Rd b/man/ivar2index.Rd index 2e0b116..6e59756 100755 --- a/man/ivar2index.Rd +++ b/man/ivar2index.Rd @@ -15,7 +15,8 @@ does not necessarily have to contain counts.} for additive log-ratio transformation. The argument will also accept feature name(s) instead of the index position(s). Set to "iqlr" to use inter-quartile log-ratio transformation. -Ignore to use centered log-ratio transformation.} +Ignore to use centered log-ratio transformation. If the data +are already transformed, skip it with \code{ivar = NA}.} } \value{ A numeric vector of indices for the reference set. diff --git a/man/pra.Rd b/man/pra.Rd index 605762b..c4b2819 100644 --- a/man/pra.Rd +++ b/man/pra.Rd @@ -4,8 +4,13 @@ \alias{pra} \title{Principal Ratio Analysis} \usage{ -pra(counts, ndim = 3, nclust = 2 * round(sqrt(ncol(counts))), - nsearch = 3, ndenom = 4) +pra( + counts, + ndim = 3, + nclust = 2 * round(sqrt(ncol(counts))), + nsearch = 3, + ndenom = 4 +) } \arguments{ \item{counts}{A data.frame or matrix. A "count matrix" with diff --git a/man/propd.Rd b/man/propd.Rd index 09e6290..d111544 100755 --- a/man/propd.Rd +++ b/man/propd.Rd @@ -62,7 +62,8 @@ whether to calculate a moderated F-statistic.} for additive log-ratio transformation. The argument will also accept feature name(s) instead of the index position(s). Set to "iqlr" to use inter-quartile log-ratio transformation. -Ignore to use centered log-ratio transformation.} +Ignore to use centered log-ratio transformation. If the data +are already transformed, skip it with \code{ivar = NA}.} } \value{ Returns a \code{propr} object. diff --git a/man/propr.Rd b/man/propr.Rd index eb41736..e3d573a 100755 --- a/man/propr.Rd +++ b/man/propr.Rd @@ -17,8 +17,15 @@ \usage{ \S4method{show}{propr}(object) -propr(counts, metric = c("rho", "phi", "phs", "cor", "vlr"), - ivar = "clr", select, symmetrize = FALSE, alpha, p = 100) +propr( + counts, + metric = c("rho", "phi", "phs", "cor", "vlr"), + ivar = "clr", + select, + symmetrize = FALSE, + alpha, + p = 100 +) phit(counts, ...) @@ -50,7 +57,8 @@ to calculate. Choose from "rho", "phi", or "phs".} for additive log-ratio transformation. The argument will also accept feature name(s) instead of the index position(s). Set to "iqlr" to use inter-quartile log-ratio transformation. -Ignore to use centered log-ratio transformation.} +Ignore to use centered log-ratio transformation. If the data +are already transformed, skip it with \code{ivar = NA}.} \item{select}{Optional. Use this to subset the final proportionality matrix without altering the result. diff --git a/man/proprPerb.Rd b/man/proprPerb.Rd index 38c8bb3..917ac57 100755 --- a/man/proprPerb.Rd +++ b/man/proprPerb.Rd @@ -15,7 +15,8 @@ does not necessarily have to contain counts.} for additive log-ratio transformation. The argument will also accept feature name(s) instead of the index position(s). Set to "iqlr" to use inter-quartile log-ratio transformation. -Ignore to use centered log-ratio transformation.} +Ignore to use centered log-ratio transformation. If the data +are already transformed, skip it with \code{ivar = NA}.} } \value{ Returns proportionality matrix. diff --git a/man/visualize.Rd b/man/visualize.Rd index dee33d8..73bf1d6 100755 --- a/man/visualize.Rd +++ b/man/visualize.Rd @@ -1,6 +1,5 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/3-proprViz.R, R/6-propdViz.R -\docType{methods} \name{visualize} \alias{visualize} \alias{plot,propr,missing-method} @@ -41,8 +40,17 @@ cytescape(object, col1, col2, prompt = TRUE, d3 = FALSE) propd2propr(object, ivar) -\S4method{plot}{propd,missing}(x, y, cutoff = 1000, col1, col2, propr, - prompt = TRUE, d3 = FALSE, plotSkip = FALSE) +\S4method{plot}{propd,missing}( + x, + y, + cutoff = 1000, + col1, + col2, + propr, + prompt = TRUE, + d3 = FALSE, + plotSkip = FALSE +) geyser(object, cutoff = 1000, k = 5, prompt = TRUE, plotly = FALSE) @@ -52,8 +60,7 @@ gemini(object, cutoff = 1000, k = 5, prompt = TRUE, plotly = FALSE) decomposed(object, cutoff = 1000) -parallel(object, cutoff = 1000, include = NA, or = TRUE, - plotly = FALSE) +parallel(object, cutoff = 1000, include = NA, or = TRUE, plotly = FALSE) } \arguments{ \item{x}{A \code{propr} or \code{propd} object.} @@ -92,7 +99,8 @@ to color \code{red} or \code{blue}, respectively.} for additive log-ratio transformation. The argument will also accept feature name(s) instead of the index position(s). Set to "iqlr" to use inter-quartile log-ratio transformation. -Ignore to use centered log-ratio transformation.} +Ignore to use centered log-ratio transformation. If the data +are already transformed, skip it with \code{ivar = NA}.} \item{cutoff}{For \code{updateCutoffs}, a numeric vector. this argument provides the FDR cutoffs to test. diff --git a/tests/testthat/test-multiomic.R b/tests/testthat/test-multiomic.R new file mode 100644 index 0000000..fc95cfd --- /dev/null +++ b/tests/testthat/test-multiomic.R @@ -0,0 +1,32 @@ +library(propr) +data(iris) +met.rel <- iris[,1:2] +mic.rel <- iris[,3:4] + +# Analyze multi-omics via back-end +clr <- function(x) sweep(log(x), 1, rowMeans(log(x)), "-") +REL <- cbind(clr(met.rel), clr(mic.rel)) +pr.r <- propr:::lr2rho(as.matrix(REL)) +colnames(pr.r) <- colnames(REL) +rownames(pr.r) <- colnames(REL) + +# Analyze multi-omics with wrapper +ivarNA <- propr(REL, ivar = NA) + +test_that("ivar = NA works as expected", { + + expect_equal( + NA, + ivarNA@ivar + ) + + expect_equal( + pr.r, + ivarNA@matrix + ) + + expect_equal( + ivarNA@counts, + ivarNA@logratio + ) +})