diff --git a/R/msqrob.R b/R/msqrob.R index 2247aa7..e48fc2b 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -259,12 +259,12 @@ msqrobLmer <- function(y, formula, data, rowdata = NULL, - tol =1e-6, + tol = 1e-6, robust = TRUE, ridge = FALSE, maxitRob = 1, doQR = TRUE, - featureGroups=NULL, + featureGroups = NULL, lmerArgs = list(control = lmerControl(calc.derivs = FALSE))){ #Get the featureGroups variable @@ -287,12 +287,12 @@ msqrobLmer <- function(y, if(is.null(rowdata)){ models <- bplapply(y, FUN = .ridge_msqrobLmer, - "formula" = formula, - "coldata" = data, - "doQR" = doQR, - "robust"=robust, - "maxitRob" = maxitRob, - "tol" =tol) + formula = formula, + coldata = data, + doQR = doQR, + robust = robust, + maxitRob = maxitRob, + tol = tol) } else{ models <- bpmapply(FUN = .ridge_msqrobLmer, y, rowdata, diff --git a/data/cptac.rda b/data/cptac.rda new file mode 100644 index 0000000..29c7c7f Binary files /dev/null and b/data/cptac.rda differ diff --git a/inst/scripts/make_data-cptac.R b/inst/scripts/make_data-cptac.R new file mode 100644 index 0000000..5282b93 --- /dev/null +++ b/inst/scripts/make_data-cptac.R @@ -0,0 +1,37 @@ + + +library("MsDataHub") +peptideFile <- cptac_peptides.txt() +peptides <- read.delim(peptideFile) + +## Subset the data to make a small dataset +allProteins <- unique(peptides$Proteins) +upsProteins <- grep("ups", allProteins, value = TRUE) +yeastProteins <- grep("ups", allProteins, value = TRUE, invert = TRUE) +set.seed(1234) +keepProteins <- c(sample(upsProteins, 10), sample(yeastProteins, 50)) +peptides <- peptides[peptides$Proteins %in% keepProteins, ] + +## Create sample annotations +quantCols <- grep("Intensity[.]", names(peptides), value = TRUE) +coldata <- data.frame(quantCols = quantCols) +coldata$lab <- rep(rep(paste0("lab", 1:3), each = 3), 5) +coldata$condition <- gsub("Intensity..(.)_.*", "\\1", quantCols) +concentrations <- c(A = 0.25, B = 0.74, C = 2.22, D = 6.67, E = 20) +coldata$concentration <- concentrations[coldata$condition] + +## Remove unnecessary feature annotations +keepAnnot <- c( + "Sequence", "Proteins", "Charges", "PEP", "Reverse", + "Potential.contaminant", quantCols +) +peptides <- peptides[, keepAnnot] + +## Create the QFeatures object +cptac <- readQFeatures( + peptides, coldata, name = "peptides", fnames = "Sequence" +) + +## Store data +save(cptac, file = file.path("data/cptac.rda"), + compress = "xz", compression_level = 9) diff --git a/tests/testthat/test_msqrob.R b/tests/testthat/test_msqrob.R new file mode 100644 index 0000000..227b577 --- /dev/null +++ b/tests/testthat/test_msqrob.R @@ -0,0 +1,163 @@ +data(cptac) +cptac <- zeroIsNA(cptac, "peptides") +cptac <- logTransform(cptac, "peptides", "peptides_log") +cptac <- normalize(cptac, "peptides_log", "peptides_norm", "center.median") + +test_that("msqrobLmer", { + y <- assay(cptac[["peptides_norm"]]) + coldata <- colData(cptac) + rowdata <- rowData(cptac)[["peptides"]] + ## Default usage + formula <- ~ condition + (1|lab) + out <- msqrobLmer(y, formula, coldata) + expect_identical(length(out), nrow(y)) + expect_true(all(sapply(out, class) == "StatModel")) + expect_identical( + as.list(table(sapply(out, function(x) x@type))), + list(fitError = 104L, lmer = 252L) + ) + expect_equal( + out[[1]]@params, + list( + coefficients = c( + `(Intercept)` = -1.43239555985396, conditionB = 0.303816995791296, conditionC = 0.400933593097784, conditionD = -0.716021732465039, `(Intercept)lablab2` = 0.277916037800659, `(Intercept)lablab3` = -0.277916037800659), + vcovUnscaled = new( + "dsCMatrix", + i = c(0L, 0L, 1L, 0L, 1L, 2L, 0L, 1L, 2L, 3L, 0L, 1L, 2L, 3L, 4L, 0L, 1L, 2L, 3L, 4L, 5L), + p = c(0L, 1L, 3L, 6L, 10L, 15L, 21L), + Dim = c(6L, 6L), + Dimnames = list( + c("(Intercept)", "conditionB", "conditionC", "conditionD", "(Intercept)lablab2", "(Intercept)lablab3"), + c("(Intercept)", "conditionB", "conditionC", "conditionD", "(Intercept)lablab2", "(Intercept)lablab3") + ), + x = c(0.972657631232327, -0.668226165745811, 1.16822616574581, -0.668226165745811, 0.668226165745811, 1.28592310873781, -0.668226165745811, 0.668226165745811, 0.651459442390789, 1.24184307090335, -0.304431465486516, 6.48495793590282e-17, -0.117696942992002, 0.0167667233550216, 0.422128408478519, -0.304431465486516, -1.3687424862664e-16, 0.117696942992002, -0.0167667233550217, 0.186734522494514, 0.422128408478519), + uplo = "U", factors = list() + ), + sigma = 0.643133756500329, + df.residual = 4.38661228982676, + w = c(0.748249658021019, 0.748249658021019, 1, 1, 1, 1, 0.353602245028265, 0.397011079589108, 1) + ) + ) + expect_equal( ## test spike-in coefs + out[[grep("ups", rowdata$Proteins)[[1]]]]@params$coefficients[["conditionB"]], + log2(0.74) - log2(0.25), + tolerance = 0.3 + ) + expect_equal( ## test constant coefs + out[[grep("YEAST", rowdata$Proteins)[[1]]]]@params$coefficients[["conditionB"]], + 0, + tolerance = 0.5 + ) + ## Test with rowData annotations + formula <- ~ condition + (1|lab) + (1 | Sequence) + featureGroups <- rowdata$Proteins + out <- msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups + ) + expect_identical( + as.list(table(sapply(out, function(x) x@type))), + list(fitError = 14L, lmer = 46L) + ) + expect_equal( ## test spike-in coefs + out[[grep("ups", names(out))[[1]]]]@params$coefficients[["conditionB"]], + log2(0.74) - log2(0.25), + tolerance = 0.3 + ) + expect_equal( ## test constant coefs + out[[grep("YEAST", names(out))[[1]]]]@params$coefficients[["conditionB"]], + 0, + tolerance = 0.3 + ) + ## Disable robust + out <- msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups, robust = FALSE + ) + expect_identical( + as.list(table(sapply(out, function(x) x@type))), + list(fitError = 14L, lmer = 46L) + ) + expect_equal( ## test spike-in coefs + out[[grep("ups", names(out))[[1]]]]@params$coefficients[["conditionB"]], + log2(0.74) - log2(0.25), + tolerance = 0.3 + ) + expect_equal( ## test constant coefs + out[[grep("YEAST", names(out))[[1]]]]@params$coefficients[["conditionB"]], + 0, + tolerance = 0.3 + ) + expect_true(all(sapply(out, function(x) all(x@params$w == 1)))) + ## Enable ridge + out <- msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups, ridge = TRUE + ) + expect_identical( + as.list(table(sapply(out, function(x) x@type))), + list(fitError = 14L, lmer = 46L) + ) + expect_equal( ## test spike-in coefs + out[[grep("ups", names(out))[[1]]]]@params$coefficients[["ridgeconditionB"]], + log2(0.74) - log2(0.25), + tolerance = 0.3 + ) + expect_equal( ## test constant coefs + out[[grep("YEAST", names(out))[[1]]]]@params$coefficients[["ridgeconditionB"]], + 0, + tolerance = 0.1 ## decrease tolerance for non DA with ridge + ) + ## Disable QR + ## no impact if no ridge + expect_identical( + msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups + ), + msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups, doQR = FALSE + ) + ) + ## disable QR when ridge + out <- msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups, ridge = TRUE, doQR = FALSE + ) + expect_identical( + as.list(table(sapply(out, function(x) x@type))), + list(fitError = 14L, lmer = 46L) + ) + expect_equal( ## test spike-in coefs + out[[grep("ups", names(out))[[1]]]]@params$coefficients[["ridgeconditionB"]], + log2(0.74) - log2(0.25), + tolerance = 0.3 + ) + expect_equal( ## test constant coefs + out[[grep("YEAST", names(out))[[1]]]]@params$coefficients[["ridgeconditionB"]], + 0, + tolerance = 0.1 ## decrease tolerance for non DA with ridge + ) + ## Test M-estimation parameters + out <- msqrobLmer( + y = y, formula = formula, data = coldata, rowdata = rowdata, + featureGroups = featureGroups, maxitRob = 10, tol = 1E-2 + ) + expect_identical( + as.list(table(sapply(out, function(x) x@type))), + list(fitError = 14L, lmer = 46L) + ) + expect_equal( ## test spike-in coefs + out[[grep("ups", names(out))[[1]]]]@params$coefficients[["conditionB"]], + log2(0.74) - log2(0.25), + tolerance = 0.3 + ) + expect_equal( ## test constant coefs + out[[grep("YEAST", names(out))[[1]]]]@params$coefficients[["conditionB"]], + 0, + tolerance = 0.3 + ) + ## Test lmer parameters + stop("lmerArgs are not used...") +})