Hi, I have a question regarding outlier replacement when including a continuous variable in DESeq2 design. I am currently analyzing Nanostring data using RUVseq normalization and DESeq2 following the method by Bhattacharya et al.
As the normalization factor from RUVseq is included in the design, outlier detection and replacement will not be performed by DESeq as described in the vignette. We do, however, have several genes with outliers, which has large impact on the output. This seem to be due to high counts in single samples for multiple genes (the issue seems to be gene specific, as the high counts coming from different samples). I have been trying different options for filtering these genes, and wanted to know if anyone has some experience on the optimal way to handle this. Especially with some high expression genes repeatedly showing up as the most significant gene in our dataset when running the analysis.
Following the suggested approach by posted here by Michal Love ([DESeq2] Cook's distance flagging for continuous variables), we use max Cook score for each gene to set a Cook cutoff based on your choice (after plotting data for visualization). Using qf(0.99, m,m-p) for cutoff (which is used by DESeq), 57 out of 770 genes are discarded as outliers (we have 126 samples, and 4 groups + normalization factor in our design). While we can be less stringent with the cutoff to discard fewer outliers, the main issue with this approach is that we are filtering the information from these genes altogether. This despite only counts for the respective genes in 1/126 sample seems to be off (at least for several of these genes).
Thus, my question is whether other approaches can be more suitable if we don't want to set the p.adj to NA, like one suggested here (DESeq2: Outliers, covariates and number of replicates)? Eg first running DESeq without continuous variable (in our case the normalization factor), then insert replaceCounts from the resulting dds object (where outlier counts, 35 in our data, have been identified and replaced by replaceOutliers in DESeq). Then this new object is used for running DESeq with the full design including the continuous variable. I know it is generally advised against using any kind of modified counts in the dds for analyses, but is this also the case when replacing count outliers (but keeping all the other original counts)? There might be an issue with this approach, as false “outliers” could be detected (that otherwise would be corrected by including the normalization factor to not be “outliers”).
The generation of the plots is at least similar to how we can assess outliers with a simplified design without the continuous variable, but will keep genes that has one large outlier in the final results rather than discarding the information altogether.
I would highly appreciate any feedback/experience with handling such issues.
Thanks, Eivind
Thanks for the reply Michael,
The design is as follow:
Each of TIL_score and Dist_met have two groups, hot/cold and yes/no. Altogether the dataset consist of 126 samples, with the following distribution considering both factors in the design (besides W_1):
The relatively large number of samples, is also the main reason why I don’t want to discard genes but rather replace counts with one outlier only (as exemplified with the PLA2G2A gene below – normal counts left, replaced counts right).
Making an assay containing gene weights I get an error due to dimnames not being identical in the new “weights” matrix:
I assume identical dimnames are not necessary for assigning individual weights as long as they are in the corresponding position to genes/cooks matrices, right (eg setting withDimnames=FALSE would be fine)? For applying manual cooks cutoff using weights, I used the following code:
Why is it preferred to first run DESeq without W_1, and generate the weights of the corresponding cooks? Including W_1 gives overall higher cooks values, and in our data genes above cooks >10 increases from 27 to 46, thus, it is mainly to have a less strict approach when "flagging" genes as outliers? I performed the analysis by including "weights" as you suggested, and also by replacingCounts as described before, and in general it seems to give quite similar outputs, but with more drastic differences on the padj with the replaceCounts for some of the "trouble" genes, liked PLA2G2A highlighted below. For the results below I made a grouped design (a variable with 4 categories, based on TIL_score and dist_met) to more easily compare two subgroups, No:hot vs Yes:hot. This to illustrate the initial issue, as PLA2G2A is identified as upregulated when comparing Yes:hot with No:hot in DESeq, but assessing the counts we can see this is mainly due to one value, and (if anything) the trend is more the other way around when this count is removed.
While both methods identifies and reduced impact of assumed outlier, through reduction of p-values, it does seem by replacing the original counts is the only way of reducing the impact of outlier on the results like PLA2G2A? If you are not strongly advising against analyzing on "modified" counts in this case, it seems more suitable for our analysis purpose. Also for visualization of difference in counts between groups as shown above, the replaceCounts approach is the only way (besides ignoring outliers if boxplots etc are used)?
Again, thank you so much for your time.
Eivind
Yes, that code is fine (you can safely use
withDimnames=FALSEthere). It's just complaining that the empty weights matrix had no row/colnames.Good question. I was just thinking to be less extreme about calling outliers. RUV factors are fit to the data, and so they are helping to explain variability (lowering dispersion estimates). When there are fewer design terms, the estimated dispersion increases, and we are less likely to call things "outliers".
I'd prefer the weights-based method with these sample sizes. We've used weights in a couple projects, I just haven't had time to automate this (and anyway it takes some careful looking by eye to convincingly identify "outliers" I think). Did PLA2G2A have a sample with Cook's distance > 10? It seems the LFC didn't change with the use of the weights matrix?
Thanks for a detailed explanation.
When it comes to the odd LFC for PLA2G2A after applying weights - you are totally correct that PLA2G2A is having a Cook's distance less than 10, and followingly not getting adjusted weights. I initially ran the code setting the weights based on the design including W_1, and assessed the Cook's distance based on this (almost 15). Without W_1, it drops to 9, explaining why the LFC of PLA2G2A doesn't change in the previous post. This is the output with reducing the Cook's distance cutoff to 8:
I will adjust the cutoff accordingly, to also include these gene I initially identified as outliers, and go ahead with the weight adjustment.
Again, thank you so much for your time.
That seems like the best solution. Glad to see it worked.