R包文檔參考：
https://www.metaboanalyst.ca/docs/RTutorial.xhtml#Workflow:%20Case%20Study%201
R包函數(shù)參考：
https://rdrr.io/github/xia-lab/MetaboAnalystR3.0/man/

MetaboAnalystR包與網(wǎng)頁同步，可以用于分析原始組學數(shù)據(jù)和批次分析，建立了一個從原始數(shù)據(jù)開始分析的優(yōu)化流程遣铝。下面的教程旨在通過提供使用R包的幾個最常見任務的逐步指導來補充我們的基于web的功能。

1.1 MetaboAnalystR 3簡介

MetaboAnalystR 3包含MetaboAnalyst web服務器下的R函數(shù)和庫莉擒，包括代謝組學數(shù)據(jù)分析酿炸、可視化和功能解釋。該包與MetaboAnalyst web服務器同步涨冀。在安裝和加載包后填硕，用戶將能夠使用從MetaboAnalyst網(wǎng)站下載的相應的R命令歷史從他們的本地計算機上重現(xiàn)相同的結果，從而實現(xiàn)最大的靈活性和可再現(xiàn)性鹿鳖。

version3 旨在通過實現(xiàn)峰值選擇的快速參數(shù)優(yōu)化算法來改進當前的全球代謝組學工作流程扁眯，并從12種成熟的方法中自動識別最適合的批量效應校正方法。此外翅帜，還增加了通過mummichog2 (PMID: 23861661)直接從m/z峰進行功能解釋的更多支持姻檀，并添加了一種新的基于路徑的方法來集成多組學數(shù)據(jù)。為了演示這一新功能涝滴，我們在本教程中對臨床IBD樣本進行端到端代謝組學數(shù)據(jù)分析绣版。在這個R包的小插圖中提供了更多的案例研究。在這里歼疮，我們更愿意提供一個全面的MetaboAnalystR 3.2安裝入門教程(更新于2021年12月16日)杂抽。

1.2 MetaboAnalystR 安裝

略

2 IBD案例研究

炎癥性腸道疾病，包括克羅恩病和潰瘍性結腸炎韩脏，已經影響了世界各地數(shù)百萬人缩麸。Jason L.等人對微生物組在IBD發(fā)病機制中的作用進行了長時間多組學研究。本文以這一新型參數(shù)優(yōu)化為例赡矢，介紹了面部樣本的代謝組學研究匙睹。

2.1 原始數(shù)據(jù)處理

原始數(shù)據(jù)處理是所有代謝組學研究的首要環(huán)節(jié)愚屁。MetaboAnalystR 3.2支持”.mzXML”、“.mzML "和".CDF”格式痕檬。原始格式(例如“.raw "和" . RAW")將很快得到支持霎槐。您可以使用 ProteoWizard (PMID: 23051804) 將原始數(shù)據(jù)轉換為支持的格式。優(yōu)先使用中心化數(shù)據(jù)梦谜。

2.1.1 加載 MetaboAnalystR包

# Load the MetaboAnalystR package
library("MetaboAnalystR")

2.1.2 下載IBD數(shù)據(jù)集QC數(shù)據(jù)

下面的步驟將使用質量控制(QC)示例丘跌。下載MS數(shù)據(jù)進行參數(shù)優(yōu)化。達到快速學習的目的唁桩，避免整個樣本(超過600個樣本)長時間運行闭树。在這里，我們只為每個CD組和非IBD組提供5個樣本荒澡。如果您想重復和驗證我們手稿中的結果报辱，請訪問IBD MultiOmics數(shù)據(jù)庫，下載完整批數(shù)據(jù)并使用此管道運行它們单山。注：數(shù)據(jù)因網(wǎng)絡無法下載

## 1 設置數(shù)據(jù)存放目錄
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數(shù)據(jù)
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數(shù)據(jù)
out <- unzip(temp, exdir = data_folder_Sample)
# Date files for normal processing example are deposited below
out

2.1.3 數(shù)據(jù)核查

在運行數(shù)據(jù)分析之前米奸，可以使用 PerformDataInspect 檢查一般數(shù)據(jù)結構和信息昼接。如果有一些非常重要的污染物，它將被直接發(fā)現(xiàn)悴晰。

# 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軸是信號強度

2.1.4 ROI（regions of interest 目標區(qū)域） 提取

QC 樣品數(shù)據(jù) ROI 提取現(xiàn)在使用 PerformROIExtraction 或 PerformDataTrimming 進行铡溪。數(shù)據(jù) ROI 提取有3種模式漂辐。默認是標準模擬方法(“ssm_trim”)。在此模式下棕硫，三維MS將先從m/z維度進行裁剪髓涯，再從RT維度進行裁剪。rt.idx可以用來調整RT數(shù)據(jù)保留的范圍饲帅。另外兩種模式被命名為 “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 初始化參數(shù)設定

已嵌入儀器特定的初始參數(shù)赦邻。目前，最流行的平臺包括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”模式均已配置用于進一步處理。此外篓跛，具體參數(shù)也可以單獨設置膝捞，如果你對參數(shù)有自己的看法，此參數(shù)優(yōu)化步驟可跳過愧沟。

# Initial platform specific parameters
param_initial <- SetPeakParam(platform = "UPLC-Q/E") 

2.1.6 參數(shù)優(yōu)化

采用“DoE”策略進行參數(shù)優(yōu)化蔬咬。初始參數(shù)來自上述優(yōu)化的參數(shù)或來自嵌入式 SetPeakParam 函數(shù)的內部默認參數(shù)。

# 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 函數(shù)讀取原始 MS 數(shù)據(jù)文件央渣，并將其保存為 OnDiskMSnExp 對象计盒。用SetPlotParam 函數(shù)集“ Plot = T ”可以輸出兩個圖——總離子色譜圖(TIC)提供整個色譜的概況，基峰色譜圖( BPC )是基于最豐富信號的光譜更清晰的輪廓芽丹。這些圖對下游參數(shù)的設置很有用北启。對于希望查看感興趣的峰值的用戶凡桥，可以使用PlotEIC 函數(shù)生成提取離子色譜圖(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 函數(shù)是從XCMS R函數(shù)中更新的峰值處理流程尤误，它自動執(zhí)行峰值檢測伙单、對齊和分組佩伤。該函數(shù)還生成兩個診斷圖凉翻，包括不同樣本中峰值總強度的統(tǒng)計信息理逊、保留時間調整圖和顯示數(shù)據(jù)清洗和統(tǒng)計分析之前總體樣本聚類的PCA圖弟疆。

# Peak Profiling with optimized parameters
mSet <- PerformPeakProfiling(rawData,param_optimized$best_parameters,
                             plotSettings = SetPlotParam(Plot = T))

2.1.8 峰注釋

PerformPeakAnnotation 函數(shù)使用 CAMERA (PMID: 22111785) 的方法注釋同位素峰和加合物峰泳猬。它將結果輸出為CSV文件(" annotated_peaklist.csv ")批钠，并將帶注釋的峰值保存到 mSet 對象。最后得封，峰值列表被格式化為 MetaboAnalystR 的正確格式埋心，并使用 FormatPeakList 函數(shù)根據(jù)用戶的規(guī)范進行過濾。該函數(shù)允許加合物的過濾(即除去除[M+H]+/[M-H]-以外的所有加合物)和同位素的過濾(即除去除單同位素峰以外的所有同位素)忙上。濾波峰值的目的是去除退化信號并減小文件大小拷呆。

# Setting the Annotation Parameters. 設定注釋參數(shù)
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行兩列

2.2 數(shù)據(jù)處理和統(tǒng)計分析

這一步可以通過 MetaboAnalyst 網(wǎng)站上的“統(tǒng)計分析”模塊輕松完成。網(wǎng)站生成的運行R代碼可以很容易地重復相同的結果。這個簡短的章節(jié)是為了展示如何產生最多的結果茬斧。同樣的方法可以實現(xiàn)更多的統(tǒng)計效用腰懂。

在相同的工作目錄中，我們將讀入過濾后的峰值表项秉。然后绣溜，我們將用一個小值(原始數(shù)據(jù)中最小正值的一半)替換缺失的值。接下來伙狐，我們將過濾掉非信息信號的變量涮毫，然后將樣本歸一化到它們的中值并執(zhí)行對數(shù)轉換。最后贷屎，我們將進行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"). 輸入峰值表咒吐，進行統(tǒng)計分析
#第一步
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 進行數(shù)據(jù)統(tǒng)計

# 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)

#樣本標準化和數(shù)值標準化繪圖
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)

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) 属划。然而恬叹，這些模塊在使用前需要對代謝物進行識別，這仍然是非靶向代謝組學中的一個重要挑戰(zhàn)同眯。相比之下绽昼，mummichog算法(Li et al. 2013, http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003123) 繞過了途徑分析之前代謝物識別的瓶頸，利用先驗途徑和網(wǎng)絡知識须蜗，根據(jù)質譜峰值直接推斷生物活性硅确。因此，我們在一個名為 “MS Peaks to Pathways” 的新模塊中實現(xiàn)了R中的mummichog算法明肮。該模塊的知識庫包括來自原始Python實現(xiàn)的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搁凸，模型參數(shù)（positive or negative）
所有輸入的文件必須為。txt格式狠毯。如果輸入是一個三列信息护糖，那么mummichog和GSEA算法(以及它們的組合)都可以應用。如果只提供p值(或按p值排序)嚼松，則只應用mummichog算法嫡良。如果只提供t-scores(或按t-scores排名)，則只應用GSEA算法献酗。
如果p值還沒有計算出來寝受，用戶可以使用 MetaboAnalyst 網(wǎng)站的“統(tǒng)計分析”模塊或這個R包上傳他們的原始峰值表，處理數(shù)據(jù)罕偎，進行 t檢驗或 fold-change 分析很澄，然后將這些結果上傳到模塊。使用該表颜及，用戶還需要使用 UpdateInstrumentParameters函數(shù)指定MS儀器的類型甩苛、離子模式(正或負)。目前俏站，MetaboAnalystR僅支持從高分辨率MS儀器(如Orbitrap)或傅里葉變換(FT)-MS儀器(原始mummichog實現(xiàn)推薦)獲得的峰的處理讯蒲。
在數(shù)據(jù)上傳后，用戶可以選擇生物庫肄扎，從中使用 SetMass 執(zhí)行非靶向通路分析墨林。PathLib 函數(shù)。然后用戶可以使用PerformPSEA對他們的數(shù)據(jù)執(zhí)行mummichog算法反浓。首先萌丈，用戶將使用setpeakenrichment method設置用于mummichog的算法。其次雷则，用戶將使用SetMummichogPval函數(shù)設置p值截斷辆雾，以區(qū)分顯著富集和非顯著富集的m/z特征。最后月劈，使用PerformPSEA來計算通路活性度迂。
這個模塊的輸出首先包含一個結果表，該表標識用戶上傳數(shù)據(jù)中豐富的top 路徑猜揪，可以在名為 “mummichog_pathway_enrichment.csv” 的工作目錄中找到惭墓。該表格包含每個路徑的總命中數(shù)、原始p值(Fisher’s或Hypergeometric)而姐、EASE評分和調整后的p值(用于排列)腊凶。在您的工作目錄中還可以找到第二個表，名為“mummichog_matched_compound_all.csv”，其中包含用戶上傳的m/z特征列表中所有匹配的代謝物钧萍。
在本教程中褐缠，我們將首先直接使用從上述IBD數(shù)據(jù)的非靶向代謝組學中獲得的示例峰值列表數(shù)據(jù)(質量精度5.0，負離子模式)风瘦。這是由原始數(shù)據(jù)處理管道生成的原始數(shù)據(jù)表队魏。

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 現(xiàn)在有一個參數(shù)万搔，讓用戶選擇是使用MS峰值到路徑算法的版本1 還是版本2胡桨。在版本2中，用戶現(xiàn)在可以上傳保留時間信息來執(zhí)行路徑分析瞬雹。請注意昧谊，只有當用戶數(shù)據(jù)包含保留時間(“rt”或“r.t”)列時，才應該使用版本2酗捌。保留時間用于將路徑分析從“化合物”空間移動到“經驗化合物”空間(詳細信息見“如何計算經驗化合物?”)揽浙。保留時間的加入將增加潛在化合物匹配的可信度和穩(wěn)健性。另一個區(qū)別是意敛，化合物直接從用戶選擇的路徑庫中移除馅巷，而不是在排列過程中從潛在的化合物命中移除。

# 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 '。
現(xiàn)在我們將使用另一個從小鼠肺泡巨噬細胞的非靶向代謝組學中獲得的峰值列表數(shù)據(jù)冰啃，使用Orbitrap LC-MS (C18負離子模式和HILIC正離子模式)邓夕。這是一個包含m.z、mode阎毅、p.value和t.score的四列表的示例焚刚。我們將執(zhí)行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峰值到路徑算法強制主離子在創(chuàng)建經驗化合物時存在(上面的步驟3)矿咕。用戶可以在UpdateInstrumentParameters 函數(shù)中通過設置“force_primary_ion”參數(shù)為“no”來禁用此功能。此外狼钮，在默認情況下痴腌，用于分割Empirical compound 的保留時間窗口計算為用戶數(shù)據(jù)中的最大保留時間乘以保留時間分數(shù)(默認為0.02)。用戶可以更改保留時間分數(shù)(rt_frac)或設置保留時間窗口(rt_tol -以秒為單位)燃领。代碼細節(jié)如下所示士聪。

# 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 批次效應校正

批次效應校正介紹
批次效應校正分析模塊對色譜峰篩選步驟生成的峰值數(shù)據(jù)表執(zhí)行自動或特定的校正。內置多種校正方法猛蔽。它們包括 combat, WaveICA, EigenMS, QC_RLSC, ANCOVA, RUV_random, RUV_2, RUVseq 系列(RUV_s, RUV_r和RUV_g)剥悟， NOMIS, CCMN灵寺。文中詳細闡述了不同方法的局限性和數(shù)學機理。這里提供了4個案例区岗，向用戶展示了批量效應和信號漂移校正的逐步工作流程略板。

數(shù)據(jù)準備
FormatPeakList函數(shù)生成的峰值表需要手工制作，以補充下面的相應信息慈缔。
第一列樣品名稱;
第二列樣品組/類信息叮称，如有QC樣品，請標記為QC;
第三列批信息藐鹤，請在同一批內提供一致的字符;
第4列訂單信息瓤檐，請?zhí)峁┱麄€樣品的注塑訂單信息;其他柱應提供特征強度。如有內部標準娱节，請用“IS”標記該特性挠蛉。
注:請?zhí)峁┍M可能多的資料。如有遺漏肄满，請留空谴古。

3.1.1 示例1—— IBD基準數(shù)據(jù)

這個Benchmark 數(shù)據(jù)是一個超過600個樣本的大批量數(shù)據(jù)。數(shù)據(jù)文件大小超過70mb稠歉。這部分是為了重復我們發(fā)表的結果而設計的掰担。此部分的運行可能需要30分鐘以上才能完成。強烈建議您使用自己的數(shù)據(jù)怒炸，而不是這個大文件恩敌，以避免消耗耐心。
數(shù)據(jù)下載

# 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)

數(shù)據(jù)準備

# 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")

數(shù)據(jù)過濾和標準化

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)

#初始化 mSet 對象
mSet <- InitDataObjects("pktable", "utils", FALSE)

#數(shù)據(jù)讀取
## 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

數(shù)據(jù)可視化——分組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)

3.1.2 示例2——信號漂移校正

除了批次效應横媚，信號漂移也是代謝組學研究面臨的一個重要問題纠炮。QC-RLSC 既可以校正批效應，也可以校正信號漂移(mmidd: 30253838)灯蝴。在這里恢口，我們展示了一個例子，用戶糾正信號漂移與模擬數(shù)據(jù)表穷躁。需要注意的一點是耕肩，對于case 1和case 2的批量效果修正，只需要批量信息问潭。但對于信號漂移校正猿诸，必須提供第4列或第4行注入順序。而批處理信息是可選的狡忙，應該添加在第3列或第3行梳虽。如果缺少，請空列或行或標記為一個批次灾茁。您可以按照以下示例準備數(shù)據(jù)表窜觉。

數(shù)據(jù)下載和準備

# 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)

#讀取數(shù)據(jù)
## 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)

3.2 化合物名稱匹配

MetaboAnalyst 現(xiàn)在為用戶提供微服務谷炸，使用我們全面的內部代謝物數(shù)據(jù)庫(20萬種化合物)執(zhí)行化合物名稱映射。在R中禀挫，要使用這個微服務旬陡，用戶首先必須安裝 httr 包。請求必須是POST請求语婴，包含化合物列表和輸入化合物類型描孟，并被發(fā)送到api.metaboanalyst.ca/mapcompounds。下面的代碼將逐步展示如何使用MetaboAnalyst API執(zhí)行復合名稱映射砰左。其他編程語言請參考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網(wǎng)站同步進行全面的代謝組學分析。生物標記分析菜职，富集分析，元分析旗闽，途徑分析酬核，綜合途徑分析，功率分析模塊适室，時間序列或雙因素設計嫡意，網(wǎng)絡瀏覽器模塊，MS峰到路徑捣辆，批量效應校正等都可以在網(wǎng)站上輕松實現(xiàn)蔬螟。更重要的是，可以使用MetaboAnalystR重復各種分析結果汽畴，包括生成圖形旧巾，以便進行進一步的高級編輯。MetaboAnalyst網(wǎng)站還為用戶提供了實時的R碼來完成這一過程忍些。因此鲁猩，在原始數(shù)據(jù)處理后，強烈建議使用網(wǎng)站進行進一步處理罢坝。為高級用戶提供一種簡便的本地分析方法廓握。我們將不同模塊的所有小插圖總結如下，作為R端深度數(shù)據(jù)挖掘的指導嘁酿。

4.1 Biomaker分析

4.2 富集分析

4.3 Meta分析

4.4 Pathway分析

4.5 整合Pathway分析

4.6 功效分析模塊

4.7 時間序列分析和雙因素設計

4.8 網(wǎng)絡分析模塊