线性泛基因组相关论文通常会获得基因存在缺失变异矩阵,接下来会使用这个矩阵构建进化树,今天的推文介绍一下使用iqtree软件利用基因存在缺失变异矩阵(0/1)矩阵构建进化树的代码
iqtree软件可以直接使用conda进行安装
如果是0/1矩阵作为输入数据,iqtree需要用到的格式是phy这个格式
http://www.iqtree.org/doc/Tutorial
我们那到的基因存在缺失变异矩阵通常的格式是 行是基因,列是样本的一个表格
这里我们用R语言把这个表格转换成iqtree需要的phy格式输入文件
R语言代码
library(tidyverse)
read_tsv("2024.data/20240123/fig1_pangenome/pres_abs.tsv") %>%
select(-c(2,3)) %>%
column_to_rownames("Orthogroup") %>%
t() %>%
as.data.frame() %>%
unite("newcol",everything(),sep="") %>%
rownames_to_column() %>%
write_tsv("2024.data/20240123/fig1_pangenome/pra.phy",col_names = FALSE)
输出文件pra.phy需要手动修改,在第一行添加两个数字,第一个数字是多少个样本,第二个数字是多少个位点,中间用制表符分隔。
这里有一个小知识点,R语言里把数据框所有列合并成一列,可以用tidyr包中的unite函数。把一列拆分成很多列可以用separate函数,参考这个链接 https://tidyr.tidyverse.org/reference/unite.html
这里的示例数据集来源于论文 Aspergillus fumigatus pan-genome analysis identifies genetic variants associated with human infection,可以直接到论文中去下载
iqtree构建进化树的代码
iqtree2 -s pra.phy -T 24 -m GTR2+FO
这里为了加快运行速度,随便选择了一个模型,没有设置其他额外参数,如果是自己的真实数据,具体参数设置需要参考iqtree的文档。
ggtree展示进化树的代码
library(ggtree)
tree<-read.tree("2024.data/20240123/fig1_pangenome/pra.phy.treefile")
ggtree(tree,layout = "circular")
这个进化树怎么美化再单独出推文进行介绍吧
欢迎大家关注我的公众号
小明的数据分析笔记本
声明:文中观点不代表本站立场。本文传送门:https://eyangzhen.com/392284.html