翻译 - LCMS data preprocessing and analysis with xcms
说明：最近在看一些代谢组学数据的处理，其中XCMS这个R包是非常好用的，但是目前对这个包还不熟悉。因此，到期官网查看相关文档，发现有几篇学习的文档。所以便翻译，边理解。也与大家共享。由于专业知识欠缺，有错误的地方请指正。
实际上，大家最好看原文才能更好的理解这个包的用法。
原文地址：http://bioconductor.org/packages/release/bioc/vignettes/xcms/inst/doc/xcms.html
LCMS data preprocessing and analysis with xcms
Time：30 April 2018
Package：xcms 3.2.0
Package: xcms
Authors: Johannes Rainer
Modified: 2018-04-30 13:21:36
Compiled: Mon Apr 30 18:47:52 2018
1. 介绍（Introduction）
本文档描述了基于xcms（version > 3）的LCMS实验数据的导入、探索、预处理和分析。实例和基本工作流程改编自Colin A. Smith的构建的基于xcms的原始LC/MS预处理和分析过程。
新的用户界面和方法使用XCMSnExp对象（代替旧的xcmsSet对象）作为预处理结果的容器。为了支持依赖于xcmsSet对象的包（packages）和管道（pipelines），可以使用as方法将XCMSnExp转换为xcmsSet对象，比如xset < - as（x，"xcmsSet")，其中x是XCMSnExp对象。
2. 数据导入（Data import）
xcms支持对（AIA/ANDI）NetCDF、mzML/mzXML和mzData格式的文件中的LC / MS数据进行分析。为了演示目的，我们将分析参考文献[1]中的一部分数据，这个数据是小鼠敲除脂肪酸酰胺水解酶（Fatty Acid Amide Hydrolase, FAAH）基因后的代谢结果。原始数据文件（采用NetCDF格式）由faahKO数据包提供。该数据集包括来自6个基因敲除小鼠和6个野生型小鼠的脊髓样品，共12个文件。每个文件包含以centroid 模式采集，且在正离子模式（positive ion mode）下分子量在200-600 m/z范围内，保留时间在2500-4500秒的数据。
下面我们加载所有必需的包，定位faahKO包中的原始CDF文件，并构建描述实验体系的phenodata数据框。
library(xcms)
library(faahKO)
library(RColorBrewer)
library(pander)
library(magrittr)
## Get the full path to the CDF files
cdfs <- dir(system.file("cdf", package = "faahKO"), full.names = TRUE,
recursive = TRUE)
## Create a phenodata data.frame
pd <- data.frame(sample_name = sub(basename(cdfs), pattern = ".CDF",
replacement = "", fixed = TRUE),
sample_group = c(rep("KO", 6), rep("WT", 6)),
stringsAsFactors = FALSE)
随后，我们使用MSnbase包中的readMSData方法加载原始数据，并保存为OnDiskMSnExp对象。虽然MSnbase软件包最初是为蛋白质组学数据处理而开发的，但其许多功能（包括原始数据导入和数据表示）可以在代谢组学数据分析中进行共享。此外，MSnbase可用于centroid 采集模式的质谱数据（请参阅MSnbase软件包中的相应的简介）。
raw_data <- readMSData(files = cdfs, pdata = new("NAnnotatedDataFrame", pd), mode = "onDisk")
OnDiskMSnExp对象包含有关光谱数（number of spectra）、保留时间（retention times）、测量的总离子流（measured total ion current）等等数据的总体信息，但不包含完整的原始数据（比如，每个测量光谱的m/z值和强度值（m/z and intensity values））。因此，OnDiskMSnExp对象占用的内存相当小，也使其成为代表大型代谢组学实验数据的理想对象，同时仍允许执行简单的质量控制，数据检查和探索以及数据子集设置等等一系列操作。后续，根据不同的需求，需要从原始数据文件中导入m/z和强度值（m/z and intensity values）。因此，进行数据导入操作后，不应更改原始数据文件的保存位置。
3. 初步数据检查（Initial data inspection）
OnDiskMSnExp对象按照光谱（spectrum）组织MS数据，并提供intensity, mz 和 rtime方法来访问文件中的原始数据（测量的强度值，相应的m / z和保留时间值（the measured intensity values, the corresponding m/z and retention time values））。此外，spectra 方法可用于返回封装在Spectrum类中的所有数据。下面我们将从OnDiskMSnExp对象中提取保留时间值。
head(rtime(raw_data))
## F01.S0001 F01.S0002 F01.S0003 F01.S0004 F01.S0005 F01.S0006
## 2501.378 2502.943 2504.508 2506.073 2507.638 2509.203
即使实验包含多个文件/样本，所有数据都作为一维向量返回（rtime的数字向量和mz和intensity的数值向量列表，每个包含来自一个光谱的值）。fromFile函数返回一个数值向量，该向量提供全部值到原始文件的映射。下面我们使用fromFile索引按照来源文件组织mz值。
mzs <- mz(raw_data)
## Split the list by file
mzs_by_file <- split(mzs, f = fromFile(raw_data))
length(mzs_by_file)
## [1] 12
作为数据的第一次评估，我们在下面为我们实验中的每一个文件绘制基峰色谱图（base peak chromatogram, BPC) ）。我们使用chromatogram方法并设置aggregationFun参数为"max"，以便为每个光谱（spectrum ）返回最大强度（maximal intensity），并根据原始数据创建BPC。如果为了总离子流图（total ion chromatogram），我们可以将aggregationFun设置为"sum"。
## Get the base peak chromatograms. This reads data from the files.
bpis <- chromatogram(raw_data, aggregationFun = "max")
## Define colors for the two groups
group_colors <- brewer.pal(3, "Set1")[1:2]
names(group_colors) <- c("KO", "WT")
## Plot all chromatograms.
plot(bpis, col = group_colors[raw_data$sample_group])
chromatogram方法返回了一个Chromatograms对象，该对象合并各个Chromatogram对象（实际上包含色谱数据）到一个二维数组中：1）列代表样品；2）行表示m/z、保留时间范围等（m/z and/or retention time ranges）（可选）。下面我们提取第一个样品的色谱数据并访问其保留时间和强度值（retention time and intensity values）。
bpi_1 <- bpis[1, 1]
head(rtime(bpi_1))
## F01.S0001 F01.S0002 F01.S0003 F01.S0004 F01.S0005 F01.S0006
## 2501.378 2502.943 2504.508 2506.073 2507.638 2509.203
head(intensity(bpi_1))
## F01.S0001 F01.S0002 F01.S0003 F01.S0004 F01.S0005 F01.S0006
## 43888 43960 43392 42632 42200 42288
chromatogram方法还支持从MS数据的m/z-rt切片（m/z-rt slice of the MS data）中提取数据。在下一节中，我们将使用此方法为选定的峰创建提取离子色谱图（extracted ion chromatogram, EIC）。
请注意：chromatogram方法从每个文件中读取原始数据以计算色谱。另一方面，bpi和tic方法不从原始文件中读取任何数据，而是使用输入文件的头定义（header definition）中提供的相应信息（可能与实际数据不同）。
下面我们创建表示每个文件的总离子流分布的箱线图。这些图对于发现有问题或运行失败的MS样本非常有用。
## Get the total ion current by file
tc <- split(tic(raw_data), f = fromFile(raw_data))
boxplot(tc, col = group_colors[raw_data$sample_group],
ylab = "intensity", main = "Total ion current")
4. 色谱峰检测（Chromatographic peak detection）
接下来，我们使用 centWave 算法[2]进行色谱峰检测。然而，在进行峰值检测之前，建议在对数据进行可视化查看，例如：提取内标或已知化合物的离子色谱图，用于评估和调整峰检测设置，因为默认设置不适合大多数LCMS实验。
centWave的两个最关键参数是peakwidth（峰宽）和ppm（centroids检测模式下，化合物分子量m/z值（m/z values of centroids）的最大预期偏差。而实验中的ppm值远大于制造商指定的ppm。为了评估特定色谱的峰宽，我们绘制了一个峰的提取离子色谱图（extracted ion chromatogram, EIC）。
## Define the rt and m/z range of the peak area
rtr <- c(2700, 2900)
mzr <- c(334.9, 335.1)
## extract the chromatogram
chr_raw <- chromatogram(raw_data, mz = mzr, rt = rtr)
plot(chr_raw, col = group_colors[chr_raw$sample_group])
注意：如果在某个扫描中（例如：在特定的保留时间内），在相应的 mz 范围内没有测量到信号，则由chromatogram方法提取的Chromatogram对象中会包含一些NA值。而这个现象可以通过上述图中未绘制为连续线的线条反映出来。
上面图中的峰值宽度约为50秒。应该设置peakwidth参数，以适应数据集中的预期宽度。我们将其设置为20,80，以适应当前的示例数据集。
而对于ppm参数，我们提取对应于上述峰的完整MS数据（强度，保留时间和m/z值（intensity, retention time and m/z values））。为此，我们首先按保留时间过滤原始数据对象（raw_data），然后按照m/z过滤，最后将筛选的数据对象按照type ="XIC"进行绘制，以生成下图。我们使用管道命令（%>%）来说明相应的工作流程。
raw_data %>%
filterRt(rt = rtr) %>%
filterMz(mz = mzr) %>%
plot(type = "XIC")
图解：对上述12个图中的每一个图：1）上半部分：横坐标是保留时间（retention time），纵坐标是对应的强度（intensity values）；2）下半部分：横坐标是保留时间（retention time），纵坐标是对应的 m/z ；3）图中的每一个点根据强度（intensity values）进行着色。
在目前的数据中，在m/z值中实际上没有变化。通常人们会看到m/z值（图下半部分）散布在该化合物的实际m/z值附近。建议对许多化合物（在样品中已知的内部标准或化合物）的m/z值的范围进行检查，并根据这些值定义 centWave 的ppm参数。
下面我们使用findchrompeak方法来执行色谱峰检测。提交的参数对象，将决定使用哪种算法，并设置该算法的相关参数设置。请注意，我们将参数``noise设置为1000`，从而只考虑在峰值检测步骤中值大于1000的信号，以此略微加快分析速度。
# define parameter object
cwp <- CentWaveParam(peakwidth = c(30, 80), noise = 1000)
# perform the chromatographic peak detection
xdata <- findChromPeaks(raw_data, param = cwp)
结果作为XCMSnExp对象返回，该对象通过存储LC/GC-MS 预处理结果来扩展OnDiskMSnExp对象。这也意味着，所有用于数据子集、数据过滤或者访问（原始）数据的方法都是从OnDiskMSnExp对象继承的。色谱峰检测的结果可以用chromPeaks 方法进行访问。
head(chromPeaks(xdata))
## mz mzmin mzmax rt rtmin rtmax into intb
## [1,] 236.1 236.1 236.1 2520.158 2501.378 2553.022 301289.53 268211.58
## [2,] 362.9 362.9 362.9 2596.840 2587.451 2604.665 19916.60 19674.79
## [3,] 302.0 302.0 302.0 2617.185 2587.451 2648.484 699458.52 687162.27
## [4,] 316.2 316.2 316.2 2668.828 2661.003 2676.653 51310.09 50854.60
## [5,] 370.1 370.1 370.1 2673.523 2643.789 2706.387 458247.09 423667.40
## [6,] 360.0 360.0 360.0 2682.913 2635.964 2718.907 9034381.61 8518838.59
## maxo sn sample is_filled
## [1,] 12957 13 1 0
## [2,] 1453 12 1 0
## [3,] 30552 73 1 0
## [4,] 4267 17 1 0
## [5,] 25672 16 1 0
## [6,] 317568 20 1 0
返回的matrix 为每个确定的色谱峰提供了m/z、保留时间范围、积分后的信号强度(“into”)、最大信号强度(“maxo”)。“样本”列包含了该项目中鉴定出峰的所有样品。
下面我们使用来自该表的数据来计算一些参数的总值。
summary_fun <- function(z) {
c(peak_count = nrow(z), rt = quantile(z[, "rtmax"] - z[, "rtmin"]))
}
T <- lapply(split.data.frame(chromPeaks(xdata),
f = chromPeaks(xdata)[, "sample"]),
FUN = summary_fun)
T <- do.call(rbind, T)
rownames(T) <- basename(fileNames(xdata))
pandoc.table(T,
caption = paste0("Summary statistics on identified chromatographic",
" peaks. Shown are number of identified peaks per",
" sample and widths/duration of chromatographic ",
"peaks."))
已确定的色谱峰的汇总统计数据
显示了每个样品的鉴定峰数和色谱峰的宽度/保留时间范围。
| peak_count | rt.0% | rt.25% | rt.50% | rt.75% | rt.100% | |
|---|---|---|---|---|---|---|
| ko15.CDF | 374 | 14.08 | 50.47 | 57.9 | 67.29 | 169 |
| ko16.CDF | 489 | 14.08 | 48.51 | 59.47 | 67.29 | 169 |
| ko18.CDF | 342 | 14.09 | 51.64 | 61.03 | 68.86 | 223.8 |
| ko19.CDF | 259 | 14.08 | 53.21 | 64.16 | 71.99 | 173.7 |
| ko21.CDF | 222 | 14.09 | 56.34 | 65.73 | 73.16 | 164.3 |
| ko22.CDF | 255 | 14.08 | 51.64 | 62.6 | 71.99 | 134.6 |
| wt15.CDF | 362 | 14.08 | 50.08 | 57.9 | 67.29 | 139.3 |
| wt16.CDF | 352 | 14.09 | 53.21 | 61.03 | 68.86 | 133 |
| wt18.CDF | 302 | 14.08 | 53.21 | 61.03 | 70.42 | 192.5 |
| wt19.CDF | 274 | 14.08 | 56.34 | 64.16 | 75.12 | 184.7 |
| wt21.CDF | 229 | 15.65 | 54.77 | 64.16 | 73.55 | 136.2 |
| wt22.CDF | 337 | 14.09 | 48.51 | 61.03 | 70.42 | 134.6 |
我们还可以用plotChromPeaks 函数在m/z - retention time空间中绘制出确定的色谱峰的位置。下面我们画出第三个样本的色谱峰。
plotChromPeaks(xdata, file = 3)
为了获得峰值检测的全局概览，我们可以在保留时间轴（retention time axis）上绘制每个文件的确定峰值的频率。这允许在MS运行中识别时间周期，在这个过程中，可以识别出更多数量的峰值，并评估这是否在整个文件中是一致的。
plotChromPeakImage(xdata)
图解：在频率中，用黄色白色代表高频率。每一行显示一个文件的峰值频率。
接下来，我们将之前已识别色谱峰作为示例峰突出显示。根据一系列已知峰或内标的峰进行作图，并进行评估，将有助于确保峰值检测的参数设置是合适的，并正确地识别出预期的峰值。
plot(chr_raw, col = group_colors[chr_raw$sample_group], lwd = 2)
highlightChromPeaks(xdata, border = group_colors[chr_raw$sample_group],
lty = 3, rt = rtr, mz = mzr)
图解：红色和蓝色分别代表敲除和野生型样本。这些矩形表示每个样本中确定的色谱峰。
请注意，我们还可以通过在chromPeaks方法中使用mz和rt参数提供相应的m/z和保留时间范围来专门提取所选区域的已识别色谱峰。
pander(chromPeaks(xdata, mz = mzr, rt = rtr),
caption = paste("Identified chromatographic peaks in a selected ",
"m/z and retention time range."))
在选定的m/z和保留时间范围内确定的色谱峰
| mz | mzmin | mzmax | rt | rtmin | rtmax | into | intb | maxo | sn | sample | is_filled |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 335 | 335 | 335 | 2783 | 2756 | 2817 | 421286 | 392692 | 16856 | 26 | 1 | 0 |
| 335 | 335 | 335 | 2788 | 2756 | 2821 | 1529680 | 1490861 | 58736 | 77 | 2 | 0 |
| 335 | 335 | 335 | 2788 | 2758 | 2822 | 221355 | 204694 | 8158 | 19 | 3 | 0 |
| 335 | 335 | 335 | 2786 | 2756 | 2821 | 299084 | 281522 | 9154 | 22 | 4 | 0 |
| 335 | 335 | 335 | 2799 | 2758 | 2838 | 174587 | 160226 | 6295 | 13 | 5 | 0 |
| 335 | 335 | 335 | 2788 | 2756 | 2822 | 954335 | 933983 | 35856 | 76 | 8 | 0 |
| 335 | 335 | 335 | 2791 | 2758 | 2827 | 1668431 | 1635820 | 54640 | 94 | 9 | 0 |
| 335 | 335 | 335 | 2786 | 2756 | 2810 | 644912 | 632013 | 20672 | 54 | 10 | 0 |
| 335 | 335 | 335 | 2794 | 2769 | 2832 | 931316 | 904196 | 27200 | 49 | 11 | 0 |
最后，我们还绘制了每个样本的峰强度分布。这允许研究样品之间的峰值信号的系统差异是否存在。
## Extract a list of per-sample peak intensities (in log2 scale)
ints <- split(log2(chromPeaks(xdata)[, "into"]),
f = chromPeaks(xdata)[, "sample"])
boxplot(ints, varwidth = TRUE, col = group_colors[xdata$sample_group],
ylab = expression(log[2]~intensity), main = "Peak intensities")
grid(nx = NA, ny = NULL)
5. 校准（Alignment）
分析物在色谱中洗脱的时间可以在样品（甚至化合物）之间变化。在上一节中作为示例显示的提取离子色谱图中已经可以观察到这种差异。而校准步骤，也称为保留时间校正，旨在沿保留时间轴移动信号，以对齐同一个实验中不同样品之间的信号，从而纠正这种偏差。
目前，存在大量的对齐算法（参见[3]），而其中一些算法也在xcms中可以实现。在xcms中执行对齐/保留时间校正（alignment/retention time correction）的方法是adjustRtime，它根据提供的参数类，可以使用不同的对齐算法。在下面的例子中，我们使用 obiwarp 方法[4]来对齐样本。我们使用参数binSize = 0.6，它在0.6的mz区间中创建翘曲函数（We use a binSize = 0.6 which creates warping functions in mz bins of 0.6.）。此外，建议修改每个实验的设置，并评估保留时间校正是否确实正确校准内部对照或已知化合物。
xdata <- adjustRtime(xdata, param = ObiwarpParam(binSize = 0.6))
adjustRtime方法除了计算每个光谱的调整保留时间外，还会调整所报告的色谱峰的保留时间。
为了提取校准后的保留时间，我们可以使用adjustRtime方法，或者只使用rtime方法从XCMSnExp对象返回默认调整的保留时间（如果存在）。
## Extract adjusted retention times
head(adjustedRtime(xdata))
## F01.S0001 F01.S0002 F01.S0003 F01.S0004 F01.S0005 F01.S0006
## 2501.378 2502.958 2504.538 2506.118 2507.699 2509.280
## Or simply use the rtime method
head(rtime(xdata))
## F01.S0001 F01.S0002 F01.S0003 F01.S0004 F01.S0005 F01.S0006
## 2501.378 2502.958 2504.538 2506.118 2507.699 2509.280
通过rtime(xdata, adjusted = FALSE)命令，可以将原始的保留时间（Raw retention times）从包含校准数据的XCMSnExp对象提取出来。
为了评估校准的影响，我们在调整后的数据上绘制基峰色谱图（base peak chromatogram, BPC) ）。此外，我们使用plotAdjustedRtime函数绘制每个样本调整后的保留时间与原始的保留时间之间的差异。
## Get the base peak chromatograms.
bpis_adj <- chromatogram(xdata, aggregationFun = "max")
par(mfrow = c(2, 1), mar = c(4.5, 4.2, 1, 0.5))
plot(bpis_adj, col = group_colors[bpis_adj$sample_group])
## Plot also the difference of adjusted to raw retention time.
plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group])
图解：上面的图是校准后的基峰色谱图（顶部）；下面的图是调整和原始保留时间之间的差异（底部）
调整后和原始保留时间之间的差异太大可能表明样品表现不佳或对齐。
或者，我们可以使用峰组对齐方法（peak groups alignment method ），通过对齐先前识别的钩峰（hook peaks）（大多数/所有样品中存在的色谱峰）来调整保留时间。理想情况下，这些钩峰（hook peaks）应该跨越保留时间范围的大部分时间。下面我们首先使用dropAdjustedRtime方法恢复原始保留时间（也是已识别峰的保留时间）。请注意，每个预处理步骤都有一个drop *方法，允许从XCMSnExp对象中删除相应处理的结果。
## Does the object have adjusted retention times?
hasAdjustedRtime(xdata)
## [1] TRUE
## Drop the alignment results.
xdata <- dropAdjustedRtime(xdata)
## Does the object have adjusted retention times?
hasAdjustedRtime(xdata)
## [1] FALSE
注意：XCMSnExp对象保持原始数据以及调整后的保留时间，在大多数情况下，子集将去除调整后的保留时间。因此，将原始保留时间替换为调整后的保留时间，有时是可能是有用的。而这个操作可以通过applyAdjustedRtime来完成。
如上所述，峰组方法（peak groups method ）需要大多数样品中存在的峰组（特征）来进行比对。因此，我们必须执行第一次对应运行（correspondence run）来识别这些峰值（有关所用算法的详细信息将在下一节中介绍）。我们在这里再次使用默认设置，但强烈建议为每个数据集调整参数。样本组的定义（比如，将实验中的单个样品分配给样品组）对于PeakDensityParam方法是强制性的。如果实验中没有样本组，sampleGroups 方法则应将每个文件设置一个单一的价值（例如，rep(1, length(fileNames(xdata))）。
## Correspondence: group peaks across samples.
pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,
minFraction = 0.8)
xdata <- groupChromPeaks(xdata, param = pdp)
## Now the retention time correction.
pgp <- PeakGroupsParam(minFraction = 0.85)
## Get the peak groups that would be used for alignment.
xdata <- adjustRtime(xdata, param = pgp)
另请注意：我们可以在校准之前对对象使用adjustRtimePeakGroups方法，以评估将执行校准的特征（峰组，peak groups）。这对于测试峰组算法的不同设置非常有用。此外，可以手动选择或定义某些峰组（即每个样品的保留时间），并使用peakGroupsMatrix方法将此矩阵提供给PeakGroupsParam类。
下图使用plotAdjustedRtime函数绘制原始和调整的保留时间之间的差异，如果使用峰组方法（peak groups method）进行对齐，则还会突出显示调整中使用的峰组。
## Plot the difference of adjusted to raw retention time.
plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group],
peakGroupsCol = "grey", peakGroupsPch = 1)
最后，我们评估了校准对测试峰值的影响。
par(mfrow = c(2, 1))
## Plot the raw data
plot(chr_raw, col = group_colors[chr_raw$sample_group])
## Extract the chromatogram from the adjusted object
chr_adj <- chromatogram(xdata, rt = rtr, mz = mzr)
plot(chr_adj, col = group_colors[chr_raw$sample_group])
6. 组装（Correspondence）
代谢组学预处理的最后一步是 correspondence，将匹配样品之间检测到的色谱峰（根据设置，如果合适，也包括样本中色谱峰匹配）。在xcms中执行 correspondence的方法是groupChromPeaks。我们将使用峰值密度法（peak density method）[5]对色谱峰进行分组。而该算法合并色谱法主要依赖于沿着 mz 维度的小切片中沿保留时间轴的峰密度（The algorithm combines chromatographic peaks depending on the density of peaks along the retention time axis within small slices along the mz dimension.）。为了说明这一点，我们在下面绘制每一个样本中具有多个色谱峰的mz切片的色谱图。下面，我们将minFraction参数设置为0.4，这样如果色谱峰存在于每个样品组中至少40％的样品中，则被归类为一个特征。使用sampleGroups参数指定样本组分配。
## Define the mz slice.
mzr <- c(305.05, 305.15)
## Extract and plot the chromatograms
chr_mzr <- chromatogram(xdata, mz = mzr, rt = c(2500, 4000))
par(mfrow = c(3, 1), mar = c(1, 4, 1, 0.5))
cols <- group_colors[chr_mzr$sample_group]
plot(chr_mzr, col = cols, xaxt = "n", xlab = "")
## Highlight the detected peaks in that region.
highlightChromPeaks(xdata, mz = mzr, col = cols, type = "point", pch = 16)
## Define the parameters for the peak density method
pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,
minFraction = 0.4, bw = 30)
par(mar = c(4, 4, 1, 0.5))
plotChromPeakDensity(xdata, mz = mzr, col = cols, param = pdp,
pch = 16, xlim = c(2500, 4000))
## Use a different bw
pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,
minFraction = 0.4, bw = 20)
plotChromPeakDensity(xdata, mz = mzr, col = cols, param = pdp,
pch = 16, xlim = c(2500, 4000))
图解：1）上图：具有多个色谱峰的mz切片的色谱图；2）中图和下图：针对不同值的bw参数，在其保留时间（x轴）和实验样本中的指数（y轴）上鉴定的色谱峰；3）黑线表示峰值密度估计值，峰的分组（基于提供的设置）由灰色矩形表示。（Upper panel: chromatogram for an mz slice with multiple chromatographic peaks. Middle and lower panel: identified chromatographic peaks at their retention time (x-axis) and index within samples of the experiments (y-axis) for different values of the bw parameter. The black line represents the peak density estimate. Grouping of peaks (based on the provided settings) is indicated by grey rectangles.）
上图中上面部分的图显示了每个样品的提取离子色谱图，其中突出显示了检测到的峰。中间和下部图显示了不同样品中每个检测到的峰的保留时间。黑色实线表示沿保留时间检测到的峰的密度分布。结合到特征（峰组）中的峰用灰色矩形表示。在PeakDensityParam方法中设置不同的bw参数值：中间图是bw = 30，下面的图是bw = 20。使用默认的bw参数值，两个相邻的色谱峰将被分组到相同的特征中，而当bw为20时，它们将被分组为单独的特征。这种分组情况主要取决于密度函数的参数以及使用PeakDensityParam方法传递给算法的其他参数。
## Perform the correspondence
pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,
minFraction = 0.4, bw = 20)
xdata <- groupChromPeaks(xdata, param = pdp)
可以使用featureDefinitions方法提取 correspondence 的结果，该方法返回一个DataFrame，其中包含了定义的特征（即，mz 和 rt 范围；并且在列peakidx中，返回每个特征的在chromPeaks矩阵中的色谱峰的索引）。
## Extract the feature definitions
featureDefinitions(xdata)
## DataFrame with 312 rows and 10 columns
## mzmed mzmin mzmax rtmed
## <numeric> <numeric> <numeric> <numeric>
## FT001 205 205 205 2790.91667807044
## FT002 206 206 206 2790.76510706472
## FT003 207.100006103516 207.100006103516 207.100006103516 2720.76064413194
## FT004 236.100006103516 236.100006103516 236.100006103516 2526.73236656529
## FT005 241.100006103516 241.100006103516 241.199996948242 3678.87576330383
## ... ... ... ... ...
## FT308 594.200012207031 594.200012207031 594.200012207031 3405.34666745362
## FT309 594.299987792969 594.299987792969 594.299987792969 3613.40135845759
## FT310 595.200012207031 595.200012207031 595.299987792969 3004.51480319748
## FT311 596.299987792969 596.299987792969 596.400024414062 3818.52945007536
## FT312 597.400024414062 597.299987792969 597.400024414062 3814.92705405312
## rtmin rtmax npeaks KO WT
## <numeric> <numeric> <numeric> <numeric> <numeric>
## FT001 2790.25914950764 2791.33708124005 12 6 6
## FT002 2787.95417597022 2794.07524052365 11 4 6
## FT003 2718.2665593745 2723.95355653805 8 4 4
## FT004 2523.80327395574 2527.55945827299 5 2 3
## FT005 3664.1493202833 3683.06286789441 9 4 2
## ... ... ... ... ... ...
## FT308 3395.29146342443 3410.28858323004 3 0 3
## FT309 3608.97879037478 3619.97980830287 6 3 2
## FT310 3000.07972689414 3019.16896250021 7 2 4
## FT311 3813.57966000225 3826.51868541398 7 4 2
## FT312 3809.07197917863 3818.5741287851 5 1 3
## peakidx
## <list>
## FT001 c(43, 416, 912, 1242, 1500, 1727, 1978, 2350, 2705, 3000, 3270, 3508)
## FT002 c(25, 1232, 1244, 1495, 1724, 1968, 2341, 2693, 2996, 3267, 3501)
## FT003 c(399, 881, 1482, 1707, 1962, 2327, 2673, 3482)
## FT004 c(1, 376, 1943, 2304, 3235)
## FT005 c(253, 1122, 1397, 1406, 1880, 1887, 2540, 3171, 3184)
## ... ...
## FT308 c(2797, 3095, 3638)
## FT309 c(217, 1605, 1865, 3390, 3396, 3698)
## FT310 c(63, 495, 2026, 2735, 3293, 3538, 3546)
## FT311 c(306, 1145, 1147, 1656, 1910, 2578, 3768)
## FT312 c(1150, 2917, 2919, 3440, 3765)
featureValues方法返回一个matrix，其中行是 features 和列是样本。可以使用value参数定义此矩阵的内容。通过设置value = "into"，将返回一个矩阵，其中包含与样本中特征相对应的峰值的积分信号（integrated signal of the peaks）。chromPeaks矩阵的任何列名都可以传递给value。下面我们提取每个特征/样本的积分峰强度（integrated peak intensity per feature/sample）。
## Extract the into column for each feature.
head(featureValues(xdata, value = "into"))
## ko15.CDF ko16.CDF ko18.CDF ko19.CDF ko21.CDF ko22.CDF wt15.CDF
## FT001 2024095.6 1832137.0 1756373.2 1275015.4 1558397.4 1274955.4 2220917.7
## FT002 229446.1 NA NA 148895.1 169782.6 181300.6 271403.0
## FT003 NA 470042.7 349545.6 NA 356775.9 240606.3 370510.3
## FT004 301289.5 250053.3 NA NA NA NA 293360.4
## FT005 784970.1 NA 1377891.5 659602.9 NA 528679.9 NA
## FT006 681973.4 266391.1 103640.3 NA NA 132666.7 682289.2
## wt16.CDF wt18.CDF wt19.CDF wt21.CDF wt22.CDF
## FT001 1717218.1 1871280.3 1551566.88 1713396.7 1438088.9
## FT002 261254.3 237429.1 225250.02 244356.5 188864.7
## FT003 388938.3 312991.1 NA NA 251502.7
## FT004 190338.9 NA NA 305942.5 NA
## FT005 1698777.3 NA 704491.40 NA NA
## FT006 NA 186512.1 25193.37 NA 1462255.2
此特征矩阵（feature matrix）包含NA，表明对应样本在特征（feature ）的m/z-rt区域中未检测到色谱峰。虽然在许多情况下，相应区域可能确实没有峰值信号，但也可能存在信号，但峰值检测算法未能检测到色谱峰。xcms提供了一个fillChromPeaks方法，根据原始文件，为此类缺失值填充强度数据（fill in intensity data）。填充的峰（filled in peaks）将添加到chromPeaks矩阵中，并在"is_filled"列中标记为1。下面我们执行这样的缺失峰值填充（filling-in of missing peaks.）。
## Filling missing peaks using default settings. Alternatively we could
## pass a FillChromPeaksParam object to the method.
xdata <- fillChromPeaks(xdata)
head(featureValues(xdata))
## ko15.CDF ko16.CDF ko18.CDF ko19.CDF ko21.CDF ko22.CDF wt15.CDF
## FT001 43 416 912 1242 1500 1727 1978
## FT002 25 3892 3986 1232 1495 1724 1968
## FT003 3798 399 881 4069 1482 1707 1962
## FT004 1 376 3987 4070 4185 4310 1943
## FT005 253 3893 1122 1397 4186 1880 4415
## FT006 46 424 916 4071 4187 1728 1984
## wt16.CDF wt18.CDF wt19.CDF wt21.CDF wt22.CDF
## FT001 2350 2705 3000 3270 3508
## FT002 2341 2693 2996 3267 3501
## FT003 2327 2673 4742 4858 3482
## FT004 2304 4627 4743 3235 4990
## FT005 2540 4628 3184 4859 4991
## FT006 4527 2706 3004 4860 3519
对于未在样品中检测到峰的特征（features，该方法提取特征的mz-rt区域中的所有强度（intensities），整合信号并将填充的峰（filled-in peak）添加到chromPeaks矩阵。如果一个特征的mz-rt区域的确没有信号被检测/可用，将不会加峰。并且对于这些特征，即使在其他特征填充缺失的峰值数据之后，在featureValues矩阵中依然报告为NA。
下面我们比较填充缺失值前后的缺失值的数量。我们可以使用featureValues方法中的filled 参数来定义是否也应该返回填充的峰值。
## Missing values before filling in peaks
apply(featureValues(xdata, filled = FALSE), MARGIN = 2,
FUN = function(z) sum(is.na(z)))
## ko15.CDF ko16.CDF ko18.CDF ko19.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF
## 98 97 84 118 133 111 117 105
## wt18.CDF wt19.CDF wt21.CDF wt22.CDF
## 121 120 139 113
## Missing values after filling in peaks
apply(featureValues(xdata), MARGIN = 2,
FUN = function(z) sum(is.na(z)))
## ko15.CDF ko16.CDF ko18.CDF ko19.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF
## 4 3 1 2 8 6 5 5
## wt18.CDF wt19.CDF wt21.CDF wt22.CDF
## 6 4 7 8
最后，我们进行主成分分析（principal component analysis），以评估本实验中样品的分组。请注意，我们没有执行任何数据规范化（data normalization ），因此分组可能（并且将）也会受到技术偏差的影响。
## Extract the features and log2 transform them
ft_ints <- log2(featureValues(xdata, value = "into"))
## Perform the PCA omitting all features with an NA in any of the
## samples. Also, the intensities are mean centered.
pc <- prcomp(t(na.omit(ft_ints)), center = TRUE)
## Plot the PCA
cols <- group_colors[xdata$sample_group]
pcSummary <- summary(pc)
plot(pc$x[, 1], pc$x[,2], pch = 21, main = "",
xlab = paste0("PC1: ", format(pcSummary$importance[2, 1] * 100,
digits = 3), " % variance"),
ylab = paste0("PC2: ", format(pcSummary$importance[2, 2] * 100,
digits = 3), " % variance"),
col = "darkgrey", bg = cols, cex = 2)
grid()
text(pc$x[, 1], pc$x[,2], labels = xdata$sample_name, col = "darkgrey",
pos = 3, cex = 2)
在 PC2 上，我们可以看到我们预期的KO（敲除）和WT（野生）样品之间的分离。在PC1上，样品基于它们的 ID 产生分离，ID <= 18的样品和 ID > 18的样品产生分离。这种分离可能是由技术偏差（例如在不同天/周进行的测量）或由于小鼠分析中的生物学特性引起的（比如，性别，年龄，同伴等）。
7. 进一步的数据处理和分析（Further data processing and analysis）
对特征的信号强度（features’ signal intensities）进行归一化是需要的，但是目前xcms尚未有相关方法的支持（某些方法可能会在不久的将来添加）。另外，对于例如具有显着强度/丰度差异的特征（features with significant different intensities/abundances）的鉴定，建议使用其他R包中封装的函数提供的功能，例如Bioconductor中优秀的limma包。为了支持依赖旧的xcmsSet对象的结果的其他包，可以使用xset <- as(x, "xcmsSet")命令强制将XCMSnExp对象转换为xcmsSet对象，其中x代表XCMSnExp对象。
8. 其他细节和说明（Additional details and notes）
有关与旧版本相比的新数据对象和更改/改进的详细说明，请参阅 new_functionality 介绍。
9. 评估过程历史（Evaluating the process history）
XCMSnExp对象允许捕获所有执行的预处理步骤以及@processHistory中使用的参数类。存储参数类（parameter class）也可确保最高程度的分析文档，并且将来能够替换或者重复全部或部分分析。所有执行的预处理步骤可以使用processHistory方法进行提取。
processHistory(xdata)
## [[1]]
## Object of class "XProcessHistory"
## type: Peak detection
## date: Mon Apr 30 18:48:18 2018
## info:
## fileIndex: 1,2,3,4,5,6,7,8,9,10,11,12
## Parameter class: CentWaveParam
## MS level(s) 1
##
## [[2]]
## Object of class "XProcessHistory"
## type: Peak grouping
## date: Mon Apr 30 18:49:15 2018
## info:
## fileIndex: 1,2,3,4,5,6,7,8,9,10,11,12
## Parameter class: PeakDensityParam
##
## [[3]]
## Object of class "XProcessHistory"
## type: Retention time correction
## date: Mon Apr 30 18:49:17 2018
## info:
## fileIndex: 1,2,3,4,5,6,7,8,9,10,11,12
## Parameter class: PeakGroupsParam
## MS level(s) 1
##
## [[4]]
## Object of class "XProcessHistory"
## type: Peak grouping
## date: Mon Apr 30 18:49:29 2018
## info:
## fileIndex: 1,2,3,4,5,6,7,8,9,10,11,12
## Parameter class: PeakDensityParam
##
## [[5]]
## Object of class "XProcessHistory"
## type: Missing peak filling
## date: Mon Apr 30 18:49:31 2018
## info:
## fileIndex: 1,2,3,4,5,6,7,8,9,10,11,12
## Parameter class: FillChromPeaksParam
还可以通过指定其类型（type）来提取特定处理步骤。而可用类型（Available types）可以通过processHistoryTypes函数列出。下面我们提取校准/保留时间调整步骤的参数类（the parameter class for the alignment/retention time adjustment step）。
ph <- processHistory(xdata, type = "Retention time correction")
ph
## [[1]]
## Object of class "XProcessHistory"
## type: Retention time correction
## date: Mon Apr 30 18:49:17 2018
## info:
## fileIndex: 1,2,3,4,5,6,7,8,9,10,11,12
## Parameter class: PeakGroupsParam
## MS level(s) 1
我们还可以提取此预处理步骤中使用的参数类。
## Access the parameter
processParam(ph[[1]])
## Object of class: PeakGroupsParam
## Parameters:
## minFraction: 0.85
## extraPeaks: 1
## smooth: loess
## span: 0.2
## family: gaussian
## number of peak groups: 50
10. 子集和过滤（Subsetting and filtering）
可以使用[方法或许多filter* 方法之一，对XCMSnEx对象进行子集化/过滤。所有这些方法都旨在确保返回对象中的数据是一致的。这意味着，例如，如果通过选择特定光谱（selecting specific spectra）对对象进行子集化（通过使用[方法]），则所有鉴定的色谱峰将被去除。如果对象被子集化，并仅包含来自所选文件的数据（使用filterFile方法），则移除 Correspondence 对应结果（即，识别的特征（identified features））。这是因为 correspondence 结果依赖于执行分析的文件，而在文件的子集上 correspondence 将会导致不同的结果。
作为例外，在任何子集化方法中将keepAdjustedRtime参数设置为TRUE，就能在子集化对象中强制保持调整的保留时间。
下面我们将结果对象分组为文件2和4的数据。
subs <- filterFile(xdata, file = c(2, 4))
## Do we have identified chromatographic peaks?
hasChromPeaks(subs)
## [1] TRUE
峰值检测分别在每个文件上单独执行，因此子集化对象（subsetted object ）包含来自两个文件的所有已识别的色谱峰。 但是，我们使用的是基于可用特征（available features）的保留时间调整（alignment步骤）。但是，所有的特征（All features）已删除，当然也包括调整后的保留时间（这是因为alignment步骤是基于从在所有文件上的色谱峰上确定的特征）。
## Do we still have features?
hasFeatures(subs)
## [1] FALSE
## Do we still have adjusted retention times?
hasAdjustedRtime(subs)
## [1] FALSE
但是，我们可以使用keepAdjustedRtime参数来强制保留调整后的保留时间。
subs <- filterFile(xdata, keepAdjustedRtime = TRUE)
hasAdjustedRtime(subs)
## [1] TRUE
filterRt方法可用于将对象子集化为特定保留时间范围内的光谱。
subs <- filterRt(xdata, rt = c(3000, 3500))
range(rtime(subs))
## [1] 3000.011 3499.909
按保留时间过滤不会改变/影响调整的保留时间（同样，如果调整的保留时间存在，则会根据调整的保留时间执行过滤）。
hasAdjustedRtime(subs)
## [1] TRUE
此外，我们可以获得在指定的保留时间范围内所有确定的色谱峰：
hasChromPeaks(subs)
## [1] TRUE
range(chromPeaks(subs)[, "rt"])
## [1] 3000.011 3499.909
在R中对任何对象进行子集的最自然方式是[。将[用于XCMSnExp对象子集上，它只会保留选定的光谱。因此，在[中使用的索引i是一个在1和和光谱总数（跨所有文件）之间的整数。下面，我们使用[和filterFile来对xdata进行子集化，以保留一个文件中的所有光谱。
## Extract all data from the 3rd file.
one_file <- filterFile(xdata, file = 3)
one_file_2 <- xdata[fromFile(xdata) == 3]
## Is the content the same?
all.equal(spectra(one_file), spectra(one_file_2))
## [1] TRUE
虽然两个对象的光谱含量相同，但one_file包含已识别的色谱峰，而one_file_2则不包含。因此，在大多数情况下，在用这两个过滤函数进行进行子集化的时候，优先于使用[。
## Subsetting with filterFile preserves chromatographic peaks
head(chromPeaks(one_file))
## mz mzmin mzmax rt rtmin rtmax into intb maxo sn
## [1,] 402.1 402.1 402.1 2515.464 2501.379 2548.328 213455.2 185881.2 8038 11
## [2,] 339.0 339.0 339.0 2515.464 2501.379 2540.503 162264.4 158324.9 5638 14
## [3,] 354.0 354.0 354.0 2515.464 2501.379 2532.678 173130.4 165213.5 6334 11
## [4,] 352.9 352.9 352.9 2517.029 2501.379 2534.243 525993.0 522760.1 20488 55
## [5,] 316.0 316.0 316.0 2517.029 2501.379 2535.808 741708.6 726925.4 27816 47
## [6,] 333.0 333.0 333.0 2515.464 2501.379 2545.198 968861.6 951705.4 33992 49
## sample is_filled
## [1,] 1 0
## [2,] 1 0
## [3,] 1 0
## [4,] 1 0
## [5,] 1 0
## [6,] 1 0
## Subsetting with [ not
head(chromPeaks(one_file_2))
## Warning in .local(object, ...): No chromatographic peaks available.
## NULL
注意：[也支持keepAdjustedRtime参数。
下面我们将对象子集化到光谱20:30。
subs <- xdata[20:30, keepAdjustedRtime = TRUE]
## Warning in xdata[20:30, keepAdjustedRtime = TRUE]: Removed preprocessing
## results
hasAdjustedRtime(subs)
## [1] TRUE
## Access adjusted retention times:
rtime(subs)
## F01.S0020 F01.S0021 F01.S0022 F01.S0023 F01.S0024 F01.S0025 F01.S0026
## 2538.347 2539.912 2541.477 2543.042 2544.607 2546.172 2547.737
## F01.S0027 F01.S0028 F01.S0029 F01.S0030
## 2549.302 2550.867 2552.432 2553.997
## Access raw retention times:
rtime(subs, adjusted = FALSE)
## F01.S0020 F01.S0021 F01.S0022 F01.S0023 F01.S0024 F01.S0025 F01.S0026
## 2531.112 2532.677 2534.242 2535.807 2537.372 2538.937 2540.502
## F01.S0027 F01.S0028 F01.S0029 F01.S0030
## 2542.067 2543.632 2545.197 2546.762
与MSnExp和OnDiskMSnExp对象一样，[[可用于从XCMSnExp对象中提取单个光谱对象。光谱的保留时间对应于调整的保留时间（如果存在）。
## Extract a single spectrum
xdata[[14]]
## Object of class "Spectrum1"
## Retention time: 42:9
## MSn level: 1
## Total ion count: 445
## Polarity: NA
最后，我们还可以使用split方法，该方法允许根据提供的因子f（factor f）拆分XCMSnExp。下面我们按文件拆分xdata。使用keepAdjustedRtime = TRUE可确保不会删除已调整的保留时间。
x_list <- split(xdata, f = fromFile(xdata), keepAdjustedRtime = TRUE)
lengths(x_list)
## 1 2 3 4 5 6 7 8 9 10 11 12
## 1278 1278 1278 1278 1278 1278 1278 1278 1278 1278 1278 1278
lapply(x_list, hasAdjustedRtime)
## $`1`
## [1] TRUE
##
## $`2`
## [1] TRUE
##
## $`3`
## [1] TRUE
##
## $`4`
## [1] TRUE
##
## $`5`
## [1] TRUE
##
## $`6`
## [1] TRUE
##
## $`7`
## [1] TRUE
##
## $`8`
## [1] TRUE
##
## $`9`
## [1] TRUE
##
## $`10`
## [1] TRUE
##
## $`11`
## [1] TRUE
##
## $`12`
## [1] TRUE
注意：对于该操作，还有一个专用的splitByFile方法，它在内部使用filterFile方法，因此，例如，不能去除已识别的色谱峰。该方法尚不支持keepAdjustedRtime参数，因此默认情况下会删除调整后的保留时间。
xdata_by_file <- splitByFile(xdata, f = factor(1:length(fileNames(xdata))))
lapply(xdata_by_file, hasChromPeaks)
## $`1`
## [1] TRUE
##
## $`2`
## [1] TRUE
##
## $`3`
## [1] TRUE
##
## $`4`
## [1] TRUE
##
## $`5`
## [1] TRUE
##
## $`6`
## [1] TRUE
##
## $`7`
## [1] TRUE
##
## $`8`
## [1] TRUE
##
## $`9`
## [1] TRUE
##
## $`10`
## [1] TRUE
##
## $`11`
## [1] TRUE
##
## $`12`
## [1] TRUE
11. 并行处理（Parallel processing）
xcms中的大多数方法都支持并行处理。并行处理由来自 Bioconductor 的BiocParallel 包处理和配置，并可以在R会话中进行全局定义。
基于Unix的系统（Linux，macOS）支持基于multicore的并行处理。要在全局配置它，我们需要register 方法注册参数类。另请注意，下面使用bpstart初始化并行进程。
register(bpstart(MulticoreParam(2)))
Windows仅支持基于 socket 的并行处理：
register(bpstart(SnowParam(2)))
请注意，基于multicore的并行处理可能是错误的或在macOS上失败。如果是这样，可以使用DoparParam代替（需要foreach包）。
有关其他选项和详细信息，请参阅BiocParallel软件包中的介绍。
12. 参考文献（References）
- Saghatelian A, Trauger SA, Want EJ, Hawkins EG, Siuzdak G, Cravatt BF: Assignment of endogenous substrates to enzymes by global metabolite profiling. Biochemistry 2004, 43:14332–9
- Tautenhahn R, Böttcher C, Neumann S: Highly sensitive feature detection for high resolution LC/MS. BMC Bioinformatics 2008, 9:504.
- Smith R, Ventura D, Prince JT: LC-MS alignment in theory and practice: a comprehensive algorithmic review. Briefings in bioinformatics 2013, 16:bbt080–117.
- Prince JT, Marcotte EM: Chromatographic alignment of ESI-LC-MS proteomics data sets by ordered bijective interpolated warping. Analytical chemistry 2006, 78:6140–6152.
- Smith CA, Want EJ, O’Maille G, Abagyan R, Siuzdak G: XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification. Analytical chemistry 2006, 78:779–787.
更新时间：2019-05-23 10:54:23
本文由 石九流 创作，如果您觉得本文不错，请随意赞赏
采用 知识共享署名4.0 国际许可协议进行许可
本站文章除注明转载/出处外，均为本站原创或翻译，转载前请务必署名
原文链接：https://blog.computsystmed.com/archives/lcms-data-preprocessing-and-analysis-with-xcms
最后更新：2019-05-23 10:54:23
