Adjusting for covariates such as duration of diabetes, age, sex and BMI

Question

[sukeshinik5]

Hello,

In RNA sequencing data, I want to adjust for the covariates such as age, BMI, sex etc. I have created colData with these covariates and factorised (for binary) and scaled (for continuous) them. My concern is I have normalised the data and then applying adjustment. However, during adjustment using both voom and removeBatchEffect, will they double adjust the data or is it the right way to perform this. Please see below the functions used.

design_matrix <- model.matrix(~Age + Gender + BMI + 
                                Batch + Condition, data=colData)
v <- voom(raw_counts, design_matrix, plot=TRUE)
fit <- lmFit(v, design_matrix)
adjusted_counts <- removeBatchEffect(v$E, covariates=phenodata[, c("Age", "BMI", "Batch")]).

This is because I wish to use this matrix for correlation studies and not for differential gene expression analysis.

Is there any other way to adjust for covariates?

Thank you,

[Gordon Smyth]

Cross-posted to Bioconductor https://support.bioconductor.org/p/9158852/

Answer 1

[dariober]

To obtain batch-corrected expression values, you could use re-start from the raw matrix of counts and apply ComBat_seq (from the sva package) to obtain batch normalised counts. Then convert these counts to cpm with edger::cpm.

[sukeshinik5]

Please correct me if i am wrong, do i separately adjust for biological covariates using limma and batch effect using ComBat? what shall be performed first?

[dariober]

You apply ComBat_seq on the raw counts (your raw_counts object, probably). You can include batch effects to be removed and, optionally, biological effects to preserve. For differential expression you still work as in your post (voom on raw_counts then lmFit). Note that in this way the workflow to correct expression values is going to be completely unrelated to the differential expression workflow so there's no point in asking which one comes first.

[sukeshinik5]

Sure, thank you so much.