Hi, I have posted a similar question few months ago applying voom + limma to a block factor design in RNA-seq experiment, but now I have a doubt if the no differences observed between cases and controls are real or caused by my design or by my misuse of block
and duplicateCorrelation
My metadata looks like this one:
sample_ID location sex polyp_type age status files 1 INTP_103 proximal 1 1 31 ctrl INTP_1031 2 INTP_103 distal 1 1 31 ctrl INTP_1032 3 INTP_103 rectum 1 1 31 ctrl INTP_1033 4 INTP_105 proximal 2 3 66 disease INTP_1051 5 INTP_105 distal 2 3 66 disease INTP_1052 6 INTP_105 rectum 2 3 66 disease INTP_1053 7 INTP_106 proximal 1 2 64 disease INTP_1061 8 INTP_106 distal 1 2 64 disease INTP_1062 9 INTP_106 rectum 1 2 64 disease INTP_1063 10 INTP_111 proximal 2 1 49 ctrl INTP_1111
Basically I have samples from controls and individuals with a disease, each of the individuals (both control and diseases) have 3 samples from the same tissue but different location. I want to test differences between locations, between cases and control, and between interactions (disease status and location), here is my code
y <- DGEList(txi_longScaled$counts) # the counts come from a kallisto experiment y <- calcNormFactors(y)
age <- metadata$age sex <- factor(metadata$sex)
Treat <- factor(paste(metadata$status, metadata$location, sep="."))
design <- model.matrix(~0+Treat + sex + age) design2 <- model.matrix(~1+status + sex + age) colnames(design)[seq_len(nlevels(Treat))] <- levels(Treat) colnames(design2) <- c("Ctrl", "Disease", "sex", "age") # Apply voom v <- voom(y, design) v2 <- voom(y, design2) #Then we estimate the correlation between measurements made on the same subject corfit <- duplicateCorrelation(v, design, block=metadata$sample_ID) corfit2 <- duplicateCorrelation(v2, design2, block=metadata$sample_ID) # Apply voom taking in account block effect ### remove correlation from the models v <- voom(y, design, block = metadata$sample_ID, correlation = corfit$consensus) v2 <- voom(y, design2, block = metadata$sample_ID, correlation = corfit2$consensus) vfit <- lmFit(v, design, block = metadata$sample_ID, correlation = corfit$consensus) vfit2 <- lmFit(v2, design2, block = metadata$sample_ID, correlation = corfit2$consensus) ## Specifigy contrasts cm <- makeContrasts( CtrlvsCasesProximal = ctrl.proximal - disease.proximal, CtrlvsCasesDistal = ctrl.distal - disease.distal, CtrlvsCasesRectum = ctrl.rectum - disease.rectum, Ctrl.proximalvsCtrl.rectum = ctrl.proximal - ctrl.rectum, Ctrl.distalvsCtrl.proximal = ctrl.distal - ctrl.proximal, Ctrl.distalvsCtrl.rectum = ctrl.distal - ctrl.rectum, Dis.proximalvsDis.distal = disease.proximal - disease.distal, Dis.proximalvsDis.rectum = disease.proximal - disease.rectum, Dis.distalvsDis.rectum = disease.distal - disease.rectum, levels=design) ## Apply limma eBayes to voom objects vfit <- contrasts.fit(vfit, cm) efit <- eBayes(vfit) efit2 <- eBayes(vfit2, robust = TRUE)
Are these designs correct? Or should I use just one single contrast and then extract the coefficients, I mean, something like this:?
design3 <- model.matrix(~0+location + status + sex + age + metadata$sample_ID)
v3 <- voom(y, design3)
vfit3 <- lmFit(v3, design3)
efit3 <- eBayes(efit3, robust = TRUE)
What is
Treat
? Don't rely on information defined in previous posts.Post updated @Aaron Lun