**10**wrote:

Hi everyone,

I'm using edgeR for differential expression genes analysis. I make 4 groups that get differential expression genes each and then compare the number of genes. But there are no common genes between 4 groups. the result is not ideal. Please, Could you help me to check on the R codes?

I have 4 groups(18H, 18L, 30H, 30L).

each group consists of 10 samples(18H_1,..., 18H_10, 18L_1, ..., 18L_10, 30H_1,..., 30H_10, 30L_1, ..., 30L_10).

I made 4 groups to get genes that differentially express between two groups(18H18L, 18H30H, 18L30L, 30H30L).

and I made venn diagram to compare the number of genes between 4 groups.

I made several scripts. But there are no common genes at all scripts.

```
#1
genecounts <- read.table("18H18L.txt", row.names = 1)
colnames(genecounts) <- c(paste("18H",1:10,sep = "_"), paste("18L",1:10,sep="_"))
condition <-c(rep("H", 10), rep("L", 10))
dge <- DGEList(counts = genecounts, group = condition)
countsPerMillion <- cpm(dge)
countCheck <- countsPerMillion > 1
keep <- which(rowSums(countCheck)>=2)
dge <- dge[keep,]
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition + 0, data = dge$samples)
colnames(design) = c("High", "Low")
disp <- estimateGLMCommonDisp(dge, design)
disp <- estimateGLMTrendedDisp(disp, design)
disp <- estimateGLMTagwiseDisp(disp, design)
fit <- glmFit(disp, design)
lrt <- glmLRT(fit)
res <- as.data.frame(lrt$table)
res <- cbind(res, FDR = p.adjust(res$PValue, method = "BH"))
fc1 <- res[res[,1] >= 1 | res[,1] <= -1,]
fc1p <- fc1[fc1$PValue <= 0.05,]
fc1q <- fc1[fc1$FDR <= 0.05,]
#2
genecounts <- read.table("18H18L.txt", row.names = 1)
colnames(genecounts) <- c(paste("18H",1:10,sep = "_"),paste("18L",1:10,sep="_"))
condition <-c(rep("H", 10), rep("L", 10))
dge <- DGEList(counts = genecounts, group = condition)
countsPerMillion <- cpm(dge)
countCheck <- countsPerMillion > 0.5
keep <- which(rowSums(countCheck)>=2)
dge <- dge[keep,]
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition + 0, data = dge$samples)
colnames(design) = c("High", "Low")
disp <- estimateGLMCommonDisp(dge, design)
disp <- estimateGLMTrendedDisp(disp, design)
disp <- estimateGLMTagwiseDisp(disp, design)
fit <- glmFit(disp, design)
highlow <- makeContrasts(High-Low,levels = design)
lrt <- glmLRT(fit, contrast =highlow)
res <- as.data.frame(lrt$table)
res <- cbind(res, FDR = p.adjust(res$PValue, method = "BH"))
fc1 <- res[res[,1] >= 1 | res[,1] <= -1,]
fc1p <- fc1[fc1$PValue <= 0.05,]
fc1q <- fc1[fc1$FDR <= 0.05,]
write.table(fc1p,"fc1p_18H18L_1.txt", quote = F, sep = "\t", row.names = T, col.names = T)
#3
genecounts <- read.table("18H18L.txt", row.names = 1)
colnames(genecounts) <- c(paste("18H",1:10,sep = "_"),paste("18L",1:10,sep="_"))
condition <-c(rep("H", 10), rep("L", 10))
dge <- DGEList(counts = genecounts, group = condition)
countsPerMillion <- cpm(dge)
countCheck <- countsPerMillion > 0.5
keep <- which(rowSums(countCheck)>=2)
dge <- dge[keep,]
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition + 0, data = dge$samples)
colnames(design) = c("High", "Low")
disp <- estimateGLMCommonDisp(dge, design)
disp <- estimateGLMTrendedDisp(disp, design)
disp <- estimateGLMTagwiseDisp(disp, design)
fit <- glmFit(disp, design)
lrt <- glmLRT(fit)
res <- as.data.frame(lrt$table)
res <- cbind(res, FDR = p.adjust(res$PValue, method = "BH"))
fc1 <- res[res[,1] >= 1 | res[,1] <= -1,]
fc1p <- fc1[fc1$PValue <= 0.05,]
fc1q <- fc1[fc1$FDR <= 0.05,]
#2-1
library(gplots)
data1 <- read.table("fc1p_18H18L_1.txt", header = T)
data2 <- read.table("fc1p_18H30H_1.txt", header = T)
data3 <- read.table("fc1p_18L30L_1.txt", header = T)
data4 <- read.table("fc1p_30H30L_1.txt", header = T)
test1 <- data1[,1]
test2 <- data2[,1]
test3 <- data3[,1]
test4 <- data4[,1]
venn(list(A = test1, B = test2, C = test3, D = test4))
A group = 105, B group = 1778, C group = 1665, D group = 47
```

There are no common genes between groups.

If there is any mistake or question, please let me know. I look forward to hearing from you.

Thanks!

**70k**• written 18 months ago by GuriGuri1565 •

**10**

Do all your groups belong to the same experiment? Why are you fitting several models with subsets of the data? Fit one model with all data and use

`makeContrasts`

to perform differential analyses for the contrasts of interest.Then, to test for differences in H18 vs L18:

In addition, from what I understood from

`2-1`

, sections`1`

,`2`

and`3`

are wrong or incomplete: sections`2`

and`3`

are (as far as I can tell) identical, and there should be a section`4`

.26kI’m so sorry for late response. Thank you~

10isn't it possible that there simply aren't any overlapping genes?

4.6k