R语言实现质谱峰的可视化
今天给大家介绍一下如何把我们在质谱中的原始数据进行质谱峰的可视化展示，质谱峰代表了在某一个质荷比值的峰强度。那么如何在R语言中重建这些峰呢。我们今天就来给大家介绍下基于MALDIquant包和MALDIquantForeign包的实现过程。这两个包关系是MALDIquantForeign包主要为MALDIquant包提供数据的导入以及导出，从而实现原始数据的分析质控。那么我们还是老套路：安装，函数介绍，实例。
首先我们看下包的安装，无需多言CRAN资源直接：
install.packages("MALDIquant")
install.packages("MALDIquantForeign")
library("MALDIquant")
library("MALDIquantForeign")接下来我们介绍下数据的导入导出：
1. import()主要数据的导入基于MALDIquantForeign包。
主要的参数：
Type主要是所能识别的文件类型，其实也可以不用设置。
txt importTxt
tab importTab
csv importCsv
fid importBrukerFlex
ciphergen importCiphergenXml
mzXML importMzXml
mzML importMzMl
imzML importImzMl
analyze importAnalyze
cdf importCdf
msd importMsdremoveEmptySpectra 是否移除空的数据，默认是TRUE。
Centroided此参数需要设为FALSE，因为后面不识别此质心算法出来后的数据。
实例:
exampleDirectory <- system.file("exampledata",package="MALDIquantForeign")
## import mzXML files
s <- import(exampleDirectory, type="mzXML")
2. export()主要数据的导出基于MALDIquantForeign包。
其主要的参数是type，也就是所支持的文件类型。
tab exportTab
csv exportCsv
imzML exportImzMl
msd exportMsd
mzML exportMzMl实例：
export(s[[1]], file="spectrum.csv")
通过以上步骤导入的数据，可以在MALDIquant中进行数据的质控以及可视化：
首先是利用基础函数进行两个判断：
1. any(sapply(s,isEmpty)) #判断是否存在空的峰
2. table(sapply(s,length))#查看每个文件峰的数量以及总文件数，同时判定每个文件是否峰数一致。
3. all(sapply(s,isRegular))#判断是否经过质心算法处理。TRUE代表没有处理过。接下来是数据的峰可视化，为了展示数据可视化的可预见性，我们使用包自带的数据。
data(fiedler2009subset)
fiedler2009subset[1:2]
plot(fiedler2009subset[[1]])
plot(fiedler2009subset[[16]])
色谱峰的修饰： 1. transformIntensity减少峰的波动，客服方差均值的相互影响。通俗讲就是标准化，其对应的方法我们就不介绍了。
实例：
spectra <- transformIntensity(fiedler2009subset,method="sqrt")2. smoothIntensity主要是对峰进行平滑操作，使峰更加显著。
其中主要的涉及两个平滑曲线的滤波器：SavitzkyGolay 和MovingAverage。主要的参数设置是halfWindowSize，此值可以根据IPO中计算的FWHM（半峰宽度）进行选择，规则就是要小于FWHM。
实例:
spectra <- smoothIntensity(spectra,method="SavitzkyGolay",halfWindowSize=10)3. estimateBaseline 基线的评估。并且提供了以下的几种方法，我们默认就好了。这里还有一个隐藏的参数，那就是迭代次数iterations。
实例：
baseline <- estimateBaseline(spectra[[16]],method="SNIP",iterations=100)
plot(spectra[[16]])
lines(baseline,col="red", lwd=2)
4. removeBaseline 主要是去除每一个样本的基线，使得峰更加统一。具体的算法，我们不再去细讲。
实例：
spectra <- removeBaseline(spectra,method="SNIP",iterations=100)
plot(spectra[[16]])
plot(spectra[[1]])
从上面的图我们可以看到，他们的强度值其实都不一样的，为了客服小批次的影响，需要我们进一步的校正，那就用到另一个函数calibrateIntensity：
参数不再去多说，我们看下实例：
spectra <- calibrateIntensity(spectra,method="TIC")
plot(spectra[[16]])
plot(spectra[[1]])
接下来就噪音的消除，我们需要用到alignSpectra，此函数可以确定我们的噪音峰的强度。我们直接看下实例：
spectra<-alignSpectra(spectra,halfWindowSize=20,SNR=2,tolerance=0.002,warpingMethod="lowess")
samples <- factor(sapply(spectra,function(x)metaData(x)$sampleName))#获取样本信息
avgSpectra <- averageMassSpectra(spectra, labels=samples,method="mean")
#每组数据的均值
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，进行接下来的峰的发现。函数我们就不展开了，直接进行实例：
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)
然后就是把所有的峰值进行保存：
peaks <- binPeaks(peaks, tolerance=0.002)# tolerance=5e-6代表误差5ppm至此我们就得到了对应的MZ和相对应的峰强度。最后数据的导出：
export(peaks[[1]], file="spectrum.csv",force=TRUE)
欢迎学习交流！
腾讯云开发者
