deseq2 filter the low counts
2
3
Entering edit mode
aristotele_m ▴ 40
@aristotele_m-6821
Last seen 7.0 years ago
Italy

Dear all,

How can I filter the  counts with low count in Deseq2?  Any suggestion on how to do?

dds<-DESEq(dt)

count<-counts(dds,normalize=TRUE)

filter<-rowsum(count)> 10

thanks so much!!

deseq2 lowcount • 49k views
ADD COMMENT
17
Entering edit mode
@mikelove
Last seen 18 hours ago
United States

hi,

If you want to filter, you can do so before running DESeq:

dds <- estimateSizeFactors(dds)
idx <- rowSums( counts(dds, normalized=TRUE) >= 5 ) >= 3

This would say, e.g. filter out genes where there are less than 3 samples with normalized counts greater than or equal to 5.

then:

dds <- dds[idx,]
dds <- DESeq(dds)

However, you typically don't need to pre-filter because independent filtering occurs within results() to save you from multiple test correction on genes with no power (see ?results and the vignette section about independent filtering, or the paper). The main reason to pre-filter would be to increase speed. Designs with many samples and many interaction terms can slow down on genes which have very few reads.

ADD COMMENT
0
Entering edit mode

Does it mean that DESeq2 has no problem estimating dispersions for low expressed features, as opposed to voom?

ADD REPLY
0
Entering edit mode
In our DESeq2 paper we discuss a case where estimation of dispersion is difficult for genes with very, very low average counts. See the methods. However it doesn't really effect the outcome because these genes have almost no power for detecting differential expression.
ADD REPLY
0
Entering edit mode

Sorry, right now I faced a point, please help me to be cleared if I am wrong. I noticed that in differential expression analysis by DESeq2, the distribution of read counts of differentially expressed genes is in favour of more highly expressed genes. I mean, likely DESeq2 has a threshold for ignoring too low expressed genes before differential expression analysis. Actually I was expected the genes with too low read counts or zeros are the reason of differential expression but box plot shows that the DE genes are among the genes with higher reads counts.

ADD REPLY
1
Entering edit mode

More highly expressed genes have higher power for detecting DE.

And yes we do have an internal filter that optimizes this.

Like all important aspects of the method it is discussed in the paper and in the vignette.

ADD REPLY
0
Entering edit mode

Thank you, you alright

ADD REPLY
0
Entering edit mode

Sorry, by using these lines

dds <- DESeq(dds, minReplicatesForReplace=Inf)
res <- results(dds, cooksCutoff=FALSE, independentFiltering=FALSE)

will I prevent internal filtering in DESeq2 to remove any genes in differential expression? I am right??

ADD REPLY
1
Entering edit mode

Yes this will turn off independent filtering (on the mean of counts), as well as outlier replacement and outlier-based gene filtering.

ADD REPLY
0
Entering edit mode

Hi Michael, thanks for your posts - they are really helpful! I was wondering though, isn't there any issue with using `estimateSizeFactors(dds)` twice? Because the DESeq function is going to do this again, no?

 

ADD REPLY
1
Entering edit mode

DESeq() does not re-estimate size factors. It will print this message also when you run it.

ADD REPLY
0
Entering edit mode

Hi Michael, would like to have an update on your explanation on this strip of code:

idx <- rowSums( counts(dds, normalized=TRUE) >= 5 ) >= 3

You said that this means that it would filter out genes where there are less than 3 samples with normalized counts greater than or equal to 5. But in turn, isn't it the opposite since we have the ">=" symbol .So it filters out genes that are more than 3 samples ,right?

Hope to hear from you regarding this.

ADD REPLY
2
Entering edit mode

Better to call it keep:

keep <- rowSums( counts(dds) >= X ) >= Y
dds <- dds[keep,]

This requires genes to have Y or more samples with counts of X or more. It therefore filters out genes that have less than Y samples with counts of X or more.

ADD REPLY
0
Entering edit mode

Hi Micheal,

I would like to know how I can choose the X and Y numbers here? How do I have to check it in my data to know how I have to put the cut-off and keep the samples?

Looking forward to hearing from you. thanks!

ADD REPLY
2
Entering edit mode

X, a good value is 10

Y, choose something like, the smallest group sample size

ADD REPLY
0
Entering edit mode

Thank you fo your comment. Just to make it clear for myself regarding the smallest group sample size; you mean if I have 20 individuals as cases and 24 individuals as controls, take the number of 20? (is it related to the number of cases or is it something else) ? keep <- rowSums( counts(dds) >= 10 ) >= 20

Am I right? thanks! looking forward to hearing from you

ADD REPLY
0
Entering edit mode

Yes, 20 (smallest sample size of the groups).

ADD REPLY
0
Entering edit mode

Hi Michael Love ,

What if there are control and disease groups, do you have any comments on filtering based on groups? For example, the X and Y cutoffs are applied to each group instead of all samples?

Thanks.

ADD REPLY
0
Entering edit mode
pkachroo ▴ 10
@pkachroo-11576
Last seen 3.5 years ago

I have a similar question: In an experiment with 5 strains in triplicates, I have a gene with the following normalized counts:

              Replicates

Strain-1: 0,0,0

Strain-2: 0,0,0

Strain-3: 1.6,1.3,0

Strain-4: 0,0, 2.6

Strain-5: 105,102,101

After running DESeq2, this gene is flagged and given "NA" for pvalue and adjusted.value, which makes sense. However, when I rerun the analysis with only first two replicates per strain (highlighted bold) and compare strain 5 and 4, this gene comes up as differentially expressed: baseMean=9.5 and log2FoldChange=3.3. I am wondering why is this gene not being flagged? and more importantly, how is deseq2 able to compute a fold change when the normalized counts for this gene in strain-4 are zeros. 

Appreciate your help.

Priyanka Kachroo

 

ADD COMMENT
1
Entering edit mode

The question about calculating fold changes when strain 4 has zeros has been answered on the site a few times but it's difficult to find the post. The short answer is that the DESeq2 statistical model (see paper) uses a prior distribution on the fold changes, and returns posterior estimates. So the posterior is a balance of the likelihood (which would give an infinite fold change) and the prior which is calculated based on the range of fold changes from the most DE genes.

Regarding the NA's:

If you read in the help page ?results about NA values in the pvalue column:

By default, results assigns a p-value of NA to genes containing count outliers, as identified using Cook’s distance. See the cooksCutoff argument for control of this behavior.

Then if you read more:

cooksCutoff - theshold on Cook’s distance, such that if one or more samples for a row have a distance higher, the p-value for the row is set to NA. The default cutoff is the .99 quantile of the F(p, m-p) distribution, where p is the number of coefficients being fitted and m is the number of samples. Set to Inf or FALSE to disable the resetting of p-values to NA. Note: this test excludes the Cook’s distance of samples belonging to experimental groups with only 2 samples.
ADD REPLY
0
Entering edit mode

Hi Michael, I have another question on a similar note.

In my project, I'm comparing 10 AML samples with TP53 mutation with 6 AML samples with wild-type TP53 - so there's a little sample imbalance. I would like to filter out genes with counts of 5 or more in less than half the amount of samples per group. Since there's two groups, TP53 MUT and TP53 WT, there'd be two filters:

keep1 <- rowSums( counts(dds1) >= 5 ) >= 5    # Where dds1 is a deseq dataset only with the TP53 MUT samples (5 is half of my N of 10) 
keep2 <- rowSums( counts(dds2) >= 5 ) >= 3    # Where dds2 is a deseq dataset only with the TP53 WT samples (3 is half of my N of 6)

Is it possible to perform such filtering, with one sample threshold per group? And if so, how would the code for manipulation of my dds (that so far contains all genes and all samples) would look like? - breaking it into 2 dds datasets (one per group), applying the two filters, merging the 2 dds datasets into a final one, removing any duplicates, and then running Deseq2? Or is there a way of applying this filter with one threshold per group without the need to create a dds per group?

Thank you so much in advance,

ADD REPLY
0
Entering edit mode

You can't filter the groups separately, you need to use a non-specific filter that doesn't make use of information about which samples are in which design.

Hence we say to use the smallest group size in the count filter. And then you would filter the whole matrix without making use of the sample design at all.

For more on why you can't use the sample information see this paper:

https://www.pnas.org/doi/10.1073/pnas.0914005107

ADD REPLY

Login before adding your answer.

Traffic: 427 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6