Hi all,
Hope someone can help me.
I have an RNA-Seq dataset I am analysing using edgeR. The data relate to a drug treatment at different concentrations (5uM, 10uM and 20uM) vs untreated (UT). I have three replicates for each concentration and three for UT.
My initial goal is simply to look for DE genes between UT vs 5uM and I've gone about this in two ways. In the first I only ingested the count data for the 3 5uM samples and the 3 UT samples, performed the analysis and got ~2K upregulated and ~2K downregulated genes. In the second, I ingested all data for all concentrations, performed the analysis with an altered design matrix and contrast to focus on the UT vs 5uM comparison and got 34 up and 69 downregulated genes. I'm obviously confused as to why I get vastly different numbers of DE genes and obviously must be doing something wrong.
In case my problem lies in my design matrix and contrasts, these are as follows:
First way: UT.0 DRUG.5 1 1 0 2 1 0 3 1 0 4 0 1 5 0 1 6 0 1 Second way: UT.0 DRUG.10 DRUG.20 DRUG.5 1 1 0 0 0 2 1 0 0 0 3 1 0 0 0 4 0 0 0 1 5 0 0 0 1 6 0 0 0 1 7 0 1 0 0 8 0 1 0 0 9 0 1 0 0 10 0 0 1 0 11 0 0 1 0 12 0 0 1 0
I used these as follows:
UTvs5 <- makeContrasts(DRUG.5-UT.0, levels=design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, contrast=UTvs5)
In the second way, I assumed that as the contrast should ignore the 10uM and 20uM data, the results for both methods should be the same in terms of DE.
I would really appreciate it if someone could tell me if I am doing something wrong here or if including the other datasets affects other parts of the standard edgeR pipeline (normalisation, dispersion estimates etc), or might it be something else!
Many thanks in advance.
EDIT - including additional analysis code as requested
At this point the count data have been ingested into countdata and I have a phenotype file called phendat, which contains
SampleID Conc Treatment 1 UT_1 0 UT 2 UT_2 0 UT 3 UT_3 0 UT 4 DRUG_5_31 5 DRUG 5 DRUG_5_32 5 DRUG 6 DRUG_5_33 5 DRUG 7 DRUG_10_34 10 DRUG 8 DRUG_10_35 10 DRUG 9 DRUG_10_36 10 DRUG 10 DRUG_20_37 20 DRUG 11 DRUG_20_38 20 DRUG 12 DRUG_20_39 20 DRUG
Retain genes if they are expressed at a counts-per-million (CPM) above 0.5 in at least three samples
myCPM <- cpm(countdata)
thresh <- myCPM > 0.5
keep <- rowSums(thresh) >= 3
counts.keep <- countdata[keep,]
y <- DGEList(counts.keep)
normalise and produce an MDS plot
mypalette.conc <- brewer.pal(4,"Set1")
col.conc <- mypalette.conc[phendat$Conc]
y <- calcNormFactors(y)
plotMDS(y, col=col.conc, dim.plot=c(1,2), xlim=c(-1.5,1.5), ylim=c(-1.5,1), cex=0.8)
I admit, I should have looked at this more closely before proceeding. Many thanks for setting me straight. To me it looks like there are definite complications with this data. I have checked with the person who generated this data and they are confident that due to the procedures and equipment used it was very unlikely there was a sample swap between DRUG_5_32 and DRUG_10_36. Also I would certainly not expect the 20uM data to be anywhere near the UT. I assume this points to some other hidden batch effect. Please let me know what you think, and many thanks again for your help.
