ggtree美化进化树

# 选择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 
# 统计序列数量，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

# 运行环境R3.4.1，ggtree版本1.8.1
# 安装ggtree包，末安装者将FALSE改为TRUE即可
if (FALSE){
  source("https://bioconductor.org/biocLite.R")
  biocLite(c("ggtree"))
}
 
# 设置工作目录：选择 Session - Set working directory - To source file location,
# 我们的脚本ggtree.r位于测试数据文件夹根本目录，执行上面的操作可调置工作目录为脚本所在文件夹
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)

# 画圈图并保存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_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表中的交集:样本可能由于实验或测序量不足而舍弃掉，每次分析都要筛选数据
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对应的数据
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)
## 添加数字矩阵
## 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()

# 树+ 组均值热图
## 有时样本过多也无法展示和阅读，需要求各组均值展示：需要将分组信息添加至样品相对丰度表，再分类汇总
## 提取实验设计中的分组信息
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
## 去除非数据列并转置
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()