R包文档参考：
https://www.metaboanalyst.ca/docs/RTutorial.xhtml#Workflow:%20Case%20Study%201
R包函数参考：
https://rdrr.io/github/xia-lab/MetaboAnalystR3.0/man/
image.png
MetaboAnalystR包与网页同步，可以用于分析原始组学数据和批次分析，建立了一个从原始数据开始分析的优化流程。下面的教程旨在通过提供使用R包的几个最常见任务的逐步指导来补充我们的基于web的功能。
1.1 MetaboAnalystR 3简介
MetaboAnalystR 3包含MetaboAnalyst web服务器下的R函数和库，包括代谢组学数据分析、可视化和功能解释。该包与MetaboAnalyst web服务器同步。在安装和加载包后，用户将能够使用从MetaboAnalyst网站下载的相应的R命令历史从他们的本地计算机上重现相同的结果，从而实现最大的灵活性和可再现性。
version3 旨在通过实现峰值选择的快速参数优化算法来改进当前的全球代谢组学工作流程，并从12种成熟的方法中自动识别最适合的批量效应校正方法。此外，还增加了通过mummichog2 (PMID: 23861661)直接从m/z峰进行功能解释的更多支持，并添加了一种新的基于路径的方法来集成多组学数据。为了演示这一新功能，我们在本教程中对临床IBD样本进行端到端代谢组学数据分析。在这个R包的小插图中提供了更多的案例研究。在这里，我们更愿意提供一个全面的MetaboAnalystR 3.2安装入门教程(更新于2021年12月16日)。
1.2 MetaboAnalystR 安装
略
2 IBD案例研究
炎症性肠道疾病，包括克罗恩病和溃疡性结肠炎，已经影响了世界各地数百万人。Jason L.等人对微生物组在IBD发病机制中的作用进行了长时间多组学研究。本文以这一新型参数优化为例，介绍了面部样本的代谢组学研究。
2.1 原始数据处理
原始数据处理是所有代谢组学研究的首要环节。MetaboAnalystR 3.2支持”.mzXML”、“.mzML "和".CDF”格式。原始格式(例如“.raw "和" . RAW")将很快得到支持。您可以使用 ProteoWizard (PMID: 23051804) 将原始数据转换为支持的格式。优先使用中心化数据。
2.1.1 加载 MetaboAnalystR包
# Load the MetaboAnalystR package
library("MetaboAnalystR")
2.1.2 下载IBD数据集QC数据
下面的步骤将使用质量控制(QC)示例。下载MS数据进行参数优化。达到快速学习的目的，避免整个样本(超过600个样本)长时间运行。在这里，我们只为每个CD组和非IBD组提供5个样本。如果您想重复和验证我们手稿中的结果，请访问IBD MultiOmics数据库，下载完整批数据并使用此管道运行它们。注：数据因网络无法下载
## 1 设置数据存放目录
data_folder_Sample <- "~/Data_IBD"
data_folder_QC <- "~/QC_IBD"
# Use Google API for data downloading here.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".zip")
# Please authorize your google account to access the data
dl <- drive_download(
as_id("10DBpPEWy2cZyXvmlLOIfpwqQYwFplKYK"), path = temp, overwrite = TRUE)
# 2 下载、放置QC数据
out <- unzip(temp, exdir = data_folder_QC)
# Date files for parameters optimization are deposited below
out
# Now, download the small example data for comparison between CD vs. nonIBD
temp <- tempfile(fileext = ".zip")
dl <- drive_download(
as_id("1-wlFUkzEwWX1afRWLJlY_KEJs7BfsZim"), path = temp, overwrite = TRUE)
# 3 下载、放置CD及nonIBD数据
out <- unzip(temp, exdir = data_folder_Sample)
# Date files for normal processing example are deposited below
out
2.1.3 数据核查
在运行数据分析之前，可以使用 PerformDataInspect 检查一般数据结构和信息。如果有一些非常重要的污染物，它将被直接发现。
# Inspect the MS data via a 3D image. "res" are used to specify the resolution for the MS data.
PerformDataInspect(data_folder_QC,res = 50)
X轴是分子量荷质比，Y轴是保留时间，Z轴是信号强度
image.png
2.1.4 ROI（regions of interest 目标区域） 提取
QC 样品数据 ROI 提取现在使用 PerformROIExtraction 或 PerformDataTrimming 进行。数据 ROI 提取有3种模式。默认是标准模拟方法(“ssm_trim”)。在此模式下，三维MS将先从m/z维度进行裁剪，再从RT维度进行裁剪。rt.idx可以用来调整RT数据保留的范围。另外两种模式被命名为 “rt_specific” 和 “mz_specific”，可用于提取(带正值)或去除(带负值)光谱的某些部分。修剪后的MS文件可以用write = T 保存，修剪后的文件将保存为“.mzML”。色谱图可用plot = T绘制。
# QC samples are trimmed with "ssm_trim" strategy.
# RT 保留20%
raw_data <- PerformROIExtraction(data_folder_QC, rt.idx = 0.2, rmConts = FALSE)
2.1.5 初始化参数设定
已嵌入仪器特定的初始参数。目前，最流行的平台包括UPLC-Q/E、UPLC-Q/TOF、UPLC-T/TOF、UPLC-Ion_trap、UPLC-Orbitrap、UPLC-G2S、HPLC-Q/TOF、HPLC-Ion_Trap、HPLC-Orbitrap和HPLC-S/Q。如果这里没有列出您的平台，则可以将“general”选项应用于此设置。“centWave”和“matchedFIlter”模式均已配置用于进一步处理。此外，具体参数也可以单独设置，如果你对参数有自己的看法，此参数优化步骤可跳过。
# Initial platform specific parameters
param_initial <- SetPeakParam(platform = "UPLC-Q/E")
2.1.6 参数优化
采用“DoE”策略进行参数优化。初始参数来自上述优化的参数或来自嵌入式 SetPeakParam 函数的内部默认参数。
# Select relative core according to your work platform
param_optimized <- PerformParamsOptimization(raw_data, param = param_initial, ncore = 8)
# Following results was generated at the Ubuntu 18.04.4
2.1.7 峰值分析
ImportRawMSData 函数读取原始 MS 数据文件，并将其保存为 OnDiskMSnExp 对象。用SetPlotParam 函数集“ Plot = T ”可以输出两个图——总离子色谱图(TIC)提供整个色谱的概况，基峰色谱图( BPC )是基于最丰富信号的光谱更清晰的轮廓。这些图对下游参数的设置很有用。对于希望查看感兴趣的峰值的用户，可以使用PlotEIC 函数生成提取离子色谱图(EIC)。
# Import raw MS data. The "SetPlotParam" parameters can be used to determine plot or not.
rawData <- ImportRawMSData(NULL, data_folder_Sample, plotSettings = SetPlotParam(Plot=FALSE))
PerformPeakProfiling 函数是从XCMS R函数中更新的峰值处理流程，它自动执行峰值检测、对齐和分组。该函数还生成两个诊断图，包括不同样本中峰值总强度的统计信息、保留时间调整图和显示数据清洗和统计分析之前总体样本聚类的PCA图。
# Peak Profiling with optimized parameters
mSet <- PerformPeakProfiling(rawData,param_optimized$best_parameters,
plotSettings = SetPlotParam(Plot = T))
2.1.8 峰注释
PerformPeakAnnotation 函数使用 CAMERA (PMID: 22111785) 的方法注释同位素峰和加合物峰。它将结果输出为CSV文件(" annotated_peaklist.csv ")，并将带注释的峰值保存到 mSet 对象。最后，峰值列表被格式化为 MetaboAnalystR 的正确格式，并使用 FormatPeakList 函数根据用户的规范进行过滤。该函数允许加合物的过滤(即除去除[M+H]+/[M-H]-以外的所有加合物)和同位素的过滤(即除去除单同位素峰以外的所有同位素)。滤波峰值的目的是去除退化信号并减小文件大小。
# Setting the Annotation Parameters. 设定注释参数
annParams <- SetAnnotationParam(polarity = "negative", mz_abs_add = 0.005)
# Perform peak annotation. 进行峰注释
annotPeaks <- PerformPeakAnnotation(mSet, annParams)
## Format and filter the peak list for MetaboAnalystR 过滤峰
maPeaks <- FormatPeakList(annotPeaks, annParams, filtIso =F, filtAdducts = FALSE, missPercent = 1)
前8行两列
image.png
2.2 数据处理和统计分析
这一步可以通过 MetaboAnalyst 网站上的“统计分析”模块轻松完成。网站生成的运行R代码可以很容易地重复相同的结果。这个简短的章节是为了展示如何产生最多的结果。同样的方法可以实现更多的统计效用。
在相同的工作目录中，我们将读入过滤后的峰值表。然后，我们将用一个小值(原始数据中最小正值的一半)替换缺失的值。接下来，我们将过滤掉非信息信号的变量，然后将样本归一化到它们的中值并执行对数转换。最后，我们将进行t检验分析，以识别差异富集的特征，将FDR调整的p值阈值设置为0.25。
2.2.1 mSet初始化用于进一步分析
# First step is to create the mSet Object, specifying that the data to be uploaded 指明上传对象
# is a peak table ("pktable") and that statistical analysis will be performed ("stat"). 输入峰值表，进行统计分析
#第一步
mSet<-InitDataObjects("pktable", "stat", FALSE)
# Second step is to read in the filtered peak list, please set the path right first
#第二步
mSet<-Read.TextData(mSet, "metaboanalyst_input.csv", "colu", "disc")
2.2.2 进行数据统计
# The third step is to perform data processing using MetaboAnalystR (filtering/normalization)
# Perform data processing - Data checking
#第三步
mSet<-SanityCheckData(mSet)
# Perform data processing - Minimum Value Replacing
#替换最小值
mSet<-ReplaceMin(mSet);
# Perform data processing - Variable Filtering and Normalization 变量过滤和标准化
mSet<-FilterVariable(mSet, "iqr", "F", 25)
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "MedianNorm", "LogNorm", "NULL", ratio=FALSE, ratioNum=20)
#样本标准化和数值标准化绘图
mSet <- PlotNormSummary(mSet, "norm_0_", "png", 72, width=NA)
mSet <- PlotSampleNormSummary(mSet, "snorm_0_", "png", 72, width=NA)
image.png
# The fourth step is to perform fold-change analysis
#第四步进行 Fold change 分析
mSet <- FC.Anal.unpaired(mSet, 2.0, 0)
mSet <- PlotFC(mSet, "fc_0_", "png", 72, width=NA)
2.2.3 单因素t-test分析
# The fifth step is to perform t-test analysis
# 第五步
mSet <- Ttests.Anal(mSet, F, 0.05, FALSE, TRUE)
mSet <- PlotTT(mSet, "tt_0_", "png", 72, width=NA)
image.png
2.2.4 多变量PCA分析
# The sixth step is to perform PCA
# 第六步
mSet <- PCA.Anal(mSet)
mSet <- PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet <- PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA, 1,2,0.95,1,0)
mSet <- PlotPCALoading(mSet, "pca_loading_0_", "png", 72, width=NA, 1,2);
mSet <- PlotPCABiplot(mSet, "pca_biplot_0_", "png", 72, width=NA, 1,2)
mSet <- PlotPCA3DScoreImg(mSet, "pca_score3d_0_", "png", 72, width=NA, 1,2,3, 40)
2.2.5 多变量PLS-DA分析
# The seventh step is to perform PLS-DA
# 第七步
mSet <- PLSR.Anal(mSet, reg=TRUE)
mSet <- PlotPLSPairSummary(mSet, "pls_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPLS2DScore(mSet, "pls_score2d_0_", "png", 72, width=NA, 1,2,0.95,1,0)
mSet <- PlotPLS3DScoreImg(mSet, "pls_score3d_0_", "png", 72, width=NA, 1,2,3, 40)
mSet <- PlotPLSLoading(mSet, "pls_loading_0_", "png", 72, width=NA, 1, 2);
mSet <- PLSDA.CV(mSet, "L",5, "Q2")
mSet <- PlotPLS.Classification(mSet, "pls_cv_0_", "png", 72, width=NA)
mSet <- PlotPLS.Imp(mSet, "pls_imp_0_", "png", 72, width=NA, "vip", "Comp. 1", 15,FALSE)
2.3 从质谱峰到通路
以前版本的MetaboAnalyst包含两个功能分析模块，代谢途径分析(MetPA) (Xia et al. 2010, https://academic.oup.com/bioinformatics/article/26/18/2342/208464) 和代谢物集富集分析(MSEA) (Xia et al. 2010, https://academic.oup.com/nar/article/38/suppl_2/W71/1101310) 。然而，这些模块在使用前需要对代谢物进行识别，这仍然是非靶向代谢组学中的一个重要挑战。相比之下，mummichog算法(Li et al. 2013, http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003123) 绕过了途径分析之前代谢物识别的瓶颈，利用先验途径和网络知识，根据质谱峰值直接推断生物活性。因此，我们在一个名为 “MS Peaks to Pathways” 的新模块中实现了R中的mummichog算法。该模块的知识库包括来自原始Python实现的5个基因组尺度代谢模型，这些模型要么是手动管理的，要么是从 BioCyc 下载的，以及来自KEGG代谢途径的21种生物的扩展库。
为了使用此模块，用户必须使用Read.PeakListData 上传一个 .txt 的表格，包含以下信息：
1、三列信息 m/z值，p-values，t-scores值或者 fold-change values
2、两列信息 m/z值，p-values 或者 t-scores值
3、一列信息 排序后的p-values 或者 t-scores值
4、四列信息 m/z值，p-values，t-scores值或者 fold-change values，模型参数（positive or negative）
所有输入的文件必须为。txt格式。如果输入是一个三列信息，那么mummichog和GSEA算法(以及它们的组合)都可以应用。如果只提供p值(或按p值排序)，则只应用mummichog算法。如果只提供t-scores(或按t-scores排名)，则只应用GSEA算法。
如果p值还没有计算出来，用户可以使用 MetaboAnalyst 网站的“统计分析”模块或这个R包上传他们的原始峰值表，处理数据，进行 t检验或 fold-change 分析，然后将这些结果上传到模块。使用该表，用户还需要使用 UpdateInstrumentParameters函数指定MS仪器的类型、离子模式(正或负)。目前，MetaboAnalystR仅支持从高分辨率MS仪器(如Orbitrap)或傅里叶变换(FT)-MS仪器(原始mummichog实现推荐)获得的峰的处理。
在数据上传后，用户可以选择生物库，从中使用 SetMass 执行非靶向通路分析。PathLib 函数。然后用户可以使用PerformPSEA对他们的数据执行mummichog算法。首先，用户将使用setpeakenrichment method设置用于mummichog的算法。其次，用户将使用SetMummichogPval函数设置p值截断，以区分显著富集和非显著富集的m/z特征。最后，使用PerformPSEA来计算通路活性。
这个模块的输出首先包含一个结果表，该表标识用户上传数据中丰富的top 路径，可以在名为 “mummichog_pathway_enrichment.csv” 的工作目录中找到。该表格包含每个路径的总命中数、原始p值(Fisher’s或Hypergeometric)、EASE评分和调整后的p值(用于排列)。在您的工作目录中还可以找到第二个表，名为“mummichog_matched_compound_all.csv”，其中包含用户上传的m/z特征列表中所有匹配的代谢物。
在本教程中，我们将首先直接使用从上述IBD数据的非靶向代谢组学中获得的示例峰值列表数据(质量精度5.0，负离子模式)。这是由原始数据处理管道生成的原始数据表。
2.3.1 Mummichog Version 1
对MS峰值进行T检验，以进行通路分析
# Perform t-test
library(MetaboAnalystR);
mSet<-InitDataObjects("pktable", "stat", FALSE);
mSet<-Read.TextData(mSet, "metaboanalyst_input.csv", "colu", "disc")
mSet<-SanityCheckData(mSet)
mSet<-ReplaceMin(mSet);
mSet<-FilterVariable(mSet, "iqr", "F", 25)
feature.nm.vec <- c("")# used to remove certain feature, i.e. c("157.0357")
smpl.nm.vec <- c("")# used to remove certain samples, i.e. c("samples 1")
grp.nm.vec <- c("") # used to remove certain groups, i.e. c("UC")
mSet <- UpdateData(mSet)
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "MedianNorm", "LogNorm", "NULL", ratio=FALSE, ratioNum=20)
mSet<-Ttests.Anal(mSet, F, 0.25, FALSE, TRUE)
# [1] "A total of 1329 significant features were found."
cmSet<-Convert2Mummichog(mSet, rt=TRUE)
mSet<-InitDataObjects("mass_all", "mummichog", FALSE)
SetPeakFormat("mprt")
mSet<-UpdateInstrumentParameters(mSet, 5, "negative");
mSet<-Read.PeakListData(mSet, "mummichog_input_your_date.txt");
mSet<-SanityCheckMummichogData(mSet)
mSet<-SetPeakEnrichMethod(mSet, "integ", version="v2")
# Please set the appropriate p value according to your data
mSet<-SetMummichogPval(mSet, 0.15)
mSet<-PerformPSEA(mSet, "hsa_mfn", "current", permNum =100)
mSet<-PlotPSEAIntegPaths(mSet, "peaks_to_paths_", "png"144)
2.3.2 Mummichog Version 2
SetPeakEnrichMethod 现在有一个参数，让用户选择是使用MS峰值到路径算法的版本1 还是版本2。在版本2中，用户现在可以上传保留时间信息来执行路径分析。请注意，只有当用户数据包含保留时间(“rt”或“r.t”)列时，才应该使用版本2。保留时间用于将路径分析从“化合物”空间移动到“经验化合物”空间(详细信息见“如何计算经验化合物?”)。保留时间的加入将增加潜在化合物匹配的可信度和稳健性。另一个区别是，化合物直接从用户选择的路径库中移除，而不是在排列过程中从潜在的化合物命中移除。
# Perform t-test
mSet<-Ttests.Anal(mSet, F, 0.25, FALSE, TRUE)
## [1] "A total of 1329 significant features were found."
mSet<-Convert2Mummichog(mSet, rt=TRUE)
mSet<-InitDataObjects("mass_all", "mummichog", FALSE)
SetPeakFormat("mprt")
mSet<-UpdateInstrumentParameters(mSet, 5, "negative");
mSet<-Read.PeakListData(mSet, "mummichog_input_your_date.txt");
mSet<-SanityCheckMummichogData(mSet)
mSet<-SetPeakEnrichMethod(mSet, "integ", version="v2")
# Please set the appropriate p value according to your data
mSet<-SetMummichogPval(mSet, 0.15)
mSet<-PerformPSEA(mSet, "hsa_mfn", "current", permNum = 100)
mSet<-PlotPSEAIntegPaths(mSet, "peaks_to_paths_", "png", 144)
2.3.3 加合物定制 & Currency 定制
加合物定制
原始质谱峰含有大量的加合物，对应于不同的质谱仪器和分析模式。“可用”面板中显示了一个全面的加合列表。使用它来定制Mummichog/GSEA算法使用的加合列表，以匹配m/z峰值与潜在化合物命中值。可供选择的加合物列表为 here: positive; here: negative ; here:mixed。
对于负离子模式，加合物列表为：M-H[-], M-2H[2-], M(C13)-H[-], M(S34)-H[-], M(Cl37)-H[-], M+Na-2H[-], M+K-2H[-], M-H2O-H[-], M+Cl[-], M+Cl37[-], M+Br[-], M+Br81[-], M+ACN-H[-], M+HCOO[-], M+CH3COO[-], and M-H+O[-].
对于正离子模式，加合物列表为：M[1+], M+H[1+], M+2H[2+], M+3H[3+], M(C13)+H[1+], M(C13)+2H[2+], M(C13)+3H[3+], M(S34)+H[1+], M(Cl37)+H[1+], M+Na[1+], M+H+Na[2+], M+K[1+], M+H2O+H[1+], M-H2O+H[1+], M-H4O2+H[1+], M-NH3+H[1+], M-CO+H[1+], M-CO2+H[1+], M-HCOOH+H[1+], M+HCOONa[1+], M-HCOONa+H[1+], M+NaCl[1+], M-C3H4O2+H[1+], M+HCOOK[1+], and M-HCOOK+H[1+].
Currency 定制
Currency 代谢物是已知存在于正常功能细胞中并参与大量代谢反应的丰富物质，如水和二氧化碳(Huss和Holme, 2007)。由于它们无处不在的特性，去除它们可以极大地改善路径分析和可视化。可用的Currency 代谢物列表可以在这里找到。默认情况下，MS从峰值到路径模块将这些代谢产物视为货币: ' C00001 '， ' C00080 '， ' C00007 '， ' C00006 '， ' C00005 '， ' C00003 '， ' C00004 '， ' C00002 '， ' C00013 '， ' C00008 '， ' C00009 '， ' C00011 '， ' G11113 '， ' H2O '， ' H+ '， ' Oxygen '， ' NADP+ '， ' NADPH '， ' NAD+ '， ' NADH '， ' ATP '， ' Pyrophosphate '， ' ADP '， '正磷酸盐'，' CO2 '。
现在我们将使用另一个从小鼠肺泡巨噬细胞的非靶向代谢组学中获得的峰值列表数据，使用Orbitrap LC-MS (C18负离子模式和HILIC正离子模式)。这是一个包含m.z、mode、p.value和t.score的四列表的示例。我们将执行Currency 和加合自定义。
# Create objects for storing processed data from the MS peaks to pathways module
mSet<-InitDataObjects("mass_all", "mummichog", FALSE)
# Set peak formart - contains m/z features, p-values and t-scores
SetPeakFormat("mpt")
mSet<-UpdateInstrumentParameters(mSet, 5.0, "mixed");
mSet<-Read.PeakListData(mSet, "https://www.metaboanalyst.ca/MetaboAnalyst/resources/data/mummichog_mixed.txt");
mSet<-SanityCheckMummichogData(mSet)
# Customize currency
curr.vec <-c("Water (C00001)""Proton (C00080)","Oxygen (C00007)","NADPH (C00005)","NADP (C00006)","NADH (C00004)","NAD (C00003)","Adenosine triphosphate (C00002)","Pyrophosphate (C00013)","Phosphate (C00009)","Carbon dioxide (C00011)","Hydrogen (C00282)","Hydrogen peroxide (C00027)","Sodium (C01330)")
# Map selected currency to user's data
mSet<-Setup.MapData(mSet, curr.vec);
mSet<-PerformCurrencyMapping(mSet)
# Now customize adducts
add.vec <-c("M [1+]","M+H [1+]","M+2H [2+]","M+3H [3+]","M+Na [1+]","M+H+Na [2+]","M+K [1+]","M+H2O+H [1+]","M-H2O+H [1+]","M-H4O2+H [1+]","M(C13)+H [1+]","M(C13)+2H [2+]","M(C13)+3H [3+]","M(Cl37)+H [1+]","M-NH3+H [1+]","M-CO+H [1+]","M-CO2+H [1+]","M-HCOOH+H [1+]","M+HCOONa [1+]","M-HCOONa+H [1+]","M+NaCl [1+]","M-C3H4O2+H [1+]","M+HCOOK [1+]","M-HCOOK+H [1+]","M-H [1-]","M-2H [2-]","M-H2O-H [1-]","M-H+O [1-]","M+K-2H [1-]","M+Na-2H [1- ]","M+Cl [1-]","M+Cl37 [1-]","M+HCOO [1-]","M+CH3COO [1-]")
# Set up the selected adducts
mSet<-Setup.AdductData(mSet, add.vec);
mSet<-PerformAdductMapping(mSet, "mixed")
# Perform mummichog algorithm using selected currency and adducts, using Version1 of the mummichog algorithm
mSet<-SetPeakEnrichMethod(mSet, "mum", "v1")
mSet<-SetMummichogPval(mSet, 1.0E-5)
mSet<-PerformPSEA(mSet, "hsa_mfn", "current", 100)
mSet<-PlotPeaks2Paths(mSet, "peaks_to_paths_0_", "png", 72, width=NA) # Image is not shown below to avoid to be overwhelmed.
2.3.4 定制经验化合物格式
默认情况下，V2 版本的MS峰值到路径算法强制主离子在创建经验化合物时存在(上面的步骤3)。用户可以在UpdateInstrumentParameters 函数中通过设置“force_primary_ion”参数为“no”来禁用此功能。此外，在默认情况下，用于分割Empirical compound 的保留时间窗口计算为用户数据中的最大保留时间乘以保留时间分数(默认为0.02)。用户可以更改保留时间分数(rt_frac)或设置保留时间窗口(rt_tol -以秒为单位)。代码细节如下所示。
# Disable force primary ion
mSet <-UpdateInstrumentParameters(mSet, instrumentOpt, msModeOpt, force_primary_ion ="no", rt_frac =0.02,rt_tol =NA)
# Change retention time fraction when calculating the retention time window
mSet <-UpdateInstrumentParameters(mSet, instrumentOpt, msModeOpt,force_primary_ion ="yes", rt_frac =0.025, rt_tol =NA)
# Set the retention time window (in seconds)
mSet <-UpdateInstrumentParameters(mSet, instrumentOpt, msModeOpt, force_primary_ion ="yes", rt_frac =0.02, rt_tol =25)
3 特色使用程序
3.1 批次效应校正
批次效应校正介绍
批次效应校正分析模块对色谱峰筛选步骤生成的峰值数据表执行自动或特定的校正。内置多种校正方法。它们包括 combat, WaveICA, EigenMS, QC_RLSC, ANCOVA, RUV_random, RUV_2, RUVseq 系列(RUV_s, RUV_r和RUV_g)， NOMIS, CCMN。文中详细阐述了不同方法的局限性和数学机理。这里提供了4个案例，向用户展示了批量效应和信号漂移校正的逐步工作流程。
数据准备
FormatPeakList函数生成的峰值表需要手工制作，以补充下面的相应信息。
第一列样品名称;
第二列样品组/类信息，如有QC样品，请标记为QC;
第三列批信息，请在同一批内提供一致的字符;
第4列订单信息，请提供整个样品的注塑订单信息;其他柱应提供特征强度。如有内部标准，请用“IS”标记该特性。
注:请提供尽可能多的资料。如有遗漏，请留空。
image.png
3.1.1 示例1—— IBD基准数据
这个Benchmark 数据是一个超过600个样本的大批量数据。数据文件大小超过70mb。这部分是为了重复我们发表的结果而设计的。此部分的运行可能需要30分钟以上才能完成。强烈建议您使用自己的数据，而不是这个大文件，以避免消耗耐心。
数据下载
# Use Google API for data downloading peak feature data generated by FormatPeakList here.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".csv")
# Please authorize your google account to access the data
dl1 <- drive_download(
as_id("1wEh2P81J_xFWJs5y4mq98-FsjxJ5wmBO"), path = temp, overwrite = TRUE)
# Use Google API for data downloading meta data here.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".csv")
# Please authorize your google account to access the data
dl2 <- drive_download(
as_id("1KaBnSNRrirVPvpRxIubGCqpjX8asNeVA"), path = temp, overwrite = TRUE)
数据准备
# Load the MetaboAnalystR package
library("MetaboAnalystR")
# Data preparation - read data in & transpose.
# This is a reference example for user to prepare their data.
# Please prepare your data table according to your data format.
MetaboAna_Data <- t(read.csv(dl1$local_path,header = T));
colnames(MetaboAna_Data) <- MetaboAna_Data[1,];
MetaboAna_Data <- MetaboAna_Data[-1,];
MetaboAna_Data <- MetaboAna_Data[order(rownames(MetaboAna_Data)),];
meta_data <- read.csv(dl2$local_path);
meta_data <- meta_data[order(meta_data[,1]),c(1,2,4)];
Prepared_Data <- cbind(meta_data,MetaboAna_Data)[,-4];
write.csv(Prepared_Data,file = "IBD_BC_correction.csv",row.names = F)
datapath <- paste0(getwd(),"/IBD_BC_correction.csv")
数据过滤和标准化
mSet<-InitDataObjects("pktable", "stat", FALSE)
mSet<-Read.TextData(mSet, dl1$local_path, "col", "disc")
mSet<-SanityCheckData(mSet)
mSet<-ReplaceMin(mSet);
mSet<-FilterVariable(mSet, "iqr", "F", 25)
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "MedianNorm", "LogNorm", "NULL", ratio=FALSE, ratioNum=20)
### Data Orgnization
normalized_set <- mSet[["dataSet"]][["norm"]]
ordered_normalized_set <- normalized_set[order(row.names(normalized_set)), ]
### import metadata
meta_data <- read.csv(dl2$local_path)
new_normalized_set <- cbind(meta_data[,2:4], ordered_normalized_set);
write.csv(new_normalized_set,file = "new_normalized_set.csv")
#perform PCA
mSet <- PCA.Anal(mSet)
mSet <- PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet <- PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA,style = 2, 1,2,0.95,0,0)
rm(mSet)
image.png
#初始化 mSet 对象
mSet <- InitDataObjects("pktable", "utils", FALSE)
#数据读取
## we set samples in "row" according to the table format. If your samples are in column, set it as "col".
mSet <- Read.BatchDataTB(mSet, "new_normalized_set.csv", "row")
#自动校正
mSet <- PerformBatchCorrection(mSet)
getwd()
# 展示距离和批次信息
mSet$dataSet$interbatch_dis
image.png
数据可视化——分组PCA图绘制
## remove the batch column in the corrected result
data_corrected <- read.csv("/home/xialab/Documents/MetaboAnalyst_batch_data.csv")
data_corrected_new <- data_corrected[,-3]
write.csv(data_corrected_new,file = "MetaboAnalyst_batch_data_stats.csv",row.names = F)
mSet <- InitDataObjects("pktable", "stat", FALSE)
mSet <- Read.TextData(mSet, "MetaboAnalyst_batch_data_stats.csv", "row", "disc")
mSet <- SanityCheckData(mSet)
mSet <- ReplaceMin(mSet);
mSet <- FilterVariable(mSet, "iqr", "F", 25)
mSet <- PreparePrenormData(mSet)
mSet <- Normalization(mSet, "NULL", "NULL", "NULL", ratio=FALSE, ratioNum=20) # No need to normalize again
mSet <- PCA.Anal(mSet)
mSet <- PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet <- PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA, style = 2, 1,2,0.95,0,0)
image.png
3.1.2 示例2——信号漂移校正
除了批次效应，信号漂移也是代谢组学研究面临的一个重要问题。QC-RLSC 既可以校正批效应，也可以校正信号漂移(mmidd: 30253838)。在这里，我们展示了一个例子，用户纠正信号漂移与模拟数据表。需要注意的一点是，对于case 1和case 2的批量效果修正，只需要批量信息。但对于信号漂移校正，必须提供第4列或第4行注入顺序。而批处理信息是可选的，应该添加在第3列或第3行。如果缺少，请空列或行或标记为一个批次。您可以按照以下示例准备数据表。
数据下载和准备
# Use Google API for data downloading peak feature.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".csv")
# Please authorize your google account to access the data
dl3 <- drive_download(
as_id("1cPs3vZhB1gVCOV3ER9BquMAFSjSvmDYe"), path = temp, overwrite = TRUE)
#初始化 mSet 对象
mSet <- InitDataObjects("pktable", "utils", FALSE)
#读取数据
## we set samples in "row" according to the table format. If your samples are in column, set it as "col".
mSet <- Read.SignalDriftData(mSet, dl3$local_path, "row")
#信号漂移校正
mSet <- PerformSignalDriftCorrection(mSet)
image.png
3.2 化合物名称匹配
MetaboAnalyst 现在为用户提供微服务，使用我们全面的内部代谢物数据库(20万种化合物)执行化合物名称映射。在R中，要使用这个微服务，用户首先必须安装 httr 包。请求必须是POST请求，包含化合物列表和输入化合物类型，并被发送到api.metaboanalyst.ca/mapcompounds。下面的代码将逐步展示如何使用MetaboAnalyst API执行复合名称映射。其他编程语言请参考MetaboAnalyst的api页面。
# First create a list containing a vector of the compounds to be queried (separated by a semi-colon)
# and another character vector containing the compound id type.
# The items in the list MUST be queryList and inputType
# Valid input types are: "name", "hmdb", "kegg", "pubchem", "chebi", "metlin"
name.vec <- c("1,3-Diaminopropane;2-Ketobutyric acid;2-Hydroxybutyric acid;2-Methoxyestrone")
toSend = list(queryList = name.vec, inputType = "name")
library(httr)
# The MetaboAnalyst API url
call <- "http://api.xialab.ca/mapcompounds"
# Use httr::POST to send the request to the MetaboAnalyst API
# The response will be saved in query_results
query_results <- httr::POST(call, body = toSend, encode = "json")
# Check if response is ok (TRUE)
# 200 is ok! 401 means an error has occured on the user's end.
query_results$status_code==200
# Parse the response into a table
# Will show mapping to "hmdb_id", "kegg_id", "pubchem_id", "chebi_id", "metlin_id", "smiles"
query_results_text <- content(query_results, "text", encoding = "UTF-8")
query_results_json <- rjson::fromJSON(query_results_text, flatten = TRUE)
query_results_table <- t(rbind.data.frame(query_results_json))
rownames(query_results_table) <- query_results_table[,1]
4 其他分析
MetaboAnalystR已设计与MetaboAnalyst网站同步进行全面的代谢组学分析。生物标记分析，富集分析，元分析，途径分析，综合途径分析，功率分析模块，时间序列或双因素设计，网络浏览器模块，MS峰到路径，批量效应校正等都可以在网站上轻松实现。更重要的是，可以使用MetaboAnalystR重复各种分析结果，包括生成图形，以便进行进一步的高级编辑。MetaboAnalyst网站还为用户提供了实时的R码来完成这一过程。因此，在原始数据处理后，强烈建议使用网站进行进一步处理。为高级用户提供一种简便的本地分析方法。我们将不同模块的所有小插图总结如下，作为R端深度数据挖掘的指导。
网友评论