I have a design of 6 conditions with three replicates each, with expression levels for ~51000 transcription products. I'm using Limma's voom transformation to make my data approximately Gaussian, which by ocular pat-down seems to be the case:

<img alt="" src="http://i.imgur.com/euA1mso.png" style="height:300px; width:421px" />

The mean-variance estimation is visible here - kinda overdispersed in the mid range but it doesn't look critically bad to me:

<img alt="" src="http://i.imgur.com/2Gxw28k.png" style="height:300px; width:463px" />

I fit the linear model thusly:

```
> labels
 [1] A A A C C C G G G L L L T T T U U U
Levels: A C G L T U
design <- model.matrix(~labels)
fit2 <- eBayes( lmFit(voomCounts, design) )
```

However, this yields the following p-value distribution from the pairwise t-tests (looks like a misspecified model)

<img alt="" src="http://i.imgur.com/yyOsYWN.png" style="height:300px; width:421px" />

The `fit2$F.p.vals` has only 10-15 genes *above* 0.01, the rest are absolutely tiny; and when I plot some of the genes called for differential expression they look very unremarkable. What am I doing wrong?

Thanks so much :)
I do not know why are you doing voom counts

but if you want differentially expressed genes from 6 samples, lets say 3 replicates each of cases and controls

`df <-` "is your dataframe containing 6 coloumns (3 for each ) of normalized expression data"

then just go for

```
library(limma)

groups<-as.factor(c(rep("Cases",3),rep("Control",3)))
design<-model.matrix(~0+groups)
colnames(design)=levels(groups)
fit<-lmFit(df, design)
cont.matrix<-makeContrasts(Cases-Control,
levels=design)
fit2<-contrasts.fit(fit, cont.matrix)
ebfit<-eBayes(fit2)
topTable(ebfit, coef=1)
topTable(ebfit, number=Inf, p=0.05, adjust.method="none", coef=1)

## topTable would be your list of DEG
```