Hello!
I'm trying to analyse an RNA-seq count data matrix with limma. From an input table with about 15000 entries, limma gets 1600 differentially expressed genes. That is OK, but these genes have an adjusted P-value about ten times higher than the required 0.05 threshold. Am I doing something wrong? Should I analyse RNA-seq data with limma in a different way? I'm using the following code:
raw_filter is the count data data frame
replicates <- 3
groups <- c(rep("A", replicates), rep("B", replicates))
DGE_List <- DGEList(counts=raw_filter, group=groups)
DGE_List <- calcNormFactors(DGE_List)
normalized_counts <- as.matrix(DGE_List)
design <- model.matrix(~0+ groups)
log2_cpm <- voom(DGE_List, design, plot=TRUE)
fit <- lmFit(log2_cpm, design)
cont_matrix <- makeContrasts(c("groupsB-groupsA"), levels=design)
fit2 <- contrasts.fit(fit, cont_matrix)
fit2 <- eBayes(fit2)
all_limma <- decideTests(fit2, p.value = 0.05, lfc = lfc)
all_limma <- topTable(fit2, adjust.method="BH", number=nrow(fit2))
all_limma <- all_limma[order(rownames(all_limma)),]
expres_diff <- subset(all_limma, all_limma$P.Value < 0.05)
Thank you in advance for your help
I'm trying to analyse an RNA-seq count data matrix with limma. From an input table with about 15000 entries, limma gets 1600 differentially expressed genes. That is OK, but these genes have an adjusted P-value about ten times higher than the required 0.05 threshold. Am I doing something wrong? Should I analyse RNA-seq data with limma in a different way? I'm using the following code:
raw_filter is the count data data frame
replicates <- 3
groups <- c(rep("A", replicates), rep("B", replicates))
DGE_List <- DGEList(counts=raw_filter, group=groups)
DGE_List <- calcNormFactors(DGE_List)
normalized_counts <- as.matrix(DGE_List)
design <- model.matrix(~0+ groups)
log2_cpm <- voom(DGE_List, design, plot=TRUE)
fit <- lmFit(log2_cpm, design)
cont_matrix <- makeContrasts(c("groupsB-groupsA"), levels=design)
fit2 <- contrasts.fit(fit, cont_matrix)
fit2 <- eBayes(fit2)
all_limma <- decideTests(fit2, p.value = 0.05, lfc = lfc)
all_limma <- topTable(fit2, adjust.method="BH", number=nrow(fit2))
all_limma <- all_limma[order(rownames(all_limma)),]
expres_diff <- subset(all_limma, all_limma$P.Value < 0.05)
Thank you in advance for your help