DESeq and edgeR up/down regulation

[murphycj]

Hi,

I' am relatively new to bioinformatics.

I' am assisting in analyzing 7 RNA-seq samples. Two conditions: (1) experimental - 6 samples, and (2) control - 1 sample.

The input for DESeq consisted of read counts. See below for the R code for DESeq:

---------------------------------------------------------------
library(DESeq)
readCounts<-read.csv("readCounts_IndividualDigEx",row.names=1,header=TRUE)

#which data are from MPSNT and Control
group<-c(rep("MPNST",4),"Control",rep("MPNST",2))

#annotate the data
cds <- newCountDataSet(readCounts,group)

#calculate library sizes and associated affect on DE testing
cds <- estimateSizeFactors( cds )

cds <- estimateVarianceFunctions( cds )

#test for DE
res <- nbinomTest( cds, "MPNST", "Control" )

#there is no expression for a gene accross all samples... need to remove those genes
whichIsNA<-which(is.na(res[,8]))
res<-res[-whichIsNA,]

#get list of DE genes, then display some of up and down-regulated genes
resSig <- res[ res$pval < .05, ]
#up-regulated
head(resSig[order(resSig$foldChange,-resSig$baseMean),])
#down-regulated
head(resSig[order(-resSig$foldChange,-resSig$baseMean),])
---------------------------------------------------------------

See the following for example data:

a3 a4 a5 a6 a7 a8 a9
MKL2 1521 991 1058 286 1965 1376 1026
OOSP1 0 0 1 3 0 1 0
CHL1 12 11 14 5154 3565 20 4


My problem is the following. The list of DE genes severely under represents up-regulated genes (low expression in control, high in experiment). For example, with a pval<0.05 (FDR corrected): 215 genes were DE, 211 were down-regulated and 4 were up regulated.

The issue does not seem related to the biology of the samples, but perhaps an error in setting up the tool. For example, I added two mock genes to test for DE:

a3 a4 a5 a6 a7 a8 a9
test1 10000 10000 10000 10000 1000 10000 10000
test2 5000 5000 5000 5000 1000 5000 5000

The mock genes were concatenated onto the other genes. The following are the p-values for each mock gene: test1: 0.25, test2: 0.30, and the adjusted p-values for both are 1.

I understand that the tool confers much less confidence to expression measurements from conditions in which there is only one replicate, but the above outcome seems rather incorrect.

Moreover, here is a sample gene which DESeq found to be significantly up-regluated (p-value = 0.0016):

a3 a4 a5 a6 a7 a8 a9
PAX7 4346 1 0 19 0 773 86


Any advice would be greatly appreciated!

Charlie

[mbblack]

Do you mean you used RPKM values as input for DEseq? If so, you need to redo that analysis using raw count data. Tools like DEseq and DEGseq are expecting to have raw or simple count data, not normalized data.

P.S. I use Bedtools' coverageBed function to get my raw counts from my mapping run BAM files (usually using a UCSC RefGene BED file as the reference).

[murphycj]

My apologies, the input data for DESeq were the read counts. I incorrectly said RPKM.

[Simon Anders]

You might have a problem with normalization. For each of your treatment library, plot the libraries individual fold w.r.t. the control library against the mean and compare with the size factors:

plot( ( counts(cds)[,1] + counts(cds)[,5] )/2, counts(cds)[,5] / counts(cds)[,1], log="xy", pch="." )
abline( h = sizeFactors(cds)[5] / sizeFactors(cds)[1] )

The horizontal line should go through the middle of the bulk of the genes, especially at the right-hand side of the plot. If it does not, the 'estimateSizefactor' function has not worked well. Maybe try this and let us know.

[murphycj]

I did what you suggested, and the plot seems to suggest that the 'estimateSizeFactors' works correctly. The horizontal line goes through the middle of the bulk of the genes. The plot looks very similar across all library comparisons (e.g. 1v5, 2v5, 3v5, etc..).

See the attached .pdf for an example.

Thanks for your reply.

[murphycj]

Maybe the following two plots will help.

(1) I have attached a histogram of the p-values. The distribution seems to indicate an error in the statistical test.

(2) The MvA plot.

Hope they help.

Additionally, the mean fold change is -0.169 with a standard deviation of 1.5.

[Simon Anders]

This seems all fine to me.

The reason for the strong asymmetry in detection power depending on the direction of change is due to your extremely unbalanced design. To get the same precission of expression strength estimation, a gene has to be much stronger expressed in control than in treatment. Hence, if a gene is low in control and high in treatment, the fold change cannot be estimated as well than if it the other way around. DESeq takes this into account and so is more reluctant to call down-regulated genes significant.

[murphycj]

Seems you are right. To test if the asymmetry is due to the unbalanced design, I replicated the the control expression data set (i.e. I inputted into DESeq 8 samples: 6 experimental and 2 control). I of course randomly modified the second control expression values by multiplication by a random number, which had a mean of 1 and low standard deviation.

What results is many fold more DE genes (due to the similarity of the two control data sets), but more importantly, roughly equal amounts of up and down-regulated genes.

Thanks for your input on this issue. The reason for the biased DESeq results make sense at least.

Our collaborators stated that they are in the process of adding more controls to the data.