Hi everybody, I'm using EdgeR to find differential expressed genes induced by diet. I have two group (5 mice in standard diet (SD) and 5 mice in Diet1 (KD)). In both groups I have males and females, respectively 2F+3M and 3F+2F

Analysis is carried out using the following code

library(edgeR)
library(ggplot2)

#make DGElist
DG <- DGEList(matrix, group= group)
DG$samples$name= x$description
DG$samples$sex= x$Sex


#filtering
keep <- filterByExpr(d) 
d <- d[keep,]
d$samples$lib.size <- colSums(d$counts)


#TMM normalization
d <- calcNormFactors(d, method = "TMM")
logCPM <- cpm(d, log=TRUE , prior.count=1)


#MDS plot (Fig1)
plotMDS(logCPM, labels=as.factor(d$samples$Experiment), col=rep(1:4, each=5))

As possible to see in MDS plot (Fig1), clustering is sex driven, indeed I lost observed clusterization when I run

logCPM_batch <- removeBatchEffect(logCPM, batch=DG$samples$sex)
plotMDS(logCPM_batch, labels=as.factor(d$samples$name), col=rep(1:4, each=5))

See Fig2

To investigate the impact of diet (removing sex effect) I have designed the matrix as follow and proceed with following code

design.mat <- model.matrix(~d$samples$sex+ d$samples$group)
d <- estimateDisp(d, design.mat) 
fit <- glmFit(d, design.mat)
lrt12 <- glmLRT(fit, coef = 2)
lrt12$table$FDR= p.adjust(lrt12$table$PValue, method = "BH")

However, looking the pvalue histogram I have a hill-shaped p-value distributions (see Fig3)

x= lrt12$table
ggplot(data = x, aes(x=PValue)) + geom_histogram()+ theme_bw(base_size = 22)+ geom_vline(xintercept = 0.05, col="red")

What I'm doing wrong?? What I can do to be sure about results?? please can you help me?

thanks a lot in advance

Giuseppe

You have tested for sex differences (using coef=2) instead of testing for group differences.