Hi Everyone, I am working on rna seq data, In my data, I have many genes which shows zero read in wild type and in treatment their is a value, so I am not able to calculate fold change and p value.
Many peoples suggested to me to give some number in place of zero, so I assign lowest detectable value after zero and hence calculate fold change. By doing so, I am getting many genes in up-regulated category, I want to limit the number of genes, even I am keeping 4 fold up, there are approx 500 genes.
Can anyone tell me some statistical measure to reduce number of genes which has some rationale too.
I am new to this field,please help.
Are you working with biological sample replicates and all those wildtypes samples are zero all and the corresponding treatments aren't? This may happen but should be rare.
Point is that you want to kown the variance in the treatment and wild type group to have proper estimates of changes and their significance. Working with replicates imo is very important.
Another point is sequencing depth. It should be high enough to capture genes with low expression rates.
I recommend to use either the DESeq2 or EdgeR R-package for analysis of differential genexpression from RNA - seq. The packages should also contain normalization schemes that account for all zero genes.
I recommend the DESeq2 package over EdgeR, since this package automatically filters your data the best way possible prior to testing. You can also choose a fold-change cut-off, before testing, to increase power (fewer tests performed and less stringent FDR calculation).
From a typical testing view point, your problem with informative zero values has similaritiies to a problem with informative missing values. You may consult the theory and methodology as is discussed in Pesarin and Salmaso (2010) "Permutation Tests for Complex Data -Theory, Applications and Software", Wiley, Chicester, pages 232-251.
You should have several samples per group. If not all the value in one group are zero, then you should get a p-value. Note that you should use a Poisson-(or negative binomial) model to get the fold-changes, because you have count data!
As others have said RNAseq tools like edgeR and DESeq2 can deal with zeros. You don't mention how big the difference between zero and non-zero counts is? If it small (
edgeR and DESeq2 has an option called "prior.count", which give you the option to add a value to all genes. This is often set to 1, since log2(1) = 0, so that 0 represent "no expression". However, the higher the value, the more you shrink the fold-changes.
EDIT: Please see Christian Coles answer below.
@Sindre if you read the manual in glmFit() prior.count has no effect on 'the likelihood ratio tests or p-values' so low expr genes will have high FDR/p-value regardless of prior.count. It is essentially there to avoid inflated logFCs caused by low count genes.
My point still stands; don't add arbitrary values to zero counts.
- Badges
- Science topic
- Similar topics
- Filing