GraPhlAn繪製的超高顏值物種樹Cladogram

GraphLan繪製教程

我們經常在文章中看到這樣的圖

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進化樹的繪製方法，如下文：

• GraPhlAn：最美進化樹或層級分類樹

今天就帶大家根據特徵表(OTU table)、和物種註釋(Taxonomy)，繪製另一類高顏值的物種樹(Cladogram，也稱進化分支圖)。並提供相關測試數據、代碼，讓你準備好輸入文件，方便一步步生成繪圖所需文件。並可按需求組合數據和樣式，達到出版要求的圖片。

代碼和數據下載鏈接：

https://github.com/YongxinLiu/Note/tree/master/R/format2graphlan

format2graphlan.Rmd # 完整代碼文件，包括R和Bash兩種語言，需要在Linux中運行

format2graphlan.html # 代碼完整運行的報告，方便閱讀，也確保代碼有效和可重複

如果鏈接失效，“宏基因組”公眾號後臺回覆“graphlan”關鍵字獲取最新數據和代碼下載鏈接。

輸入文件

文件夾內要準備至少兩個文件：OTU表和物種註釋

<code># 從現在項目中複製文件，準備起始數據
cd ~
/github/
Note/R/format2graphlan
cp ~
/ehbio/
amplicon/
22
Pipeline/result/otutab.txt ./
cp ~
/ehbio/
amplicon/
22
Pipeline/result/taxonomy.txt .//<code>

OTU表<code>otutab.txt/<code>格式如下：行名為特徵OTU/ASV，列名為樣本名，可以為原始值或標準化的小數均可。

<code>



/<code>

物種註釋<code>taxonomy.txt/<code>：包括OTUID和7級註釋，末知的補Unassigned

<code>
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/<code>

首選我們要對原始數據進行篩選，因為結果過少或過多都不美觀。如根據丰度進行篩選Top 150的特徵進行展示。

1. 特徵表求均值並按丰度篩選

輸入文件：OTU表+物種註釋

可以指定丰度或數量篩選，兩個條件選擇共有部分

輸出文件：OTU對應均值，篩選後的OTU表+物種註釋

<code># 參數設置
# 按丰度篩選，如
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))/<code>

2. 繪製樹骨架文件

輸入文件為篩選後的taxonomy文件：filtered_taxonomy.txt

本處主要篩選了門、綱、目、科、屬和OTU作為樹枝，按科添加標籤，並對應門著色。由於Unassigned末分類的較多，重名會引著色混亂(每個標籤是獨立著色的，名稱必須唯一，不唯一時後出現的名稱會覆蓋之前的顏色值。)，本文去除了在科水平無註釋的分類單元。

<code># 讀取篩選後的文件，不設置行名
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=
""
)/<code>

此時生成了兩個文件

樹骨架

tree1_backbone.txt

是一點相連的各級物種分類名稱，添加p, c等為減少不同級別的不規範重名引起顏色混亂

<code>
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
/<code>

樹顏色和標籤

tree2_label_color.txt

科水平的標籤、標籤旋轉角度和與門對應的顏色。

<code>
Thermomonosporaceae
 
annotation
 *
 
Thermomonosporaceae
 
annotation_rotation
 90
Thermomonosporaceae
 
annotation_background_color
 
#F58D8D

Comamonadaceae
 
annotation
 *
Comamonadaceae
 
annotation_rotation
 90
Comamonadaceae
 
annotation_background_color
 
#85F29B
/<code>

基本樹繪圖

繪製樹，還需要一些參數文件，見cfg目錄，是我提前編寫好的樣本，可以調整更多樣式。

<code>cfg/global.cfg/<code>設置了圖型的基本樣式，配色等，

以下部分以bash中操作，需要在Linux上的Rstudio或Rstudio server中操作。或自己使用終端連接服務器執行。

一定要提前安裝過graphlan這個軟件，安裝方法<code>conda install graphlan/<code>

<code>
rm
 -rf track*







/<code>

現在用以上代碼為大家寫出了一套註釋方案，這要是手動編寫和優化出這方案，也可能要花上幾天至幾周。

我們需要從樹文件中獲得節點名稱，並添加註釋數據。

如獲得結點的丰度，在下面很多註釋都會基於丰度信息

<code>
 獲得最終出圖的結點ID

cut -f 6 -d 
'.'
 tree1_backbone.txt > tree1_backbone.id


/<code>

形狀標籤有無

樣式1. 如篩選丰度，用紫色方塊標出大於千分之5的結點

<code>
 環1，篩選千分之五的結果註釋為方塊，cfg/ring1.cfg中的m代表紫色，R代表方塊

cat cfg/ring1.cfg '$2>0.5'
 tree1_backbone.mean | cut -f 1 | sed 
's/$/\tring_shape\t1\tR/'
) > track1

 



/<code>

樣式2. 如篩選丰度，用第二環位置橙色倒三角標出小於千分之5的結點

註釋：ring2.cfg為第二環，顏色y為yellow橙色，註釋track中也為2

<code>
 環1篩選千分之五的結果註釋為方塊，cfg/ring1.cfg中的m代表紫色，R代表方塊 

cat cfg/ring2.cfg '$2<=0.5'
 tree1_backbone.mean | cut -f 1 | sed 
's/$/\tring_shape\t2\tv/'
) > track2




/<code>

熱圖展示丰度

添加所有樣品均值作為熱圖，作為第3環。

本質上熱圖即環形條帶的透明度

<code>
 環3用綠色不同的透明度展示丰度

cat cfg/heat3.cfg 's/\t/\tring_alpha\t3\t/g'
 tree1_backbone.mean) > track3




/<code>

我們可以用同樣原理，添加每個組，或每個樣品的丰度熱圖。

柱狀圖顯示丰度

<code>
 環4用藍色柱狀圖展示丰度

cat cfg/bar4.cfg 's/\t/\tring_height\t4\t/g'
 tree1_backbone.mean) > track4




/<code>

附錄1. 顏色

顏色有三種設置方法

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_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

Reference

http://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

生物科普: 腸道細菌 人體上的生命 生命大躍進 細胞暗戰 人體奧秘

寫在後面

學習16S擴增子、宏基因組科研思路和分析實戰，關注“宏基因組”