library(xcms)
library(ggplot2)
library(ggrepel)
dda_data <- readMSData("PestMix1_DDA.mzML", mode = "onDisk")
檢出peak
官方推薦參數(shù)
cwp <- CentWaveParam(snthresh = 5, noise = 100, ppm = 10,
peakwidth = c(3, 30))
dda_data <- findChromPeaks(dda_data, param = cwp)
查看檢出的features
plotChromPeaks(dda_data)
上圖中一條線（實(shí)際為矩形）代表一個(gè)離子救崔，線的長度代表峰寬，寬度代表質(zhì)荷比m/z的范圍捏顺。
圖形不夠美觀六孵，也不能帶標(biāo)記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()
可以觀察到許多features具有非常相似的質(zhì)荷比mz和非常接近的保留時(shí)間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()
提取EIC圖
dda_data |>
filterMz(mzr)|>
filterRt(rtr)|>
plot(type = "XIC")
abline(v = data$rt[i],col = "red")
abline(h = data$mz[i],col = "blue")
證明他們確實(shí)是同一個(gè)features
設(shè)置合并參數(shù)將其合并。
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函數(shù)的參數(shù)有4個(gè)挪凑，分別為：expandRt, expandMz , ppm和 minProp孕索。
expandRt將保留時(shí)間窗口向兩邊擴(kuò)展多少秒，
expandMz躏碳，將每個(gè)色譜峰的m/z范圍擴(kuò)大(兩邊)搞旭。
minProp, 數(shù)值在0-1之間，表示合并兩個(gè)峰所需要達(dá)到的強(qiáng)度比例。 默認(rèn)(minProp = 0.75)表示只有當(dāng)信號(hào)的中間部分大于兩個(gè)峰值的“maxo”(峰值頂點(diǎn)的最大強(qiáng)度)中最小值的75%時(shí)肄渗，峰才會(huì)合并镇眷。
比較合并前和合并后的峰形
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")
合并前后對(duì)比
提取上述色譜峰的二級(jí)質(zhì)譜
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
共檢測(cè)到158個(gè)二級(jí)質(zhì)譜，有些色譜峰沒有二級(jí)質(zhì)譜恳啥，如CP02偏灿，有些色譜峰又有好幾個(gè)二級(jí)質(zhì)譜，如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]])
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")
整合同一個(gè)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]])
將整合后的碎片信息光譜與整合前做對(duì)比
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]])
比對(duì)評(píng)分
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進(jìn)行LC-MS數(shù)據(jù)預(yù)處理和分析 ? xcms (sneumann.github.io)
Chromatographic peak detection using the centWave method — findChromPeaks-centWave ? xcms (sneumann.github.io)
centWave Parameters ? msbrowser (tkimhofer.github.io)
