Dear all,

I have 6 samples which are three High marbling Hanwoo(Korean cattle) and three

Low marbling and each sample have 4 different stages(6mon, 12mon, 24mon, 30mon). 

I'd like to find DEG from each month and intersected genes which are located in center of Venn-diagram 

of all months.

Initially, I have analyzed with cufflinks and cuffdiff package to detect DEGs and now, try to do with limma.

Here is my question

I find more than 238 DEGs with Cuffdiff(24mon, q-value < 0.05), However, when I used limma package,

it gave me only 1 DEG(adj.P.Val < 0.05). I know these two tools have different protocol, and way to find DEGs

such as statistical test, normalization, expression value etc.

I'm just wondering there are some wrong codes in my scripts, if not, why only 1 gene came from.

My R scripts are below, any help should be appreciated. thanks. 

library(edgeR)
library(limma)

data <- read.table("24mon_H_L.txt", sep = "\t",header = T, row.names = 1)
data <- data[,-1]
head(data)
dim(data)

#filtering
isexpr <- rowSums(cpm(data) > 1) >= 3
table(isexpr)
data <- data[isexpr,]
dim(data)
genes <- rownames(data)
colnames(data)

#makes model matrix
group <- factor(c(0,0,0,1,1,1))
group
design <- model.matrix(~0 + group)
design
colnames(design) <- c("High", "Low")

# calculate normalization factors
nf <- calcNormFactors(data, method = "TMM")
nf

# normalize the read count with voom
y <- voom(data, design, lib.size = colSums(data) * nf)
y

# extract normalized read counts
counts.voom <- y$E
head(counts.voom)
dim(counts.voom)

# save normalized expression
write.table(counts.voom, "counts.voom.txt", row.names = T, quote = F, sep = "\t")

# linear model
fit <- lmFit(y, design)
fit

# contrast matrix
cont.matrix <- makeContrasts(High-Low, levels = design)
cont.matrix

fit <- contrasts.fit(fit, cont.matrix)
fit

# Bayes
fit <- eBayes(fit)
dim(fit)
fit

# topTable
colnames(fit$coefficients)
mycoef = "High - Low"
topTable(fit, coef = mycoef)

limma.res <- topTable(fit, coef = mycoef, n = dim(fit)[1])
deg <- limma.res[which(limma.res$adj.P.Val <= 0.05),]
deg

I don't see anything wrong with your code, so it's probably not that. Instead it is probably due to the fact that cuffdiff doesn't do a very good job of controlling type I errors, whereas limma-voom does. In addition (and here I am on thin ice because I have only used cufflinks et al once, and vowed never to repeat that experience), the main goal for the tophat-cufflinks-cuffdiff pipeline is to detect differential transcript abundance, not differential gene expression. Those are not the same thing.

Presumably you can tell cuffdiff you want DEG, and maybe that's what you did. Anyway, even Lior Pachter has stopped using cufflinks/cuffdiff and is now recommending that people use kallisto and sleuth for such things, so it's probably not a super good idea to rely on software that isn't even being used by the guy whose lab developed it.