removeBatchEffect with non-linear model fit
1
0
Entering edit mode
@2289c15f
Last seen 10 weeks ago
Germany

Hello, I am attempting to use limma's removeBatchEffect for visualization purposes (heatmat & PCA) while fitting non-linear models (splines) to my expression data in DESeq2. Given that my design is balanced, would this approach work within the framework of this function? The results _looks_ correct, but I want to be sure as I've never seen it used like that.

My code:

mat <- assay(vsd)
mm <- model.matrix(~ ns(age_scaled, df = 3), colData(vsd))
mat <- limma::removeBatchEffect(mat, batch = vsd$genotype, design = mm)
assay(vsd) <- mat

Thanks a lot in advance.

sessionInfo( )
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] splines   stats4    stats     graphics  grDevices utils     datasets  methods  
[9] base     

other attached packages:
 [1] edgeR_3.42.4                limma_3.56.2               
 [3] ReactomePA_1.44.0           DOSE_3.26.2                
 [5] enrichplot_1.20.3           clusterProfiler_4.8.3      
 [7] shiny_1.8.0                 xlsx_0.6.5                 
 [9] readxl_1.4.3                EnhancedVolcano_1.18.0     
[11] ggrepel_0.9.3               RColorBrewer_1.1-3         
[13] pheatmap_1.0.12             DESeq2_1.40.2              
[15] SummarizedExperiment_1.30.2 Biobase_2.60.0             
[17] MatrixGenerics_1.12.3       matrixStats_1.0.0          
[19] GenomicRanges_1.52.0        GenomeInfoDb_1.36.4        
[21] IRanges_2.34.1              S4Vectors_0.38.1           
[23] BiocGenerics_0.46.0         lubridate_1.9.3            
[25] forcats_1.0.0               stringr_1.5.1              
[27] dplyr_1.1.2                 purrr_1.0.2                
[29] readr_2.1.4                 tidyr_1.3.0                
[31] tibble_3.2.1                ggplot2_3.4.4              
[33] tidyverse_2.0.0             BiocManager_1.30.22        

loaded via a namespace (and not attached):
  [1] rstudioapi_0.15.0       jsonlite_1.8.7          magrittr_2.0.3         
  [4] farver_2.1.1            rmarkdown_2.25          fs_1.6.3               
  [7] zlibbioc_1.46.0         vctrs_0.6.3             memoise_2.0.1          
 [10] RCurl_1.98-1.12         ggtree_3.8.2            htmltools_0.5.7        
 [13] S4Arrays_1.0.5          cellranger_1.1.0        gridGraphics_0.5-1     
 [16] plyr_1.8.9              cachem_1.0.8            igraph_1.5.1           
 [19] mime_0.12               lifecycle_1.0.4         pkgconfig_2.0.3        
 [22] Matrix_1.6-3            R6_2.5.1                fastmap_1.1.1          
 [25] gson_0.1.0              GenomeInfoDbData_1.2.10 digest_0.6.33          
 [28] aplot_0.2.2             colorspace_2.1-0        patchwork_1.1.3        
 [31] AnnotationDbi_1.62.2    RSQLite_2.3.1           timechange_0.2.0       
 [34] fansi_1.0.4             httr_1.4.7              polyclip_1.10-6        
 [37] abind_1.4-5             compiler_4.3.1          bit64_4.0.5            
 [40] withr_2.5.2             downloader_0.4          graphite_1.46.0        
 [43] BiocParallel_1.34.2     viridis_0.6.4           DBI_1.1.3              
 [46] ggforce_0.4.1           MASS_7.3-60             rappdirs_0.3.3         
 [49] DelayedArray_0.26.7     HDO.db_0.99.1           tools_4.3.1            
 [52] ape_5.7-1               scatterpie_0.2.1        httpuv_1.6.11          
 [55] glue_1.6.2              promises_1.2.1          nlme_3.1-162           
 [58] GOSemSim_2.26.1         grid_4.3.1              shadowtext_0.1.2       
 [61] reshape2_1.4.4          fgsea_1.26.0            generics_0.1.3         
 [64] gtable_0.3.4            tzdb_0.4.0              hms_1.1.3              
 [67] data.table_1.14.8       tidygraph_1.2.3         utf8_1.2.3             
 [70] XVector_0.40.0          pillar_1.9.0            yulab.utils_0.1.0      
 [73] later_1.3.1             rJava_1.0-6             tweenr_2.0.2           
 [76] treeio_1.24.3           lattice_0.21-8          bit_4.0.5              
 [79] tidyselect_1.2.0        GO.db_3.17.0            locfit_1.5-9.8         
 [82] Biostrings_2.68.1       reactome.db_1.84.0      knitr_1.45             
 [85] gridExtra_2.3           xfun_0.40               graphlayouts_1.0.1     
 [88] stringi_1.7.12          lazyeval_0.2.2          ggfun_0.1.3            
 [91] yaml_2.3.7              xlsxjars_0.6.1          evaluate_0.23          
 [94] codetools_0.2-19        ggraph_2.1.0            qvalue_2.32.0          
 [97] graph_1.78.0            ggplotify_0.1.2         cli_3.6.1              
[100] xtable_1.8-4            munsell_0.5.0           Rcpp_1.0.11            
[103] png_0.1-8               parallel_4.3.1          ellipsis_0.3.2         
[106] blob_1.2.4              bitops_1.0-7            viridisLite_0.4.2      
[109] tidytree_0.4.5          scales_1.2.1            crayon_1.5.2           
[112] rlang_1.1.1             cowplot_1.1.1           fastmatch_1.1-4        
[115] KEGGREST_1.40.1
limma DESeq2 • 294 views
ADD COMMENT
2
Entering edit mode
@gordon-smyth
Last seen 5 hours ago
WEHI, Melbourne, Australia

Assuming that the full model is splines + genotype, then you are using removeBatchEffect correctly.

ADD COMMENT

Login before adding your answer.

Traffic: 508 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6