Phylophlan（三）將新物種插入進(jìn)化樹(shù)

導(dǎo)讀

Prodigal是原核生物基因預(yù)測(cè)軟件，常被用于原核生物de novo組裝分析中。Prodigal能預(yù)測(cè)由de novo組裝得到的新物種的基因組草圖（bin）或contigs或scaffold中有哪些基因序列，并同時(shí)將這些序列翻譯成蛋白。Phylophlan能通過(guò)將由基因組草圖預(yù)測(cè)和翻譯得到的基因蛋白序列數(shù)據(jù)與Phylophlan自帶的3000+種微生物的蛋白數(shù)據(jù)做比對(duì)分析新物種與3000+微生物的進(jìn)化關(guān)系。利用ggtree可視化結(jié)果可觀察新物種在已知物種在進(jìn)化樹(shù)中的位置。我這里采用的輸入數(shù)據(jù)是宏基因組分箱得到的基因組草圖，獲取方法請(qǐng)?jiān)冢?a href="" target="_blank">宏基因組分箱（二）Metabat2分箱實(shí)戰(zhàn)。接下來(lái)是將新物種插入進(jìn)化樹(shù)的具體方法：（1）Prodigal蛋白預(yù)測(cè)；（2）Phyloplan進(jìn)化分析；（3）ggtree可視化1、2、3、4。

1. Prodigal蛋白預(yù)測(cè)

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（基因組草圖）進(jìn)行基因預(yù)測(cè)，獲得蛋白序列信息。
# 其實(shí)我這里只有兩個(gè)bin：all.1.fasta和all.5.fasta

    # prodigal部分參數(shù)：
    -a: 輸出的蛋白序列文件名
    -d: 輸出的核酸序列文件名
    -f: 輸出文件格式(gbk, gff, or sco)，默認(rèn)gbk
    -g: 指定密碼子，原核為第11套，默認(rèn)是11
    -o: 輸出文件，默認(rèn)為屏幕輸出，標(biāo)準(zhǔn)輸出
    -p: 選擇方式，是單菌還是meta樣品
    -s: 將所有帶有得分的潛在基因?qū)懭氲街付ㄎ募?    -i: 指定FASTA/Genbank輸入文件(默認(rèn)標(biāo)準(zhǔn)輸入)

二、phylophlan進(jìn)化分析

將包含基因預(yù)測(cè)得到的faa(FASTA Amino Acid file)文件的整個(gè)“all”文件夾拷貝到phylophlan軟件的input文件夾中，然后進(jìn)行系統(tǒng)發(fā)生分析。這一步會(huì)將3000+種已知微生物與我的兩個(gè)基因組草圖放在一起進(jìn)行進(jìn)化分析，因此速度非常慢。

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
# 進(jìn)化分析

ll output/all/
# 查看結(jié)果文件，如下：

    總用量 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

三、配色和標(biāo)簽

利用nsegata-phylophlan/data/ppafull.tax.txt制作個(gè)性化畫(huà)圖所需的配色和標(biāo)簽文件。ppafull.tax.txt內(nèi)含Phylophlan自帶的3000+物種的分類(lèi)注釋信息。該文件內(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的第一列全部?jī)?nèi)容和第二列的門(mén)注釋信息，接著添加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畫(huà)樹(shù)狀圖：長(zhǎng)方型、加對(duì)齊線(xiàn)

關(guān)鍵參數(shù)：
%<+% map # 引入注釋文件
layout="rectangular" # 長(zhǎng)方型樹(shù)狀圖
align=TRUE # 添加對(duì)齊虛線(xiàn)
col=Phylum # 以Phylum給虛線(xiàn)著色
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")

# 長(zhǎng)方形（線(xiàn)型）
gra=ggtree(tree, layout="rectangular", size=0.1) %<+% map +
# 樹(shù)型、線(xiàn)粗細(xì)、末端顏色 + 注釋信息
geom_tiplab(aes(label=NA, col=Phylum), hjust=-0.5, align=TRUE, linesize=0.1)+
# 注釋、顏色、高度、對(duì)其、虛點(diǎn)大小
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()

打開(kāi)繪圖結(jié)果tree_rectangular_line.pdf，如下：

圖片.png

五、ggtree畫(huà)樹(shù)狀圖：長(zhǎng)方型、末端著色

關(guān)鍵參數(shù)：
ggtree(aes(col=Phylum)) # 末端“枝”顏色
geom_tippoint(aes(color=Phylum)) # 末端“點(diǎn)”顏色

gra=ggtree(tree, aes(col=Phylum), layout="rectangular", size=0.1) %<+% map +
# 樹(shù)型、線(xiàn)粗細(xì)、末端顏色 + 注釋信息
geom_tippoint(aes(color=Phylum), size=0.1)+
# 端點(diǎn)顏色、大小
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()

打開(kāi)繪圖結(jié)果tree_rectangular_point.pdf，如下：

圖片.png

六、ggtree畫(huà)樹(shù)狀圖：圓型、加對(duì)齊線(xiàn)

關(guān)鍵參數(shù)：
layout="circular" # 畫(huà)圓型樹(shù)狀圖

# 圓形（線(xiàn)型）
gra=ggtree(tree, layout="circular", size=0.1) %<+% map +
# 樹(shù)型、線(xiàn)粗細(xì)、末端顏色 + 注釋信息
geom_tiplab(aes(label=NA, col=Phylum), hjust=2, align=TRUE, linesize=0.1)+
# 注釋、顏色、高度、對(duì)其、虛點(diǎn)大小
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()

打開(kāi)繪圖結(jié)果tree_circular_line.pdf，如下：

圖片.png

七、ggtree畫(huà)樹(shù)狀圖：圓型、末端著色

關(guān)鍵參數(shù);
ggtree(aes(col=Phylum)) # 末端“枝”顏色
geom_tippoint(aes(color=Phylum)) # 末端“點(diǎn)”顏色

gra=ggtree(tree, aes(col=Phylum), layout="circular", size=0.1) %<+% map +
# 樹(shù)型、線(xiàn)粗細(xì)、末端顏色 + 注釋信息
geom_tippoint(aes(color=Phylum), size=0.1)+
# 端點(diǎn)顏色、大小
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()

打開(kāi)繪圖結(jié)果tree_circular_point.pdf，如下：

圖片.png

八、ggtree畫(huà)樹(shù)狀圖：長(zhǎng)方型、加對(duì)齊線(xiàn)、加注釋

1 準(zhǔn)備注釋文件

# 制作注釋文件
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

# 繪圖準(zhǔn)備
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 +
# 樹(shù)型、線(xiàn)粗細(xì)、末端顏色 + 注釋信息
geom_tiplab(aes(label=Species, col=Phylum), align=TRUE, size=0.18, linesize=0.05)+
# 注釋、顏色、高度、對(duì)其、虛點(diǎn)大小
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 +
# 樹(shù)型、線(xiàn)粗細(xì)、末端顏色 + 注釋信息
geom_tiplab(aes(col=Phylum), align=TRUE, size=2.5) +
# 注釋、顏色、高度、對(duì)其、虛點(diǎn)大小
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進(jìn)行基因家族基因進(jìn)化樹(shù)構(gòu)建
用ggtree來(lái)繪制進(jìn)化樹(shù)
在R中繪制進(jìn)化樹(shù)：ggtree學(xué)習(xí)筆記

\color{green}{????原創(chuàng)文章，碼字不易，轉(zhuǎn)載請(qǐng)注明出處????}

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容