Microarray data analysis with and without contrast

Question

[1769mkc]

I have 6 samples of control and test condition normalized microarray intensities which Im using a starting point of my analysis.

This is my basic analysis code

library(readxl)
library(dplyr)
library(tidyverse)
library(limma)  # t

# Create a data frame with the sample names and group names
metdata <- data.frame(
  SampleName = c("C14WEx_MBr_1", "C14WEx_MBr_2", "C14WEx_MBr_4", "T814WEx_MBr_1", "T814WEx_MBr_2", "T814WEx_MBr_4"),
  GroupName = c("Control", "Control", "Control", "TRPM8KO", "TRPM8KO", "TRPM8KO")
)

metdata$GroupName <- as.factor(metdata$GroupName)
mrna_df <- read_excel("data.xlsx", 
                      skip = 24)
mrna_norm <- mrna_df %>% select(Transcript_ID,c(8:13))

mrna_norm_df <- mrna_norm %>% column_to_rownames('Transcript_ID')

# Extract the SampleName column from the metadata dataframe
sample_names <- metdata$SampleName

# Rename the columns of mrna_norm_df with the SampleName values
colnames(mrna_norm_df) <- sample_names

# Print the updated dataframe
head(mrna_norm_df)


class(metdata$GroupName)
design_matrix <- model.matrix(~GroupName, metdata) ## with intercept 
design_matrix
levels(metdata$GroupName)



#https://alexslemonade.github.io/refinebio-examples/02-microarray/differential-expression_microarray_01_2-groups.html
fit <- lmFit(mrna_norm_df, design_matrix)
fit <- eBayes(fit)
tempOutput <- topTable(fit, coef="GroupNameTRPM8KO",n=Inf,
                       adjust.method = 'BH')


length(which(tempOutput$adj.P.Val < 0.05))
#5
#####################################################

Running the fit with contrast

# with contrast# 
contrasts <- makeContrasts(Control - TRPM8KO, levels = metdata$GroupName)
#contrasts <- relevel(contrasts, ref = "Control")
contrasts <- makeContrasts(Control - TRPM8KO, levels = metdata$GroupName, ref = "Control")

contrasts
fit <- lmFit(mrna_norm_df, design = design_matrix)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)

tempOutput1 <- topTable(fit2, coef="Control - TRPM8KO",n=Inf,
                        adjust.method = 'BH')

length(which(tempOutput1$adj.P.Val < 0.05))

#19814

So when I run limma without contrast as i get only 5 genes with the 0.05 threshold but when Im running with contrast with same design then im getting 19814. Am i doing something wrong in both the case which is making this sharp change in the list of genes?

Anysuggestion or help would be really helpful

Answer 1

[Gordon Smyth]

The second approach is incorrect. makeContrasts is used to construct contrasts from the coefficients defined by the design matrix. You are trying to compare Control to TRPM8KO but the design matrix doesn't have coefficients by either of those names, so the contrast doesn't match the design matrix.

The use of design matrices with and without contrasts is discussed in detail in the limma User's Guide.

[1769mkc]

so the first one what i have done so far is correct?

" You are trying to compare Control to TRPM8KO but the design matrix doesn't have coefficients by either of those names," this i tried to fix after going through different answers but i was not able to fix that issue. What is it im doing it incorrectly ?

[ATpoint]

so the first one what i have done so far is correct?

Yes, that tests for differences between the two groups.

[1769mkc]

why can;t i see your answer which i do see in my inbox but not here? regarding the contrast

[ATpoint]

I removed it because my wording was not good and I was busy with other things so I did not have time to properly rephrase it. Just stick with coef=2 with the intercept design, that is correct for DE.

[1769mkc]

i tried what you suggested now i get same result as my first method as well as contrast