Hello,

I'm working on a RNA-seq differential expression analysis on human samples, and have been trying to understand the effect of the various packages that allows the weighting of samples or features - e.g. voom's sample quality weights, or edgeR-robust's implementation of observation weights.

I'm doing this because I'm particularly concerned about the effect of extreme outliers in my analysis - i.e. samples with counts for a gene that are 2 orders of magnitude different to the average across the entire sample set. Based on my experience of previous differential analyses on this dataset, I'm worried about a large proportion of the significant genes being called as significant due to the effect of 2-3 of these hugely outlier samples. (The dataset consists of 48 samples, divided into 4 groups.)

I'm also wondering about batch effects in my data, so I've tried out the sva package as well, which - from my understanding - uses a similar sample weighting approach to account for confounding variables.

My question is: given that the weights (of all types) are generally used in the fitting of the GLM to call differential genes, is it reasonable to include the weights when plotting a PCA or MDS, to see the effect of the weights on the sample group separation?

To my mind, if you're doing the MDS/PCA on the unweighted count data, then you weight the data before doing the actual differential analysis, it should still be valid to include the weights at the MDS/PCA stage.

Thanks in advance for any suggestions and insights into the statistics!

If you include the weights in the PCA or MDS, you by definition cannot see the effect of the weights because the plot now includes the weights! Anyway, the weights aren't used to modify the data - instead they are used as part of the linear regression to decrease the impact of the outliers on the model fit. There is a pretty good explanation here (http://www.itl.nist.gov/div898/handbook/pmd/section1/pmd143.htm and http://www.itl.nist.gov/div898/handbook/pmd/section4/pmd432.htm).

Instead, I generally make the PCA plot and then plot the weights as well (right above or below each sample), so I can see which samples are being down-weighted, and by how much.

Apologies for the delayed response. The voomWithQualityWeights() function combines observational level weights (i.e. a distinct weight for each and every gene expression value) from voom with sample-specific weights by multiplying them together. Running this function with plot=TRUE provides a summary of the abundance related trend in variability estimated by voom in one panel and sample-specific variability in a second panel. Sample weights are often (but not always) related to how well samples cluster in the PCA / MDS plot (i.e. lower weights are assigned to samples that cluster less well).

In practice, trying an analysis with different weighting schemes (i.e. in the case of RNA-seq voom only versus voomWithQualityWeights) and seeing if you get more differentially expressed genes is a good way to empirically determine whether the extra effort of modelling sample-level variability was worthwhile. I hope this helps.