library(xcms)
library(ggplot2)
library(ggrepel)
dda_data <- readMSData("PestMix1_DDA.mzML", mode = "onDisk")

检出peak
官方推荐参数

cwp <- CentWaveParam(snthresh = 5, noise = 100, ppm = 10,
                     peakwidth = c(3, 30))
dda_data <- findChromPeaks(dda_data, param = cwp)

参数设置.png

查看检出的features

plotChromPeaks(dda_data)

plotChromPeaks

上图中一条线（实际为矩形）代表一个离子，线的长度代表峰宽，宽度代表质荷比m/z的范围。
图形不够美观，也不能带标记label。
使用ggplot2绘图。

data <- chromPeaks(dda_data)|> data.frame()
ggplot(data)+ 
  geom_rect(aes(xmin = rtmin,xmax = rtmax ,ymin = mzmin ,ymax = mzmax),color = "black",alpha = 0)+
  geom_text_repel(aes(x = rt,y = mz,label = row.names(data)),size = 3)+
  xlab("Retention time")+
  theme_bw()

geom_rect

可以观察到许多features具有非常相似的质荷比mz和非常接近的保留时间rt。如CP083，CP084。

i = 83
mz = c(data$mzmin[i],data$mzmax[i])
mz
rt = c(data$rtmin[i],data$rtmax[i])
rt
mzr <- mz + c(-0.01, 0.01)
rtr <- rt + c(-5, 5)
chromPeaks(dda_data, mz = mz, rt = rt)

#           mz    mzmin    mzmax      rt   rtmin   rtmax     into     intb     maxo  sn sample
#CP83 217.144 217.1437 217.1443 508.258 504.478 513.617 473.4116 468.0113 184.1682 116      1
#CP84 217.144 217.1437 217.1443 508.678 508.678 513.617 207.2879 204.4068 163.8278 116      1

看看他们是否是同一种feature。
可视化

ggplot(data)+ 
  geom_rect(aes(xmin = rtmin,xmax = rtmax ,ymin = mzmin ,ymax = mzmax),color = "black",alpha = 0)+
  geom_text_repel(aes(x = rt,y = mz,label = row.names(data)),size = 3)+
  xlab("Retention time")+
  xlim(rtr)+
  ylim(mzr)+
  theme_bw()

geom_rect

提取EIC图

dda_data |> 
  filterMz(mzr)|>
  filterRt(rtr)|>
  plot(type = "XIC")
abline(v = data$rt[i],col = "red")
abline(h = data$mz[i],col = "blue")

XIC

证明他们确实是同一个features
设置合并参数将其合并。

mpp <- MergeNeighboringPeaksParam(expandRt = 2, expandMz = 0.01)
xdata_pp <- refineChromPeaks(dda_data, mpp)

# Evaluating 5 peaks in file PestMix1_DDA.mzML for merging ... OK

#Merging reduced 99 chromPeaks to 99.
#Warning message:
#In serialize(data, node$con) :
#  'package:stats' may not be available when loading

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

比较合并前和合并后的峰形

chr_raw <- chromatogram(dda_data, mz = mzr, rt = rtr)
plot(chr_raw)
abline(v = data$rt[i],col = "red")
chr_raw <- chromatogram(xdata_pp, mz = mzr, rt = rtr)
plot(chr_raw)
abline(v = data$rt[i],col = "red")

chromatogram

合并前后对比

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

dda_spectra <- chromPeakSpectra(dda_data)
mcols(dda_spectra) %>% nrow()
#[1] 158
mcols(dda_spectra)$peak_id |>table()

#CP01 CP04 CP05 CP06 CP08 CP11 CP12 CP13 CP14 CP18 CP22 CP25 CP26 
#   3    2    2    2    2    2    2    4    4    1    5    5    6 
#CP33 CP34 CP35 CP36 CP41 CP42 CP44 CP46 CP47 CP48 CP50 CP51 CP52 
#   2    5    5    1    3    5    2    4    3    3    1    3    3 
#CP53 CP57 CP59 CP60 CP61 CP63 CP64 CP65 CP66 CP67 CP69 CP71 CP72 
#   5    5    2    2    2    4    5    1    4    3    3    4    3 
#CP73 CP81 CP82 CP85 CP88 CP89 CP90 CP91 CP93 CP94 CP95 CP98 CP99 
 #  1    4    3    2    2    3    3    3    6    4    1    2    1

共检测到158个二级质谱，有些色谱峰没有二级质谱，如CP02，有些色谱峰又有好几个二级质谱，如CP42。

i = 93
id = rownames(data)
id[i]
ex_spectra_S1 <- dda_spectra[mcols(dda_spectra)$peak_id == id[i]]
ex_spectra_S1
plot(ex_spectra_S1[[1]])

MS2

e = c()
l = length(ex_spectra_S1)
for (i in 1:l){
  tar <- ex_spectra_S1[[i]]
  for (t in 1:l) {
    ex_spectra <- ex_spectra_S1[[t]]
    d <-compareSpectra(tar,ex_spectra, binSize = 0.1,fun = "dotproduct")
    e <-append(e, d)
  }
}
e
matrix <- matrix(e,l,l)
names(ex_spectra_S1)
row.names(matrix)<- colnames(matrix) <- names(ex_spectra_S1)
plot(hclust(as.dist(matrix)))
#data1
pheatmap::pheatmap(matrix,
                   display_numbers = T,
                   cluster_rows = F,
                   cluster_cols = F,
                   number_format = "%.2f")

pheatmap

整合同一个feature的多张图谱

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

参考资料：

XCMS官方说明文档
LC-MS/MS data analysis with xcms (bioconductor.org)
使用xcms进行LC-MS数据预处理和分析 • xcms (sneumann.github.io)
Chromatographic peak detection using the centWave method — findChromPeaks-centWave • xcms (sneumann.github.io)
centWave Parameters • msbrowser (tkimhofer.github.io)