【原】五彩進化樹與熱圖更配-ggtree美顏進化樹(宏基因組擴增子)

宏基因組 2020-10-09

展開全文

Y叔創(chuàng)建的R包——ggtree,，進化樹美化神器。

軟件原文G Yu, DK Smith, H Zhu, Y Guan, TTY Lam. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution. doi:10.1111/2041-210X.12628. 軟件發(fā)表剛滿二年,，引用已經206次(Google scholar, 截止19年1月6日),，比Nature文章的平均引用水平還高一倍。

研究基因功能的人建個樹,，需要找近緣物種,、外類群十幾至幾十個物種，費N天的勁才能做個樹,。而宏基因組領域的人不用去收集其它物種,，因為研究的對象本身就是成百上千的物種，為了方便閱讀或展示主要信息,，我們反而會挑選結果中前100以內的物種去分析并展示,。這是我們絕對的優(yōu)勢,。

文中使用所有文件都是擴增子分析常用結果,，后回復“擴增子”獲取

篩選合適數(shù)量的物種

正如上面所説，擴增子結果有成千上萬的OTU,，全畫是看清的,。需要挑重點，怎么挑,，才能數(shù)據(jù)量即不多,，也不少呢？

我們在《擴增子分析解讀6進化樹》中構建了進化樹result/rep_seqs.tree文件,。這是基于全部的OTU,，用于計算Alpha和Beta多樣性的。我們需要篩選高豐度的OTU,，構建進化樹用于展示,，比如豐度大于0.5%或0.1%。《擴增子分析解讀7篩選進化樹》中采用0.1%豐度篩選獲得104條序列,，其實還是有點多,，因為過百的序列展示過于密集，人類讀起來有困難，用于圈圖360度展示還可以接受,，但是矩形不太適合,。本文為了展示多種組合圖形，采用OTU豐度大于0.5%的篩選閾值,。

# 選擇OTU表中豐度大于0.5%的OTU
filter_otus_from_otu_table.py --min_count_fraction 0.005 -i result/otu_table4.biom -o temp/otu_table_k5.biom
# 獲得對應的fasta序列
filter_fasta.py -f result/rep_seqs.fa -o temp/rep_seqs_k5.fa -b temp/otu_table_k5.biom 
# 統(tǒng)計序列數(shù)量,，28條；30條以字比較大,，容易看清標簽
grep -c '>' temp/rep_seqs_k5.fa # 28
# clusto多序列比對
clustalo -i temp/rep_seqs_k5.fa -o temp/rep_seqs_k5_clus.fa --seqtype=DNA --full --force --threads=30
# 使用qiime的腳本調用fastree建樹
make_phylogeny.py -i temp/rep_seqs_k5_clus.fa -o temp/rep_seqs_k5.tree
# 格式轉換為R ggtree可用的樹,，之前加載報錯，后來刪除'符號后正常
sed "s/'//g" temp/rep_seqs_k5.tree > result/rep_seqs_k5.tree # remove '
# 獲得序列ID
grep '>' temp/rep_seqs_k5_clus.fa|sed 's/>//g' > temp/rep_seqs_k5_clus.id
# 獲得這些序列的物種注釋,，用于樹上著色顯示不同分類信息
awk 'BEGIN{OFS="\t";FS="\t"} NR==FNR {a[$1]=$0} NR>FNR {print a[$1]}' result/rep_seqs_tax_assignments.txt temp/rep_seqs_k5_clus.id|sed 's/; /\t/g'|cut -f 1-5 |sed 's/;/\t/g' |cut -f 1-5 > result/rep_seqs_k5.tax

OTU按分類門級別上色

在Rstudio中,，設置工作目錄為下載測序數(shù)據(jù)的目錄。有三種方法可選： Ctrl+Shift+H,，或在Session菜單中選擇,，或使用setwd()命令設置工作目錄。

# 運行環(huán)境R3.4.1,，ggtree版本1.8.1
# 安裝ggtree包,，末安裝者將FALSE改為TRUE即可
if (FALSE){
  source("https:///biocLite.R")
  biocLite(c("ggtree"))
}

# 設置工作目錄：選擇 Session - Set working directory - To source file location,
# 我們的腳本ggtree.r位于測試數(shù)據(jù)文件夾根本目錄，執(zhí)行上面的操作可調置工作目錄為腳本所在文件夾
rm(list=ls()) 
# 設置工作文件夾進入result,，我們使用的大部分文件在此目錄
setwd("result")
# 加載ggtree包
library("ggtree")

# 讀入分析相關文件
# 讀取樹文件
tree <- read.tree("rep_seqs_k5.tree")
# 讀取樹物種注釋信息
tax <- read.table("rep_seqs_k5.tax", row.names=1)
# 物種注釋等級標簽,，共七級，但細菌末分類物種太多,，一般只能在門,、綱、目水平比較確定
colnames(tax) = c("kingdom","phylum","class","order")

# 按門水平建樹并上色
## 給每個OTU按門分類分組,，此處可以更改為其它分類級別,，如綱、目等,，即phylum替換為order或class即可
groupInfo <- split(row.names(tax), tax$phylum) # OTU and phylum for group
## 將分組信息添加到樹中
tree <- groupOTU(tree, groupInfo)
# 畫樹,，按組上色
ggtree(tree, aes(color=group))+  theme(legend.position = "right")+geom_tiplab(size=3)

按門水著色的矩形進化樹。我們可以看到我們基于V5-V7 rDNA區(qū)構建的樹與分類學門水平相關性較好,，如上面一大枝藍色為Proteobacteria(變形菌門),，下面紅色的一大枝為Actinobacteria(放線菌門)；也有一些異常結果,，如最下面的藍色變形菌門與黃色擬桿菌門聚在一起,，導致分類與進化關系不一致原因很多，如V5-V7甚至rDNA并不能代表菌的真實分類,、相同rDNA也可能有不同功能或新物種等原因,。

畫圈圖

# 畫圈圖并保存PDF
pdf(file="ggtree_circle_color.pdf", width=9, height=5)
## tiplab2保證標簽自動角度,，默認無圖例，要顯示需要+theme
ggtree(tree, layout="fan", ladderize = FALSE, branch.length = "none",aes(color=group))+geom_tiplab2(size=3)+ theme(legend.position = "right")
dev.off()

樹+樣品豐度熱圖

我們最常用的是高豐度,、或差異OTU的樹圖,，結合各樣品的表達熱圖，來展示各樣品間的重復性好壞,、以及組間的差異,。

# 樹+豐度熱圖
# 思路：矩形樹右端添加每個樣品的表達豐度。
## 讀取OTU表
otu_table = read.delim("otu_table.txt", row.names= 1,  header=T, sep="\t")
## 讀取實驗設計
design = read.table("design.txt", header=T, row.names= 1, sep="\t")
## 取實驗設計和OTU表中的交集:樣本可能由于實驗或測序量不足而舍棄掉,，每次分析都要篩選數(shù)據(jù)
idx=intersect(rownames(design),colnames(otu_table))
sub_design=design[idx,]
## 按實驗設計的樣品順序重排列
otu_table=otu_table[,idx]
## 將OTU表count轉換為百分比
norm = t(t(otu_table)/colSums(otu_table,na=T)) * 100 # normalization to total 100
## 篩選樹中OTU對應的數(shù)據(jù)
tax_per = norm[rownames(tax),]

## 保存樹圖于變量,，align調置樹OTU文字對齊，linesize設置虛線精細
p = ggtree(tree, aes(color=group))+  theme(legend.position = "right")+geom_tiplab(size=3, align=TRUE, linesize=.5)
p
pdf(file="ggtree_heat_sample.pdf", width=9, height=5)
## 添加數(shù)字矩陣
## offset設置兩者間距,，用于解決圖重疊問題,；width設置熱圖相對樹圖的寬度，解決熱圖和樹圖大小關系,；font.size設置熱圖文字大小,，解決文字過大重疊；colnames_angle調整熱圖標簽角度,，解決文字重疊問題,；hjust調整熱圖標簽位置，解決文字與熱圖重疊問題,。
gheatmap(p, tax_per, offset = .15, width=3, font.size=3, colnames_angle=-45, hjust=-.1)
dev.off()

現(xiàn)在我們看到了所有高豐度OTU的進化樹,，結合所有樣品的相對豐度熱圖。樹圖為了讓標簽同時作為右側熱圖的標簽,，采用了添加虛線左對齊的方式,；熱圖展示樣品中每個OTU的相對豐度百分比，我們可以看到各組內樣品間的重復情況,，也能看到各組間是否有差別,，比如OTU_1在WT中整體偏高,，而OTU_111卻在WT中整體偏低,。樣品名為了顯示全采用傾斜45度角，右側還有門分類學,、相對豐度百分比的圖例,。

樹+組均值熱圖

# 樹+ 組均值熱圖
## 有時樣本過多也無法展示和閱讀，需要求各組均值展示：需要將分組信息添加至樣品相對豐度表,，再分類匯總
## 提取實驗設計中的分組信息
sampFile = as.data.frame(sub_design$genotype,row.names = row.names(sub_design))
colnames(sampFile)[1] = "group"
## OTU表轉置,，讓樣品名為行
mat_t = t(tax_per)
## 合并分組信息至豐度矩陣，并去除樣品名列
mat_t2 = merge(sampFile, mat_t, by="row.names")[,-1]
## 按組求均值
mat_mean = aggregate(mat_t2[,-1], by=mat_t2[1], FUN=mean) # mean
## 去除非數(shù)據(jù)列并轉置
mat_mean_final = do.call(rbind, mat_mean)[-1,]
## 重命名列名為組名
colnames(mat_mean_final) = mat_mean$group

## 按組均值熱圖
pdf(file="ggtree_heat_group.pdf", width=7, height=5)
gheatmap(p, mat_mean_final, offset = .05, width=1, font.size=3, hjust=-.1)
dev.off()