Hi,

I have multiple rna-seq samples with treated and untreated with a drug. Something which looks like:

tissue    treatment     gender 
patient1    treated    female
patient1    untreated    female
patient2    treated    female
patient2    untreated    female
patient3    treated    female
patient3    untreated    female
patient4    treated    male
patient4    untreated    male
patient5    treated    male
patient5    untreated    male

I am interested to see the DE genes between treated and untreated considering the gender. There might be other factors that add noise to data such as the time of collection, age etc. 

So far what I have done is the following. Is it the correct approach ?  Which would be the best way to get rid of unknown noise and find DE genes ?

library(edgeR)
rawdata <- read.delim("selected_raw_count_matrix.txt",header=T,row.names=1)
y <- DGEList(counts=rawdata) 
y <- calcNormFactors(y)
gender=c("female","female","female","female","female","female","male","male","male","male")
treat=c("treated","untreated","treated","untreated","treated","untreated","treated","untreated","treated","untreated")
design <- model.matrix(~gender+treat)
y <- estimateGLMCommonDisp(y, design, verbose=TRUE)
y <- estimateGLMTrendedDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit)
topTags(lrt)​

You current design assumes that the treatment effect will be the same between males and females. To distinguish between the sexes, I would use a design like this:

> design <- model.matrix(~0 + gender + gender:treat)
> colnames(design)
[1] "genderfemale"                "gendermale"                 
[3] "genderfemale:treatuntreated" "gendermale:treatuntreated"  

The third coefficient represents the log-fold change of untreated over treated in females, while the fourth coefficient represents the corresponding log-fold change in males. To test for treatment-induced DE for each sex, all you need to do is to set coef=3 in glmLRT for females and coef=4 for males.

As for the other factors; I would check on a MDS plot to see whether there's good separation between treated and untreated samples. Similarly, I would have a look at the dispersion values with plotBCV to see if the dispersions are small (at or below 0.05 for samples collected in controlled laboratory conditions) and have a decreasing trend with respect to abundance. If everything looks good, you may not need to include other blocking factors. Otherwise, you might have to look at packages like RUVSeq.

--

Edit: Actually, your design would be even better with a patient blocking effect:

> patient <- factor(rep(1:5, each=2))
> design <- model.matrix(~ 0 + patient + gender:treat)
> design <- design[,-(8:9)] # get back to full rank
> colnames(design)
[1] "patient1"                  "patient2"                 
[3] "patient3"                  "patient4"                 
[5] "patient5"                  "genderfemale:treattreated"
[7] "gendermale:treattreated"  

The first five coefficients represent patient-specific effects on gene expression. The sixth coefficient represents the effect of treatment (over untreated) in females, while the seventh coefficient represents the treatment effect in males. This model ensures that systematic differences in the absolute expression levels between patients do not inflate the variances, so long as the log-fold change between treatment and untreated is consistent. Of course, you have fewer residual d.f. with this larger model. This results in less reliable dispersion estimates and reduced scope for detecting hidden factors of variation, though that's an acceptable price to pay from my perspective.

This experiment has what I have been calling a multi-level design. Designs like this are covered in Section 3.5 of the edgeR User's Guide. The analysis is equivalent to Aaron's edited model.