Comparison of different statistical tests on different amount of biological replicates I am testing differentially expressed genes using different statistical tests on different amount of biological replicates (Wilcoxon test, t-test and negative binomial test). Figure below is what I got.
.
My questions are:
1: Why the No. of DEG decreases when replicates reach to 500?
2: Why Negative binomial test detect less DEG than other two tests?
Does this trend make sense?
 A: I'm kind of guessing here:

1: Why the No. of DEG decreases when replicates reach to 500?

The difference between 100 and 500 is not much and it may be explained by random sampling. If you try different samples for N=100 and N=500, do you still see the same picture?

2: Why Negative binomial test detect less DEG than other two tests?

Negative binomial allows extra variation in addition to random sampling from a Poisson distribution and this should make it more conservative compared to t and Wilcox tests, especially if the data is close to normally distributed and underdispersed compared to a Poisson.
Here's a simple simulation of 1000 genes, 10 replicates per condition, true difference is 50 vs 52. Data comes from a normal distribution with standard deviation 5:
library(data.table)
library(MASS)

N <- 1000
nb <- rep(NA, N)
wx <- rep(NA, N)
tt <- rep(NA, N)
alldat <- list()
for(i in 1:N) {
    if(i %% 100 == 0) {
        print(i)
    }
    n <- 10
    set.seed(i)
    dat <- data.table(
        cond= rep(c('A', 'B'), each= n),
        cnt= round(c(rnorm(n= n, mean= 50, sd= 5), rnorm(n= n, mean= 52, sd= 5))),
        i= i
    )
    smry <- summary(glm.nb(cnt ~ cond, data = dat))
    nb[i] <- smry$coefficients[2, 4]
    wx[i] <- wilcox.test(dat[cond == 'A']$cnt, dat[cond == 'B']$cnt)$p.value
    tt[i] <- t.test(dat[cond == 'A']$cnt, dat[cond == 'B']$cnt)$p.value
    dat[, nb:= nb[i]]
    dat[, wx :=wx[i]]
    dat[, tt :=tt[i]]
    alldat[[length(alldat) + 1]] <- dat
}
alldat <- rbindlist(alldat)
sum(nb < 0.01) # 3 genes with p < 0.01
sum(wx < 0.01) # 32
sum(tt < 0.01) # 36

If you replace sd= 5 with sd= 8, thus making the data more dispersed, negbin model is more powerful than the others:
sum(nb < 0.01)
[1] 30
sum(wx < 0.01)
[1] 17
sum(tt < 0.01)
[1] 18

Maybe this helps to explain your results?
