導讀
Prodigal是原核生物基因預測軟件侄旬，常被用于原核生物de novo組裝分析中鸽捻。Prodigal能預測由de novo組裝得到的新物種的基因組草圖（bin）或contigs或scaffold中有哪些基因序列糊余，并同時將這些序列翻譯成蛋白。Phylophlan能通過將由基因組草圖預測和翻譯得到的基因蛋白序列數(shù)據(jù)與Phylophlan自帶的3000+種微生物的蛋白數(shù)據(jù)做比對分析新物種與3000+微生物的進化關系。利用ggtree可視化結果可觀察新物種在已知物種在進化樹中的位置豫缨。我這里采用的輸入數(shù)據(jù)是宏基因組分箱得到的基因組草圖，獲取方法請在：宏基因組分箱（二）Metabat2分箱實戰(zhàn)端朵。接下來是將新物種插入進化樹的具體方法：（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部分參數(shù)：
-a: 輸出的蛋白序列文件名
-d: 輸出的核酸序列文件名
-f: 輸出文件格式(gbk, gff, or sco)乘凸，默認gbk
-g: 指定密碼子厕诡，原核為第11套，默認是11
-o: 輸出文件营勤，默認為屏幕輸出灵嫌，標準輸出
-p: 選擇方式信柿，是單菌還是meta樣品
-s: 將所有帶有得分的潛在基因?qū)懭氲街付ㄎ募? -i: 指定FASTA/Genbank輸入文件(默認標準輸入)
二、phylophlan進化分析
將包含基因預測得到的faa(FASTA Amino Acid file)文件的整個“all”文件夾拷貝到phylophlan軟件的input文件夾中醒第，然后進行系統(tǒng)發(fā)生分析渔嚷。這一步會將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內(nèi)含Phylophlan自帶的3000+物種的分類注釋信息。該文件內(nèi)容如下：
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的第一列全部內(nèi)容和第二列的門注釋信息司恳，接著添加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畫樹狀圖：長方型扔傅、加對齊線
關鍵參數(shù)：
%<+% 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畫樹狀圖：長方型倘屹、末端著色
關鍵參數(shù)：
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畫樹狀圖：圓型、加對齊線
關鍵參數(shù)：
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畫樹狀圖：圓型、末端著色
關鍵參數(shù);
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學習筆記
