导读
Prodigal是原核生物基因预测软件，常被用于原核生物de novo组装分析中。Prodigal能预测由de novo组装得到的新物种的基因组草图（bin）或contigs或scaffold中有哪些基因序列，并同时将这些序列翻译成蛋白。Phylophlan能通过将由基因组草图预测和翻译得到的基因蛋白序列数据与Phylophlan自带的3000+种微生物的蛋白数据做比对分析新物种与3000+微生物的进化关系。利用ggtree可视化结果可观察新物种在已知物种在进化树中的位置。我这里采用的输入数据是宏基因组分箱得到的基因组草图，获取方法请在：宏基因组分箱（二）Metabat2分箱实战。接下来是将新物种插入进化树的具体方法：（1）Prodigal蛋白预测；（2）Phyloplan进化分析；（3）ggtree可视化1、2、3、4。
1. Prodigal蛋白预测
for I in bin/all/*; do
BASE=${I##*/}
SAMPLE=${BASE%.*}
prodigal -a bin_function/all/$SAMPLE.faa -d bin_function/all/$SAMPLE.fna -f gff -g 11 -o bin_function/all/$SAMPLE.out -p single -s bin_function/all/$SAMPLE.pot -i $I &
done
# 将bin/all文件夹里的完成度大于90%的bin（基因组草图）进行基因预测，获得蛋白序列信息。
# 其实我这里只有两个bin：all.1.fasta和all.5.fasta
# prodigal部分参数：
-a: 输出的蛋白序列文件名
-d: 输出的核酸序列文件名
-f: 输出文件格式(gbk, gff, or sco)，默认gbk
-g: 指定密码子，原核为第11套，默认是11
-o: 输出文件，默认为屏幕输出，标准输出
-p: 选择方式，是单菌还是meta样品
-s: 将所有带有得分的潜在基因写入到指定文件中
-i: 指定FASTA/Genbank输入文件(默认标准输入)
二、phylophlan进化分析
将包含基因预测得到的faa(FASTA Amino Acid file)文件的整个“all”文件夹拷贝到phylophlan软件的input文件夹中，然后进行系统发生分析。这一步会将3000+种已知微生物与我的两个基因组草图放在一起进行进化分析，因此速度非常慢。
ll input/all/
# 查看输入文件，如下：
总用量 1608
drwxrwxr-x 2 cheng WST 4096 9月 30 18:01 ./
drwxrwxr-x 8 cheng WST 4096 10月 24 14:26 ../
-rw-rw-r-- 1 cheng WST 885979 9月 30 18:01 all.1.faa
-rw-rw-r-- 1 cheng WST 748679 9月 30 18:01 all.5.faa
phylophlan.py -i -t all --nproc 16
# 进化分析
ll output/all/
# 查看结果文件，如下：
总用量 7552
drwxrwxr-x 2 cheng WST 4096 10月 17 10:14 ./
drwxrwxr-x 9 cheng WST 4096 10月 24 14:28 ../
-rw-rw-r-- 1 cheng WST 107348 10月 8 11:26 all.tree.int.nwk
-rw-rw-r-- 1 cheng WST 219 10月 8 11:30 imputed_conf_high-conf.txt
三、配色和标签
利用nsegata-phylophlan/data/ppafull.tax.txt制作个性化画图所需的配色和标签文件。ppafull.tax.txt内含Phylophlan自带的3000+物种的分类注释信息。该文件内容如下：
less -S ppafull.tax.txt
# 查看文件：
t643692001 d__Bacteria.p__Acidobacteria.c__Acidobacteria.o__Acidobacteriales.f__Acidobacteriaceae.g__Acidobacterium.s__capsulatum.t__ATCC
t648276601 d__Bacteria.p__Acidobacteria.c__Acidobacteria.o__Acidobacteriales.f__Acidobacteriaceae.g__Granulicella.s__sp_.t__MP5ACTX8
t649633002 d__Bacteria.p__Acidobacteria.c__Acidobacteria.o__Acidobacteriales.f__Acidobacteriaceae.g__Granulicella.s__sp_.t__MP5ACTX9
t649633100 d__Bacteria.p__Acidobacteria.c__Acidobacteria.o__Acidobacteriales.f__Acidobacteriaceae.g__Terriglobus.s__saanensis.t__SP1PR4
t637000001 d__Bacteria.p__Acidobacteria.c__Acidobacteria.o__Acidobacteriales.f__Korebacteraceae.g__Korebacter.s__versatilis.t__Ellin345
t639633060 d__Bacteria.p__Acidobacteria.c__Solibacteres.o__Solibacterales.f__Solibacteraceae.g__Solibacter.s__usitatus.t__Ellin6076
...
提取ppafull.tax.txt的第一列全部内容和第二列的门注释信息，接着添加Bin信息，可得如下文件：
id Phylum
all.1 Bin
all.5 Bin
t643692001 Acidobacteria
t648276601 Acidobacteria
t649633002 Acidobacteria
t649633100 Acidobacteria
t637000001 Acidobacteria
t639633060 Acidobacteria
t644736322 Actinobacteria
t639633001 Actinobacteria
t643886017 Actinobacteria
...
四、ggtree画树状图：长方型、加对齐线
关键参数：
%<+% map # 引入注释文件
layout="rectangular" # 长方型树状图
align=TRUE # 添加对齐虚线
col=Phylum # 以Phylum给虚线着色
legend.position="bottom" # legend至于底部
legend.box="horizontal" # legend水平放置
library(ggplot2)
library(ggtree)
tree=read.tree("all.tree.int.nwk")
data=fortify(tree)
map=read.table("mapping.txt", header=T, sep="\t")
# 长方形（线型）
gra=ggtree(tree, layout="rectangular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=NA, col=Phylum), hjust=-0.5, align=TRUE, linesize=0.1)+
# 注释、颜色、高度、对其、虚点大小
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)
pdf("tree_rectangular_line.pdf")
gra
dev.off()
打开绘图结果tree_rectangular_line.pdf，如下：
图片.png
五、ggtree画树状图：长方型、末端着色
关键参数：
ggtree(aes(col=Phylum)) # 末端“枝”颜色
geom_tippoint(aes(color=Phylum)) # 末端“点”颜色
gra=ggtree(tree, aes(col=Phylum), layout="rectangular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tippoint(aes(color=Phylum), size=0.1)+
# 端点颜色、大小
geom_tiplab(aes(label=NA), size=0.8)+
# 注释
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)
pdf("tree_rectangular_point.pdf")
gra
dev.off()
打开绘图结果tree_rectangular_point.pdf，如下：
图片.png
六、ggtree画树状图：圆型、加对齐线
关键参数：
layout="circular" # 画圆型树状图
# 圆形（线型）
gra=ggtree(tree, layout="circular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=NA, col=Phylum), hjust=2, align=TRUE, linesize=0.1)+
# 注释、颜色、高度、对其、虚点大小
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)
pdf("tree_circular_line.pdf")
gra
dev.off()
打开绘图结果tree_circular_line.pdf，如下：
图片.png
七、ggtree画树状图：圆型、末端着色
关键参数;
ggtree(aes(col=Phylum)) # 末端“枝”颜色
geom_tippoint(aes(color=Phylum)) # 末端“点”颜色
gra=ggtree(tree, aes(col=Phylum), layout="circular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tippoint(aes(color=Phylum), size=0.1)+
# 端点颜色、大小
geom_tiplab(aes(label=NA, col=NA), size=0.8)+
# 注释、注释的颜色
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)
pdf("tree_circular_point.pdf")
gra
dev.off()
打开绘图结果tree_circular_point.pdf，如下：
图片.png
八、ggtree画树状图：长方型、加对齐线、加注释
1 准备注释文件
# 制作注释文件
cd /home/cheng/huty/softwares/phylophlan/data
cat ppafull.tax.txt | sed 's/.p__/\t/g' | sed 's/.c__/\t/g' | sed 's/.g__/\t/g' | awk 'BEGIN{OFS="\t"}{print $1, $3, $5}' > tax.txt
# 绘图准备
touch tax_bin.txt
echo -e 'id\tPhylum\tSpecies' > tax_bin.txt
cat bin_taxonomy.txt | sed '1d' | sed 's/.p__/\t/' | sed 's/.c__/\t/' | awk '{printf("%s\t%s\tBin\n", $1, $3)}' >> tax_bin.txt
cat /home/cheng/huty/softwares/phylophlan/data/tax.txt >> tax_bin.txt
2 绘图
# 绘制insert tree
library(ggplot2)
library(ggtree)
tree=read.tree("wm.insert.tree.int.nwk")
data=fortify(tree)
map=read.table("tax_bin.txt", header=T, sep="\t")
gra=ggtree(tree, layout="rectangular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=Species, col=Phylum), align=TRUE, size=0.18, linesize=0.05)+
# 注释、颜色、高度、对其、虚点大小
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5))) +
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)
ggsave(gra, filename="tree_rectangular.pdf", height=48, width=10)
九、ggtree + gheatmap组合图
library(ggplot2) # 加载ggplot2
library(ggtree) # 加载ggtree
tree=read.tree(list.files()[grepl("denovo.tree.nwk", list.files())]) # 读取nwk文件
map=read.table("map.txt", header=T, sep="\t", comment.char="")
data=fortify(tree)
abun = read.table("/home/cheng/client2/wangmeng/X101SC19061144-Z01-J041/Bin_all/Bin_quant/bin_abundance_table.txt", header=T, sep="\t", row.names=1)
abun = abun[, order(colnames(abun), decreasing=F)]
p = ggtree(tree, layout="rectangular", ladderize=FALSE, branch.length="none", size=0.8, aes(color = Phylum)) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(col=Phylum), align=TRUE, size=2.5) +
# 注释、颜色、高度、对其、虚点大小
theme(legend.position = 'none')
ggsave(p, file="bin_pheatmap.pdf", width=8, height=8)
p2 = gheatmap(p, abun, offset = 0.6, width = 0.5,
font.size=3, low = "green", high = "red", color = "white",
colnames_level = colnames(abun),
colnames_angle=90, hjust=.90)
ggsave(p2, file="bin_pheatmap2.pdf", width=14)
参考：
使用Y叔神包ggtree进行基因家族基因进化树构建
用ggtree来绘制进化树
在R中绘制进化树：ggtree学习笔记
