本来应该这是一个很正常的学习过程,之前总结了一篇博文Bioconductor的质谱蛋白组学数据分析,对蛋白组学定量那块比较感兴趣,正好看到一个R包-MSstats,其可用来对DDA,SRM和DIA的结果进行蛋白差异分析,这R包发表于2014年,那时来说还是很不错的(还在不断更新维护),并且其还支持Maxquant查库结果文件作为输入(主要我有些此类测试文件),非常有兴趣的想尝试下看看结果,然后就入坑了。。。
从其官网http://msstats.org/可看出,其现在的功能还是非常全面的,当然我只暂时用到其一小部分功能
正常流程一般是从Biocondutor上其R包的使用说明看起,然后拿自己的测试数据走一遍流程。但是!!从第一个函数MaxQtoMSstatsFormat
(用于将Maxquant查库结果文件转化为其标准输入格式)起就遇到了一个槛,由于其需要一个annotation文件(感觉类似于甲基化芯片分析的Champ包那个注释文件),但其官网文档并没有给出完整的格式,翻了一些资料并看了这函数的源码才整理出一份正确的格式(我的是Label-free分级测试文件),如下:
MSstate_annotation
本以为后续应该可以比较正常的进行下去了,结果在第二个对输入数据进行处理的函数dataProcess
就开始不停的报错,我本以为是我测试数据的问题,结果换了好几批数据后才觉得,可能是我的数据跟其认为的maxquant的正常输出结果不一致,因此决定从其源码开始找寻报错的原因,结果从其代码中发现了几处'BUG'(至少对我的测试数据来说)
其实对R包这种开源的软件来说,如果遇到报错信息,首先想到的是查阅官网网站或文档寻找答案;再者去网上寻找是否有跟自己遇到相同问题的人,并且看看是否有对应的解决办法;最后如果还是不行的话,那么就看R包对应函数的源码吧(至少能解决不少问题)。其实有些R包并不复杂,而且看源码的过程也是一种学习的过程,等以后自己写R包的时候也能用上一些技巧嘛
下面则是我看了MSstats
包的几个重要函数后的随笔,记录了个人理解下的其运行的原理(主要其发表的文章中并未提起原理部分,只是粗略讲了软件的功能以及结果展现形式,所以还是得从代码入手),但有些地方还是不太理解
下面这段摘抄自其官网文档,主要讲了这个函数做了哪些事情:
Data pre-processing and quality control of MS runs of the original raw data into quantitative data for model fitting and group comparison. Log transformation is automatically applied and additional variables are created in columns for model fitting and group comparison process.
一个样本,循环,第一个run1开始,对应去重后的肽段,统计每个样本在不同run下的有abundance的肽段个数,将同个样本中最大数目肽段的run对应的FACTION赋予对应的run1的肽段,循环至最后一个run,然后停止。
其原理我认为是这样的:在MS分级中,会将样本等量分成N份,对应的就是N个run,从而将不同特性的离子区分开;那么在某个run里鉴定到肽段数目一般相比其他run来说应该是最大的
代码中有个问题,其只循环了样本1,而其他样本则不管了,这里会导致一些肽段万一没有在样本1出现,那么其对应的FRACTION则是NA,会导致代码报错;但如果再循环每个样本的话,后面的样本会影响前面样本的结果,个人觉得应该每个样本各自确定自己RUN对应的FRACTION即可。最终我觉得这步骤主要还是为了确定FRACTION而已
我这里有个疑问,为啥不根据结果的RUN的那列数据来确定FRACTION呢,至少命名上还是可以看出不同RUN到底是属于哪个FRACTION的。PS.如果做了5次分级,6个样本,那么其实了跑了30次LC/MS,那么有30个RUN,但是FRACTION还只是5次哦
对于同一肽段在不同FRACTION下有不同的丰度(abundance),那么该如何确定哪个abundance应该赋予这个肽段,这个也是我之前一直想知道的,作者则是做了以下处理:
经过上述步骤过滤完后,则肽段只有在某一FRACTION下才有丰度值,其实主要是针对分级的Label-free数据而言的(PS. 真心觉得Label-free还是分级的好,不分级鉴定数目还是少了点);而言也做了Logarithm transformation,接下来则是标准化了,作者主要给了以下几个选项:
normalize.quantiles
函数,原理参照分位数标准化除了标准化外,作者还用了Tukey's median polish(TMP)算法确定一个cutoff,根据这个阈值过滤掉低丰度值的肽段(即将这些肽段的丰度值变为NA),涉及的默认参数是summaryMethod
,censoredInt
和MBimpute
此外,作者对于那些Protein Group只有一个肽段的(也就是说这个Protein Group只是由一个肽段定性出来的,可信度比较低。。。),并没有舍去,只是给了提醒,让用户考虑是否要舍去
最后这步是这个函数比较重要的一个统计分析,其以一个Protein Group为一个单位,依次进行分析(也就是一个Protein Group一次分析),主要做了以下处理:
summaryMethod
参数的模型不同,其结果有略有差异,主要分TMP和linear两大类remove50missing
参数,则会过滤掉那些缺失值超过50%的run最终的结果为一个processedquant列表用于后续分析
这个函数主要是用于可视化的,展示一个或者N个Protein Group下的肽段丰度分布,主要有箱线图和折线图等,用一下就知道了,主要输入对象是上面的processedquant
MSstate_boxplot
这个函数则比较重要了,Finding differentially abundant proteins across conditions,也就是差异蛋白分析,也是我一直比较想了解的分析过程。虽然使用很简单,但是既然上面几个步骤的源码都看了,也不差这一个函数了,顺便也看了一遍,归纳下主要是以下几个步骤:
使用的话,可以分为两步,先建立比较组矩阵(类似NGS那些差异分析软件的分组矩阵),然后就用groupComparison
进行分析了,如下:
comparison <- matrix(c(-1,1),nrow=1)
row.names(comparison) <- "B-A"
comparison_result <- groupComparison(contrast.matrix = comparison, data = processedquant)
head(comparison_result$ComparisonResult)
接下来我想通过源码了解作者是如何在Protein层面上进行差异分析(结果我关键的算法步骤还是没看懂,统计方面实在太差,以后再补上了),其他步骤大概是这样的:
折腾了一周,每天断断续续看了其代码,虽然代码比较长,但是很多部分是一些数据判定检查过程(可跳过),所以代码思路非常严谨,从其代码中也学到了不少小技巧,也算一个额外的收获。并且还了解了作者对于蛋白组定量差异分析的一些思路,可能不是最优的,但也是一个参考。虽然这个R包对我的数据来说无法正常使用(因为必须先修改其函数中的部分代码才行),但理解其思路才是最主要的!