16s分析之进化树+差异分析(二)
上回我们说到,R语言16s进化树的构建,单单一个优秀的进化树对于我们不是做进化的人来讲,还是不够,有些华而不实,但是Y叔这个ggtree功能当然不止如此,Y叔也经常宣称其为真正的’支持图形语法’;相比功能和可塑性都会提升,就下来我们就进一步深入了解ggtree
这里我要做的就是将进化树连同OTU丰度信息展示出来,ggtree可以做,我们就来一一展示,套用上次推送的代码做前面处理:(如果没有看,请点击:16s分析之进化树+差异分析(一))
p<-ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group)) %<+% tax1+
geom_tiplab2(size=1.5)+ theme(legend.position= "right")
#查看节点;这里简单说几句,geom_label2是Y树参照ggplot的geom_label在ggtree中重新写的,主要特点就是可以使用subset,直接进行数据取舍,这里的意思就是不显示进化树叶端点
p+geom_label2(aes(subset=!isTip,label=node),size=2)
##加上阴影,这里我们可以根据自己的需求选择合适的端点来添加高亮阴影,就16s数据,我挑选了这么几个大类
p1<-p+
geom_hilight(node=323, fill="gray",alpha=0.5) +
geom_hilight(node=320, fill="pink",alpha=0.5) +
geom_hilight(node=313,fill="beige", alpha=0.5) +
geom_hilight(node=298,fill="yellow", alpha=0.5)+
geom_hilight(node=255, fill="blue",alpha=0.5) +
geom_hilight(node=237, fill="red",alpha=0.5) +
geom_hilight(node=194,fill="green", alpha=0.5)+
geom_hilight(node=269,fill="lightblue", alpha=0.5)
p1#这就是我们上次得到的图形,这里我们再次展示一下:
#这里我仅仅添加了OTU名称信息
添加数据,今天的重头戏在于添加数据,表示差异,这里表示:我接下来的代码参考了宏基因组的推送,进化树,上期我已经转发过来了,大家可以自行查看(抄作业从不脸红)
###下面我们对进化树添加数据
## 读取OTU表
otu_table =read.delim("otu_table.txt", row.names= 1, header=T, sep="\t")
## 分组文件我这里只选择我这次做的三个组
design =read.table("map_lxdjhg_CK.txt", header=T, row.names= 1,sep="\t")
## 这里我进行筛选
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
## 筛选树中OTU对应的数据
tax_tomato =norm[rownames(tax),]
head(tax_tomato)
# 树+ 组均值热图
## 提取实验设计中的分组信息
sampFile =as.data.frame(sub_design$SampleType,row.names = row.names(sub_design))
colnames(sampFile)[1]= "group"
## OTU表转置,让样品名为行
mat_t = t(tax_tomato)
## 合并分组信息至丰度矩阵,并去除样品名列
mat_t2 =merge(sampFile, mat_t, by="row.names")[,-1]
## 按组求均值
mat_mean =aggregate(mat_t2[,-1], by=mat_t2[1], FUN=mean)
## 去除非数据列并转置
mean_tomato =do.call(rbind, mat_mean)[-1,]
## 重命名列名为组名
colnames(mean_tomato)= mat_mean$group
wt2<-sqrt(mean_tomato)#
p4 <- gheatmap(p3, wt2, low="blue",high="red",offset = 15, width=0.2, font.size=1, hjust=-.1)
plot(p4)
#这个颜色似乎很重口味,换一个呗
p( p1, wt2,offset = 15,width=0.2, low="white", high="black",colnames_position = "top", font.size=2)
# plot
plot(p5)#这个颜色似乎就清爽许多
当然这还不够,达不到我们的要求
###这里调整热图和标签位置
p1<-ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group)) %<+% tax1+
geom_tiplab2(size=1.5, offset = 5)+theme(legend.position = "right")+
geom_hilight(node=323, fill="gray",alpha=0.5) +
geom_hilight(node=320, fill="pink",alpha=0.5) +
geom_hilight(node=313,fill="beige", alpha=0.5) +
geom_hilight(node=298,fill="yellow", alpha=0.5)+
geom_hilight(node=255, fill="blue",alpha=0.5) +
geom_hilight(node=237, fill="red",alpha=0.5) +
geom_hilight(node=194,fill="green", alpha=0.5)+
geom_hilight(node=269,fill="lightblue", alpha=0.5)
p1
p5 <- gheatmap( p1, wt2, offset = 0,width=0.2,low="white", high="black", colnames_position ="top", font.size=1)
# plot
plot(p5)
##下面我们来对热图添加分类信息
p5+geom_cladelabel(node=323,"Galagoidea", offset=10,angle = 350, fontsize=3, barsize=1,offset.text=1.5,color="gray") +
geom_cladelabel(320, "Lemuroidea",offset=9,angle = 10, fontsize=3, barsize=1, offset.text=1.5,color="pink") +
geom_cladelabel(313, "Tarsioidea",offset=11,angle = 50, fontsize=3, barsize=1, offset.text=1.5,color="beige") +
geom_cladelabel(298,"Ceboidea",offset=9.5,angle = 80, fontsize=3, barsize=1,offset.text=1.5,color="yellow") +
geom_cladelabel(255, "Hominoidea",offset=9.5,angle = 305, fontsize=3, barsize=1, offset.text=1.5,color="blue") +
geom_cladelabel(237,"Cercopithecoidea",offset=10.5,angle = 360, fontsize=3, barsize=1,offset.text=1.5,color="red") +
geom_cladelabel(194,"Cercopithecoidea",offset=11.5,angle = 70, fontsize=3, barsize=1,offset.text=1.5,color="green")+
geom_cladelabel(269,"Cercopithecoidea",offset=9,angle = 320, fontsize=3, barsize=1,offset.text=1.5,color="lightblue")
这里角度我没有调整好,所以自己调整,参数如下:
angle:就是设置分类label的角度
offset:距离中心的距离
fontsize:字体大小
barsize:标签柱子宽度
到这里,我们的圈图就到这里了,似乎在圈图在ggtree中只可以添加热图,柱状图,饼图目前也没有看到添加方法;
这里提到Y树ggtree的facet_plot,可以添加各种图形同树图一块展示,但是在layout为fan情况下,要想展示,就必须将其坐标转化为极坐标,而且需要相应调整,因为这种树图和数据结合采用圈图方式的确是逼格满满:
这里展示iTOL的主页:的确是漂亮
于是我问Y树,可不可以做成这种形式,答案在此:
这张代表我的心情:
我的耳旁响起:一千个伤心的理由,一千个伤心的理由,
怎么办呢,等Y叔的扩展,还是自己写一个,我写的了吗?不试试谁知道呢?
学习永无止境,分享永不停歇!