R语言实现质谱峰的可视化

2019-07-31 14:36:56 浏览数 (1)

今天给大家介绍一下如何把我们在质谱中的原始数据进行质谱峰的可视化展示,质谱峰代表了在某一个质荷比值的峰强度。那么如何在R语言中重建这些峰呢。我们今天就来给大家介绍下基于MALDIquant包和MALDIquantForeign包的实现过程。这两个包关系是MALDIquantForeign包主要为MALDIquant包提供数据的导入以及导出,从而实现原始数据的分析质控。那么我们还是老套路:安装,函数介绍,实例。

首先我们看下包的安装,无需多言CRAN资源直接:

代码语言:javascript复制
install.packages("MALDIquant")
install.packages("MALDIquantForeign")
 
library("MALDIquant")
library("MALDIquantForeign")

接下来我们介绍下数据的导入导出:

1. import()主要数据的导入基于MALDIquantForeign包。

主要的参数:

Type主要是所能识别的文件类型,其实也可以不用设置。

代码语言:javascript复制
txt  importTxt
tab  importTab
csv  importCsv
fid  importBrukerFlex
ciphergen   importCiphergenXml
mzXML   importMzXml
mzML  importMzMl
imzML   importImzMl
analyze  importAnalyze
cdf   importCdf
msd   importMsd

removeEmptySpectra 是否移除空的数据,默认是TRUE。

Centroided此参数需要设为FALSE,因为后面不识别此质心算法出来后的数据。

实例:

代码语言:javascript复制
exampleDirectory <- system.file("exampledata",package="MALDIquantForeign")
 
## import mzXML files
s <- import(exampleDirectory, type="mzXML")
 

2. export()主要数据的导出基于MALDIquantForeign包。

其主要的参数是type,也就是所支持的文件类型。

代码语言:javascript复制
tab  exportTab
csv  exportCsv
imzML   exportImzMl
msd  exportMsd
mzML  exportMzMl

实例:

代码语言:javascript复制
export(s[[1]], file="spectrum.csv")

通过以上步骤导入的数据,可以在MALDIquant中进行数据的质控以及可视化:

首先是利用基础函数进行两个判断:

代码语言:javascript复制
1.  any(sapply(s,isEmpty)) #判断是否存在空的峰
2.  table(sapply(s,length))#查看每个文件峰的数量以及总文件数,同时判定每个文件是否峰数一致。
3.  all(sapply(s,isRegular))#判断是否经过质心算法处理。TRUE代表没有处理过。

接下来是数据的峰可视化,为了展示数据可视化的可预见性,我们使用包自带的数据。

代码语言:javascript复制
data(fiedler2009subset)
fiedler2009subset[1:2]
代码语言:javascript复制
plot(fiedler2009subset[[1]])
代码语言:javascript复制
plot(fiedler2009subset[[16]])

色谱峰的修饰: 1. transformIntensity减少峰的波动,客服方差均值的相互影响。通俗讲就是标准化,其对应的方法我们就不介绍了。

实例:

代码语言:javascript复制
spectra <- transformIntensity(fiedler2009subset,method="sqrt")

2. smoothIntensity主要是对峰进行平滑操作,使峰更加显著。

其中主要的涉及两个平滑曲线的滤波器:SavitzkyGolay 和MovingAverage。主要的参数设置是halfWindowSize,此值可以根据IPO中计算的FWHM(半峰宽度)进行选择,规则就是要小于FWHM。

实例:

代码语言:javascript复制
spectra <- smoothIntensity(spectra,method="SavitzkyGolay",halfWindowSize=10)

3. estimateBaseline 基线的评估。并且提供了以下的几种方法,我们默认就好了。这里还有一个隐藏的参数,那就是迭代次数iterations。

实例:

代码语言:javascript复制
baseline <- estimateBaseline(spectra[[16]],method="SNIP",iterations=100)
plot(spectra[[16]])
lines(baseline,col="red", lwd=2)

4. removeBaseline 主要是去除每一个样本的基线,使得峰更加统一。具体的算法,我们不再去细讲。

实例:

代码语言:javascript复制
spectra <- removeBaseline(spectra,method="SNIP",iterations=100)
 plot(spectra[[16]])
代码语言:javascript复制
plot(spectra[[1]])

从上面的图我们可以看到,他们的强度值其实都不一样的,为了客服小批次的影响,需要我们进一步的校正,那就用到另一个函数calibrateIntensity:

参数不再去多说,我们看下实例:

代码语言:javascript复制
spectra <- calibrateIntensity(spectra,method="TIC")
plot(spectra[[16]])
代码语言:javascript复制
plot(spectra[[1]])

接下来就噪音的消除,我们需要用到alignSpectra,此函数可以确定我们的噪音峰的强度。我们直接看下实例:

代码语言:javascript复制
spectra<-alignSpectra(spectra,halfWindowSize=20,SNR=2,tolerance=0.002,warpingMethod="lowess")

samples <- factor(sapply(spectra,function(x)metaData(x)$sampleName))#获取样本信息
代码语言:javascript复制
avgSpectra <- averageMassSpectra(spectra, labels=samples,method="mean")
#每组数据的均值
代码语言:javascript复制
 
noise <- estimateNoise(avgSpectra[[1]])#噪音评估
 
plot(avgSpectra[[1]], xlim=c(4000, 5000), ylim=c(0, 0.002))
lines(noise, col="red")
lines(noise[,1], noise[, 2]*2, col="blue")

从图中我们可以看出当信噪比为2时比较靠谱,能有效去除噪音峰。因此我们选择SNR=2,进行接下来的峰的发现。函数我们就不展开了,直接进行实例:

代码语言:javascript复制
peaks <- detectPeaks(avgSpectra, method="MAD",halfWindowSize=20,SNR=2)


plot(avgSpectra[[1]], xlim=c(4000, 5000), ylim=c(0, 0.002))
points(peaks[[1]], col="red", pch=4)

然后就是把所有的峰值进行保存:

代码语言:javascript复制
peaks <- binPeaks(peaks, tolerance=0.002)# tolerance=5e-6代表误差5ppm

至此我们就得到了对应的MZ和相对应的峰强度。最后数据的导出:

代码语言:javascript复制
export(peaks[[1]], file="spectrum.csv",force=TRUE)

欢迎学习交流!

0 人点赞