Hi Lima and DESeq2 community,

I'm trying to study treatment effects on two different Phenotypes (say, Disease and WT). 7 patients were collected for the Disease group while 6 normal people were collected for WT. For each of these 13 subjects, RNA-seq data were generated after 3 different Treatments which were a Vehicle control, a Drug A treatment and a Drug B treatment.

From the PCA plot, these 39 samples are clearly separated into two groups according to Phenotype. In the following analysis I want to explore:

1) difference between WT Vehicle and Disease Vehicle (baseline difference, as expected and indicated by PCA)

2) effect of Drug A on WT (off-target effect of A)

3) effect of Drug B on WT (off-target effect of B)

4) difference between WT with Drug A and Disease with Drug A (effect of A specific to Disease)

5) difference between WT with Drug B and Disease with Drug B  (effect of B specific to Disease)

I applied Limma, as it can answer all the 5 questions within one model (Thanks for @Gordon Smyth suggestion). I combined Phenotype and Treatment and applied a design blocking for Individuals (whose level is 1 to 13 here) in the experiments:

combined.cond <- factor(paste(group$Phenotype,group$Treatment,sep="."))

design <- model.matrix(~0+combined.cond)

colnames(design) <- levels(combined.cond)

v <- voom(dge,design)

corfit <- duplicateCorrelation(v,design,block=group$Individual)

fit <- lmFit(v,design,block=group$Individual,correlation=corfit$consensus)

cm <- makeContrasts(Disease.Veh_vs_WT.Veh = disease.Vehicle - wt.Vehicle, WT.A = wt.A, WT.B = wt.B, Disease.A_vs_WT.A = disease.A - wt.A, Disease.B_vs_WT.B = disease.B - wt.B, levels=design)

fit2 <- contrasts.fit(fit, cm)

fit2 <- eBayes(fit2)

By using: 

topTable(fit2, coef="Disease.Veh_vs_WT.Veh",number=30000,p.value=0.05)

I found 618 genes were differently expressed between WT Vehicle and Disease Vehicle.

topTable(fit2, coef="Disease.B_vs_WT.B",number=30000,p.value=0.05)

I found 1704 genes were differently expressed between WT with Drug B and Disease with Drug B  (effect of B specific to Disease). But then I found something "wired":

topTable(fit2, coef="WT.B",number=30000,p.value=0.05)

This should tell me how many genes are differently expressed when treating WT with Drug B (off-target effect of B). This number is 19820, more than half of genes were changed. Actually Drug B was previously reported as a therapy for this specific disease. I tried to summary the logFC of these 19820 genes and get

> summary(resx$logFC)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 

 -6.720  -2.502   2.578   1.603   5.094  15.410 

Seems most changes were quite dramatic, which made no sense according to publications out there. In addition, a very closely number was found for Drug A effect on WT. 

How can I explain this? Am I applying a wrong model? Any suggestion is really appreciated. Thanks in advance.

Update: (Sorry @Michael Love, I'm back again...) 

I also tried to nest both Individual and Treatment under Phenotype (i.e. repeated measurements on individuals).  I used DESeq2 with the following design, although I will need another model on a subset of full data to answer question (1):

design <- ~Phenotype+Phenotype:Individual.nested+Phenotype:Treatment.

It returns about 8700 genes differently expressed for Drug B effect on WT. Still a number seems to be extremely high. But most of them have small fold changes. Apply a naive cutoff at log2FC=1 will leave around 1000 genes:

> summary(resy$log2FoldChange)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 

-6.07000 -0.35700  0.08636 -0.01365  0.33110  5.02800 

Does this log2FC=1 cutoff make sense here? I'm going to divide genes with adjusted p < 0.05 into two groups: log2FC >=1 and log2FC<1, and run a GO analysis on them. 

In your code, you drop the wt.B coefficient alone. This is a gibberish contrast, as you're using a design matrix without an intercept; each coefficient represents the average log-expression of the named group, so your contrast named WT.B is effectively testing whether log-expression is equal to zero. Obviously, many genes will have absolute log-expression different from zero, which is why you get loads and loads of DE genes. If you want to test for DE upon treatment with drug B, your contrast should look something like:

con <- makeContrasts(wt.B - wt.vehicle, levels=design)

As an unrelated aside, you could also consider testing for a differential drug response between WT and disease:

con <- makeContrasts((wt.B - wt.vehicle) - (disease.B - disease.vehicle), levels=design)

This might give you more relevant results than a direct comparison between the WT/disease treated samples, as it would correct for any disease-specific differences in baseline expression before treatment (i.e., DE that's already there in the vehicle).