使用iqtree软件利用基因存在缺失变异矩阵（0/1）矩阵构建进化树

线性泛基因组相关论文通常会获得基因存在缺失变异矩阵，接下来会使用这个矩阵构建进化树，今天的推文介绍一下使用iqtree软件利用基因存在缺失变异矩阵（0/1）矩阵构建进化树的代码

iqtree软件可以直接使用conda进行安装

如果是0/1矩阵作为输入数据，iqtree需要用到的格式是phy这个格式

http://www.iqtree.org/doc/Tutorial

image.png

我们那到的基因存在缺失变异矩阵通常的格式是 行是基因，列是样本的一个表格

image.png

这里我们用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")

这个进化树怎么美化再单独出推文进行介绍吧