LCMS data preprocessing and analysis with xcms
Package: xcms
!官網(wǎng)
1. 介紹
該文件描述了使用xcms >= 3版本的LCMS實驗的數(shù)據(jù)導(dǎo)入、探索、預(yù)處理和分析。這些例子和基本工作流程改編自Colin A. Smith的LC/MS預(yù)處理和分析小冊子。
新的用戶界面和方法使用XCMSnExp對象（而不是舊的xcmsSet對象）作為預(yù)處理結(jié)果的容器。為了支持依賴xcmsSet對象的包和管道，可以使用as方法將XCMSnExp轉(zhuǎn)換成xcmsSet對象（即xset <- as(x, "xcmsSet")），x是一個XCMSnExp對象。
2. 數(shù)據(jù)導(dǎo)入（讀數(shù)據(jù)）
xcms支持分析來自（AIA/ANDI）NetCDF、mzXML和mzML格式文件的LC/MS數(shù)據(jù)。對于實際的數(shù)據(jù)導(dǎo)入，使用了Bioconductor的mzR。為了演示，我們將分析[1]中的一個數(shù)據(jù)子集，其中研究了敲除小鼠脂肪酸酰胺水解酶（FAAH）基因的代謝后果。原始數(shù)據(jù)文件（NetCDF格式）由faahKO數(shù)據(jù)包提供。該數(shù)據(jù)集由6只敲除型和6只野生型小鼠的脊髓樣本組成。每個文件都包含在正離子模式下獲得的200-600 m/z和2500-4500秒的中心模式的數(shù)據(jù)。為了加快這個小冊子的處理速度，我們將把分析限制在只有8個文件和2500-3500秒的保留時間范圍內(nèi)。
下面我們加載所有需要的軟件包，在faahKO軟件包中找到原始CDF文件，并建立一個描述實驗設(shè)置的phenodata數(shù)據(jù)框。請注意，對于真正的實驗，建議定義一個文件（表），包含原始數(shù)據(jù)文件的文件名以及每個文件的樣本描述作為附加列。這樣的文件可以用例如read.table作為變量pd導(dǎo)入（而不是像下面的例子那樣在R中定義），文件名可以傳遞給下面的readMSData函數(shù)，例如files = paste0(MZML_PATH, "/", pd$mzML_file) 其中MZML_PATH是文件所在目錄的路徑，"mzML_file "是phenodata文件中包含文件名的列名。
初步接觸和摸索代謝或者蛋白組的數(shù)據(jù)，還是萌萌的。mzXML這個格式我是搞定了的，下機的‘.raw’，被我轉(zhuǎn)成mzXML之后讀入進行分析了，這里還是尊重官網(wǎng)教程走一遍，后面有機會歇一歇怎么轉(zhuǎn)數(shù)據(jù)吧。
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方法將原始數(shù)據(jù)加載為OnDiskMSnExp對象。MSnbase提供了處理質(zhì)譜數(shù)據(jù)的基礎(chǔ)結(jié)構(gòu)和基礎(chǔ)設(shè)施。此外，MSnbase還可以用來對剖面模式的質(zhì)譜數(shù)據(jù)進行中心化處理（見MSnbase包中的相應(yīng)小冊子）。
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.
接下來我們將數(shù)據(jù)集限制在2500到3500秒的保留時間范圍內(nèi)。這只是為了減少這個小插曲的處理時間。也不知道這個閾值怎么確定的
raw_data <- filterRt(raw_data, c(2500, 3500))
由此產(chǎn)生的OnDiskMSnExp對象包含關(guān)于光譜數(shù)量、保留時間、測量的總離子電流等一般信息，但不包含完整的原始數(shù)據(jù)（即每個測量光譜的m/z和強度值）。因此，它的內(nèi)存占用相當(dāng)小，使其成為代表大型代謝組學(xué)實驗的理想對象，同時允許執(zhí)行簡單的質(zhì)量控制、數(shù)據(jù)檢查和探索以及數(shù)據(jù)子集操作。m/z和強度值是按要求從原始數(shù)據(jù)文件中導(dǎo)入的，因此在初始數(shù)據(jù)導(dǎo)入后，原始數(shù)據(jù)文件的位置不應(yīng)改變。
3.初始數(shù)據(jù)檢查（相當(dāng)于質(zhì)控了）
OnDiskMSnExp按光譜組織MS數(shù)據(jù)，并提供強度、mz和rtime方法來訪問文件中的原始數(shù)據(jù)（測量的強度值、相應(yīng)的m/z和保留時間值）。此外，光譜方法可用于返回封裝在光譜對象中的所有數(shù)據(jù)。下面我們從對象中提取保留時間值。
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
所有的數(shù)據(jù)都以一維向量的形式返回（一個用于rtime的數(shù)字向量和一個用于mz和強度的數(shù)字向量列表，每個向量包含一個光譜的值），即使實驗由多個文件/樣品組成。fromFile函數(shù)返回一個整數(shù)向量，提供對原始文件的數(shù)值的映射。下面我們使用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
作為對數(shù)據(jù)的第一次評估，我們在下面為實驗中的每個文件繪制了基礎(chǔ)峰色譜圖（BPC）。我們使用色譜方法，并將aggregationFun設(shè)置為 "max"，以返回每個光譜的最大強度，從而從原始數(shù)據(jù)中創(chuàng)建BPC。為了創(chuàng)建一個總的離子色譜圖，我們可以將aggregationFun設(shè)置為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對象，該對象將各個色譜對象（實際上包含色譜數(shù)據(jù)）組織成一個二維數(shù)組：列代表樣品，行（可選）代表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
色譜方法也支持從質(zhì)譜數(shù)據(jù)的m/z-rt切片中提取色譜數(shù)據(jù)。在下一節(jié)中，我們將使用這種方法為一個選定的峰創(chuàng)建一個提取的離子色譜圖（EIC）。
請注意，chromatogram會從每個文件中讀取原始數(shù)據(jù)來計算色譜圖。另一方面，bpi和tic方法不從原始文件中讀取任何數(shù)據(jù)，而是使用輸入文件的頭定義中提供的相應(yīng)信息（可能與實際數(shù)據(jù)不同）。
下面我們創(chuàng)建了代表每個文件的總離子電流分布的圖表。這樣的圖對于發(fā)現(xiàn)有問題或失敗的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")
另外，我們可以根據(jù)樣品基峰色譜的相似性對其進行分組。這也有助于發(fā)現(xiàn)實驗中可能有問題的樣品，或者一般來說，對實驗中的樣品分組有一個初步的了解。由于樣品之間的保留時間并不完全相同，我們使用bin函數(shù)沿保留時間軸將強度分組在固定的時間范圍（bin）內(nèi)。在本例中，我們使用的bin大小為1秒，默認為0.5秒。聚類是在分檔后的基峰色譜的成對相關(guān)性上使用完全聯(lián)系的層次結(jié)構(gòu)聚類法進行的。
## 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))
這些樣本以成對的方式聚類，樣本指數(shù)的KO和WT樣本具有最相似的BPC。
4.色譜峰檢測
接下來我們使用centWave算法[2]進行色譜峰值檢測。然而，在運行峰值檢測之前，強烈建議目視檢查，例如內(nèi)標或已知化合物的提取離子色譜圖，以評估和調(diào)整峰值檢測設(shè)置，因為默認設(shè)置將不適合大多數(shù)LCMS實驗。centWave最關(guān)鍵的兩個參數(shù)是峰寬（色譜峰寬的預(yù)期范圍）和ppm（對應(yīng)于一個色譜峰的中心點的m/z值的最大預(yù)期偏差；這通常比制造商指定的ppm大很多）參數(shù)。為了評估典型的色譜峰寬，我們繪制一個峰的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])
請注意，如果在某次掃描中（即對于特定的保留時間）沒有在相應(yīng)的mz范圍內(nèi)測量到信號，則通過色譜法提取的色譜對象包含一個NA值。這反映在上面的圖中，這些線沒有被畫成連續(xù)線。
上面的峰的寬度約為50秒。峰寬參數(shù)的設(shè)置應(yīng)適應(yīng)數(shù)據(jù)集中預(yù)期的峰寬。對于本例數(shù)據(jù)集，我們將其設(shè)置為20,80。
對于ppm參數(shù)，我們提取與上述峰值相對應(yīng)的完整MS數(shù)據(jù)（強度、保留時間和m/z值）。為此，我們首先按保留時間過濾原始對象，然后按m/z過濾，最后用type = "XIC "繪制對象，產(chǎn)生下面的圖。我們使用管道（%>%）命令更好地說明了相應(yīng)的工作流程。還要注意的是，在這種類型的繪圖中，如果存在已識別的色譜峰，將默認顯示。
raw_data %>%
filterRt(rt = rtr) %>%
filterMz(mz = mzr) %>%
plot(type = "XIC")
對于每個圖：上圖：色譜圖繪制強度值與保留時間的關(guān)系，下圖m/z與保留時間的關(guān)系。各個數(shù)據(jù)點根據(jù)強度的不同而著色。
在目前的數(shù)據(jù)中，m/z值實際上沒有變化。通常情況下，人們會看到m/z值（下圖）在化合物的真實m/z值周圍散開。centWave算法的第一步是根據(jù)連續(xù)掃描的m/z值的差異來定義所謂的感興趣區(qū)域（ROI）。詳細地說，如果連續(xù)掃描的m/z值與ROI的平均m/z值之間的差異小于用戶定義的ppm參數(shù)，則來自連續(xù)掃描的m/z值被納入ROI。因此，對ppm的合理選擇可以是作為色譜峰一部分的相鄰掃描/譜段的數(shù)據(jù)點的最大m/z差。建議檢查許多化合物（內(nèi)部標準或已知存在于樣品中的化合物）的m/z值范圍，并根據(jù)這些定義centWave的ppm參數(shù)。
請注意，我們也可以對提取的離子色譜圖進行峰值檢測。這可以幫助評估不同的峰值檢測設(shè)置。請注意，提取的離子色譜圖上的峰值檢測將不考慮ppm參數(shù)，對背景信號的估計與完整數(shù)據(jù)集上的峰值檢測不同；因此，snthresh的值將有不同的后果。下面我們用findChromPeaks函數(shù)對提取的離子色譜圖進行峰值檢測。提交的參數(shù)對象定義了將使用的算法，并允許定義該算法的設(shè)置。我們使用默認設(shè)置的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矩陣平行的還有一個數(shù)據(jù)框chromPeakData，可以為每個色譜峰添加任意的注釋。下面我們提取這個數(shù)據(jù)框，默認情況下，它只包含鑒定峰的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參數(shù)為各個色譜線著色。也可以為已識別的色譜峰指定顏色，peakCol為前景/邊框顏色，peakBg為背景/填充顏色。必須為chromPeaks列出的每個色譜峰提供一種顏色。下面我們定義一種顏色來表示樣品所在的樣品組，并使用峰的 "樣品 "欄中的樣品信息來為每個色譜峰分配正確的顏色。更多的色譜峰高亮選項將在下文進一步描述。
sample_colors <- group_colors[xchr$sample_group]
plot(xchr, col = sample_colors,
peakBg = sample_colors[chromPeaks(xchr)[, "column"]])
紅色和藍色分別代表KO和野生型樣品。確定的色譜峰的峰面積在樣品組顏色中突出顯示。
最后我們對完整的數(shù)據(jù)集進行色譜峰檢測。請注意，我們將參數(shù)prefilter設(shè)置為c(6, 5000)，噪聲設(shè)置為5000，以減少本小節(jié)的運行時間。通過這個設(shè)置，我們在峰值檢測步驟中只考慮數(shù)值大于5000的信號。
cwp <- CentWaveParam(peakwidth = c(20, 80), noise = 5000,
prefilter = c(6, 5000))
xdata <- findChromPeaks(raw_data, param = cwp)
結(jié)果以XCMSnExp對象形式返回，該對象通過存儲LC/GC-MS預(yù)處理結(jié)果擴展了OnDiskMSnExp對象。這也意味著所有用于子集和過濾數(shù)據(jù)或訪問（原始）數(shù)據(jù)的方法都是從OnDiskMSnExp對象繼承的，因此可以重復(fù)使用。還要注意的是，通過調(diào)用參數(shù)add = TRUE的findChromPeaks，可以在xdata對象上執(zhí)行額外的幾輪峰值檢測（例如，在MS級別>1的數(shù)據(jù)上）。
色譜峰檢測的結(jié)果可以通過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函數(shù)提取。這個數(shù)據(jù)框也可以用來為每個檢測到的峰添加/儲存任意的注釋。
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
峰值檢測并不總是完美的，會導(dǎo)致峰值檢測的假象，如重疊的峰值或人為分割的峰值。refineChromPeaks功能允許通過移除不符合特定標準的識別峰或合并人為分割的色譜峰來完善峰值檢測結(jié)果。通過參數(shù)對象CleanPeaksParam和FilterIntensityParam，可以分別移除保留時間范圍內(nèi)的峰或強度低于閾值的峰（更多細節(jié)和例子見各自的幫助頁面）。通過MergeNeighboringPeaksParam，可以合并色譜峰。下面我們對色譜峰檢測結(jié)果進行后處理，如果重疊的色譜峰之間的信號低于較小的色譜峰最大強度的75%，則在每個文件的4秒窗口中合并這些色譜峰。請參閱MergeNeighboringPeaksParam幫助頁面，了解設(shè)置和方法的詳細說明。
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個全區(qū)峰和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)
還要注意的是，可以對提取的離子色譜圖進行峰的細化。例如，這可以用來微調(diào)參數(shù)的設(shè)置。為了說明這一點，我們在下面對提取的離子色譜圖chr_1進行峰值細化，減少minProp參數(shù)以強制連接兩個峰。
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)
[to be continued]
