Question: I'm analysing a microarray (single-channel) study with 2 experimental groups, each one with only 2 samples.
I wanted to compare the two groups in order to find DE genes, so, I ran limma analysis on this dataset.
The results were weird...lots of genes w/ very low p-values (lowest p-value: 6.144e-213). I made a volcano plot to see what was happening with the data...got this **weird** volcano plot: http://s33.postimg.org/mxnx39j67/weird_volcano.png

Do you guys know what kind of stuff would lead to that?

>print(annot)
GEO_ID group
1 GSM321605 uninfected
2 GSM321606 uninfected
3 GSM321607 infected
4 GSM321608 infected

>design <- model.matrix(~0 +annot[, "group"])
>head(design)
annot[, "group"]infected annot[, "group"]uninfected
1 0 1
2 0 1
3 1 0
4 1 0

>colnames(design) <- c("infected", "uninfected")
>cm <- makeContrasts(InfvsUninf = infected-uninfected, levels=design)
>print(cm)
Contrasts
Levels InfvsUninf
infected 1
uninfected -1

>fit <- lmFit(exprs, design)
>fit2 <- contrasts.fit(fit, cm)
>fit2 <- eBayes(fit2)
>results <- topTable(fit2, "InfvsUninf", number=Inf)
>head(results)
ID logFC AveExpr t P.Value adj.P.Val
IDO1 7.496098 8.077363 31.43580 6.144487e-213 8.017326e-209
TNFAIP6 6.945901 7.551685 29.12848 1.359307e-183 8.868116e-180
RSAD2 6.532472 8.066029 27.39471 6.405851e-163 2.786118e-159
IFI27 6.230206 7.727492 26.12712 1.461824e-148 4.768471e-145
IFITM1 6.128964 6.775055 25.70255 6.766299e-144 1.765733e-140
INHBA 5.915550 7.068516 24.80758 2.674626e-134 5.816421e-131

Answer: I think one of the reasons that you are getting weird p-values is because the **filtering of the dataset is not proper** or the standard error across all the probes in your dataset is large. Limma calculates **moderated t-statistic** which uses posterior variance rather than sample variance in t-test. Or, in the layman language it is is the ratio of the M-value to its standard error similar to T-statistic. However, in this case the standard errors are moderated across all genes/ probes in your `exprs`. You may find more about it [here][1]. I would recommend you to use all the probes or use `varFilter` to filter the low expressed probes/genes.

Code:

> dim(exprs)
[1] 54675 4
> design <- model.matrix(~ group + 0)
> colnames(design) <- levels(group)
> head(design)
uninfected infected
1 1 0
2 1 0
3 0 1
4 0 1
> cont.matrix <- makeContrasts(InfvsUninf = infected-uninfected, levels=design)
> print(cont.matrix)
Contrasts
Levels InfvsUninf
infected 1
uninfected -1
> fit1 <- lmFit(exprs, design)
> fit2 <- contrasts.fit(fit1, cont.matrix)
> fit2 <- eBayes(fit2)
Warning message:
Zero sample variances detected, have been offset
> results <- topTable(fit2, "InfvsUninf", number=nrow(exprs.dat))
> head(results)
logFC AveExpr t P.Value adj.P.Val B
210029_at 7.496098 8.077363 53.69020 5.996081e-09 0.0001206758 9.926007
202411_at 6.230206 7.727492 52.40836 6.882552e-09 0.0001206758 9.869859
201601_x_at 5.956054 7.549371 52.02982 7.173197e-09 0.0001206758 9.852669
222838_at 5.748147 7.996823 49.34098 9.709206e-09 0.0001206758 9.721920
206025_s_at 6.945901 7.551685 48.24558 1.103574e-08 0.0001206758 9.663952
201108_s_at 5.564837 7.035012 45.45633 1.549903e-08 0.0001401571 9.502414

Now, filtering the low variable probes and hence, the value of moderate t-statistic (t column) changes. P-value is dependent upon t value and therefore, it become lower than the previous one.

> library(genefilter)
> exprs = varFilter(exprs)
> dim(exprs)
[1] 27337 4
> design <- model.matrix(~ group + 0)
> colnames(design) <- levels(group)
> head(design)
uninfected infected
1 1 0
2 1 0
3 0 1
4 0 1
> cont.matrix <- makeContrasts(InfvsUninf = infected-uninfected, levels=design)
> print(cont.matrix)
Contrasts
Levels InfvsUninf
infected 1
uninfected -1
> fit1 <- lmFit(exprs, design)
> fit2 <- contrasts.fit(fit1, cont.matrix)
> fit2 <- eBayes(fit2)
> results <- topTable(fit2, "InfvsUninf", number=nrow(exprs.dat))
> head(results)
logFC AveExpr t P.Value adj.P.Val B
210029_at 7.496098 8.077363 33.85290 3.055145e-50 8.351851e-46 103.28349
206025_s_at 6.945901 7.551685 31.34093 1.140235e-47 1.558530e-43 97.62137
213797_at 6.532472 8.066029 29.35202 1.661339e-45 1.513868e-41 92.82696
202411_at 6.230206 7.727492 28.24293 3.021780e-44 2.065160e-40 90.02318
204698_at 6.131048 7.734954 27.51300 2.145762e-43 1.173174e-39 88.12390
214022_s_at 6.128964 6.775055 27.39214 2.980586e-43 1.358005e-39 87.80514

In case you are interested in top ~ 5000 genes/probes only then use `varFilter(exprs, var.cutoff = 0.9)`. Try playing with various different values of var.cutoff to get desired number of genes.

[1]: https://stat.ethz.ch/pipermail/bioconductor/2004-September/006132.html