Hello everyone

I`m working with RNA-seq gene expression data derived from multiple tissue samples collected from different patients. My primary goal is to identify differentially expressed genes (DEGs) while minimizing the confounding effects of tissue of origin on the results.

A brief overview of my approach:

I`ve converted the raw count expression data into CPM values using edgeR to normalize between samples, accounting for library size. Filtering was applied using filterByExpr to retain relevant genes. Normalization was conducted using the TMM method:

keep <- filterByExpr(counts)
counts <- counts[keep, , keep.lib.sizes=FALSE]
counts <- normLibSizes(counts, method = "TMM")
counts_cpm <- cpm(counts, log = TRUE)

To account for patient and tissue variability, I applied a mixed linear model using the variancePartition package. My formula models the contribution of both patient and tissue to the gene expression variation:

form <- ~ (1 | Tissue) + (1 | Patient)
vp_modelFit <- fitVarPartModel(counts_cpm, form, df)
vp_modelFit_res <- residuals(vp_modelFit)

My understanding is that the residuals from this model should, in theory, represent gene expression values devoid of tissue- and patient-specific effects, potentially revealing the intrinsic cancer-related signals.

Question:

1. Is this approach statistically sound for achieving my aim? Specifically, does this methodology appropriately remove the unwanted variation from tissue and patient sources while preserving biologically relevant signals?
2. Any recommendations for improving the robustness of this approach, especially in terms of ensuring that intrinsic cancer-related signals are not inadvertently removed?

Thank you!

1. You should always use voomWithDreamWeights() on the counts before fitting the model. This produces log2 CPM and computes precision weights

1. Doing an analysis in two steps by a) computing residuals and then b) performing a second regression for statistical methods will increase the false positive rate since the covariance between the variable of interest and other variables is ignored. I _strongly_ recommend fitting a single model with ~ (1 | Tissue) + (1 | Patient) + variable and using dream()

2. These issues are not specific to the dream workflow, they are general properties of regression that apply to limma as well

I generally agree with Gabriel Hoffman although, without knowing more about the design of your experiment, it is impossible to say whether treating Tissue and Patient as random is appropriate.

No, your proposed approach is not statistically sound. The whole purpose of a statistical analysis is to evaluate systematic changes relative to variation between patients and tissue samples, so trying to remove patient and sample variation before the analysis begins doesn't make sense. Such an approach will give an unrealistic impression of significance.