I am using Limma for downstream analysis of mass spectrometry proteomic data.

I am comparing differential protein expression between healthy controls (n= 40) and two disease subtypes (disease_1 n= 59, disease_2 n= 41). I have an 'intensities' matrix with the log2 transformed LFQ expression for 1300 proteins in rows with individual samples in columns. I also have a patient variables matrix with the columns including sample name, disease subtype, batch, age & sex. I wish to adjust for batch (A-E), as well as sex and age.

I read the Limma user guide and searched the Bioconductor support pages- both of which I found very helpful.

From my understanding, it is most correct to include batch adjustment (and other covariants?) in the linear model. Therefore I have designed my code as outlined below.

My question: Is it possible to extract a fitted expression matrix from the lmfit function in limma (post batch adjustment) which can be used for PCA?

Or is the only possible to get a fitted expression matrix for PCA by independently using the remove::batch function? I have also included this code.

Thank you.

 design<-model.matrix(~patientvar$Subtype+patientvar$Batch+ patientvar$Sex +patientvar$Age)
    fit<-lmFit(intensities,design)
    fit2 <- eBayes(fit)
    topTable(fit2, coef=2)
    topTable(fit2, coef=3)

    design0 = model.matrix(~patientvar$Subtype)
    rembatch_intensities <- removeBatchEffect(intensities, batch=patientvar$Batch,
    batch2 = patientvar$Sex, covariates = patientvar$Age, design = design0)

> sessioninfo()
 setting  value
 version  R version 4.3.0 (2023-04-21 ucrt)
 os       Windows 11 x64 (build 22621)
 system   x86_64, mingw32
 ui       RStudio
 language (EN)
 collate  English_Ireland.utf8
 ctype    English_Ireland.utf8
 tz       America/Los_Angeles
 date     2023-05-26
 rstudio  2023.03.1+446 Cherry Blossom (desktop)
 pandoc   2.19.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

package     * version date (UTC) lib source
 BiocManager   1.30.20 2023-02-24 [1] CRAN (R 4.3.0)
 cli           3.6.1   2023-03-23 [1] CRAN (R 4.3.0)
 digest        0.6.31  2022-12-11 [1] CRAN (R 4.3.0)
 evaluate      0.20    2023-01-17 [1] CRAN (R 4.3.0)
 fastmap       1.1.1   2023-02-24 [1] CRAN (R 4.3.0)
 htmltools     0.5.5   2023-03-23 [1] CRAN (R 4.3.0)
 knitr         1.42    2023-01-25 [1] CRAN (R 4.3.0)
 limma       * 3.56.0  2023-04-25 [1] Bioconductor
 rlang         1.1.1   2023-04-28 [1] CRAN (R 4.3.0)
 rmarkdown     2.21    2023-03-26 [1] CRAN (R 4.3.0)
 rstudioapi    0.14    2022-08-22 [1] CRAN (R 4.3.0)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.3.0)
 xfun          0.39    2023-04-20 [1] CRAN (R 4.3.0)
 yaml          2.3.7   2023-01-23 [1] CRAN (R 4.3.0)

 [1] C:/Users/sarah/AppData/Local/R/win-library/4.3
 [2] C:/Program Files/R/R-4.3.0/library

Using removeBatchEffect() is the only way to get batch-corrected expression values for plotting purposes. removeBatchEffect() and lmFit() fit the same linear model but lmFit doesn't know which factors are batch effects or covariates and which are of primary interest.