介绍
本文档描述了 xcms（version>=3）的 LCMS 实验的数据导入、探索、预处理和分析。示例和基本工作流程改编自 Colin A.Smith 的原始 LC/MS 预处理和分析。
新版本使用 XCMSnExpe 对象（而不是旧的 xcmsSet 对象）作为预处理结果的容器。然而，为了支持依赖于 xcmsSet 对象的包和管道，可以使用 as 方法（即 xset<-as (x, "xcmsSet")）将 XCMSnExpe 转换为 xcmsSet 对象，其中 x 是 XCMSnxp 对象。
安装
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("xcms")导入数据
xcms 支持（AIA/ANDI）NetCDF、mzXML 和 mzML 格式的 LC/MS 数据。对于实际数据导入，使用 Bioconductor 的 mzR。出于演示目的，我们将分析 [1] 中的数据子集，其研究了敲除脂肪酸酰胺水解酶（FAAH）基因在小鼠中的代谢后果。原始数据文件（NetCDF 格式）随 faahKO 数据包提供。数据集包括来自 6 只敲除小鼠和 6 只野生型小鼠脊髓的样本。每个文件包含在正离子模式下以 200-600 m/z 和 2500-4500 秒采集的质心模式（centroid mode）数据。为了加快处理速度，我们仅分析 8 个文件，且限定保留时间范围为 2500 到 3500 秒。
下面我们加载所有必需的包，在 faahKO 包中找到原始 CDF 文件，并构建描述实验设置的表型（phenodata）数据框。注意，对于实际实验，建议定义一个文件（表），其中包含原始数据文件的文件名以及每个文件的样本描述，作为附加列。然后，可以使用例如 read.table 导入该文件作为变量 pd（而不是在下面的示例中在 R 中定义），并且可以将文件名传递给下面的 readMSData 函数，例如 files=paste0(MZML_PATH, "/", pd$MZML_file) ，其中 MZML_PATH 是文件所在目录的路径， "mzMLfile" 是 phenodata 文件中包含文件名那一列的列名。
library(xcms)
library(faahKO)
library(RColorBrewer)
library(pander)
library(magrittr)
library(pheatmap)
library(SummarizedExperiment)
## Get the full path to the CDF files
cdfs <- dir(system.file("cdf", package = "faahKO"), full.names = TRUE, recursive = TRUE)[c(1, 2, 5, 6, 7, 8, 11, 12)]
## Create a phenodata data.frame
pd <- data.frame(sample_name = sub(basename(cdfs), pattern = ".CDF", replacement = "", fixed = TRUE), sample_group = c(rep("KO", 4), rep("WT", 4)), stringsAsFactors = FALSE)随后，我们使用 MSnbase 包中的 readMSData 方法将原始数据加载为 OnDiskMSnExpe 对象。MSnbase 为质谱数据处理提供了基础结构和基础设施。此外，MSnbase 可用于形心剖面模式（centroid profile-mode）MS 数据（参见 MSnbase 包中的相应介绍）。
raw_data <- readMSData(files = cdfs, pdata = new("NAnnotatedDataFrame", pd), mode = "onDisk")接下来，我们将数据集限制在 2500 到 3500 秒的保留时间范围内。这仅仅是为了减少该示例的处理时间。
raw_data <- filterRt(raw_data, c(2500, 3500))生成的 OnDiskMSnExpe 对象包含关于 spectrum 数量、保留时间、测量的总离子电流等的一般信息，但不包含完整的原始数据（即每个测量光谱的 m/z 和强度值）。因此，它的内存占用相当小，是代表大型代谢组学实验的理想对象，同时允许执行简单的质量控制、数据检查和探索以及数据取子集操作。m/z 和强度值根据需要从原始数据文件中导入，因此原始数据文件的位置在初始数据导入后不应更改。
初始数据检查
OnDiskMSnExp 按 spectrum 组织 MS 数据，并提供 intensity 、 mz 和 rtime 方法以访问文件中的原始数据（测量的强度值、相应的 m/z 和保留时间值）。此外， spectra 方法可用于返回封装在 spectrum 对象中的所有数据。下面我们从对象中提取保留时间值。
head(rtime(raw_data))## F1.S0001 F1.S0002 F1.S0003 F1.S0004 F1.S0005 F1.S0006
## 2501.378 2502.943 2504.508 2506.073 2507.638 2509.203即使实验由多个文件 / 样本组成，所有数据都会作为一维矢量（one-dimensional vectors，rtime 的一个数值型向量，一个向量列表 for mz 和 intensity，每个各自都包含一个 spectrum 的值）返回。 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] 8作为数据的第一次评估，我们在实验中绘制了每个文件的基峰色谱图（base peak chromatogram, BPC）。我们使用 chromatogram 方法并将 aggregationFun 设置为 "max" ，以返回每个 spectrum 的最大强度，从而根据原始数据创建 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 <- paste0(brewer.pal(3, "Set1")[1:2], "60")
names(group_colors) <- c("KO", "WT")
## Plot all chromatograms.
png(file="xcms_1.png")
plot(bpis, col = group_colors[raw_data$sample_group])
dev.off()
chromatogram 返回了一个 MCchromatograms 对象，该对象以二维矩阵形式组织各个 Chromatogram 对象（实际上包含 chromatographic 数据）：列表示样本，行（可选）表示 m/z 和 / 或保留时间范围。下面我们提取第一个样品的 chromatogram，并获取其保留时间和强度值。
bpi_1 <- bpis[1, 1]
head(rtime(bpi_1))## F1.S0001 F1.S0002 F1.S0003 F1.S0004 F1.S0005 F1.S0006
## 2501.378 2502.943 2504.508 2506.073 2507.638 2509.203head(intensity(bpi_1))## F1.S0001 F1.S0002 F1.S0003 F1.S0004 F1.S0005 F1.S0006
## 43888 43960 43392 42632 42200 42288 chromatogram 法还支持从 MS 数据的 m/z-rt 切片中提取 chromatographic 数据。在下一节中，我们将使用此方法为选定的峰创建提取的离子色谱图（extracted ion chromatogram, EIC）。
注意， chromatogram 从每个文件读取原始数据以计算 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")
此外，我们可以根据样品的 BPC 的相似性对样品进行聚类。这也有助于在实验中发现潜在的有问题的样本，或获得实验中样本分组的初步概述。由于样本之间的保留时间不完全相同，我们使用 bin 函数将强度分组在沿保留时间轴的固定时间范围（bins）中。在本示例中，我们使用的 bin 大小为 2 秒，默认值为 0.5 秒。聚类是使用完全连锁分级聚类（complete linkage hierarchical clustering）对合并的 BPC 的成对相关性进行的。
## Bin the BPC
bpis_bin <- MSnbase::bin(bpis, binSize = 2)
## Calculate correlation on the log2 transformed base peak intensities
cormat <- cor(log2(do.call(cbind, lapply(bpis_bin, intensity))))
colnames(cormat) <- rownames(cormat) <- raw_data$sample_name
## Define which phenodata columns should be highlighted in the plot
ann <- data.frame(group = raw_data$sample_group)
rownames(ann) <- raw_data$sample_name
## Perform the cluster analysis
pheatmap(cormat, annotation = ann, annotation_color = list(group = group_colors))
样本以成对方式聚类，样本索引的 KO 和 WT 样本具有最相似的 BPC。
色谱峰检测
接下来，我们使用 centWave 算法进行色谱峰检测 [2]。然而，在运行峰值检测之前，强烈建议目视检查内参或已知化合物的提取离子色谱图，以评估和调整峰值检测设置，因为默认设置不适用于大多数 LCMS 实验。centWave 的两个最关键的参数是 peakwidth （色谱峰宽度的预期范围）和 ppm （对应于一个色谱峰的 centroids m/z 值的最大预期偏差；这通常比制造商规定的 ppm 大得多）参数。为了评估典型的色谱峰宽度，我们绘制了一个峰的 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 值）。为此，我们首先按保留时间过滤原始对象，然后按 m/z 过滤，最后用 type=“XIC” 绘制对象以生成下图。我们使用 pipe（ %>% ）命令更好地说明了相应的工作流。还应注意，在这种类型的图中，如果存在色谱峰，则默认显示色谱峰。
raw_data %>%
filterRt(rt = rtr) %>%
filterMz(mz = mzr) %>%
plot(type = "XIC")
图注：对于每个图：上面板：根据保留时间绘制强度值的色谱图；下面板：m/z 根据保留时间绘图。各个数据点根据强度进行着色。
在目前的数据中，m/z 值实际上没有变化。通常可以看到 m/z 值（下面板）散布在化合物的实际 m/z 值周围。centWave 算法的第一步基于连续扫描的 m/z 值的差异定义所谓的感兴趣区域（regions of interest, ROI）。具体而言，如果 ROI 的 m/z 和平均 m/z 之间的差小于用户定义的 ppm 参数，则来自连续扫描的 m/z 值被包括在 ROI 中。因此，ppm 的合理选择可以是来自作为色谱峰一部分的相邻扫描 / 光谱（scans/spectra）的数据点的最大 m/z 差。建议检查许多化合物（内标或样品中已知存在的化合物）的 m/z 值范围，并根据这些确定 centWave 的 ppm 参数。
注意，我们也可以对提取的离子色谱图（ion chromatogram）进行峰检测。这有助于评估不同的峰值检测设置。只需注意，提取的离子色谱图上的峰值检测不会考虑 ppm 参数，并且背景信号的估计与完整数据集上的峰值探测不同；因此， snthresh 的值将具有不同的结果。下面我们使用 findChromPeaks 函数对提取的离子色谱图进行峰值检测。提交的参数对象定义将使用的算法，并允许定义此算法的设置。我们使用带有默认设置的 centWave 算法， snthresh 除外。
xchr <- findChromPeaks(chr_raw, param = CentWaveParam(snthresh = 2))我们可以使用 chromPeaks 函数访问已识别的色谱峰。
head(chromPeaks(xchr))## rt rtmin rtmax into intb maxo sn row column
## [1,] 2781.505 2761.160 2809.674 412134.25 355516.37 16856 13 1 1
## [2,] 2786.199 2764.290 2812.803 1496244.21 1391821.33 58736 20 1 2
## [3,] 2734.556 2714.211 2765.855 21579.37 18449.43 899 4 1 3
## [4,] 2797.154 2775.245 2815.933 159058.78 150289.31 6295 12 1 3
## [5,] 2784.635 2761.160 2808.109 54947.54 37923.53 2715 2 1 4
## [6,] 2859.752 2845.668 2878.532 13895.21 13874.87 905 904 1 4与 chromPeaks 矩阵相似，还有一个数据框 chromPeakData ，允许为每个色谱峰添加任意注释。下面我们提取这个数据框，默认情况下，该数据框只包含识别峰值的 MS 水平。
chromPeakData(xchr)## DataFrame with 12 rows and 4 columns
## ms_level is_filled row column
## <integer> <logical> <integer> <integer>
## 1 1 FALSE 1 1
## 2 1 FALSE 1 2
## 3 1 FALSE 1 3
## 4 1 FALSE 1 3
## 5 1 FALSE 1 4
## ... ... ... ... ...
## 8 1 FALSE 1 4
## 9 1 FALSE 1 5
## 10 1 FALSE 1 6
## 11 1 FALSE 1 7
## 12 1 FALSE 1 8接下来，我们在提取的离子色谱图中绘制已识别的色谱峰。我们使用 col 参数为各个色谱线上色。还可以为识别的峰值指定颜色， peakCol 为前景 / 边框颜色， peakBg 为背景 / 填充颜色。必须为 chromPeaks 列出的每个色谱峰提供一种颜色。下面我们定义了一种颜色，以指示样本所在的样本组，并使用峰 “sample” 列中的样本信息为每个色谱峰分配正确的颜色。更多峰值突出显示选项将在下面进一步描述。
sample_colors <- group_colors[xchr$sample_group]
plot(xchr, col = sample_colors,
peakBg = sample_colors[chromPeaks(xchr)[, "column"]])
图注：红色和蓝色分别代表 KO 和野生型样本。已识别色谱峰的峰面积用样品组颜色突出显示。
最后，我们对整个数据集进行色谱峰检测。请注意，我们将参数 prefilter 设置为 c(6, 5000) ，将 noise 设置为 5000 ，以减少该示例的运行时间。使用此设置，我们只考虑峰值检测步骤中值大于 5000 的信号。
cwp <- CentWaveParam(peakwidth = c(20, 80), noise = 5000,
prefilter = c(6, 5000))
xdata <- findChromPeaks(raw_data, param = cwp)结果作为 XCMSnExp 对象返回，该对象通过存储 LC/GC-MS 预处理结果来扩展 OnDiskMSnExp 对象。这也意味着，所有设置和过滤数据或访问（原始）数据的方法都是从 OnDiskMSnExp 对象继承的，因此可以重新使用。还请注意，通过使用参数 add＝TRUE 调用 findChromPeaks ，可以对 xdata 对象执行额外的峰值检测（例如，在 MS level＞1 的数据上）。
色谱峰检测的结果可通过 chromPeaks 方法获得。
head(chromPeaks(xdata))## mz mzmin mzmax rt rtmin rtmax into intb maxo sn sample
## CP0001 453.2 453.2 453.2 2509.203 2501.378 2527.982 1007409.0 1007380.8 38152 38151 1
## CP0002 236.1 236.1 236.1 2518.593 2501.378 2537.372 253501.0 226896.3 12957 11 1
## CP0003 594.0 594.0 594.0 2601.535 2581.191 2637.529 161042.2 149297.3 7850 13 1
## CP0004 577.0 577.0 577.0 2604.665 2581.191 2626.574 136105.2 129195.5 6215 13 1
## CP0005 369.2 369.2 369.2 2587.451 2556.151 2631.269 483852.3 483777.1 7215 7214 1
## CP0006 369.2 369.2 369.2 2568.671 2557.716 2578.061 144624.8 144602.9 7033 7032 1返回的 matrix 提供了每个已识别色谱峰的 m/z 和保留时间范围，以及整合的信号强度（“into”）和最大峰强度（“maxo”）。“sample” 列包含确定峰值的对象 / 实验中样本的索引。
可以使用 chromPeakData 函数提取每个峰的注释。该数据框还可以用于为每个检测到的峰添加 / 存储任意注释。
chromPeakData(xdata)## DataFrame with 1707 rows and 2 columns
## ms_level is_filled
## <integer> <logical>
## CP0001 1 FALSE
## CP0002 1 FALSE
## CP0003 1 FALSE
## CP0004 1 FALSE
## CP0005 1 FALSE
## ... ... ...
## CP1703 1 FALSE
## CP1704 1 FALSE
## CP1705 1 FALSE
## CP1706 1 FALSE
## CP1707 1 FALSE峰检测不会总是完美地工作，从而导致峰检测伪影（artifacts），例如重叠峰或人为分割峰。 refineChromPeaks 函数允许通过去除未通过特定标准的已识别峰或合并人工拆分的色谱峰来优化峰检测结果。使用参数对象 CleanPeaksParam 和 FilterIntensityParam ，可以分别删除保留时间范围或强度低于阈值的峰（有关更多详细信息和示例，请参阅各自的帮助页面）。使用 MergeNeighbringPeaksParam 可以合并色谱峰。下面我们对峰检测结果进行进一步处理，如果它们之间的信号低于较小峰最大强度的 75%，合并每个文件在 4 秒窗口中重叠的峰值（merging peaks overlapping in a 4 second window per file）。有关设置和方法的详细说明，请参阅 MergeNeighborgingPeaksParam 帮助页面。
mpp <- MergeNeighboringPeaksParam(expandRt = 4)
xdata_pp <- refineChromPeaks(xdata, mpp)如下展示了合并峰的示例。
mzr_1 <- 305.1 + c(-0.01, 0.01)
chr_1 <- chromatogram(filterFile(xdata, 1), mz = mzr_1)
chr_2 <- chromatogram(filterFile(xdata_pp, 1), mz = mzr_1)
par(mfrow = c(1, 2))
png(file="xcms_7.png")
plot(chr_1)
dev.off()
png(file="xcms_8.png")
plot(chr_2)
dev.off() |  |
对于以上色谱图 centWave 检测到 3 个峰（1 个峰代表整个区域，两个较小的峰，见上图中的左侧面板）。 MergeNeighborgingPeaksParam 的峰细化（refinement ）将它们减少为单个峰（上图中的右面板）。注意，这种细化不会合并相邻峰值，因为它们之间的信号低于一定比例（参见下图）。
mzr_1 <- 496.2 + c(-0.01, 0.01)
chr_1 <- chromatogram(filterFile(xdata, 1), mz = mzr_1)
chr_2 <- chromatogram(filterFile(xdata_pp, 1), mz = mzr_1)
par(mfrow = c(1, 2))
png(file="xcms_9.png")
plot(chr_1)
dev.off()
png(file="xcms_10.png")
plot(chr_2)
dev.off() |  |
还应注意，可以对提取的离子色谱图进行峰细化。这可以例如用于微调参数的设置。为了说明这一点，我们下面对提取的离子色谱图 chr_1 进行了峰细化，减少了 minProp 参数，以强制连接两个峰。
res <- refineChromPeaks(chr_1, MergeNeighboringPeaksParam(minProp = 0.05))
chromPeaks(res)## mz mzmin mzmax rt rtmin rtmax into intb maxo sn
## CPM1 496.2 496.2 496.2 3384.012 3294.809 3412.181 45940118 NA 1128960 177
## sample row column
## CPM1 1 1 1png(file="cxms_11.png")
plot(res)
dev.off()
在继续之前，我们用峰细化的结果替换 xdata 对象。
xdata <- xdata_pp下面我们使用 chromPeaks 矩阵中的数据来计算每个文件的一些摘要。
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** 400 10.95 34.43 42.25 53.21 319.2
**ko16.CDF** 522 10.95 32.86 42.25 53.21 151.8
**ko21.CDF** 215 10.95 42.25 50.08 64.95 164.3
**ko22.CDF** 239 10.95 37.56 46.95 59.47 147.1
**wt15.CDF** 408 10.95 31.3 42.25 53.21 161.2
**wt16.CDF** 366 10.95 35.99 45.38 59.47 162.8
**wt21.CDF** 228 10.95 35.99 48.51 65.73 172.1
**wt22.CDF** 334 10.95 35.99 43.82 57.9 228.5
------------------------------------------------------------------------
Table: Summary statistics on identified chromatographic peaks. Shown are number of identified peaks per sample and widths/duration of chromatographic peaks.我们还可以使用 plotChromPeaks 函数绘制一个文件的 m/z—— 保留时间空间中已识别色谱峰的位置。下面我们绘制了第三个样品的色谱峰。
png(file="xcms_12.png")
plotChromPeaks(xdata, file = 3)
dev.off()
为了获得峰检测的全局概览，我们可以沿着保留时间轴绘制每个文件的已识别峰的频率。这允许识别 MS 运行中的时间段，在该时间段中识别出更多数量的峰，并评估这在文件之间是否一致。
png(file="xcms_13.png")
plotChromPeakImage(xdata)
dev.off()
图注：频率是彩色编码的，较高的频率由黄 - 白表示。每行显示一个文件的峰值频率。
接下来，我们重点介绍前面示例峰的已识别色谱峰。在对应于已知峰或内标的峰列表上评估这样的图有助于确保峰检测设置是适当的并且正确地识别预期峰。我们从峰检测结果对象中提取离子色谱图，其中还包含我们可以使用 chromPeaks 函数提取的该离子的已识别色谱峰。
chr_ex <- chromatogram(xdata, mz = mzr, rt = rtr)
chromPeaks(chr_ex)## mz mzmin mzmax rt rtmin rtmax into intb maxo sn
## CP0045 335 335 335 2781.505 2761.160 2809.674 412134.3 383167.4 16856 23
## CP0501 335 335 335 2786.199 2764.290 2812.803 1496244.2 1461187.3 58736 72
## CP1041 335 335 335 2797.154 2775.245 2815.933 159058.8 149229.6 6295 13
## CP2002 335 335 335 2786.199 2764.290 2812.803 932645.2 915333.8 35856 66
## CP2391 335 335 335 2792.461 2768.987 2823.760 876585.5 848569.1 27200 36
## sample row column
## CP0045 1 1 1
## CP0501 2 1 2
## CP1041 3 1 3
## CP2002 6 1 6
## CP2391 7 1 7我们还可以绘制提取的离子色谱图。已识别的色谱峰将在图中自动突出显示。下面我们用矩形突出显示色谱峰，从峰的最小到最大 rt，从 0 到峰的最大信号。
sample_colors <- group_colors[chr_ex$sample_group]
png(file="xcms_14.png")
plot(chr_ex, col = sample_colors, peakType = "rectangle", peakCol = sample_colors[chromPeaks(chr_ex)[, "sample"]], peakBg = NA)
dev.off()
图注：红色和蓝色分别代表 KO 和野生型样本。矩形表示每个样品的色谱峰。
作为上述矩形可视化的替代方案，可以用一个点表示每个峰的顶点位置（将参数 type="point" 传递给函数），或通过指定 type="polygon" 绘制实际识别的峰。为了完全省略突出显示已识别的峰（例如绘制基本峰色谱图或类似图），可以使用 type = "none" ”。下面我们使用 type="polygon" 填充每个样品中每个已识别色谱峰的峰面积。然而，在这样的图中是否仍然可以识别单个峰取决于从中绘制峰的样本数量。
png(file="xcms_15.png")
plot(chr_ex, col = group_colors[chr_raw$sample_group], lwd = 2,
peakBg = sample_colors[chromPeaks(chr_ex)[, "sample"]])
dev.off()
图注：红色和蓝色分别代表 KO 和野生型样本。已识别色谱峰的信号区域用颜色填充。
注意，我们还可以通过在 chromPeaks 方法中提供各自的 m/z 和保留时间范围以及 mz 和 rt 参数，提取选定区域的已鉴定色谱峰。
pander(chromPeaks(xdata, mz = mzr, rt = rtr), caption = paste("Identified chromatographic peaks in a selected ", "m/z and retention time range."))-----------------------------------------------------------------------------
mz mzmin mzmax rt rtmin rtmax into intb
------------ ----- ------- ------- ------ ------- ------- --------- ---------
**CP0045** 335 335 335 2782 2761 2810 412134 383167
**CP0501** 335 335 335 2786 2764 2813 1496244 1461187
**CP1041** 335 335 335 2797 2775 2816 159059 149230
**CP2002** 335 335 335 2786 2764 2813 932645 915334
**CP2391** 335 335 335 2792 2769 2824 876586 848569
-----------------------------------------------------------------------------
Table: Identified chromatographic peaks in a selected m/z and retention time range. (continued below)
----------------------------------
maxo sn sample
------------ ------- ---- --------
**CP0045** 16856 23 1
**CP0501** 58736 72 2
**CP1041** 6295 13 3
**CP2002** 35856 66 6
**CP2391** 27200 36 7
----------------------------------最后，我们还绘制了每个样本的峰强度分布。这允许调查样本之间的峰信号是否存在系统差异。
## Extract a list of per-sample peak intensities (in log2 scale)
ints <- split(log2(chromPeaks(xdata)[, "into"]), f = chromPeaks(xdata)[, "sample"])
png(file="xcms_16.png")
boxplot(ints, varwidth = TRUE, col = group_colors[xdata$sample_group], ylab = expression(log[2]~intensity), main = "Peak intensities")
grid(nx = NA, ny = NULL)
dev.off()
注意，除了上述色谱峰的识别之外，还可以使用 manualChromPeaks 函数手动定义和添加色谱峰（更多信息请参见 ?manualChronPeaks 帮助页面）。
Alignment
分析物在色谱中洗脱的时间可能因样品（甚至化合物）而异。在上一节中作为示例所示的提取离子色谱图中已经可以观察到这种差异。对齐步骤，也称为保留时间校正，旨在通过沿保留时间轴移动信号来调整这一点，以在实验中对齐不同样本之间的信号。
存在很多的对齐算法（参见 [3]），其中一些算法也在 xcms 中实现。在 xcms 中执行对齐 / 保留时间校正的方法是 adjustRtime ，其根据所提供的参数类别使用不同的对齐算法。
在下面的示例中，我们使用 obiwarp 方法 [4] 来对齐样本。我们使用 binSize＝0.6 ，它在 mz 为 0.6 的 bins 中创建扭曲函数（warping functions）。此外，这里建议修改每个实验的设置，并评估保留时间校正是否对齐了内参或已知化合物。
xdata <- adjustRtime(xdata, param = ObiwarpParam(binSize = 0.6))adjustRtime ，除了计算每个光谱的校正保留时间外，还调整了已鉴定色谱峰的报告保留时间。
为了提取校正后的保留时间，我们可以使用 adjustedRtime 方法，或者简单地使用 rtime 方法，如果存在的话，默认情况下会从 XCMSnExp 对象返回校正后的保留时间。
## Extract adjusted retention times
head(adjustedRtime(xdata))## F1.S0001 F1.S0002 F1.S0003 F1.S0004 F1.S0005 F1.S0006
## 2501.378 2502.958 2504.538 2506.118 2507.699 2509.280## Or simply use the rtime method
head(rtime(xdata))## F1.S0001 F1.S0002 F1.S0003 F1.S0004 F1.S0005 F1.S0006
## 2501.378 2502.958 2504.538 2506.118 2507.699 2509.280原始保留时间可以从包含与 rtime(xdata，adjusted=FALSE) 对齐的数据的 XCMSnExp 中提取。
为了评估对齐的影响，我们利用校正后的数据绘制 BPC。此外，我们使用 plotAdjustedRtime 函数绘制每个样本的校正后保留时间与原始保留时间的差异。对于基峰色谱图，从结果对象中提取已识别的色谱峰是没有意义的。因此，我们在 chromatogram 调用中使用参数 include="none" ，以在返回的对象中不包含色谱峰。请注意，也可以通过在 plot 调用中设置 peakType = "none" 来避免绘制它们。
## Get the base peak chromatograms.
bpis_adj <- chromatogram(xdata, aggregationFun = "max", include = "none")
par(mfrow = c(2, 1), mar = c(4.5, 4.2, 1, 0.5))
png(file="xcms_17.png")
plot(bpis_adj, col = group_colors[bpis_adj$sample_group])
dev.off()
## Plot also the difference of adjusted to raw retention time.
png(file="xcms_18.png")
plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group])
dev.off()
图注：对齐后的基峰色谱图（顶部）和沿保留时间轴的校正保留时间与原始保留时间之间的差异（底部）。
校正后的保留时间与原始保留时间之间的差异太大，可能表明样品或对齐效果不佳（poorly performing samples or alignment）。
注意： XCMSnExp 对象保存原始数据以及校正后的保留时间，在大多数情况下，子设置（subsetting）将删除校正后的保持时间。因此，一些情况下用校正后的保留时间代替原始保留时间可能是有用的。这可以通过 applyAdjustedRtime 完成。
最后，我们评估了对齐对测试峰的影响。
par(mfrow = c(2, 1))
## Plot the raw data
png(file="xcms_19.png")
plot(chr_raw, col = group_colors[chr_raw$sample_group])
dev.off()
## Extract the chromatogram from the adjusted object
chr_adj <- chromatogram(xdata, rt = rtr, mz = mzr)
png(file="xcms_20.png")
plot(chr_adj, col = group_colors[chr_raw$sample_group], peakType = "none")
dev.off()
Subset-based alignment
在一些实验中，仅根据样本的子集进行比对可能会有所帮助，例如，如果 QC 样本定期注射（injected ），或者如果实验包含空白样本。xcms 中的对齐方法允许估计样本子集（不包括空白样本的所有样本或在测量运行期间定期注入的质控样本）的保留时间漂移，并使用这些来调整整个数据集。
参数 subset （ PeakGroupsParam 或 ObiwarpParam 对象）可用于定义样本子集，整个数据集的对齐将基于该样本子集（例如，子集是 QC 样本的索引），参数 subsetAdjust 允许指定调整遗漏样本（left-out samples）的方法。目前有两个选项可用：
subsetAdjust = "previous"：根据先前子集样本（例如 QC 样本）的对齐结果调整非子集样本的保留时间。如果样本的顺序为 A1、B1、B2、A2、B3、B4，A 代表 QC 样本，B 代表研究样本，则使用subset = c(1, 4)和subsetAdjust = "previous"将导致所有 A 样本彼此对齐，而非子集样本 B1 和 B2 将根据子集样本 A1 和 B3 以及 B4 对 A2 样本的对齐结果进行调整。subsetAdjust = "average"：基于先前和后续子集样本的对齐结果的内插（interpolation）来调整非子集样本的保留时间。在上述示例中，B1 将基于子集（QC）样本 A1 和 A2 之间的校正保留时间的平均值进行调整。由于在非子集样本 B3 和 B4 之后没有子集样本，因此将仅基于 A2 的对齐结果来校正这些样本。注意，加权平均值用于计算经校正的保留时间平均值，其使用非子集样本与子集样本的索引之差的倒数作为权重。因此，如果我们有类似 A1、B1、B2、A2 的设置，在校正非子集样本 B1 时，A1 的校正保留时间将比 A2 的更大，从而使其校正保留时间更接近 A1 的保留时间，而不是 A2 的保留时间。请参见下面的示例。
这两种情况都需要对对象内的样本进行有意义的 / 正确的排序（例如按注入索引排序）。
以下示例旨在说明这些对齐选项的效果。我们假设 faahKO 数据集中的样本 1、4 和 7 是质控样本（样本池）。因此，我们希望基于这些样本执行对齐，随后基于来自相邻子集（QC）样本的结果的插值来校正剩余的样本（2、3、5、6 和 8）的保留时间。在初始峰分组后，我们使用峰组法（peak groups method）进行校准，通过子集参数传递我们希望校准所基于的样本的索引，并指定 subsetAdjust = "average" 以根据相邻子集 / QC 样本的校准结果的插值校正研究样本。
注意，对于任何子集对齐，所有参数（如 minFraction ）都是相对于子集的，而不是整个实验！
要重新执行对齐，我们可以首先使用 dropAdjustedRtime 函数删除以前的对齐结果。
xdata <- dropAdjustedRtime(xdata)
## Define the experimental layout
xdata$sample_type <- "study"
xdata$sample_type[c(1, 4, 7)] <- "QC"接下来，我们必须进行初始对应分析，因为峰群对齐方法通过对齐之前识别的 hook 峰（存在于大多数 / 所有样品中的色谱峰；下一节将介绍所用算法的详细信息）来校正保留时间。我们在这里使用默认设置，但强烈建议为每个数据集调整参数。PeakDensityParam 必须定义样本组（即，将单个样本分配给实验中的样本组）。如果实验中没有样本组，则应将每个文件的样本组设置为单个值（例如 rep(1, length(fileNames(xdata)) ）。
## Initial peak grouping. Use sample_type as grouping variable
pdp_subs <- PeakDensityParam(sampleGroups = xdata$sample_type,
minFraction = 0.9)
xdata <- groupChromPeaks(xdata, param = pdp_subs)
## Define subset-alignment options and perform the alignment
pgp_subs <- PeakGroupsParam(minFraction = 0.85,
subset = which(xdata$sample_type == "QC"),
subsetAdjust = "average", span = 0.4)
xdata <- adjustRtime(xdata, param = pgp_subs)下面，我们绘制了对齐结果，用绿色标记了作为子集一部分的样本，其他样本用灰色标注。这很好地显示了 subsetAdjust = "average" 的插值是如何工作的：样本 2 的保留时间是根据子集样本 1 和 4 的保留时间进行校正的，但是，给更接近的子集样本 1 赋予了更大的权重，这导致校正后的样本 2 的保持时间与样本 1 的保持时间更相似。另一方面，样本 3 得到校正，给第二子集样本（4）赋予更多权重。
clrs <- rep("#00000040", 8)
clrs[xdata$sample_type == "QC"] <- c("#00ce0080")
par(mfrow = c(2, 1), mar = c(4, 4.5, 1, 0.5))
png(file="xcms_21.png")
plot(chromatogram(xdata, aggregationFun = "sum"), col = clrs, peakType = "none")
dev.off()
png(file="xcms_22.png")
plotAdjustedRtime(xdata, col = clrs, peakGroupsPch = 1, peakGroupsCol = "#00ce0040")
dev.off()
图注：沿保留时间轴校正保留时间与原始保留时间之间的差异。对校准模型进行估算的样本显示为绿色，研究样本显示为灰色。
选项 subsetAdjust = "previous" 基于单个子集样本（先前的）校正非子集样本的保留时间，这在大多数情况下导致非子集样本校正后的保留时间与用于校正的子集样本高度相似。
Correspondence
代谢组学预处理的最后一步是匹配样本之间（取决于设置，如果样本相邻，也在样本内）检测到的色谱峰的对应关系。在 xcms 中执行对应的方法是 groupChromPeaks 。我们将使用峰密度法 [5] 对色谱峰进行分组。该算法结合了色谱峰，这取决于沿着 mz 维度的小切片中沿着保留时间轴的峰密度。为了说明这一点，我们绘制了每个样品中具有多个色谱峰的 mz 切片的色谱图。我们使用 0.4 的 minFraction 参数值，因此每个样本组中至少有 40% 的样本中存在的色谱峰被分组为特征（feature）。样品组分配是使用 sampleGroups 参数指定的。
## Define the mz slice.
mzr <- c(305.05, 305.15)
## Extract and plot the chromatograms
chr_mzr <- chromatogram(xdata, mz = mzr)
## Define the parameters for the peak density method
pdp <- PeakDensityParam(sampleGroups = xdata$sample_group, minFraction = 0.4, bw = 30)
png(file="xcms_23.png")
plotChromPeakDensity(chr_mzr, col = sample_colors, param = pdp, peakBg = sample_colors[chromPeaks(chr_mzr)[, "sample"]], peakCol = sample_colors[chromPeaks(chr_mzr)[, "sample"]], peakPch = 16)
dev.off()
图注：上图：具有多个色谱峰的 mz 切片的色谱图。下图：确定保留时间的色谱峰 (x 轴) 和不同 bw 参数值的实验样品内的指数 (y 轴)。黑线表示峰密度估计值。峰分组（基于提供的设置）由灰色矩形表示。
上图中的上部面板显示了每个样品的提取离子色谱图，突出显示了检测到的峰。中间和下部曲线显示了不同样品中每个检测到的峰的保留时间。黑色实线表示检测到的峰沿保留时间的密度分布。组合成特征（峰值组）的峰值用灰色矩形表示。这种类型的可视化非常适合在将示例 m/z 切片应用于完整数据集之前测试它们的对应设置。
下面，我们将基于定义设置（defined settings）对整个数据集进行对应分析。
## Perform the correspondence
pdp <- PeakDensityParam(sampleGroups = xdata$sample_group, minFraction = 0.4, bw = 30)
xdata <- groupChromPeaks(xdata, param = pdp)基于 xcms 的预处理结果可以通过 quantify 法从 SummarizedExperiment 包中汇总为 SummarizedExperiment 对象。该对象将包含作为分析矩阵的特征丰度、作为 rowData （即行注释）的特征定义（其 m/z、保留时间和其他元数据）和作为 colData （即列注释）的样本 / 表型信息。所有处理历史记录都将放入对象的元数据中（metadata）。然后，该对象可以用于任何进一步的（与 xcms 无关的）处理和分析。
下面我们使用 quantify 来生成当前分析的结果对象。参数值和任何其他附加参数将传递给 featureValues 方法，该方法在内部用于创建特征丰度矩阵。
res <- quantify(xdata, value = "into")样品注释可以通过 colData 方法访问。
colData(res)## DataFrame with 8 rows and 3 columns
## sample_name sample_group sample_type
## <character> <character> <character>
## ko15.CDF ko15 KO QC
## ko16.CDF ko16 KO study
## ko21.CDF ko21 KO study
## ko22.CDF ko22 KO QC
## wt15.CDF wt15 WT study
## wt16.CDF wt16 WT study
## wt21.CDF wt21 WT QC
## wt22.CDF wt22 WT study通过 rowData 注释特征：
rowData(res)DataFrame with 357 rows and 11 columns
mzmed mzmin mzmax rtmed rtmin rtmax npeaks
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
FT001 200.1 200.1 200.1 2902.25 2882.24 2922.27 2
FT002 205.0 205.0 205.0 2789.40 2782.36 2795.94 8
FT003 206.0 206.0 206.0 2788.64 2780.79 2793.78 7
FT004 207.1 207.1 207.1 2718.25 2713.28 2726.63 7
FT005 219.1 219.1 219.1 2518.60 2517.35 2521.09 3
... ... ... ... ... ... ... ...
FT353 595.25 595.2 595.3 3010.39 2992.58 3014.38 6
FT354 596.20 596.2 596.2 2997.60 2992.58 3002.61 2
FT355 596.30 596.3 596.3 3818.98 3811.68 3835.78 4
FT356 597.40 597.4 597.4 3821.10 3817.96 3825.14 3
FT357 599.30 599.3 599.3 4070.45 4042.10 4123.52 3
KO WT peakidx ms_level
<numeric> <numeric> <list> <integer>
FT001 2 0 463,1180 1
FT002 4 4 47,448,960,... 1
FT003 3 4 32, 435,1164,... 1
FT004 4 3 19,425,943,... 1
FT005 1 2 1140,1379,2382 1
... ... ... ... ...
FT353 2 3 67, 529,1466,... 1
FT354 0 2 1457,2452 1
FT355 2 2 332,1110,2071,... 1
FT356 1 2 838,2355,2680 1
FT357 1 2 375,1757,2137 1特征丰度可以通过 assay 方法获得。还要注意， SummarizedExperiment 支持多个这样的分析矩阵。
head(assay(res)) ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF
FT001 NA 506848.9 NA 169955.6 NA NA NA
FT002 1924712.0 1757151.0 1383416.7 1180288.2 2129885.1 1634342.0 1623589.2
FT003 213659.3 289500.7 NA 178285.7 253825.6 241844.4 240606.0
FT004 349011.5 451863.7 343897.8 208002.8 364609.8 360908.9 NA
FT005 NA NA NA 107348.5 223951.8 NA NA
FT006 286221.4 NA 164009.0 149097.6 255697.7 311296.8 366441.5
wt22.CDF
FT001 NA
FT002 1354004.93
FT003 185399.47
FT004 221937.53
FT005 84772.92
FT006 271128.02此外，还可以使用 featureDefinitions 和 featureValues 方法分别从对应分析中提取结果，前者返回带有特征定义的数据框 (即 mz 和 rt 范围，在列 peakidx 中为每个特征的 chromPeaks 矩阵中的色谱峰指数)，后者为特征丰度。
## Extract the feature definitions
featureDefinitions(xdata)DataFrame with 357 rows and 11 columns
mzmed mzmin mzmax rtmed rtmin rtmax npeaks
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
FT001 200.1 200.1 200.1 2902.25 2882.24 2922.27 2
FT002 205.0 205.0 205.0 2789.40 2782.36 2795.94 8
FT003 206.0 206.0 206.0 2788.64 2780.79 2793.78 7
FT004 207.1 207.1 207.1 2718.25 2713.28 2726.63 7
FT005 219.1 219.1 219.1 2518.60 2517.35 2521.09 3
... ... ... ... ... ... ... ...
FT353 595.25 595.2 595.3 3010.39 2992.58 3014.38 6
FT354 596.20 596.2 596.2 2997.60 2992.58 3002.61 2
FT355 596.30 596.3 596.3 3818.98 3811.68 3835.78 4
FT356 597.40 597.4 597.4 3821.10 3817.96 3825.14 3
FT357 599.30 599.3 599.3 4070.45 4042.10 4123.52 3
KO WT peakidx ms_level
<numeric> <numeric> <list> <integer>
FT001 2 0 463,1180 1
FT002 4 4 47,448,960,... 1
FT003 3 4 32, 435,1164,... 1
FT004 4 3 19,425,943,... 1
FT005 1 2 1140,1379,2382 1
... ... ... ... ...
FT353 2 3 67, 529,1466,... 1
FT354 0 2 1457,2452 1
FT355 2 2 332,1110,2071,... 1
FT356 1 2 838,2355,2680 1
FT357 1 2 375,1757,2137 1featureValues 方法返回一个矩阵，其中行是特征，列是样本。这个矩阵的内容可以使用 value 参数定义。默认值 value = "into" 返回一个矩阵，其中包含样本中某个特征对应的峰的整合信号。 chromPeaks 矩阵的任何列名都可以传递给参数 value 。下面我们提取每个特征 / 样本的整合峰值强度。
## Extract the into column for each feature.
head(featureValues(xdata, value = "into")) ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF
FT001 NA 506848.9 NA 169955.6 NA NA NA
FT002 1924712.0 1757151.0 1383416.7 1180288.2 2129885.1 1634342.0 1623589.2
FT003 213659.3 289500.7 NA 178285.7 253825.6 241844.4 240606.0
FT004 349011.5 451863.7 343897.8 208002.8 364609.8 360908.9 NA
FT005 NA NA NA 107348.5 223951.8 NA NA
FT006 286221.4 NA 164009.0 149097.6 255697.7 311296.8 366441.5
wt22.CDF
FT001 NA
FT002 1354004.93
FT003 185399.47
FT004 221937.53
FT005 84772.92
FT006 271128.02该特征矩阵包含某些样品在特征 m/z-rt 区域未检测到色谱峰（NA 值）。虽然在许多情况下，在相应的区域可能确实没有峰信号，但也可能有信号，但峰检测算法未能检测到色谱峰 (例如，因为信号太低或噪声太大)。xcms 提供 fillChromPeaks 方法来填充原始文件中这些缺失值的强度数据。填充的峰被添加到 chromPeaks 矩阵中，并在 chromPeakData 数据框的 “is_filled” 列中用 TRUE 表示。下面我们执行这样的填充。
xdata <- fillChromPeaks(xdata, param = ChromPeakAreaParam())
head(featureValues(xdata))
ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF
FT001 135315.2 506848.88 111783.48 169955.6 210166.7 141768.0 233086.7
FT002 1924712.0 1757150.96 1383416.72 1180288.2 2129885.1 1634342.0 1623589.2
FT003 213659.3 289500.67 164404.50 178285.7 253825.6 241844.4 240606.0
FT004 349011.5 451863.66 343897.76 208002.8 364609.8 360908.9 225332.1
FT005 161425.0 37819.63 82574.72 107348.5 223951.8 135617.6 191928.5
FT006 286221.4 286834.04 164008.97 149097.6 255697.7 311296.8 366441.5
wt22.CDF
FT001 142142.23
FT002 1354004.93
FT003 185399.47
FT004 221937.53
FT005 84772.92
FT006 271128.02对于样本中没有检测到峰的特征，该方法提取特征的 mz-rt 区域中的所有强度，对信号进行整合（integrates），并将填充的峰添加到 chromPeaks 矩阵中。如果特征的 mz-rt 区域没有测量 / 可用信号，则不添加峰。对于这些，即使在填写了缺失的峰数据后，也会在 featureValues 矩阵中报告 NA。
可以使用不同的选项来定义特征的 mz-rt 区域。使用上面使用的 ChromPeakAreaParam() 参数对象，使用其所有（检测到的）色谱峰的 m/z 和 rt 范围定义特征区域：区域中低的 m/z 值定义为特征所有峰的 “mzmin” 值的下四分位（25% 分位数），高 m/z 值为 “mzmax” 值的上四分位数（75% 分位数），低 rt 值作为 “rtmin” 值的下四分位（25% 分位数），高 rt 值为 “rtmax” 值的上四分位数（75% 分位数）。这确保了信号是从特定特征区域集成的。
或者，可以在 fillChromPeaks 调用中使用 FillChromPpeaksParam 参数对象，这类似于初始（old）xcms 实现的方法。
下面我们比较填写缺失值前后缺失值的数量。我们可以使用 featureValues 方法的 filled 参数来定义是否也应该返回填充的峰值（filled-in peak values）。
## Missing values before filling in peaks
apply(featureValues(xdata, filled = FALSE), MARGIN = 2, FUN = function(z) sum(is.na(z)))## ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF wt22.CDF
## 91 93 168 154 95 119 167 123## Missing values after filling in peaks
apply(featureValues(xdata), MARGIN = 2, FUN = function(z) sum(is.na(z)))## ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF wt22.CDF
## 5 4 5 5 4 7 6 2接下来，我们使用 featureSummary 函数获取每个特征的一般摘要，其中包括找到峰的样本数量或为特征分配了多个峰的样本数。同时指定样本组可以分解每个样本组的汇总统计信息。
head(featureSummary(xdata, group = xdata$sample_group)) count perc multi_count multi_perc rsd KO_count KO_perc
FT001 2 25.0 0 0 0.7039537 2 50
FT002 8 100.0 0 0 0.1936518 4 100
FT003 7 87.5 0 0 0.1717662 3 75
FT004 7 87.5 0 0 0.2609145 4 100
FT005 3 37.5 0 0 0.5385767 1 25
FT006 7 87.5 0 0 0.3016973 3 75
KO_multi_count KO_multi_perc KO_rsd WT_count WT_perc WT_multi_count
FT001 0 0 0.7039537 0 0 0
FT002 0 0 0.2178920 4 100 0
FT003 0 0 0.2501505 4 100 0
FT004 0 0 0.2957873 3 75 0
FT005 0 0 NA 2 50 0
FT006 0 0 0.3765933 4 100 0
WT_multi_perc WT_rsd
FT001 0 NA
FT002 0 0.1918936
FT003 0 0.1327983
FT004 0 0.2575039
FT005 0 0.6375539
FT006 0 0.1641781我们可以将特征值矩阵与缺失峰值的填充数据一起添加到我们的 SummarizedExperiment 对象 res 中，作为附加矩阵（assay）：
assays(res)$raw_filled <- featureValues(xdata, filled = TRUE)我们现在有两个可用的矩阵（assays），一个具有检测值的矩阵，一个包含检测值和填充值的矩阵，每个都可以通过其名称访问。
assayNames(res)## [1] "raw" "raw_filled"head(assay(res, "raw")) ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF
FT001 NA 506848.9 NA 169955.6 NA NA NA
FT002 1924712.0 1757151.0 1383416.7 1180288.2 2129885.1 1634342.0 1623589.2
FT003 213659.3 289500.7 NA 178285.7 253825.6 241844.4 240606.0
FT004 349011.5 451863.7 343897.8 208002.8 364609.8 360908.9 NA
FT005 NA NA NA 107348.5 223951.8 NA NA
FT006 286221.4 NA 164009.0 149097.6 255697.7 311296.8 366441.5
wt22.CDF
FT001 NA
FT002 1354004.93
FT003 185399.47
FT004 221937.53
FT005 84772.92
FT006 271128.02head(assay(res, "raw_filled")) ko15.CDF ko16.CDF ko21.CDF ko22.CDF wt15.CDF wt16.CDF wt21.CDF
FT001 135315.2 506848.88 111783.48 169955.6 210166.7 141768.0 233086.7
FT002 1924712.0 1757150.96 1383416.72 1180288.2 2129885.1 1634342.0 1623589.2
FT003 213659.3 289500.67 164404.50 178285.7 253825.6 241844.4 240606.0
FT004 349011.5 451863.66 343897.76 208002.8 364609.8 360908.9 225332.1
FT005 161425.0 37819.63 82574.72 107348.5 223951.8 135617.6 191928.5
FT006 286221.4 286834.04 164008.97 149097.6 255697.7 311296.8 366441.5
wt22.CDF
FT001 142142.23
FT002 1354004.93
FT003 185399.47
FT004 221937.53
FT005 84772.92
FT006 271128.02应始终通过检查提取的离子色谱图（例如已知化合物、内标或一般识别特征）来评估峰检测、比对（alignment）和对应（correspondence）的性能。 featureChromatogram 函数允许提取 featureDefinitions 中每个特征的色谱图。返回的 MCchromatograms 对象包含每个特征（每行包含一个特征的数据）和样本（每列代表一个样本的数据）的离子色谱图。下面我们提取前 4 个特征的色谱图。
feature_chroms <- featureChromatograms(xdata, features = 1:4)
feature_chroms## XChromatograms with 4 rows and 8 columns
## 1 2 3 4
## <XChromatogram> <XChromatogram> <XChromatogram> <XChromatogram>
## [1,] peaks: 0 peaks: 1 peaks: 0 peaks: 1
## [2,] peaks: 1 peaks: 1 peaks: 1 peaks: 1
## [3,] peaks: 1 peaks: 1 peaks: 0 peaks: 1
## [4,] peaks: 1 peaks: 1 peaks: 1 peaks: 1
## 5 6 7 8
## <XChromatogram> <XChromatogram> <XChromatogram> <XChromatogram>
## [1,] peaks: 0 peaks: 0 peaks: 0 peaks: 0
## [2,] peaks: 1 peaks: 1 peaks: 1 peaks: 1
## [3,] peaks: 1 peaks: 1 peaks: 1 peaks: 1
## [4,] peaks: 1 peaks: 1 peaks: 0 peaks: 1
## phenoData with 3 variables
## featureData with 5 variables
## - - - xcms preprocessing - - -
## Chromatographic peak detection:
## method: centWave
## Correspondence:
## method: chromatographic peak density
## 4 feature(s) identified.绘制提取的离子色谱图。我们再次使用每个鉴定峰的组颜色来填充该区域。
png(file="xcms_24.png")
plot(feature_chroms, col = sample_colors, peakBg = sample_colors[chromPeaks(feature_chroms)[, "sample"]])
dev.off()
为了访问第二个特征的 EIC，我们可以简单地将 feature_chros 对象取子集。
eic_2 <- feature_chroms[2, ]
chromPeaks(eic_2) mz mzmin mzmax rt rtmin rtmax into intb maxo sn
CP0055 205 205 205 2791.599 2771.048 2815.323 1924712 1850331 84280 64
CP0502 205 205 205 2795.355 2773.257 2820.612 1757151 1711473 68384 69
CP1049 205 205 205 2795.936 2773.877 2821.143 1383417 1334570 47384 54
CP1277 205 205 205 2788.642 2768.195 2812.229 1180288 1126958 48336 32
CP1567 205 205 205 2782.355 2761.949 2805.895 2129885 2054677 93312 44
CP2005 205 205 205 2787.051 2766.690 2812.111 1634342 1566379 67984 53
CP2392 205 205 205 2790.157 2763.587 2821.427 1623589 1531573 49208 28
CP2649 205 205 205 2787.032 2766.714 2812.043 1354005 1299188 55712 35
sample row column
CP0055 1 1 1
CP0502 2 1 2
CP1049 3 1 3
CP1277 4 1 4
CP1567 5 1 5
CP2005 6 1 6
CP2392 7 1 7
CP2649 8 1 8最后，我们进行了主成分分析，以评估本实验中样本的分组。注意，我们没有执行任何数据标准化，因此分组可能（也将）受到技术偏见的影响。
## Extract the features and log2 transform them
ft_ints <- log2(assay(res, "raw_filled"))
## 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)
png(file="xcms_25.png")
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)
dev.off()
我们可以在 PC2 上看到 KO 和 WT 样品之间的预期分离。在 PC1 上样本根据其 ID 进行分离，ID<=18 的样本与 ID>18 的样本被分开。这种分离可能是由于技术偏差（例如，在不同的天 / 周进行的测量）或由于所分析小鼠的生物学特性（性别、年龄、产仔等）造成的。
Further data processing and analysis
标准化特征的信号强度是必需的，但目前 xcms 还不支持（在不久的将来可能会添加一些方法）。建议将 quantify 方法返回的 SummarizedExperiment 用于任何进一步的数据处理，因为这种类型的对象在同一对象中存储特征定义（feature definitions）、样本注释以及特征丰度。为了识别例如具有显著不同强度 / 丰度的特征，建议使用其他 R 包中提供的功能，例如 Bioconductor 的 limma 包。为了支持依赖于旧 xcmsSet 结果对象的其他包，可以使用 xset <- as(x, "xcmsSet") 将新的 XCMSnExp 对象强制转换为 xcmsSet 对象，其中 x 是 XCMSnExp 对象。
参考
关注公众号 “生信之巅”。
 |  |
敬告：使用文中脚本请引用本文网址，转载请注明出处请尊重本人的劳动成果，谢谢！
