代谢组数据分析（十七）：基于structToolbox代谢组数据分析全流程讲解

    library(tidyverse)
    library(patchwork)
    library(cowplot)
    
    # BiocManager::install("structToolbox")
    # BiocManager::install(c('pmp', 'ropls', 'BiocFileCache'))
    library(structToolbox)
    library(pmp)
    library(ropls)
    library(BiocFileCache)
    library(gridExtra)
    library(openxlsx)
    
    # rm(list = ls())
    options(stringsAsFactors = F)
    options(future.globals.maxSize = 1000 * 1024^2)

    profile <- read.csv("./data/MTBLS79_metabolites.csv", row.names = 1)
    
    head(profile[, 1:6])
    
    metadata <- read.csv("./data/MTBLS79_metadata.csv", row.names = 1)

    data_frame <- as.data.frame(profile %>% t())
    sample_meta_frame <- metadata
    
    # 使用虚拟variable数据框
    variable_meta_frame <- data.frame(
      idx = 1:ncol(data_frame),
      row.names = colnames(data_frame)
    )
    
    # 构建DE数据对象
    DE <- struct::DatasetExperiment(
      data = data_frame,
      sample_meta = sample_meta_frame,
      variable_meta = variable_meta_frame,
      name = "MTBLS79",
      description = "Converted datasets into DE object")
    
    # 添加样本顺序
    DE$sample_meta$run_order <- 1:nrow(DE)
    
    # 添加样本分组标签
    Type <- as.character(DE$sample_meta$Class)
    Type[Type != 'QC'] <- 'Sample'
    DE$sample_meta$Type <- factor(Type)
    
    # 添加批次分组标签
    DE$sample_meta$batch_qc <- DE$sample_meta$Batch
    DE$sample_meta$batch_qc[DE$sample_meta$Type=='QC'] <- 'QC'
    
    # 离散型变量转换成因子变量
    DE$sample_meta$Batch <- factor(DE$sample_meta$Batch)
    DE$sample_meta$Type <- factor(DE$sample_meta$Type)
    DE$sample_meta$Class <- factor(DE$sample_meta$Class)
    DE$sample_meta$batch_qc <- factor(DE$sample_meta$batch_qc)
    
    DE

    A "DatasetExperiment" object
    ----------------------------
    name:          MTBLS79
    description:   Converted datasets into DE object
    data:          172 rows x 2488 columns
    sample_meta:   172 rows x 7 columns
    variable_meta: 2488 rows x 1 columns

    # metadata <- DE$sample_meta
    # 
    # profile <- DE$data

    # batch correction: 设置批次校正的参数
    M <- structToolbox::sb_corr(
      order_col = 'run_order',
      batch_col = 'Batch', 
      qc_col = 'Type', 
      qc_label = 'QC',
      spar_lim = c(0.6,0.8))
    
    M <- model_apply(M, DE)
    
    M

    The number of NA and <= 0 values in peaksData before QC-RSC: 18222
    A "sb_corr" object
    ------------------
    name:          Signal/batch correction for mass spectrometry data
    description:   Applies Quality Control Robust Spline (QC-RSC) method to correct for signal drift and batch
                 differences in mass spectrometry data.
    input params:  order_col, batch_col, qc_col, smooth, use_log, min_qc, qc_label, spar_lim 
    outputs:       corrected, fitted 
    predicted:     corrected
    seq_in:        data

    C <- structToolbox::feature_profile(
      run_order = 'run_order',
      qc_label = 'QC',
      qc_column = 'Type',
      colour_by = 'batch_qc',
      feature_to_plot = '200.03196',
      plot_sd = FALSE)
    
    plot_grid(
      (structToolbox::chart_plot(C, M, DE) + ylab('Peak area') + ggtitle('Before')), 
      (structToolbox::chart_plot(C, struct::predicted(M)) + ylab('Peak area') + ggtitle('After')),
      ncol = 2, align = "v"          
      )

    M2 <- structToolbox::filter_na_count(
      threshold = 3,
      factor_name = 'Batch')
    
    M2 <- structToolbox::model_apply(M2, struct::predicted(M))
    
    # QC-RSC未能校正的代谢物数目
    nc <- ncol(DE) - ncol(struct::predicted(M2))
    
    cat(paste0('Number of features removed: ', nc))

    M3 <- structToolbox::kw_rank_sum(
      alpha = 0.0001,
      mtc = 'none',
      factor_names = 'Batch',
      predicted = 'significant') +
      structToolbox::filter_by_name(
    mode = 'exclude',
    dimension = 'variable',
    seq_in = 'names', 
    names = 'seq_fcn', 
    seq_fcn = function(x){return(x[, 1])})
    
    M3 <- structToolbox::model_apply(M3, struct::predicted(M2))
    
    nc <- ncol(struct::predicted(M2)) - ncol(struct::predicted(M3))
    
    cat(paste0('Number of features removed: ', nc))

    M4 <- structToolbox::wilcox_test(
      alpha = 1e-14,
      factor_names = 'Type', 
      mtc = 'none', 
      predicted = 'significant') +
      structToolbox::filter_by_name(
    mode = 'exclude',
    dimension = 'variable',
    seq_in = 'names', 
    names = 'place_holder',
    seq_fcn = function(x){return(x$significant)})
    
    M4 <- structToolbox::model_apply(M4, struct::predicted(M3))
    
    # M4[2]@filtered$sample_meta
    
    nc <- ncol(struct::predicted(M3)) - ncol(struct::predicted(M4))
    cat(paste0('Number of features removed: ', nc))

    M5 <- structToolbox::rsd_filter(
      rsd_threshold = 20,
      factor_name = 'Type')
    
    M5 <- structToolbox::model_apply(M5, struct::predicted(M4))
    
    nc <- ncol(struct::predicted(M4)) - ncol(struct::predicted(M5))
    cat(paste0('Number of features removed: ', nc))

    M6 <- structToolbox::pqn_norm(
    qc_label = 'QC',
    factor_name = 'Type') + 
      structToolbox::knn_impute(neighbours = 5) +
      structToolbox::glog_transform(
    qc_label = 'QC',
    factor_name = 'Type')
    
    M6 <- structToolbox::model_apply(M6, struct::predicted(M5))

    # PCA
    M7 <- structToolbox::mean_centre() + 
      structToolbox::PCA(number_components = 2)
    M7 <- structToolbox::model_apply(M7, struct::predicted(M6))
    
    # plot pca scores
    C <- structToolbox::pca_scores_plot(
      factor_name = c('Sample_Rep','Class'),
      ellipse = 'none')
    
    C2 <- structToolbox::pca_scores_plot(
      factor_name = c('Batch'),
      ellipse = 'none')
    
    plot_grid(
      (structToolbox::chart_plot(C, M7[2]) + coord_fixed() + guides(colour = FALSE)), 
      (structToolbox::chart_plot(C2, M7[2]) + coord_fixed()),
      ncol = 2, align = "v"          
      )

    # PCA
    M7 <- structToolbox::mean_centre() + 
      structToolbox::PCA(number_components = 6)
    M7 <- structToolbox::model_apply(M7, struct::predicted(M6))
    
    C <- structToolbox::pca_scree_plot()
    g1 <- structToolbox::chart_plot(C, M7[2])
    
    C <- structToolbox::pca_loadings_plot()
    g2 <- structToolbox::chart_plot(C, M7[2])
    
    C <- structToolbox::pca_dstat_plot(alpha = 0.95)
    g3 <- structToolbox::chart_plot(C, M7[2])
    
    p1 <- plot_grid(plotlist = list(g1, g2), align = 'h', nrow = 1, axis = 'b')
    p2 <- plot_grid(plotlist = list(g3), nrow = 1)
    
    plot_grid(p1 ,p2, nrow = 2, align = "v")

    QC_paras <- structToolbox::filter_by_name(
      mode = "exclude",
      dimension = "sample",
      names = grep("QC", rownames(struct::predicted(M6)$sample_meta), value = T))
    DE_new <- structToolbox::model_apply(QC_paras, struct::predicted(M6))
    
    
    # QC_paras <- structToolbox::filter_smeta(
    #   mode = 'include',
    #   factor_name = 'Class',
    #   levels = c('C', 'S'))
    # DE_new <- structToolbox::model_apply(QC_paras, struct::predicted(M6))
    
    DE_new

    A "filter_by_name" object
    -------------------------
    name:          Filter by name
    description:   Filter samples/variables by row/column name, index or logicals.
    input params:  mode, dimension, names 
    outputs:       filtered 
    predicted:     filtered
    seq_in:        data

    M <- structToolbox::autoscale() + 
      structToolbox::PLSDA(factor_name = 'Class')
    
    M <- structToolbox::model_apply(M, struct::predicted(DE_new))
    
    C <- structToolbox::pls_scores_plot(factor_name = 'Class')
    
    structToolbox::chart_plot(C, M[2])

    # convert gender to numeric
    DE_filter <- struct::predicted(DE_new)
    
    DE_filter$sample_meta$Class <- as.numeric(DE_filter$sample_meta$Class)
    
    # models sequence
    M <- structToolbox::autoscale(mode = 'both') + 
      structToolbox::PLSR(
    factor_name = 'Class', 
    number_components=3)
    
    M <- structToolbox::model_apply(M, DE_filter)
    
    # some diagnostic charts
    C <- structToolbox::plsr_cook_dist()
    g1 <- structToolbox::chart_plot(C, M[2])
    
    C <- structToolbox::plsr_prediction_plot()
    g2 <- structToolbox::chart_plot(C, M[2])
    
    C <- structToolbox::plsr_qq_plot()
    g3 <- structToolbox::chart_plot(C, M[2])
    
    C <- structToolbox::plsr_residual_hist()
    g4 <- structToolbox::chart_plot(C, M[2])
    
    plot_grid(plotlist = list(g1, g2, g3, g4), 
          nrow=2, align='vh')

    # model sequence
    M <- structToolbox::kfold_xval(
      folds = 5, 
      factor_name = 'Class') * 
      (structToolbox::autoscale(mode = 'both') + 
     structToolbox::PLSR(factor_name = 'Class'))
    M <- structToolbox::run(M, DE_filter, structToolbox::r_squared())
    
    print(paste0("Training set R2: ", paste(round(M@metric.train, 4), collapse = ";")))
    
    print(paste0("Test set Q2: ", round(M@metric.test, 4)))

    # reset gender to original factor
    DE_filter$sample_meta$Class <- droplevels(struct::predicted(DE_new)$sample_meta$Class)
    # model sequence
    M <- structToolbox::permutation_test(
      number_of_permutations = 10, 
      factor_name = 'Class') *
    structToolbox::kfold_xval(
      folds = 5,
      factor_name = 'Class') *
    (structToolbox::autoscale() + 
       structToolbox::PLSDA(
         factor_name = 'Class',
         number_components = 3))
    
    M <- structToolbox::run(M, DE_filter, structToolbox::balanced_accuracy())
    
    C <- permutation_test_plot(style = 'boxplot')
    
    chart_plot(C, M) + ylab('1 - balanced accuracy')

    # prepare model
    TT <- structToolbox::filter_smeta(
      mode = 'include',
      factor_name = 'Class',
      levels = c('C', 'S')) +  
      structToolbox::ttest(
    alpha = 0.05,
    mtc = 'fdr',
    factor_names = 'Class')
    
    # apply model
    TT <- structToolbox::model_apply(TT, struct::predicted(M6))
    
    # keep the data filtered by group for later
    filtered <- struct::predicted(TT[1])
    
    # convert to data frame
    out <- structToolbox::as_data_frame(TT[2])
    
    # show first few features
    head(out)

    QC_paras <- structToolbox::filter_smeta(
      mode = 'include',
      factor_name = 'Class',
      levels = c('C', 'S'))
    DE <- structToolbox::model_apply(QC_paras, struct::predicted(M6))
    
    DE_filter <- struct::predicted(DE)
    DE_filter$sample_meta$Class_num <- as.numeric(DE_filter$sample_meta$Class)
    
    # prepare model
    set.seed(123)
    M <- structToolbox::stratified_split(
      p_train = 0.75,
      factor_name = 'Class')
    # apply to filtered data
    M <- structToolbox::model_apply(M, DE_filter)
    
    # get data from object
    train <- M$training
    train
    
    cat('\n')
    test <- M$testing
    test

    A "DatasetExperiment" object
    ----------------------------
    name:          (Training set)
    description:   A subset of the data has been selected as a training set
    data:          100 rows x 1579 columns
    sample_meta:   100 rows x 8 columns
    variable_meta: 1579 rows x 1 columns
    
    A "DatasetExperiment" object
    ----------------------------
    name:          (Testing set)
    description:   A subset of the data has been selected as a test set
    data:          34 rows x 1579 columns
    sample_meta:   34 rows x 8 columns
    variable_meta: 1579 rows x 1 columns

    # prepare model sequence
    MS <- structToolbox::grid_search_1d(
      param_to_optimise = 'number_components',
      search_values = as.numeric(c(1:6)),
      model_index = 2,
      factor_name = 'Class_num',
      max_min = 'max') *
      structToolbox::permute_sample_order(
    number_of_permutations = 10) *
      structToolbox::kfold_xval(
    folds = 5,
    factor_name = 'Class_num') * 
      (structToolbox::mean_centre(mode = 'sample_meta') + 
     structToolbox::PLSR(factor_name = 'Class_num'))
    
    # run the validation
    MS <- struct::run(MS, train, structToolbox::r_squared())
    
    # plot
    C <- structToolbox::gs_line()
    
    structToolbox::chart_plot(C, MS)

    # prepare the discriminant model
    P <- structToolbox::PLSDA(
      number_components = 3, 
      factor_name = 'Class')
    
    # apply the model
    fit <- structToolbox::model_apply(P, train)
    
    # charts
    C <- structToolbox::plsda_predicted_plot(
      factor_name = 'Class',
      style = 'boxplot')
    g1 <- structToolbox::chart_plot(C, fit)
    
    C <- structToolbox::plsda_predicted_plot(
      factor_name = 'Class',
      style = 'density')
    g2 <- structToolbox::chart_plot(C, fit) + 
      xlim(c(-2, 2))
    
    C <- structToolbox::plsda_roc_plot(factor_name = 'Class')
    g3 <- structToolbox::chart_plot(C, fit)
    
    plot_grid(g1, g2, g3,
          align = 'vh', axis = 'tblr',
          nrow = 1, labels = c('A', 'B', 'C'))

    # AUC for comparison 
    MET <- structToolbox::calculate(
      structToolbox::AUC(), 
      fit$y$Class,
      fit$yhat[, 1])
    
    MET

    A "AUC" object
    --------------
    name:          Area under ROC curve
    description:   The area under the ROC curve of a classifier is estimated using the trapezoid method.
    value:         1

    # model sequence
    MS <- structToolbox::permutation_test(
      number_of_permutations = 20,
      factor_name = 'Class_num') * 
      structToolbox::kfold_xval(
    folds = 5,
    factor_name = 'Class_num') *
      (structToolbox::mean_centre(mode = 'sample_meta') + 
     structToolbox::PLSR(factor_name = 'Class_num', number_components = 3))
    
    # run iterator
    MS <- struct::run(MS, train, structToolbox::r_squared())
    
    # chart
    C <- structToolbox::permutation_test_plot(style = 'density') 
    
    structToolbox::chart_plot(C, MS) + xlim(c(-1, 1)) + xlab('R Squared')

    # prepare the discriminant model
    P <- structToolbox::PLSDA(
      number_components = 2, 
      factor_name = 'Class')
    
    # apply the model
    P <- structToolbox::model_apply(P, train)
    
    C <- structToolbox::pls_scores_plot(components = c(1, 2), factor_name = 'Class')
    structToolbox::chart_plot(C, P)

    # prepare chart
    C <- structToolbox::pls_vip_plot(ycol = 'C')
    g1 <- structToolbox::chart_plot(C, P)
    
    C <- structToolbox::pls_regcoeff_plot(ycol = 'C')
    g2 <- structToolbox::chart_plot(C, P)
    
    plot_grid(g1, g2, align = 'hv', axis = 'tblr', nrow = 2)

    QC_paras <- structToolbox::filter_smeta(
      mode = 'include',
      factor_name = 'Class',
      levels = c('C', 'S'))
    DE <- structToolbox::model_apply(QC_paras, struct::predicted(M6))
    
    DE_filter <- struct::predicted(DE)
    DE_filter$sample_meta$Class_num <- as.numeric(DE_filter$sample_meta$Class)
    
    # model sequence and pls model (NB data already centred)
    MS <- structToolbox::PLSDA(factor_name = 'Class', number_components = 2)
    
    # apply PLS model
    MS <- structToolbox::model_apply(MS, DE_filter)
    
    # plot the data
    C <- structToolbox::pls_scores_plot(factor_name = 'Class')
    
    structToolbox::chart_plot(C, MS)
    
    # new DatasetExperiment object from the PLS scores
    DE2 <- struct::DatasetExperiment(
      data = MS$scores$data, 
      sample_meta = DE_filter$sample_meta,
      variable_meta = data.frame('LV' = c(1:2), row.names = colnames(MS$scores)),
      name = 'Illustrativate SVM dataset',
      description = 'Generated by applying PLS to the processed Metabolites data')
    
    DE2

    A "DatasetExperiment" object
    ----------------------------
    name:          Illustrativate SVM dataset
    description:   Generated by applying PLS to the processed Metabolites data
    data:          134 rows x 2 columns
    sample_meta:   134 rows x 8 columns
    variable_meta: 2 rows x 1 columns

    # SVM model
    M <- structToolbox::SVM(
      factor_name = 'Class',
      kernel = 'linear')
    
    # apply model
    M <- structToolbox::model_apply(M, DE2)
    
    # low cost
    M$cost <- 0.01
    M <- structToolbox::model_apply(M, DE2)
    C <- structToolbox::svm_plot_2d(factor_name = 'Class')
    g1 <- structToolbox::chart_plot(C, M, DE2)
    
    # medium cost
    M$cost <- 0.05
    M <- structToolbox::model_apply(M, DE2)
    C <- structToolbox::svm_plot_2d(factor_name = 'Class')
    g2 <- structToolbox::chart_plot(C, M, DE2)
    
    # high cost
    M$cost <- 100
    M <- structToolbox::model_apply(M, DE2)
    C <- structToolbox::svm_plot_2d(factor_name = 'Class')
    g3 <- structToolbox::chart_plot(C, M, DE2)
    
    # plot
    prow <- plot_grid(
      g1 + theme(legend.position = "none"),
      g2 + theme(legend.position = "none"),
      g3 + theme(legend.position = "none"),
      align = 'vh',
      labels = c("Low cost", "Medium cost", "High cost"),
      hjust = -1,
      nrow = 2
    )
    
    legend <- get_legend(
      # create some space to the left of the legend
      g1 + guides(color = guide_legend(nrow = 1)) +
      theme(legend.position = "bottom")
    )
    
    plot_grid(prow, legend, ncol=1, rel_heights = c(1, .1))

    R version 4.3.3 (2024-02-29)
    Platform: aarch64-apple-darwin20 (64-bit)
    Running under: macOS Sonoma 14.2
    
    Matrix products: default
    BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
    LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
    
    locale:
    [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    
    time zone: Asia/Shanghai
    tzcode source: internal
    
    attached base packages:
    [1] stats     graphics  grDevices utils     datasets  methods   base     
    
    other attached packages:
     [1] structToolbox_1.17.3 struct_1.14.1        cowplot_1.1.3        patchwork_1.2.0      lubridate_1.9.3      forcats_1.0.0       
     [7] stringr_1.5.1        dplyr_1.1.4          purrr_1.0.2          readr_2.1.5          tidyr_1.3.1          tibble_3.2.1        
    [13] ggplot2_3.5.1        tidyverse_2.0.0     
    
    loaded via a namespace (and not attached):
     [1] SummarizedExperiment_1.32.0 gtable_0.3.4                xfun_0.43                   Biobase_2.62.0             
     [5] lattice_0.22-6              tzdb_0.4.0                  vctrs_0.6.5                 tools_4.3.3                
     [9] bitops_1.0-7                generics_0.1.3              stats4_4.3.3                fansi_1.0.6                
    [13] pkgconfig_2.0.3             Matrix_1.6-5                S4Vectors_0.40.2            lifecycle_1.0.4            
    [17] GenomeInfoDbData_1.2.11     compiler_4.3.3              munsell_0.5.0               GenomeInfoDb_1.38.8        
    [21] RCurl_1.98-1.14             pillar_1.9.0                crayon_1.5.2                DelayedArray_0.28.0        
    [25] abind_1.4-5                 tidyselect_1.2.1            stringi_1.8.3               ggthemes_5.1.0             
    [29] grid_4.3.3                  colorspace_2.1-0            cli_3.6.3                   SparseArray_1.2.4          
    [33] magrittr_2.0.3              S4Arrays_1.2.1              utf8_1.2.4                  withr_3.0.0                
    [37] scales_1.3.0                sp_2.1-3                    timechange_0.3.0            XVector_0.42.0             
    [41] matrixStats_1.2.0           gridExtra_2.3               hms_1.1.3                   knitr_1.45                 
    [45] GenomicRanges_1.54.1        IRanges_2.36.0              rlang_1.1.4                 ontologyIndex_2.12         
    [49] glue_1.7.0                  BiocGenerics_0.48.1         rstudioapi_0.16.0           R6_2.5.1                   
    [53] MatrixGenerics_1.14.0       zlibbioc_1.48.2