我們經常在文章中看到這樣的圖
Yang Bai, Daniel B. Müller, Girish Srinivas, Ruben Garrido-Oter, Eva Potthoff, Matthias Rott, Nina Dombrowski, Philipp C. Münch, Stijn Spaepen, Mitja Remus-Emsermann, Bruno Hüttel, Alice C. McHardy, Julia A. Vorholt & Paul Schulze-Lefert. Functional overlap of the Arabidopsis leaf and root microbiota. Nature. 2015, 528: 364-369. doi:10.1038/nature16192
還有這樣的圖
Jingying Zhang, Yong-Xin Liu, Na Zhang, Bin Hu, Tao Jin, Haoran Xu, Yuan Qin, Pengxu Yan, Xiaoning Zhang, Xiaoxuan Guo, Jing Hui, Shouyun Cao, Xin Wang, Chao Wang, Hui Wang, Baoyuan Qu, Guangyi Fan, Lixing Yuan, Ruben Garrido-Oter, Chengcai Chu & Yang Bai. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology. 2019, 37: 676-684. doi:10.1038/s41587-019-0104-4
是不是很漂亮
之前公眾號已經為大家介紹了GraPhlAn進化樹的繪製方法，如下文：
今天就帶大家根據特徵表(OTU table)、和物種注釋(Taxonomy)，繪製另一類高顏值的物種樹(Cladogram，也稱進化分支圖)。並提供相關測試數據、代碼，讓你準備好輸入文件，方便一步步生成繪圖所需文件。並可按需求組合數據和樣式，達到出版要求的圖片。
代碼和數據下載連結：
https://github.com/YongxinLiu/Note/tree/master/R/format2graphlan
format2graphlan.Rmd # 完整代碼文件，包括R和Bash兩種語言，需要在Linux中運行
format2graphlan.html # 代碼完整運行的報告，方便閱讀，也確保代碼有效和可重複
如果連結失效，「宏基因組」公眾號後臺回復「graphlan」關鍵字獲取最新數據和代碼下載連結。
輸入文件文件夾內要準備至少兩個文件：OTU表和物種注釋
# 從現在項目中複製文件，準備起始數據
cd ~/github/Note/R/format2graphlan
cp ~/ehbio/amplicon/22Pipeline/result/otutab.txt ./
cp ~/ehbio/amplicon/22Pipeline/result/taxonomy.txt ./
OTU表otutab.txt格式如下：行名為特徵OTU/ASV，列名為樣本名，可以為原始值或標準化的小數均可。
#OTUID KO1 KO2 KO3
ASV_1 1113 1968 816
ASV_2 1922 1227 2355
ASV_3 568 460 899
物種注釋taxonomy.txt：包括OTUID和7級注釋，末知的補Unassigned
OTUID Kingdom Phylum Class Order Family Genus Species
ASV_1 Bacteria Actinobacteria Actinobacteria Actinomycetales Thermomonosporaceae Unassigned Unassigned
ASV_2 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Pelomonas Pelomonas_puraquae
ASV_3 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Rhizobacter Rhizobacter_bergeniae
首選我們要對原始數據進行篩選，因為結果過少或過多都不美觀。如根據豐度進行篩選Top 150的特徵進行展示。
1. 特徵表求均值並按豐度篩選輸入文件：OTU表+物種注釋
可以指定豐度或數量篩選，兩個條件選擇共有部分
輸出文件：OTU對應均值，篩選後的OTU表+物種注釋
# 參數設置
# 按豐度篩選，如0.01即代表0.01%，即萬分之一
abundance = 0.01
# 按數量篩選，如150即代表最高豐度的150個特徵
number = 150
# 讀取輸入文件
otutab = read.table("otutab.txt", sep="\t", header = TRUE, row.names = 1, stringsAsFactors = F, comment.char = "")
taxonomy = read.table("taxonomy.txt", sep="\t", header = TRUE, row.names = 1, stringsAsFactors = F, comment.char = "")
# 數據篩選
# 標準化並求均值
norm = as.data.frame(t(t(otutab)/colSums(otutab,na=T)*100))
# 豐度由大到小排序
idx = order(rowMeans(norm), decreasing = T)
norm = norm[idx,]
# 按豐度篩選
idx = rowMeans(norm) > abundance
filtered_otutab = norm[idx,]
# 按數量篩選
filtered_otutab = head(norm, number)
# 添加均值並保留4位小數
filtered_otutab = round(cbind(rowMeans(filtered_otutab), filtered_otutab), digits = 4)
colnames(filtered_otutab)[1] = "Mean"
# 對應過濾物種注釋
idx = rownames(filtered_otutab) %in% rownames(taxonomy)
filtered_otutab = filtered_otutab[idx,]
filtered_taxonomy = taxonomy[rownames(filtered_otutab),]
# 保存輸出文件
# 過濾的OTU表
write.table("OTUID\t", file="filtered_otutab.txt", append = F, sep="\t", quote=F, eol = "", row.names=F, col.names=F)
suppressWarnings(write.table(filtered_otutab, file="filtered_otutab.txt", append = T, sep="\t", quote=F, row.names=T, col.names=T))
# 過濾的物種注釋
write.table("OTUID\t", file="filtered_taxonomy.txt", append = F, sep="\t", quote=F, eol = "", row.names=F, col.names=F)
suppressWarnings(write.table(filtered_taxonomy, file="filtered_taxonomy.txt", append = T, sep="\t", quote=F, row.names=T, col.names=T))
輸入文件為篩選後的taxonomy文件：filtered_taxonomy.txt
本處主要篩選了門、綱、目、科、屬和OTU作為樹枝，按科添加標籤，並對應門著色。由於Unassigned末分類的較多，重名會引著色混亂(每個標籤是獨立著色的，名稱必須唯一，不唯一時後出現的名稱會覆蓋之前的顏色值。)，本文去除了在科水平無注釋的分類單元。
# 讀取篩選後的文件，不設置行名
tax = read.table("filtered_taxonomy.txt", sep="\t", header = TRUE, stringsAsFactors = F)
# 篩選門-屬5級+OTUID
tree = data.frame(tax[,c(3:7,1)], stringsAsFactors = F)
# head(tree)
## clarify taxonomy，解決不同級別重名問題，為可識別級別，且與Greengene格式保持一致
tree[,1] = paste("p__",tree[,1],sep = "")
tree[,2] = paste("c__",tree[,2],sep = "")
tree[,3] = paste("o__",tree[,3],sep = "")
# tree[,4] = paste("f__",tree[,4],sep = "")
tree[,5] = paste("g__",tree[,5],sep = "")
# save tree backbone, 按點分隔格式
# 解決科標籤重名問題
idx = tree[,4] %in% "Unassigned"
# 方法1. 重名標籤添加數字編號，但結果有太多Unassigned
# tree[idx,4] = paste0(tree[idx,4], 1:length(tree[idx,4]))
# 方法2. 過濾掉科末注釋的條目，數量會減少，但圖片更美觀
tree = tree[!idx,]
# 簡化一些代_的不規則科名
tree[,4] = gsub('_\\w*',"",tree[,4])
write.table (tree, file="tree1_backbone.txt", sep=".", col.names=F, row.names=F, quote=F)
# 列出現在有門、綱、目、科、屬，用於設置與門對應的背景色
Phylum = unique(tree[,1])
Class = unique(tree[,2])
Order = unique(tree[,3])
Family = unique(tree[,4])
Genus = unique(tree[,5])
# 篩選四大菌門中的科並按門著色
# 修改為目，則將tree的4列改為3列，Family改為Order
pro = tree[tree[,1]=="p__Proteobacteria",4]
act = tree[tree[,1]=="p__Actinobacteria",4]
bac = tree[tree[,1]=="p__Bacteroidetes",4]
fir = tree[tree[,1]=="p__Firmicutes",4]
# 對每個科進行標籤、文字旋轉、按門注釋背景色
# 也可調整為其它級別，如Order, Class或Genus
label_color = data.frame(stringsAsFactors = F)
for (element in Family)
{
# element
anno = data.frame(stringsAsFactors = F)
anno[1,1] = element
anno[1,2] = "annotation"
anno[1,3] = "*"
# 設置文字旋轉90度
anno[2,1] = element
anno[2,2] = "annotation_rotation"
anno[2,3] = "90"
# 設置背景色，四大門各指定一種色，其它為灰色
anno[3,1] = element
anno[3,2] = "annotation_background_color"
if (element %in% pro)
{
anno[3,3] = "#85F29B"
} else if (element %in% act)
{
anno[3,3] = "#F58D8D"
} else if (element %in% fir)
{
anno[3,3] = "#F7C875"
} else if (element %in% bac)
{
anno[3,3] = "#91DBF6"
} else {
anno[3,3] = "grey"
}
label_color = rbind(label_color,anno)
}
write.table(label_color, "tree2_label_color.txt", sep = "\t", quote = F,col.names = F,row.names = F, na="")
此時生成了兩個文件
樹骨架tree1_backbone.txt
是一點相連的各級物種分類名稱，添加p, c等為減少不同級別的不規範重名引起顏色混亂
p__Actinobacteria.c__Actinobacteria.o__Actinomycetales.Thermomonosporaceae.g__Unassigned.ASV_1
p__Proteobacteria.c__Betaproteobacteria.o__Burkholderiales.Comamonadaceae.g__Pelomonas.ASV_2
p__Proteobacteria.c__Gammaproteobacteria.o__Pseudomonadales.Pseudomonadaceae.g__Rhizobacter.ASV_3
tree2_label_color.txt
科水平的標籤、標籤旋轉角度和與門對應的顏色。
Thermomonosporaceae annotation *
Thermomonosporaceae annotation_rotation 90
Thermomonosporaceae annotation_background_color #F58D8D
Comamonadaceae annotation *
Comamonadaceae annotation_rotation 90
Comamonadaceae annotation_background_color #85F29B
繪製樹，還需要一些參數文件，見cfg目錄，是我提前編寫好的樣本，可以調整更多樣式。
cfg/global.cfg設置了圖型的基本樣式，配色等，
以下部分以bash中操作，需要在Linux上的Rstudio或Rstudio server中操作。或自己使用終端連接伺服器執行。
一定要提前安裝過graphlan這個軟體，安裝方法conda install graphlan
rm -rf track*
# 生成樹的默認參數，可手動調整更多樣式
cat cfg/global.cfg tree2_label_color.txt > track0
# 合併所有的注釋，接下來會生成更多track，使樹更複雜
cat track* > graphlan_annotate.txt
# 注釋樹
graphlan_annotate.py --annot graphlan_annotate.txt tree1_backbone.txt graphlan.xml
# 繪圖，size決定圖片大小，越大字越小
graphlan.py graphlan.xml graphlan0_tree.pdf --size 5
現在用以上代碼為大家寫出了一套注釋方案，這要是手動編寫和優化出這方案，也可能要花上幾天至幾周。
我們需要從樹文件中獲得節點名稱，並添加注釋數據。
如獲得結點的豐度，在下面很多注釋都會基於豐度信息
# 獲得最終出圖的結點ID
cut -f 6 -d '.' tree1_backbone.txt > tree1_backbone.id
# 注釋結果豐度均值
awk 'BEGIN{OFS=FS="\t"} NR==FNR{a[$1]=$2} NR>FNR {print $1,a[$1]}' filtered_otutab.txt tree1_backbone.id > tree1_backbone.mean
樣式1. 如篩選豐度，用紫色方塊標出大於千分之5的結點
# 環1，篩選千分之五的結果注釋為方塊，cfg/ring1.cfg中的m代表紫色，R代表方塊
cat cfg/ring1.cfg <(awk '$2>0.5' tree1_backbone.mean | cut -f 1 | sed 's/$/\tring_shape\t1\tR/') > track1
# 繪圖，加第一環矩形，展示豐度大於千萬的特徵
cat track* > graphlan_annotate.txt
graphlan_annotate.py --annot graphlan_annotate.txt tree1_backbone.txt graphlan.xml
graphlan.py graphlan.xml graphlan1_rectangle.pdf --size 5
樣式2. 如篩選豐度，用第二環位置橙色倒三角標出小於千分之5的結點
注釋：ring2.cfg為第二環，顏色y為yellow橙色，注釋track中也為2
# 環1篩選千分之五的結果注釋為方塊，cfg/ring1.cfg中的m代表紫色，R代表方塊
cat cfg/ring2.cfg <(awk '$2<=0.5' tree1_backbone.mean | cut -f 1 | sed 's/$/\tring_shape\t2\tv/') > track2
# 繪圖，加第一環矩形，展示豐度大於千萬的特徵
cat track* > graphlan_annotate.txt
graphlan_annotate.py --annot graphlan_annotate.txt tree1_backbone.txt graphlan.xml
graphlan.py graphlan.xml graphlan2_triangle.pdf --size 5
添加所有樣品均值作為熱圖，作為第3環。
本質上熱圖即環形條帶的透明度
# 環3用綠色不同的透明度展示豐度
cat cfg/heat3.cfg <(sed 's/\t/\tring_alpha\t3\t/g' tree1_backbone.mean) > track3
# 繪圖綠色不同的透明度的3號環
cat track* > graphlan_annotate.txt
graphlan_annotate.py --annot graphlan_annotate.txt tree1_backbone.txt graphlan.xml
graphlan.py graphlan.xml graphlan3_heatmap.pdf --size 5
我們可以用同樣原理，添加每個組，或每個樣品的豐度熱圖。
柱狀圖顯示豐度# 環4用藍色柱狀圖展示豐度
cat cfg/bar4.cfg <(sed 's/\t/\tring_height\t4\t/g' tree1_backbone.mean) > track4
# 繪圖，環4用藍色柱狀圖展示豐度
cat track* > graphlan_annotate.txt
graphlan_annotate.py --annot graphlan_annotate.txt tree1_backbone.txt graphlan.xml
graphlan.py graphlan.xml graphlan4_bar.pdf --size 5
顏色有三種設置方法
1. 顏色英文名稱
blue, green, red, cyan, magenta, yellow, black, white
2. 單個字母的縮寫
『b』 (blue), 『g』 (green), 『r』 (red), 『c』 (cyan), 『m』 (magenta), 『y』 (yellow), 『k』 (black), 『w』 (white)
3. RGB模式顏色
rrggbb, for example #FF0000 corresponds to (full) red附錄2. 形狀選擇『.』 : 點 point marker
『,』 : pixel marker
『o』 : 圈 circle marker
『v』 : 下三角 triangle_down marker
『^』 : triangle_up marker
『<』 : triangle_left marker
『>』 : triangle_right marker
『1』 : tri_down marker
『2』 : tri_up marker
『3』 : tri_left marker
『4』 : tri_right marker
『s』 : square marker
『R』 : 矩陣 rectangle marker
『p』 : pentagon marker
『*』 : star marker
『h』 : hexagon1 marker
『H』 : hexagon2 marker
『+』 : plus marker
『x』 : x marker
『D』 : diamond marker
『d』 : thin_diamond marker
『|』 : vline marker
『_』 : hline marker
Referencehttp://huttenhower.sph.harvard.edu/graphlan
Asnicar, Francesco, George Weingart, Timothy L. Tickle, Curtis Huttenhower, and Nicola Segata. 2015. 『Compact graphical representation of phylogenetic data and metadata with GraPhlAn』, PeerJ, 3: e1029.
Jingying Zhang, Yong-Xin Liu, Na Zhang, Bin Hu, Tao Jin, Haoran Xu, Yuan Qin, Pengxu Yan, Xiaoning Zhang, Xiaoxuan Guo, Jing Hui, Shouyun Cao, Xin Wang, Chao Wang, Hui Wang, Baoyuan Qu, Guangyi Fan, Lixing Yuan, Ruben Garrido-Oter, Chengcai Chu & Yang Bai. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology. 2019, 37: 676-684. doi:10.1038/s41587-019-0104-4
註：本文的代碼來自我之前發表的Nature Biotechnology中的圖5，如果參考本文代碼繪製類似圖，請引用以上兩篇文章，謝謝！
猜你喜歡10000+：菌群分析 寶寶與貓狗 梅毒狂想曲 提DNA發Nature Cell專刊 腸道指揮大腦
系列教程：微生物組入門 Biostar 微生物組 宏基因組
專業技能：學術圖表 高分文章 生信寶典 不可或缺的人
一文讀懂：宏基因組 寄生蟲益處 進化樹
必備技能：提問 搜索 Endnote
文獻閱讀 熱心腸 SemanticScholar Geenmedical
擴增子分析：圖表解讀 分析流程 統計繪圖
16S功能預測 PICRUSt FAPROTAX Bugbase Tax4Fun
在線工具：16S預測培養基 生信繪圖
科研經驗：雲筆記 雲協作 公眾號
編程模板: Shell R Perl
生物科普: 腸道細菌 人體上的生命 生命大躍進 細胞暗戰 人體奧秘
寫在後面為鼓勵讀者交流、快速解決科研困難，我們建立了「宏基因組」專業討論群，目前己有國內外5000+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職稱/年級」。PI請明示身份，另有海內外微生物相關PI群供大佬合作交流。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍未解決群內討論，問題不私聊，幫助同行。
學習16S擴增子、宏基因組科研思路和分析實戰，關注「宏基因組」
點擊閱讀原文，跳轉最新文章目錄閱讀