基于IBS矩阵在R语言中构建NJ进化树写出nwk文件-Toy模板网

这篇具有很好参考价值的文章主要介绍了基于IBS矩阵在R语言中构建NJ进化树写出nwk文件。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

构建系统发育树有很多方法，但是mega太慢，DNAman太丑。

当时，小Y还有30分钟进行工作汇报，眼瞅着来不及，在友人小湖的指点下，第一次在R语言上进行了进化树的绘制。随后下载nwk文件，在itol网站上进行美化。

一个小白的分享哈

流程：

1、使用plink进行IBS矩阵的构建；

2、基于IBS矩阵，在R语言构建nj进化树，并写出nwk文件；

3、将nwk文件在itol网站上进行美化。

1、使用plink进行IBS矩阵的构建；

plink --bfile 45_zk -recode -out 45_zk ## 转换二进制文件45_zk 为文本格式

plink --file 45_zk --cluster --matrix --noweb ##构建IBS矩阵

基于IBS矩阵在R语言中构建NJ进化树写出nwk文件

会出现如上图所示的文件。接下来，我们会使用plink.mibs 和 plink.mibs.id

2、基于IBS矩阵，在R语言构建nj进化树，并写出nwk文件；

#包的安装与读取  

install.packages("age")
install.packages("phangorn")
install.packages("seqinr")

library(ape)
library(phangorn)
library(seqinr)


##读入plink.mibs文件
dir()
m<-as.matrix(read.table("plink.mibs"))

基于IBS矩阵在R语言中构建NJ进化树写出nwk文件

这是小Y的一个数据，可以看到读入IBS矩阵后，共有个体366个


m<-as.matrix(read.table("plink.mibs"))
dimnames(m) <- list(1:366, 1:366) #366为个体数


#构建单位矩阵g
g=matrix(1,nrow=366,ncol=366,byrow=FALSE)#记得更新个体数

#计算遗传距离 D=1-IBS，
D=g-m

此时有的一些数据情况如下图

基于IBS矩阵在R语言中构建NJ进化树写出nwk文件

ID_number <- read.table("plink.mibs.id")
rownames(D) <- ID_number$V1  #补充每分支的ID

tr2 <- bionj(D) #使用bionj函数进行建树

write.tree(tr,file="tr.nwk")  #写出nwk文件，可在itol中美化树

如果直接在R 中出图，可以用以下命令

plot(tr2,cex =0.5) #输出系统发育树

plot(tr2,type="unrooted") #输出系统发育树

但是小Y半瓶子咣当当，不会调参数，就直接在itol 中进行了美化