Commit bdb15ba2 authored by Leon-Charles Tranchevent's avatar Leon-Charles Tranchevent
Browse files

Limma anlaysis: integration of co-factor and block key usage (for datasets with paired samples).

parent 231e050f
...@@ -42,7 +42,7 @@ preprocess_data_agilent_limma <- function(input_data_dir, output_data_file, ...@@ -42,7 +42,7 @@ preprocess_data_agilent_limma <- function(input_data_dir, output_data_file,
batch_data <- log2(batch$E) batch_data <- log2(batch$E)
# We change the probe ids (1 to 45015) to be the probe names instead (A_XX_PXXXX). # We change the probe ids (1 to 45015) to be the probe names instead (A_XX_PXXXX).
rownames(batch_data) <- (batch$genes)$ProbeName rownames(batch_data) <- (batch$genes)$ProbeName # nolint
remove(raw_data_input_dir, batch) remove(raw_data_input_dir, batch)
# We run the LIMMA pre-processing method on the data. # We run the LIMMA pre-processing method on the data.
......
...@@ -25,25 +25,25 @@ run_limma <- function(pheno_data, eset, limma_parameters, ...@@ -25,25 +25,25 @@ run_limma <- function(pheno_data, eset, limma_parameters,
clinical_factor <- limma_parameters[[1]] clinical_factor <- limma_parameters[[1]]
coefficients <- limma_parameters[[2]] coefficients <- limma_parameters[[2]]
has_block <- !is.null(cofactor_name) has_block <- !is.null(cofactor_name)
is_factorial_limma <- length(coefficients) == 3
# We create the design matrix (regular case - no block). # We create the design matrix (regular case - no block).
clinical <- factor(pheno_data[[clinical_factor]]) clinical <- factor(pheno_data[[clinical_factor]])
design <- stats::model.matrix(~ 0 + clinical) design <- stats::model.matrix(~ 0 + clinical)
colnames(design) <- levels(clinical) colnames(design) <- levels(clinical)
# If necessary we add a block based on a co-factor. # If necessary we add a block based on a co-factor.
# For instance, if we have paired samples. # For instance, if we have paired samples.
if (has_block) { if (has_block) {
clinical_cofactor <- factor(pheno_data[[cofactor_name]]) clinical_cofactor <- factor(pheno_data[[cofactor_name]])
design <- stats::model.matrix(~ 0 + clinical + clinical_cofactor) design <- stats::model.matrix(~ 0 + clinical + clinical_cofactor)
colnames(design) <- c(levels(clinical), clinical_cofactor) colnames(design) <- c(levels(clinical), cofactor_name)
} }
# We build the contrast (default only one). # We build the contrast (default only one).
contrast <- eval(substitute(limma::makeContrasts(coeff1 = c1, levels = design), contrast <- eval(substitute(limma::makeContrasts(coeff1 = c1, levels = design),
list(c1 = coefficients[1]))) list(c1 = coefficients[1])))
# For larger lists of coefficients, we add them all. # For larger lists of coefficients, we add them all.
if (is_factorial_limma) { if (limma_parameters$factorial) {
contrast <- eval(substitute(limma::makeContrasts(coeff1 = c1, contrast <- eval(substitute(limma::makeContrasts(coeff1 = c1,
coeff2 = c2, coeff2 = c2,
coeff3 = c3, coeff3 = c3,
...@@ -60,7 +60,7 @@ run_limma <- function(pheno_data, eset, limma_parameters, ...@@ -60,7 +60,7 @@ run_limma <- function(pheno_data, eset, limma_parameters,
# We create a block if necessary. # We create a block if necessary.
if (has_block) { if (has_block) {
corfit <- limma::duplicateCorrelation(eset, design, block = pheno_data[[block_key]]) corfit <- limma::duplicateCorrelation(eset, design, block = pheno_data[[block_key]])
lmfit_gender <- limma::lmFit(eset, design, block = pheno_data[[block_key]], lmfit <- limma::lmFit(eset, design, block = pheno_data[[block_key]],
correlation = corfit$consensus) correlation = corfit$consensus)
} }
fit <- limma::contrasts.fit(lmfit, contrast) fit <- limma::contrasts.fit(lmfit, contrast)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment