之前在公众号里已经粗略的介绍过一个怎么用snp callingc出来 的vcf文件画进化树
之前见过师兄师姐们用蛋白序列画进化树,基于氨基酸序列的差异,通过比较序列类似性,基于氨基酸相对突变率矩阵计算不同序列差异性积分作为它们的差异性量度。一般都是用Clustal进行多序列比对,转化为mega的输入格式来画树(我很少画蛋白的树,如果有什么不对的请各位小伙伴多包涵)。
用核酸序列画树也是基于SNP和Indel这些差异信息来的。根据SNP或者Indel 构建其系统进化树,可以清晰的展示群体中不同个体的相互关系。
1.构建进化树的算法
- 基于距离矩阵的方法:NJ(邻接法)
- MP(最大简约法)
- ML(最大似然法)
- 贝叶斯法
具体用选择哪种算法进行画树,和你的数据有关,一般近缘序列的话,使用MP;远源序列,一般使用NJ或者ML。
在实战的时候,我们通常都要对raw vcf文件进行初步的过滤,如果样本数很大,vcf文件也会很大,而且一些低质量的位点会影响距离的估算,而且还有很多的多态位点也会带来很大的误差。
我一般都是过滤掉低qual值和多态位点,snp和indel数据也会分开,直接用filter后的snp数据作为画树的初始数据。
1构建进化树的方法
方法一:
将vcf 格式转成phylip格式
[phylip下载地址]:(https://github.com/edgardomortiz/vcf2phylip)
这个包就是一个python的脚本,解压后就可以直接用
不过这个包调用需要python2,python3会报错
/usr/bin/python2.7 /data/home/mjchen/app/vcf2phylip-master/vcf2phylip.py --input Four_parent_F2.ref_beagle.vcf
用mega7把转化好的Four_parent_F2.ref_beagle.min4.phy转化为画树的mega文件
最后就可以直接用mega文件画进化树了,一般的画进化树的软件都可以用
方法二:
生成tfam和tped文件
这一步用plink实现或者自己写脚本实现都可以,我用的是自己写的脚本
- 缺失数据用 “0”
- 与ref一直用 “1”
- 与alt一致用 “2”
具体每一列代表什么意思,在我之前的文章里面已经介绍过了,这里就不重复了
用emmax计算IBS的kinship
/data/home/mjchen/app/emmax-beta-07Mar2010/emmax-kin -v -h -s -d 10 Four_parent_F2.ref_beagle_tr
转化为画neighbor软件的输入文件
这一步也是通过自己写脚本实现的
就是在kinship矩阵的首行加入一行样品数目,接下来的数字行都为(1-IBS)的结果
用neighbor.exe创建画树的文件
直接输入文件名,按提示操作就好
形成这个文件后,就可以直接用mega打开画树了
方法三:
实际上这个方法和第二种基本是一致的,区别在于计算kinship的方式不一样,方法二用的emmax,三用的plink
plink --file Four_parent_F2.ref_beagle_tr --maf 0.05 --map3 --noweb --cluster --distance-matrix --out Four_parent_F2.ref_beagle_tr
这里需要map和ped的输入文件,不是.tped和.tfam了
具体格式 参考链接里的介绍
[plink格式介绍]:(http://www.cog-genomics.org/plink/1.9/formats)
接下来就和方法二一样了
这是我所知道的几种画进化树的方式,有其他方法的欢迎小伙伴们补充