前言

批次校正是数据预处理常见的步骤，在整合多批次数据的时候，需要考虑到这些非生物学因素对treatment的影响。更多知识分享请到 https://zouhua.top/。

加载R包

  1. library(dplyr)
  2. library(tibble)
  3. library(Biobase)
  4. library(limma)
  5. library(ggplot2)
  7. # rm(list = ls())
  8. options(stringsAsFactors = F)
  9. options(future.globals.maxSize = 1000 * 1024^2)
  11. subgrp <- c("HC", "CP", "PDAC")
  12. grp.col <- c("#568875", "#73FAFC", "#EE853D")
  14. ExprSet <- readRDS("MS_Protein_ImputeExprSet.RDS")

Function

  1. PCAFun <- function(dataset=ExprSet){
  3.   # dataset=ExprSet
  5.   metadata <- pData(dataset)
  6.   profile <- exprs(dataset)
  8.   # pca 
  9.   pca <- prcomp(scale(t(profile), center = T, scale = T))
  10.   require(factoextra)
  11.   eig <- get_eig(pca)
  12.   # explains variable 
  13.   explains <- paste0(paste0("PC", seq(2)), "(", paste0(round(eig[1:2, 2], 2), "%"), ")")
  14.   # principal component score of each sample
  15.   score <- inner_join(pca$x %>% data.frame() %>% 
  16.                         dplyr::select(c(1:2)) %>% 
  17.                         rownames_to_column("SampleID"), 
  18.                       metadata %>% rownames_to_column("SampleID"),
  19.                       by = "SampleID") %>%
  20.     mutate(SubGroup=factor(SubGroup, levels = subgrp))
  23.   # PERMANOVA
  24.   require(vegan)
  25.   set.seed(123)
  26.   if(any(profile < 0)){
  27.     res_adonis <- adonis(vegdist(t(profile), method = "manhattan") ~ metadata$SubGroup, permutations = 999) 
  28.   }else{
  29.     res_adonis <- adonis(vegdist(t(profile), method = "bray") ~ metadata$SubGroup, permutations = 999)    
  30.   }
  31.   adn_pvalue <- res_adonis[[1]][["Pr(>F)"]][1]
  32.   adn_rsquared <- round(res_adonis[[1]][["R2"]][1],3)
  33.   signi_label <- paste(cut(adn_pvalue,breaks=c(-Inf, 0.001, 0.01, 0.05, Inf), label=c("***", "**", "*", ".")))
  34.   adn_res_format <- bquote(atop(atop("PERMANOVA",R^2==~.(adn_rsquared)),
  35.                                   atop("p-value="~.(adn_pvalue)~.(signi_label), phantom()))) 
  37.   pl <- ggplot(score, aes(x=PC1, y=PC2))+
  38.               #geom_point(aes(fill=SubGroup), size=3.5, shape=21, stroke = .8, color = "black")+
  39.               geom_point(aes(color=SubGroup, shape=Omics), size=2.5)+
  40.               stat_ellipse(aes(color=SubGroup), level = 0.95, linetype = 1, size = 1.5)+
  41.               labs(x=explains[1], y=explains[2])+
  42.               scale_color_manual(values = grp.col)+
  43.               # scale_fill_manual(name = "Condition", 
  44.               #                   values = grp.col)+
  45.               scale_shape_manual(name = "Batch", 
  46.                                  values = c(15, 16, 17))+    
  47.               annotate("text", x = max(score$PC1) - 8,
  48.                        y = min(score$PC1),
  49.                        label = adn_res_format,
  50.                        size = 6)+ 
  51.               #guides(color="none")+
  52.               theme_bw()+
  53.               theme(axis.title = element_text(size = 10, color = "black", face = "bold"),
  54.                     axis.text = element_text(size = 9, color = "black"),
  55.                     text = element_text(size = 8, color = "black", family = "serif"),
  56.                     strip.text = element_text(size = 9, color = "black", face = "bold"), 
  57.                     panel.grid = element_blank(),
  58.                     legend.title = element_text(size = 11, color = "black", family = "serif"),
  59.                     legend.text = element_text(size = 10, color = "black", family = "serif"),
  60.                     legend.position = c(0, 0),
  61.                     legend.justification = c(0, 0),
  62.                     legend.background = element_rect(color = "black", fill = "white", linetype = 2, size = 0.5))
  63.   return(pl)
  64. }
  66. RemoveBatchExprSet <- function(datset=ExprSet){
  68.   # datset=ExprSet
  70.   pheno <- pData(datset)
  71.   expr <- exprs(datset)
  72.   feature <- fData(datset)
  74.   # Remove batch effects
  75.   phen <- pheno %>% arrange(Omics)
  76.   prof <- data.frame(expr) %>% dplyr::select(rownames(phen))
  77.   if(!any(rownames(phen) == colnames(prof))){
  78.     stop("The Order of SampleID in DiscoverySet is wrong")
  79.   }
  81.   # CA19-9 should be removed batch effects
  82.   prof_CA199 <- prof[rownames(prof) == "CA199", , F]
  83.   prof_NoCA199 <- prof[rownames(prof) != "CA199", , F]
  84.   prof_rbe <- removeBatchEffect(as.matrix(prof_NoCA199), batch = phen$Omics)
  85.   prof_merge <- rbind(prof_rbe, prof_CA199)
  87.   fdf <- new("AnnotatedDataFrame", data=feature)
  88.   exprs_dis <- as.matrix(prof_merge)
  89.   adf_dis <-  new("AnnotatedDataFrame", data=phen)
  90.   experimentData <- new("MIAME",
  91.           name="Hua", lab="Lab",
  92.           title="tumor Experiment",
  93.           abstract="The Mass Spectrometry ExpressionSet with overlap samples",
  94.           url="www.zouhua.top",
  95.           other=list(notes="The Intersect Proteins and Samples"))
  96.   res <- new("ExpressionSet", 
  97.                 exprs=exprs_dis,
  98.                 phenoData=adf_dis,
  99.                 featureData=fdf,
  100.                 experimentData=experimentData)
  102.   return(res)  
  103. }
  106. RemoveBatchExprSet_Design <- function(datset=ExprSet){
  108.   # datset=ExprSet
  110.   pheno <- pData(datset)
  111.   expr <- exprs(datset)
  112.   feature <- fData(datset)
  114.   # design
  115.   Design <- model.matrix(~SubGroup, data = pheno)
  117.   # Remove batch effects
  118.   phen <- pheno %>% arrange(Omics)
  119.   prof <- data.frame(expr) %>% dplyr::select(rownames(phen))
  120.   if(!any(rownames(phen) == colnames(prof))){
  121.     stop("The Order of SampleID in DiscoverySet is wrong")
  122.   }
  124.   # CA19-9 should be removed batch effects
  125.   prof_CA199 <- prof[rownames(prof) == "CA199", , F]
  126.   prof_NoCA199 <- prof[rownames(prof) != "CA199", , F]
  127.   prof_rbe <- removeBatchEffect(as.matrix(prof_NoCA199), batch = phen$Omics, design = Design)
  128.   prof_merge <- rbind(prof_rbe, prof_CA199)
  130.   fdf <- new("AnnotatedDataFrame", data=feature)
  131.   exprs_dis <- as.matrix(prof_merge)
  132.   adf_dis <-  new("AnnotatedDataFrame", data=phen)
  133.   experimentData <- new("MIAME",
  134.           name="Hua", lab="Lab",
  135.           title="tumor Experiment",
  136.           abstract="The Mass Spectrometry ExpressionSet with overlap samples",
  137.           url="www.zouhua.top",
  138.           other=list(notes="The Intersect Proteins and Samples"))
  139.   res <- new("ExpressionSet", 
  140.                 exprs=exprs_dis,
  141.                 phenoData=adf_dis,
  142.                 featureData=fdf,
  143.                 experimentData=experimentData)
  145.   return(res)  
  146. }
  149. # Batch mean centering
  150. RemoveBatchExprSet_BMC <- function(datset=ExprSet){
  152.   # datset=ExprSet
  154.   pheno <- pData(datset)
  155.   expr <- exprs(datset)
  156.   feature <- fData(datset)
  158.   # Remove batch effects
  159.   phen <- pheno %>% arrange(Omics)
  160.   prof <- data.frame(expr) %>% dplyr::select(rownames(phen))
  161.   # CA19-9 should be removed batch effects
  162.   prof_CA199 <- prof[rownames(prof) == "CA199", , F]
  163.   prof_NoCA199 <- prof[rownames(prof) != "CA199", , F]
  165.   prof_list <- list()
  166.   batchlevels <- unique(phen$Omics)
  167.   for(i in 1:length(batchlevels)){
  168.     prof_list[[i]] <- scale(t(prof_NoCA199 %>% 
  169.                                 dplyr::select(rownames(phen[phen$Omics == batchlevels[i], ]))),
  170.                             center = TRUE, scale = FALSE)
  171.   }
  172.   prof_merge <- rbind(t(do.call(rbind, prof_list)), prof_CA199)
  173.   if(!any(rownames(phen) == colnames(prof_merge))){
  174.     stop("The Order of SampleID in DiscoverySet is wrong")
  175.   }  
  177.   fdf <- new("AnnotatedDataFrame", data=feature)
  178.   exprs_dis <- as.matrix(prof_merge)
  179.   adf_dis <-  new("AnnotatedDataFrame", data=phen)
  180.   experimentData <- new("MIAME",
  181.           name="Hua", lab="Lab",
  182.           title="tumor Experiment",
  183.           abstract="The Mass Spectrometry ExpressionSet with overlap samples",
  184.           url="www.zouhua.top",
  185.           other=list(notes="The Intersect Proteins and Samples"))
  186.   res <- new("ExpressionSet", 
  187.              exprs=exprs_dis,
  188.              phenoData=adf_dis,
  189.              featureData=fdf,
  190.              experimentData=experimentData)
  192.   return(res)  
  193. }
  196. # ComBat
  197. RemoveBatchExprSet_ComBat <- function(datset=ExprSet){
  199.   # datset=ExprSet
  201.   pheno <- pData(datset)
  202.   expr <- exprs(datset)
  203.   feature <- fData(datset)
  205.   # design
  206.   Design <- model.matrix(~SubGroup, data = pheno)  
  208.   # Remove batch effects
  209.   phen <- pheno %>% arrange(Omics)
  210.   prof <- data.frame(expr) %>% dplyr::select(rownames(phen))
  211.   # CA19-9 should be removed batch effects
  212.   prof_CA199 <- prof[rownames(prof) == "CA199", , F]
  213.   prof_NoCA199 <- prof[rownames(prof) != "CA199", , F]
  215.   prof_combat <- sva::ComBat(prof_NoCA199, batch = phen$Omics, mod = Design, par.prior = F, prior.plots = F)
  216.   prof_merge <- rbind(prof_combat, prof_CA199)
  217.   if(!any(rownames(phen) == colnames(prof_merge))){
  218.     stop("The Order of SampleID in DiscoverySet is wrong")
  219.   }  
  221.   fdf <- new("AnnotatedDataFrame", data=feature)
  222.   exprs_dis <- as.matrix(prof_merge)
  223.   adf_dis <-  new("AnnotatedDataFrame", data=phen)
  224.   experimentData <- new("MIAME",
  225.           name="Hua", lab="Lab",
  226.           title="tumor Experiment",
  227.           abstract="The Mass Spectrometry ExpressionSet with overlap samples",
  228.           url="www.zouhua.top",
  229.           other=list(notes="The Intersect Proteins and Samples"))
  230.   res <- new("ExpressionSet", 
  231.              exprs=exprs_dis,
  232.              phenoData=adf_dis,
  233.              featureData=fdf,
  234.              experimentData=experimentData)
  236.   return(res)  
  237. }

Run

  1. ExprSet_RBE <- RemoveBatchExprSet(datset=ExprSet)
  2. ExprSet_RBE_Design <- RemoveBatchExprSet_Design(datset=ExprSet)
  3. ExprSet_RBE_BMC <- RemoveBatchExprSet_BMC(datset=ExprSet)
  4. ExprSet_RBE_ComBat <- RemoveBatchExprSet_ComBat(datset=ExprSet)
  6. cowplot::plot_grid(#PCAFun(dataset=ExprSet), 
  7.                    PCAFun(dataset=ExprSet_RBE),
  8.                    PCAFun(dataset=ExprSet_RBE_Design),
  9.                    PCAFun(dataset=ExprSet_RBE_BMC),
  10.                    PCAFun(dataset=ExprSet_RBE_ComBat),
  11.                    align = "hv", 
  12.                    ncol = 2,
  13.                    labels = c(#"Origin", 
  14.                               "Batch corrected", 
  15.                               "Batch corrected with Design", 
  16.                               "Batch corrected with BMC",
  17.                               "Batch corrected with ComBat"))

systemic information

  1. sessionInfo()
  1. R version 4.0.2 (2020-06-22)
  2. Platform: x86_64-conda_cos6-linux-gnu (64-bit)
  3. Running under: CentOS Linux 8 (Core)
  5. Matrix products: default
  6. BLAS/LAPACK: /disk/share/anaconda3/lib/libopenblasp-r0.3.10.so
  8. locale:
  9.  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
  10.  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
  11.  [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
  13. attached base packages:
  14. [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
  16. other attached packages:
  17.  [1] DEP_1.12.0          gg.gap_1.4          ggpubr_0.4.0        gtsummary_1.4.0     vegan_2.5-6        
  18.  [6] lattice_0.20-45     permute_0.9-5       factoextra_1.0.7    ggrepel_0.9.1.9999  sampling_2.8       
  19. [11] data.table_1.14.2   convert_1.64.0      marray_1.68.0       limma_3.46.0        Biobase_2.50.0     
  20. [16] BiocGenerics_0.36.0 ggplot2_3.3.5       tibble_3.1.6        dplyr_1.0.7        
  22. loaded via a namespace (and not attached):
  23.   [1] utf8_1.2.2                  shinydashboard_0.7.1        reticulate_1.18            
  24.   [4] gmm_1.6-5                   tidyselect_1.1.1            htmlwidgets_1.5.4          
  25.   [7] grid_4.0.2                  BiocParallel_1.24.1         lpSolve_5.6.15             
  26.  [10] norm_1.0-9.5                munsell_0.5.0               codetools_0.2-18           
  27.  [13] preprocessCore_1.52.1       DT_0.20                     umap_0.2.7.0               
  28.  [16] withr_2.4.3                 colorspace_2.0-2            knitr_1.37                 
  29.  [19] rstudioapi_0.13             stats4_4.0.2                robustbase_0.93-6          
  30.  [22] bayesm_3.1-4                ggsignif_0.6.0              mzID_1.28.0                
  31.  [25] MatrixGenerics_1.2.1        labeling_0.4.2              GenomeInfoDbData_1.2.4     
  32.  [28] farver_2.1.0                vctrs_0.3.8                 generics_0.1.1             
  33.  [31] xfun_0.29                   R6_2.5.1                    doParallel_1.0.16          
  34.  [34] GenomeInfoDb_1.26.4         clue_0.3-60                 bitops_1.0-7               
  35.  [37] DelayedArray_0.16.3         assertthat_0.2.1            promises_1.2.0.1           
  36.  [40] scales_1.1.1                gtable_0.3.0                Cairo_1.5-12.2             
  37.  [43] affy_1.68.0                 sandwich_3.0-0              rlang_0.4.12               
  38.  [46] mzR_2.24.1                  GlobalOptions_0.1.2         splines_4.0.2              
  39.  [49] rstatix_0.7.0               impute_1.64.0               broom_0.7.9                
  40.  [52] BiocManager_1.30.16         yaml_2.2.1                  abind_1.4-5                
  41.  [55] crosstalk_1.2.0             backports_1.4.1             httpuv_1.6.4               
  42.  [58] tensorA_0.36.2              tools_4.0.2                 affyio_1.60.0              
  43.  [61] ellipsis_0.3.2              jquerylib_0.1.4             RColorBrewer_1.1-2         
  44.  [64] MSnbase_2.16.1              Rcpp_1.0.7                  plyr_1.8.6                 
  45.  [67] zlibbioc_1.36.0             purrr_0.3.4                 RCurl_1.98-1.3             
  46.  [70] openssl_1.4.6               GetoptLong_1.0.5            cowplot_1.1.0              
  47.  [73] S4Vectors_0.28.1            zoo_1.8-8                   SummarizedExperiment_1.20.0
  48.  [76] haven_2.3.1                 cluster_2.1.0               tinytex_0.36               
  49.  [79] magrittr_2.0.1              RSpectra_0.16-0             openxlsx_4.2.3             
  50.  [82] circlize_0.4.10             pcaMethods_1.80.0           mvtnorm_1.1-1              
  51.  [85] ProtGenerics_1.22.0         matrixStats_0.61.0          xtable_1.8-4               
  52.  [88] hms_1.1.1                   mime_0.12                   evaluate_0.14              
  53.  [91] XML_3.99-0.6                rio_0.5.16                  readxl_1.3.1               
  54.  [94] IRanges_2.24.1              shape_1.4.5                 compiler_4.0.2             
  55.  [97] gt_0.2.2                    ncdf4_1.17                  crayon_1.4.2               
  56. [100] htmltools_0.5.2             mgcv_1.8-38                 later_1.2.0                
  57. [103] tzdb_0.2.0                  tidyr_1.1.4                 DBI_1.1.2                  
  58. [106] ComplexHeatmap_2.6.2        MASS_7.3-54                 tmvtnorm_1.4-10            
  59. [109] broom.helpers_1.3.0         compositions_2.0-2          Matrix_1.4-0               
  60. [112] car_3.0-10                  readr_2.1.1                 cli_3.1.0                  
  61. [115] vsn_3.58.0                  imputeLCMD_2.0              GenomicRanges_1.42.0       
  62. [118] forcats_0.5.0               pkgconfig_2.0.3             foreign_0.8-81             
  63. [121] MALDIquant_1.19.3           foreach_1.5.1               bslib_0.3.1                
  64. [124] XVector_0.30.0              stringr_1.4.0               digest_0.6.29              
  65. [127] rmarkdown_2.11              cellranger_1.1.0            curl_4.3.2                 
  66. [130] shiny_1.7.1                 rjson_0.2.20                lifecycle_1.0.1            
  67. [133] nlme_3.1-150                jsonlite_1.7.2              carData_3.0-4              
  68. [136] askpass_1.1                 fansi_0.5.0                 pillar_1.6.4               
  69. [139] fastmap_1.1.0               DEoptimR_1.0-8              survival_3.2-13            
  70. [142] glue_1.6.0                  zip_2.1.1                   png_0.1-7                  
  71. [145] iterators_1.0.13            stringi_1.4.6               sass_0.4.0

Reference

参考文章如引起任何侵权问题，可以与我联系，谢谢。