diff --git a/R/anova.R b/R/anova.R index 7a284595f6073fcc999749fea92eecffbc86bc63..1917e5cf00e4710e074deb98974f236548e7b56d 100644 --- a/R/anova.R +++ b/R/anova.R @@ -116,7 +116,36 @@ batch_effect <- function(data) { return(data) } -#' build rm model between drugs accounting for batch effect +#' build lm model between drugs accounting for batch effect +#' +#' @param data a data.frame +#' @param formula (default: "ratio ~ drug + batch") the formula of the model +#' @param lower (default: TRUE) tests if "y" (the ratio) is lower than in controls +#' @param outdir set the outdir directory to a path (default extracted from fcs) +#' @return a data.frame +#' @examples +#' \dontrun{ +#' data <- anova_lm(data) +#' } +#' @importFrom grDevices dev.off pdf +#' @importFrom stats as.formula quantile +#' @export anova_lm +anova_lm <- function(data, formula = "ratio ~ drug + batch", lower = TRUE, + outdir) { + variable_name <- gsub("(.*) ~.*", "\\1", formula) + model <- lm(stats::as.formula(formula), + data = data) + model_anova <- compute_pval(model, lower = lower) + if (missing(outdir)) { + outdir <- mk_outdir(data, "/test") + } + data <- export_lm_results(data, model_anova) + save(data, model, file = paste0(outdir, "anova_rlm.Rdata")) + export_drug_table(data, model_anova, outdir) + return(data) +} + +#' build rlm model between drugs accounting for batch effect #' #' @param data a data.frame #' @param formula (default: "ratio ~ drug + batch") the formula of the model @@ -181,6 +210,27 @@ export_rlm_results <- function(data, model_anova) { return(data) } +#' @importFrom utils write.csv +export_lm_results <- function(data, model_anova) { + data$coef <- NA + data$coef_std <- NA + data$pval <- NA + data$tval <- NA + for (drug in levels(data$drug)) { + if (!(drug %in% "None")) { + data$coef[data$drug %in% drug] <- model_anova$Estimate[grepl(drug, + rownames(model_anova))] + data$coef_std[data$drug %in% drug] <- model_anova[grepl(drug, + rownames(model_anova)), 2] + data$tval[data$drug %in% drug] <- model_anova$t.value[grepl(drug, + rownames(model_anova))] + data$pval[data$drug %in% drug] <- model_anova$pval[grepl(drug, + rownames(model_anova))] + } + } + return(data) +} + export_drug_table <- function(data, model_anova, outdir, channels = c("Y1.A", "B1.A")) { drug_table <- model_anova[grepl("drug", rownames(model_anova)), ] diff --git a/R/project.R b/R/project.R index 1955d9a0d4d89feb42907b6e84d850d8b09a559b..e88168befd2ce8ec0d6ec9851af7f7ea70e052d9 100644 --- a/R/project.R +++ b/R/project.R @@ -37,13 +37,14 @@ set_analysis <- function(data_path = "data/", meta = F) { #' perform all the analysis for a set of data sets #' #' @param data_path path to folder containing the data sets +#' @param rlm_model (default: TRUE) should rlm model be use or lm ? #' @return TRUE if everythings ran correclty #' @examples #' \dontrun{ #' analysis("data/set_test") #' } #' @export analysis -analysis <- function(data_path = "data/") { +analysis <- function(data_path = "data/", rlm_model = TRUE) { if (base::file.info(data_path)$isdir) { set_folders <- list.dirs(data_path, full.names = F)[-1] } else { @@ -88,8 +89,13 @@ analysis <- function(data_path = "data/") { rm(set_data) } data$set <- as.factor(data$set) - data <- anova_rlm(data, formula = "ratio ~ drug + batch + set", - outdir = outdir_rlm) + if (rlm_model) { + data <- anova_rlm(data, formula = "ratio ~ drug + batch + set", + outdir = outdir_rlm) + } else { + data <- anova_lm(data, formula = "ratio ~ drug + batch + set", + outdir = outdir_rlm) + } print(summary(data)) for (folder in set_folders) { message(paste0("plotting for ", folder)) diff --git a/criblejurkat_final.pbs b/criblejurkat_final.pbs index 6208edec0b6f093ced47770186ae553b203275b2..6bbf09b474603278039220863aaf8df854857189 100644 --- a/criblejurkat_final.pbs +++ b/criblejurkat_final.pbs @@ -24,6 +24,6 @@ source /usr/share/lmod/lmod/init/bash ml R/3.4.3 export R_LIBS="/Xnfs/lbmcdb/common/R/x86_64-pc-linux-gnu-library/3.4/" -R -e '.libPaths(); library("criblejurkat"); analysis("data/final"); traceback()' +R -e '.libPaths(); library("criblejurkat"); analysis("data/final", rlm_model = F); traceback()' diff --git a/criblejurkat_full.pbs b/criblejurkat_full.pbs index f8df3418705d1516bffcdcad22d7c0850c990d02..4fcaf3ce45122b635a8aaa643a606a8116b5dab8 100644 --- a/criblejurkat_full.pbs +++ b/criblejurkat_full.pbs @@ -24,6 +24,6 @@ source /usr/share/lmod/lmod/init/bash ml R/3.4.3 export R_LIBS="/Xnfs/lbmcdb/common/R/x86_64-pc-linux-gnu-library/3.4/" -R -e '.libPaths(); library("criblejurkat"); analysis("data/full"); traceback()' +R -e '.libPaths(); library("criblejurkat"); analysis("data/full", rlm_model = F); traceback()'