用R语言对vcf文件进行数据挖掘.11 CNV分析

Python027

用R语言对vcf文件进行数据挖掘.11 CNV分析,第1张

目录

在之前的文章里介绍了如何通过直方图来可视化等位杂合碱基的比例来判断物种的染色体倍数性。在本文里会继续向下挖掘,介绍如何可视化染色体上的拷贝数变化(CNVs)。

和前文一样的操作,使用包自带的数据

我们需要去除过高和过低深度的数据。和前文的操作一样,提取vcf文件里的深度数据"AD"。

然后过滤出10%~90%的数据,当然此处可以根据实际情况进行微调。然后对第一种出现频率最高的碱基进行可视化。(一般情况下一个位点上会有两种碱基,具体参考前文。)

同样也可以对出现频率第二高的碱基进行同样的操作,这里节约篇幅就省略了。

为了避免复杂的基于AD比例的模型假设,程序里设计了非参数估计法来计算峰值。计算完了以后可以直接对染色体进行拆分以后可视化进行校验。

根据尺寸把染色体分割成合适的大小

然后用 freq_peak 函数计算峰值。并对数据进行处理,去掉负数和Na值。

计算到此为止,可以可视化实际数据来验证计算的正确性。

仔细想一下,峰值计算的结果其实就是CNV的结果。这里根据窗口大小把染色体分成了若干段。(那么是不是可以给每一段 CDS进行细分然后计算出每一个CDS的具体数字呢????)

当然也可以把所有样本组合到一起。

什么是R语言?应该如何开始学习/使用R语言呢?

学习R有几个月了,总算是摸着了一点门道。

写一些自己的心得和经验,方便自己进一步鼓捣R。如果有人看到我写的东西而得到了帮助,那就更好了。

什么是R?R的优点何在?

R是一个数据分析软件。简单点说,R可以看做MATLAB的“替代品”,而且具有免费开源的优势。R可以像MATLAB一样解决有关数值计算的问题,而且具有强大的数据处理,绘图功能。

R拥有大量的统计分析工具包,我的感觉是——只有我们没听说过的工具,绝对没有R没有的工具包。配合着各种各样的工具包,你可以毁灭任何关于数据和统计的问题。因为数据包的数量庞大,所以查找自己需要的数据包,可能很烦恼。

如果有以下技能,学R会很方便:

1.已经了解些高级程序语言(非常重要)

2.英语不坏

3.概率统计理论基础

4.看数据不头疼

5.看cmd or terminal 也不头疼

你需要一本适合你的R语言教材

我开始学习R的时候,找到了这个帖子

非常强大的关于R语言教材综述。我非常感谢原帖作者。你可以参考这个帖子选一本适合你的教材。

我这里在说一下我主要使用的几本教材的心得:

1. 统计建模与R软件(薛毅著):非常优秀的R语言入门教材,涵盖了所有R的基础应用&方法,示例代码也很优秀。作为一本中文的程序语言教材,绝对是最优秀的之一。但是要看懂这本书,还是需要“已经了解些高级程序语言”。PS:我亲爱的吉林大学图书馆,有两本该教材流通,我常年霸占一本。

2. R in Nutshell:从讲解内容上看,与上一本差别不大,在R语言的应用上都是比较初级的入门,但是有些R软件&语言上的特性,写得比薛毅老师的教材深刻。这本书最大的优点就是工具书,方便开始入门时候,对有些“模棱两可”的东西的查询。PS:我将这本书打印了出来,简单的从头到尾翻过,最大的用途就是像一本字典一样查询。

3. ggplot2 Elegant Graphics for Data:这是一本介绍如何使用ggplot2包,进行绘图的书。ggplot2包,非常强大的绘图工具,几乎可以操作任何图中的元素,而且是提供添加图层的方式让我们可以一步步的作图。提到ggplot2包,应该提到一个词——“潜力无穷”,每一个介绍

ggplot2的人,都会用这个形容词。这本书最大的作用也是当做一本绘图相关的工具书,书中讲解详细,细致,每个小参数的变动都会配图帮你理解。PS:这本书我也打印出来了,非常适合查询。

几个可以逐步提高R能力的网站

1.R-bloggers: 这里有关于R和数据的一切讨论,前沿的问题,基础的问题,应有尽有。可以说这些家伙们让R变得越来越强大。我RSS了这个网站,每天都看一下有什么我感兴趣的方法和话题,慢慢的积累一些知识,是一个很有意思的过程。

2.统计之都: 这是一个有大量R使用者交流的论坛,你可以上去提问题,总有好心人来帮助你的。

3.R客: 是关于R的一个博客,更新不快,偏重国内R的一些发展。

R的使用环境

如果你看见terminal or cmd就打怵的话,一定要使用Rstudio。Rstudio的优点是,集成了Rconsole、脚本编辑器、可视化的数据查询、历史命令、帮助查询等,还有的完美的脚本和console的互动。毕竟是可视化的界面,有许多按钮可以用。R 的脚本编辑器很蛋疼,就比记事本多了个颜色高亮吧,不适合编写脚本,但适合调试脚本。

最后,说一下,刚开始学习R或者其他什么语言,都有一个通病,就是一些小细节的不知道,或者是记得不清楚,往往一个蛋疼的bug就可以耗掉大量的时间,这是一个让人想砸电脑的过程。我往后,会在博客里记录一些让我蛋很疼的小细节。本文分为6个部分,分别介绍初级入门,高级入门,绘图与可视化,计量经济学,时间序列分析,金融等。

1.初级入门

《An Introduction to R》,这是官方的入门小册子。其有中文版,由丁国徽翻译,译名为《R导论》。《R4Beginners》,这本小册子有中文版应该叫《R入门》。除此之外,还可以去读刘思喆的《153分钟学会R》。这本书收集了R初学者提问频率最高的153个问题。为什么叫153分钟呢?因为最初作者写了153个问题,阅读一个问题花费1分钟时间,全局下来也就是153分钟了。有了这些基础之后,要去读一些经典书籍比较全面的入门书籍,比如《统计建模与R软件》,国外还有《R Cookbook》和《R in action》,本人没有看过,因此不便评论。

最后推荐,《R in a Nutshell》。对,“果壳里面的R”!当然,是开玩笑的,in a Nutshell是俚语,意思大致是“简单的说”。目前,我们正在翻译这本书的中文版,大概明年三月份交稿!这本书很不错,大家可以从现在开始期待,并广而告知一下!

2.高级入门

读了上述书籍之后,你就可以去高级入门阶段了。这时候要读的书有两本很经典的。《Statistics with R》和《The R book》。之所以说这两本书高级,是因为这两本书已经不再限于R基础了,而是结合了数据分析的各种常见方法来写就的,比较系统的介绍了R在线性回归、方差分析、多元统计、R绘图、时间序列分析、数据挖掘等各方面的内容,看完之后你会发现,哇,原来R能做的事情这么多,而且做起来是那么简洁。读到这里已经差不多了,剩下的估计就是你要专门攻读的某个方面内容了。下面大致说一说。

3.绘图与可视化

亚里斯多德说,“较其他感觉而言,人类更喜欢观看”。因此,绘图和可视化得到很多人的关注和重视。那么,如何学习R画图和数据可视化呢?再简单些,如何画直方图?如何往直方图上添加密度曲线呢?我想读完下面这几本书你就大致会明白了。

首先,画图入门可以读《R Graphics》,个人认为这本是比较经典的,全面介绍了R中绘图系统。该书对应的有一个网站,google之就可以了。更深入的可以读《Lattice:Multivariate Data Visualization with R》。上面这些都是比较普通的。当然,有比较文艺和优雅的——ggplot2系统,看《ggplot2:Elegant Graphics for Data Analysis》。还有数据挖掘方面的书:《Data Mining with Rattle and R》,主要是用Rattle软件,个人比较喜欢Rattle!当然,Rattle不是最好的,Rweka也很棒!再有就是交互图形的书了,著名的交互系统是ggobi,这个我已经喜欢两年多了,关于ggobi的书有《Interactive and Dynamic Graphics for Data Analysis With R and GGobi》,不过,也只是适宜入门,更多更全面的还是去ggobi的主页吧,上面有各种资料以及包的更新信息!

特别推荐一下,中文版绘图书籍有《现代统计图形》。

4.计量经济学

关于计量经济学,首先推荐一本很薄的小册子:《Econometrics In R》,做入门用。然后,是《Applied Econometrics with R》,该书对应的R包是AER,可以安装之后配合使用,效果甚佳。计量经济学中很大一部分是关于时间序列分析的,这一块内容在下面的地方说。

5.时间序列分析

时间序列书籍的书籍分两类,一种是比较普适的书籍,典型的代表是:《Time Series Analysis and Its Applications :with R examples》。该书介绍了各种时间序列分析的经典方法及实现各种经典方法的R代码,该书有中文版。如果不想买的话,建议去作者主页直接下载,英文版其实读起来很简单。时间序列分析中有一大块儿是关于金融时间序列分析的。这方面比较流行的书有两本《Analysis of financial time series》,这本书的最初是用的S-plus代码,不过新版已经以R代码为主了。这本书适合有时间序列分析基础和金融基础的人来看,因为书中关于时间序列分析的理论以及各种金融知识讲解的不是特别清楚,将极值理论计算VaR的部分就比较难看懂。另外一个比较有意思的是Rmetrics推出的《TimeSeriesFAQ》,这本书是金融时间序列入门的东西,讲的很基础,但是很难懂。对应的中文版有《金融时间序列分析常见问题集》,当然,目前还没有发出来。经济领域的时间序列有一种特殊的情况叫协整,很多人很关注这方面的理论,关心这个的可以看《Analysis of Integrated and Cointegrated Time Series with R》。最后,比较高级的一本书是关于小波分析的,看《Wavelet Methods in Statistics with R》。附加一点,关于时间序列聚类的书籍目前比较少见,是一个处女地,有志之士可以开垦之!

6.金融

金融的领域很广泛,如果是大金融的话,保险也要被纳入此间。用R做金融更多地需要掌握的是金融知识,只会数据分析技术意义寥寥。我觉得这些书对于懂金融、不同数据分析技术的人比较有用,只懂数据分析技术而不动金融知识的人看起来肯定如雾里看花,甚至有人会觉得金融分析比较低级。这方面比较经典的书籍有:《Advanced Topics in Analysis of Economic and Financial Data Using R》以及《Modelling Financial Time Series With S-plus》。金融产品定价之类的常常要用到随机微分方程,有一本叫《Simulation Inference Stochastic Differential Equations:with R examples》的书是关于这方面的内容的,有实例,内容还算详实!此外,是风险度量与管理类。比较经典的有《Simulation Techniques in Financial Risk Management》、《Modern Actuarial Risk Theory Using R》和《Quantitative Risk Management:Concepts, Techniques and Tools》。投资组合分析类和期权定价类可以分别看《Portfolio Optimization with R》和《Option Pricing and Estimation of Financial Models with R》。

7.数据挖掘

这方面的书不多,只有《Data Mining with R:learing with case studies》。不过,R中数据挖掘方面的包已经足够多了,参考包中的帮助文档就足够了。

目录

vcfR 可以直接读取vcf格式的数据。如果同时读取参照序列fasta格式的序列文件和gff格式文件的注释文件还可以获取更完整的信息(此步骤并非必须,可以只读取vcf数据)。在此处便于重复用到了 pinfsc50 包。这个包里是植物致病微生物的基因序列测序结果。包含了一个vcf文件,一个fasta文件和一个gff文件。

这里用到参照序列的数据。

当这些数据被读取到内存的时候就可以开始对染色体名字或者其它一些东西进行修改了。由于 vcfR 更擅长对的单独染色体进行分析,所以当你的基因过大或者有很多样本的时候,建议对数据进行拆分。

读取完数据以后就可以建立 chromR ,来对数据进行详细的分析。

首先对数据进行初步的可视化,

我们在上面的图里得到很多信息,比方说测序深度(DP)的峰在500,但是拖着尾巴,这个尾巴表示数据里包含着CNV信息。然后比对质量(MQ)的峰值在60,于是我们可以以60为中心对数据进行过滤。

使用 masker 可以对数据进行过滤标记。然后可视化过滤以后的数据。

是不是顺眼多了。当然我们也可以看一下SNP的分布情况。注意右下角的图。

用 chromoqc() 可以对数据进行更完整的可视化。包括外显子内含子的分布,GC含量的分布等等。

最后可以用函数 write.vcf() 把数据输出成新的vcf文件。