R语言读入比对好的fasta文件然后做NJ树并做boostrap检验

2021-01-20 16:04:37 浏览数 (1)

今天的推文内容主要参考

  • https://www.rpubs.com/michelleprem/683962
  • https://fuzzyatelin.github.io/bioanth-stats/module-24/module-24.html
首先是读入数据

今天推文用到的示例数据是参考链接2中提供的usflu.fasta,fasta文件已经比对好,R语言里读入fasta格式的数据可以使用adegenet包中的fasta2DNAbin函数

代码语言:javascript复制
#install.packages("adegenet")
library(adegenet)
dna<-fasta2DNAbin(file = "usflu.fasta")
dna
计算距离矩阵
代码语言:javascript复制
library(ape)
dd<-dist.dna(dna)

用到的是ape包中的dist.dna()函数

构建NJ树
代码语言:javascript复制
tree<-nj(dd)

用到的是ape包中的nj()函数

ggtree进行可视化
代码语言:javascript复制
library(ggtree)

ggtree(tree) 
  geom_tiplab(size=2)

image.png

做bootstrap检验
代码语言:javascript复制
bs.tree<-boot.phylo(tree,dna,
                    function(x)nj(dist.dna(x)),1000,
                    trees = TRUE)
将得到的bootstrap值合并到tree中
代码语言:javascript复制
tree$node.label<-bs.tree$BP

这一步不知道对不对,好像是有问题,暂时还不知道如何验证

结果里展示bootstrap值
代码语言:javascript复制
ggtree(tree) 
  geom_tiplab(size=2) 
  geom_nodelab(hjust=-0.2,size=2)

image.png

一下是之前关于ggtree的一些推送

0 人点赞