LCMS data preprocessing and analysis with xcms
Package: xcms
!官网
1. 介绍
该文件描述了使用xcms >= 3版本的LCMS实验的数据导入、探索、预处理和分析。这些例子和基本工作流程改编自Colin A. Smith的LC/MS预处理和分析小册子。
新的用户界面和方法使用XCMSnExp对象（而不是旧的xcmsSet对象）作为预处理结果的容器。为了支持依赖xcmsSet对象的包和管道，可以使用as方法将XCMSnExp转换成xcmsSet对象（即xset <- as(x, "xcmsSet")），x是一个XCMSnExp对象。
2. 数据导入（读数据）
xcms支持分析来自（AIA/ANDI）NetCDF、mzXML和mzML格式文件的LC/MS数据。对于实际的数据导入，使用了Bioconductor的mzR。为了演示，我们将分析[1]中的一个数据子集，其中研究了敲除小鼠脂肪酸酰胺水解酶（FAAH）基因的代谢后果。原始数据文件（NetCDF格式）由faahKO数据包提供。该数据集由6只敲除型和6只野生型小鼠的脊髓样本组成。每个文件都包含在正离子模式下获得的200-600 m/z和2500-4500秒的中心模式的数据。为了加快这个小册子的处理速度，我们将把分析限制在只有8个文件和2500-3500秒的保留时间范围内。
下面我们加载所有需要的软件包，在faahKO软件包中找到原始CDF文件，并建立一个描述实验设置的phenodata数据框。请注意，对于真正的实验，建议定义一个文件（表），包含原始数据文件的文件名以及每个文件的样本描述作为附加列。这样的文件可以用例如read.table作为变量pd导入（而不是像下面的例子那样在R中定义），文件名可以传递给下面的readMSData函数，例如files = paste0(MZML_PATH, "/", pd$mzML_file) 其中MZML_PATH是文件所在目录的路径，"mzML_file "是phenodata文件中包含文件名的列名。
初步接触和摸索代谢或者蛋白组的数据，还是萌萌的。mzXML这个格式我是搞定了的，下机的‘.raw’，被我转成mzXML之后读入进行分析了，这里还是尊重官网教程走一遍，后面有机会歇一歇怎么转数据吧。
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方法将原始数据加载为OnDiskMSnExp对象。MSnbase提供了处理质谱数据的基础结构和基础设施。此外，MSnbase还可以用来对剖面模式的质谱数据进行中心化处理（见MSnbase包中的相应小册子）。
raw_data <- readMSData(files = cdfs, pdata = new("NAnnotatedDataFrame", pd),
mode = "onDisk")
# Polarity can not be extracted from netCDF files, please set manually the polarity with the 'polarity' method.
接下来我们将数据集限制在2500到3500秒的保留时间范围内。这只是为了减少这个小插曲的处理时间。也不知道这个阈值怎么确定的
raw_data <- filterRt(raw_data, c(2500, 3500))
由此产生的OnDiskMSnExp对象包含关于光谱数量、保留时间、测量的总离子电流等一般信息，但不包含完整的原始数据（即每个测量光谱的m/z和强度值）。因此，它的内存占用相当小，使其成为代表大型代谢组学实验的理想对象，同时允许执行简单的质量控制、数据检查和探索以及数据子集操作。m/z和强度值是按要求从原始数据文件中导入的，因此在初始数据导入后，原始数据文件的位置不应改变。
3.初始数据检查（相当于质控了）
OnDiskMSnExp按光谱组织MS数据，并提供强度、mz和rtime方法来访问文件中的原始数据（测量的强度值、相应的m/z和保留时间值）。此外，光谱方法可用于返回封装在光谱对象中的所有数据。下面我们从对象中提取保留时间值。
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
所有的数据都以一维向量的形式返回（一个用于rtime的数字向量和一个用于mz和强度的数字向量列表，每个向量包含一个光谱的值），即使实验由多个文件/样品组成。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
作为对数据的第一次评估，我们在下面为实验中的每个文件绘制了基础峰色谱图（BPC）。我们使用色谱方法，并将aggregationFun设置为 "max"，以返回每个光谱的最大强度，从而从原始数据中创建BPC。为了创建一个总的离子色谱图，我们可以将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.
plot(bpis, col = group_colors[raw_data$sample_group])
# ## Plot all chromatograms. [保存图片]
# pdf("all_sample_bpis.pdf",width = 7,height = 5)
# plot(bpis, col = group_colors[raw_data$sample_group])
# dev.off()
我画的图，不是官网的图，一样呢
色谱方法返回了一个MCromatograms对象，该对象将各个色谱对象（实际上包含色谱数据）组织成一个二维数组：列代表样品，行（可选）代表m/z和/或保留时间范围。下面我们提取第一个样品的色谱图并访问其保留时间和强度值。
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.203
head(intensity(bpi_1))
## F1.S0001 F1.S0002 F1.S0003 F1.S0004 F1.S0005 F1.S0006
## 43888 43960 43392 42632 42200 42288
色谱方法也支持从质谱数据的m/z-rt切片中提取色谱数据。在下一节中，我们将使用这种方法为一个选定的峰创建一个提取的离子色谱图（EIC）。
请注意，chromatogram会从每个文件中读取原始数据来计算色谱图。另一方面，bpi和tic方法不从原始文件中读取任何数据，而是使用输入文件的头定义中提供的相应信息（可能与实际数据不同）。
下面我们创建了代表每个文件的总离子电流分布的图表。这样的图对于发现有问题或失败的MS运行非常有用。
tc <- split(tic(raw_data), f = fromFile(raw_data))
boxplot(tc, col = group_colors[raw_data$sample_group],
ylab = "intensity", main = "Total ion current")
依旧是我产生的图呢
另外，我们可以根据样品基峰色谱的相似性对其进行分组。这也有助于发现实验中可能有问题的样品，或者一般来说，对实验中的样品分组有一个初步的了解。由于样品之间的保留时间并不完全相同，我们使用bin函数沿保留时间轴将强度分组在固定的时间范围（bin）内。在本例中，我们使用的bin大小为1秒，默认为0.5秒。聚类是在分档后的基峰色谱的成对相关性上使用完全联系的层次结构聚类法进行的。
## 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))
聚类啊，样本的相关性，说实话pheatmap产生的图真的不怎么样
这些样本以成对的方式聚类，样本指数的KO和WT样本具有最相似的BPC。
4.色谱峰检测
接下来我们使用centWave算法[2]进行色谱峰值检测。然而，在运行峰值检测之前，强烈建议目视检查，例如内标或已知化合物的提取离子色谱图，以评估和调整峰值检测设置，因为默认设置将不适合大多数LCMS实验。centWave最关键的两个参数是峰宽（色谱峰宽的预期范围）和ppm（对应于一个色谱峰的中心点的m/z值的最大预期偏差；这通常比制造商指定的ppm大很多）参数。为了评估典型的色谱峰宽，我们绘制一个峰的EIC。
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])
一个离子的peak
请注意，如果在某次扫描中（即对于特定的保留时间）没有在相应的mz范围内测量到信号，则通过色谱法提取的色谱对象包含一个NA值。这反映在上面的图中，这些线没有被画成连续线。
上面的峰的宽度约为50秒。峰宽参数的设置应适应数据集中预期的峰宽。对于本例数据集，我们将其设置为20,80。
对于ppm参数，我们提取与上述峰值相对应的完整MS数据（强度、保留时间和m/z值）。为此，我们首先按保留时间过滤原始对象，然后按m/z过滤，最后用type = "XIC "绘制对象，产生下面的图。我们使用管道（%>%）命令更好地说明了相应的工作流程。还要注意的是，在这种类型的绘图中，如果存在已识别的色谱峰，将默认显示。
raw_data %>%
filterRt(rt = rtr) %>%
filterMz(mz = mzr) %>%
plot(type = "XIC")
image.png
对于每个图：上图：色谱图绘制强度值与保留时间的关系，下图m/z与保留时间的关系。各个数据点根据强度的不同而着色。
在目前的数据中，m/z值实际上没有变化。通常情况下，人们会看到m/z值（下图）在化合物的真实m/z值周围散开。centWave算法的第一步是根据连续扫描的m/z值的差异来定义所谓的感兴趣区域（ROI）。详细地说，如果连续扫描的m/z值与ROI的平均m/z值之间的差异小于用户定义的ppm参数，则来自连续扫描的m/z值被纳入ROI。因此，对ppm的合理选择可以是作为色谱峰一部分的相邻扫描/谱段的数据点的最大m/z差。建议检查许多化合物（内部标准或已知存在于样品中的化合物）的m/z值范围，并根据这些定义centWave的ppm参数。
请注意，我们也可以对提取的离子色谱图进行峰值检测。这可以帮助评估不同的峰值检测设置。请注意，提取的离子色谱图上的峰值检测将不考虑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_colors <- group_colors[xchr$sample_group]
plot(xchr, col = sample_colors,
peakBg = sample_colors[chromPeaks(xchr)[, "column"]])
image.png
红色和蓝色分别代表KO和野生型样品。确定的色谱峰的峰面积在样品组颜色中突出显示。
最后我们对完整的数据集进行色谱峰检测。请注意，我们将参数prefilter设置为c(6, 5000)，噪声设置为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级别>1的数据上）。
色谱峰检测的结果可以通过chromPeaks方法获取。
head(chromPeaks(xdata))
## mz mzmin mzmax rt rtmin rtmax into intb maxo
## CP0001 453.2 453.2 453.2 2509.203 2501.378 2527.982 1007409.0 1007380.8 38152
## CP0002 236.1 236.1 236.1 2518.593 2501.378 2537.372 253501.0 226896.3 12957
## CP0003 594.0 594.0 594.0 2601.535 2581.191 2637.529 161042.2 149297.3 7850
## CP0004 577.0 577.0 577.0 2604.665 2581.191 2626.574 136105.2 129195.5 6215
## CP0005 369.2 369.2 369.2 2587.451 2556.151 2631.269 483852.3 483777.1 7215
## CP0006 369.2 369.2 369.2 2568.671 2557.716 2578.061 144624.8 144602.9 7033
## sn sample
## CP0001 38151 1
## CP0002 11 1
## CP0003 13 1
## CP0004 13 1
## CP0005 7214 1
## CP0006 7032 1
返回的矩阵提供了每个识别的色谱峰的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
峰值检测并不总是完美的，会导致峰值检测的假象，如重叠的峰值或人为分割的峰值。refineChromPeaks功能允许通过移除不符合特定标准的识别峰或合并人为分割的色谱峰来完善峰值检测结果。通过参数对象CleanPeaksParam和FilterIntensityParam，可以分别移除保留时间范围内的峰或强度低于阈值的峰（更多细节和例子见各自的帮助页面）。通过MergeNeighboringPeaksParam，可以合并色谱峰。下面我们对色谱峰检测结果进行后处理，如果重叠的色谱峰之间的信号低于较小的色谱峰最大强度的75%，则在每个文件的4秒窗口中合并这些色谱峰。请参阅MergeNeighboringPeaksParam帮助页面，了解设置和方法的详细说明。
mpp <- MergeNeighboringPeaksParam(expandRt = 4)
xdata_pp <- refineChromPeaks(xdata, mpp)
# An example for a merged peak is given below.
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))
plot(chr_1)
plot(chr_2)
左图：处理前的数据，右图：细化后的数据。分裂的峰被合并成一个
对于上述色谱图中的第一条痕量，centWave检测到了3个峰（1个全区峰和2个较小的峰，见上图中的左面板）。用MergeNeighboringPeaksParam进行的峰细化将它们减少到一个峰（上图中的右面板）。请注意，如果相邻的峰之间的信号低于一定的比例，这种细化就不会合并它们（见下图）。
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))
plot(chr_1)
plot(chr_2)
左图：处理前的数据，右图：细化后的数据。峰值没有被合并
还要注意的是，可以对提取的离子色谱图进行峰的细化。例如，这可以用来微调参数的设置。为了说明这一点，我们在下面对提取的离子色谱图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 1255
## sample row column
## CPM1 1 1 1
plot(res)
image.png
[to be continued]
网友评论