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

yy160417

已于 2023-01-26 00:09:50 修改

阅读量2.2k

点赞数 3

分类专栏：亲缘关系文章标签： r语言

于 2023-01-26 00:02:27 首次发布

本文链接：https://blog.csdn.net/2301_76218064/article/details/128763151

版权

亲缘关系专栏收录该内容

1 篇文章 0 订阅

订阅专栏

构建系统发育树有很多方法，但是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矩阵

会出现如上图所示的文件。接下来，我们会使用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"))

这是小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

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

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 中进行了美化