导入数据

library(xcms)
library(faahKO)
library(RColorBrewer)
library(pander)
library(magrittr)
library(pheatmap)
library(SummarizedExperiment)

1.单个样本

dda_file <- system.file("TripleTOF-SWATH/PestMix1_DDA.mzML",package = "msdata")
dda_data <- readMSData(dda_file, mode = "onDisk")
table(msLevel(dda_data))

3 色谱峰检测

3.1 设置参数

使用centWave算法对centroid模式的高分辨LC-MS进行色谱峰检测。centWave算法最适用于高分辨率 centroid模式 的LC/{TOF、OrbiTrap、FTICR}-MS数据。在第一阶段，该方法确定了感兴趣区域(ROI)，这些区域代表了LC/MS连续扫描时小于ppm m/z偏差的质量轨迹。 详细地说，从单个m/z开始，如果在下一次扫描(频谱)中发现的m/z，其与平均m/z的差异小于用户定义的m/z的ppm，则合并为一个ROI。 考虑到新加入的m/z值，ROI的平均m/z值也随之更新。

cwp <- CentWaveParam(snthresh = 5, noise = 100, ppm = 14,
                     peakwidth = c(1, 30))

参数说明

snthresh -- 定义峰检测时信噪比阈值
prefilter，c(k,I)-- 在第一步(ROI检测)指定预过滤步骤。 只有当它们包含至少k个强度值>= I的峰时，质量才会被保留。
noise --设置在第一个分析步骤中考虑的质心所需的最小强度(强度<噪声的将从ROI检测中忽略)
ppm -- 定义峰检测时连续扫描的m/z最大允许偏差，以ppm（parts per million）为单位
peakwidth -- 定义峰检测时峰宽范围，类似色谱峰的峰宽，该参数指定每个峰峰宽的最大值和最小值，以秒为单位。

centWave两个最关键参数是peakwidth （色谱峰宽的预期范围）和ppm（对应于一个色谱峰的质心 m/z 值的最大预期偏差；通常远大于制造商指定的 ppm）。

3.2 检测一级质谱

dda_data <- findChromPeaks(dda_data, param = cwp)

dda_data.png

此时dda_data多了msFeatureData属性，检测到的色谱峰信息储存在
dda_data@msFeatureData[["chromPeaks"]]

dda_data@msFeatureData[["chromPeaks"]]%>% head()

#         mz    mzmin    mzmax      rt   rtmin   rtmax       into       intb      maxo  sn sample
#CP001 142.0084 142.0082 142.0085 100.967  98.508 102.562   288.7129   285.3309  125.2583  16      1
#CP002 142.9926 142.9922 142.9931 130.205 126.506 132.984  1148.2069  1141.9128  346.7006  88      1
#CP003 221.0918 221.0906 221.0925 239.257 236.657 240.897   459.3357   455.0957  189.3477  66      1
#CP004 221.0920 221.0916 221.0923 241.847 236.657 254.562   367.0682   357.8542  212.5239  15      1
#CP005 220.0985 220.0978 220.0988 241.018 237.187 247.874  1835.6351  1827.5389  585.3036  84      1
#CP006 219.0957 219.0950 219.0962 241.018 236.253 246.863 14940.2013 14880.2656 4877.1162 132      1

dda_data@msFeatureData[["chromPeaks"]]%>% nrow()
[1] 111

共检测到111个一级质谱,将检测到的离子可视化

plotChromPeaks(dda_data, xlim = c(80, 900), ylim = c(130, 650))

plotChromPeaks.png

上图中一条线代表一个离子，线的长度代表峰宽，宽度代表质荷比m/z的范围。
也可以画成散点图

data=dda_data@msFeatureData[["chromPeaks"]]%>% data.frame()
plot(data$mz,data$rt)

plot.png

或者ggplot2出图

p1=ggplot(data,aes(mz,rt))+geom_point(alpha=0.2)+theme_bw()
p1

ggplot.png

途中可以看到许多重叠的点，说明这些离子具有相似甚至相同的保留时间和质荷比。

合并

合并伪峰和重叠峰

设置合并参数

mpp <- MergeNeighboringPeaksParam(expandRt = 2,expandMz = 0.45)
xdata_pp <- refineChromPeaks(dda_data, mpp)
#Merging reduced 111 chromPeaks to 55.

expandRt将保留时间窗口向两边扩展多少秒，
expandMz，将每个色谱峰的m/z范围扩大(两边)。
minProp, 数值在0-1之间，表示峰连接所需的强度比例。 默认(minProp = 0.75)表示只有当信号的中间部分大于两个峰值的“maxo”(峰值顶点的最大强度)中最小值的75%时，峰值才会连接。
111个峰合并为55个。
合并前后对比

data=xdata_pp@msFeatureData[["chromPeaks"]]%>% data.frame()
p2=ggplot(data,aes(mz,rt))+geom_point(alpha=0.2)+theme_bw()
p1|p2

对比.png

提取上述色谱峰的二级质谱

dda_spectra <- chromPeakSpectra(dda_data)
mcols(dda_spectra) %>% nrow()
#[1] 122

共检测到122个二级质谱，有些色谱峰没有二级质谱，如CP001，有些色谱峰又有好几个二级质谱，如CP002。

ex_spectra <- dda_spectra[mcols(dda_spectra)$peak_id == "CP002"]
plot(ex_spectra[[1]])+theme_bw()

MS2.png

整合多个图谱

ex_spectrum <- combineSpectra(ex_spectra, method = consensusSpectrum, mzd = 0,
                              ppm = 20, minProp = 0.8, weighted = FALSE,
                              intensityFun = median, mzFun = median)
ex_spectrum
plot(ex_spectrum[[1]])

将整合后的碎片信息光谱与整合前做对比

par(mfrow = c(1, 3))
plot(ex_spectrum[[1]],ex_spectra[[1]])
plot(ex_spectrum[[1]],ex_spectra[[2]])
plot(ex_spectrum[[1]],ex_spectra[[3]])

质谱比对.png

比对评分

compareSpectra(ex_spectrum[[1]],ex_spectra[[1]], binSize = 0.02,
               fun = "dotproduct")
compareSpectra(ex_spectrum[[1]],ex_spectra[[2]], binSize = 0.02,
               fun = "dotproduct")
compareSpectra(ex_spectrum[[1]],ex_spectra[[3]], binSize = 0.02,
               fun = "dotproduct")

参考资料：
https://bioconductor.org/packages/release/bioc/html/xcms.html
LCMS data preprocessing and analysis with xcms (bioconductor.org)