too many de genes

[mtiwaridros]

Hi

I have processed my paired-end rna-seq data (2 conditions x 3 replicates) based on the following pipeline:
mapping with STAR>rsem-calculate-expression>de testing with ebseq

However, at the normally accepted filters (log2FC>2 AND adjustedp<0.05), I am getting close to 8000 differentially expressed genes.
Changing the filter values (for instance, reducing adjustedp to 0.001) still results in ~7000 de genes.

This is a large list, and I do not expect my samples to have so many differences.

Can anyone suggest a resolution to the issue?

Thanks.

Best
M

[Birdman]

I would say this is caused by the low power of statistics on n=3. You basically get a lot of DE genes by chance. P-values tend to be extremely low in these kind of stats so you can consider having higher thresholds. You could also try to have a higher n if possible, change the stats (for example try a glm with multiple variables, if you have others e.g. sex, age, ...) or use GO terms to have groups of genes to compare. BTW, you are not alone, most people face this problem with RNA-seq

[Bukowski]

you are not alone, most people face this problem with RNA-seq

Having run lots and lots of RNA-Seq experiments I would counter that 'most' people face this problem. I have certainly seen experiments with upwards of a thousand DE (on interfering with splicing machinery), but I've also seen very small results sets from very specific siRNA KD experiments. This high level of DE is not usual in my experience.

But you can't draw any conclusions about too few replicates, when to be honest, at least some replicates have actually been run.

I think the real question that needs asking here was 'What was the experiment you ran?' and 'how close are the replicates to each other?'.

[mtiwaridros]

Hi

Thanks for the inputs.
My replicates look similar (PCA, MA plots) and the mapping quality is also similar in each of them (70-80%).
I agree that there must be many false positives in the data, and this begets the question of how to proceed from here?
I do have some specific questions in mind, like impact on cell cycle, specific tfs etc. but wouldn't pre-selection by supplying GO terms bias the output.
Also, there might be several lowly-expressed transcripts (belonging especially to, say, miRNA regulators) that might be filtered out if I increase the thresholds.
Surfing around, it seems wgcna might be a viable option.
Any suggestions.

Thanks again.
Best
M

[Bukowski]

Also, there might be several lowly-expressed transcripts (belonging especially to, say, miRNA regulators) that might be filtered out if I increase the thresholds.
Surfing around, it seems wgcna might be a viable option.

It's true that in some experiments I have a lot of stuff expressed at very low levels that has significant p-values and reasonable foldchanges, but I tend to apply some kind of minimum FPKM filter in both conditions - which generally cuts a swathe of data out, as I just don't trust data hovering just above background. I've heard of, but not used WGCNA, so I'm afraid I can't advise on that.

[dariober]

It's true that in some experiments I have a lot of stuff expressed at very low levels that has significant p-values and reasonable foldchanges, but I tend to apply some kind of minimum FPKM filter in both conditions - which generally cuts a swathe of data out, as I just don't trust data hovering just above background. I've heard of, but not used WGCNA, so I'm afraid I can't advise on that.

Hello,

Loosely related to this topic... If after differential expression analysis there are too many genes differentially expressed at FDR < x, one could rank them by logFC. However, this tends to select genes towards the low expressed range. On the other hand, ranking by FDR will typically skew the selection towards genes highly expressed even if they have small logFC.

An alternative I found quite nice is to apply an "intensity dependent" filter where genes are ranked according to how extreme they are in terms of logFC relative to genes with similar expression level. On an MAplot it means choosing genes towards to periphery of the cloud. Interesting genes would be those with FDR < x and |z-score| > y (say z-score > 2, where the z-score is the measure of how far a gene is from the cloud).

This helps in filtering genes when lots of them are DE. Then the reason why one gets 1000s of genes DE depends on the experimental design, of course.

I wrote a couple of R functions (attached) to implement the idea. Here's an example usage:

source('intensityFilter.R')

## Simulate some data
set.seed(1234)

## Expression level
logCPM<- sort(rnorm(n= 10000, mean= 5, sd= 1))

## Log fold change
set.seed(12345)
logFC<- sapply(length(logCPM):1, function(i) rnorm(n= 1, mean= 0, sd= i^(1/1.5)))

smoothScatter(x= logCPM, y= logFC, nrpoints= 1000)

## z-score
z<- localZ(logCPM, logFC, nbins= 20)

## Highlight genes on the periphery of the cloud
points(logCPM, logFC, col= ifelse(abs(z) > 1.5, 'red', NA), pch= '.', cex= 0.5)

(NB: The idea is not mine, I got it in part from this post by simonandrews http://seqanswers.com/forums/showpos...8&postcount=34)

Any comments much appreciated!

Dario