之前我便说过,16s出图主要是用R,自这两年来Y叔的ggtree出来,好评不断,引用不断攀升
这里当然借用Y树的包来做我们的进化树,当然Y叔公众号讲的更加深入和详细,有需要请关注:biobabble
这里Y树曾经详细介绍过ggtree,里面还顺带介绍了ggplot,有兴趣的亲查看(反正我收获颇丰):
这里我参考了宏基因组公众号的代码:(参考别人的嘛,也不难为情,工具这种事情,总是不多看别人怎么做,自己学着搞)
这是刘永鑫的博客:http://blog.sciencenet.cn/blog-3334560-1079159.html
现在我们开始:
这里我的R包很老了,Y叔的包还是很新的,不多说,只能更新R包了
##################首先我更新了一下R,已经不能安装ggtree###################
#由于R太老了,所以我跟新了R
#安装更新工具
install.packages("installr")
library(installr)
updateR()
#安装ggtree
source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")
#更新的R包有好多依赖包没有,手动安装,还有好多包需要升级,慢慢弄吧
source("https://bioconductor.org/biocLite.R")
biocLite(c("treeio"))
install.packages("tidyr")
##################首先我更新了一下R,已经不能安装ggtree###################
下面开始我们的分析,在画进化树之前首先我们的前处理在qiime中:
# 选择OTU表中丰度大于0.1%的OTU,足够我们分析的了
filter_otus_from_otu_table.py--min_count_fraction 0.001 -i otu_table.biom -o hg_tree/otu_table_k1.biom
# 获得对应的fasta序列rep_seqs.fa
filter_fasta.py-f rep_seqs.fa -o hg_tree/rep_seqs_k1.fa -b hg_tree/otu_table_k1.biom
# 统计序列数量,正则表达式,如果看不懂,百度搜索,一大堆,这里我简单说一下,单引号中即为所搜索的字符,-c就是计算找到 '搜寻字符串' 的次数,输出含有>的次数
grep -c '>'hg_tree/rep_seqs_k1.fa
这里我们安装clustalo
#qiime已经集成了安装多序列比对的依赖包了,我下载的版本:clustalo-1.2.4-Ubuntu-x86_64,使用qiime1.9.2亲测可直接安装此包,无需依赖
sudo apt-get installclustalo
# clusto多序列比对
clustalo -ihg_tree/rep_seqs_k1.fa -o hg_tree/rep_seqs_k1_clus.fa --full --force--threads=5
# 使用qiime的脚本调用fastree建树
make_phylogeny.py-i hg_tree/rep_seqs_k1_clus.fa -o hg_tree/rep_seqs_k1.tree
# 格式转换,这里我们要提一下,为什么要做格式转换,这就是我们使用qiime的make_phylogeny.py命令做的树文件,首先它是一个无根树,其次,最重要的是,每个OTU名称都被单引号夹起来了,当然我们不需要这种形式,所以下面开始去除
到这里,可能就有很多人都有满肚子疑问了,tree的文件到底是一种什么格式,我们这种格式得给的通俗的解释啊!别着急,下面我来做这个工作:
Newick格式,就是我们使用qiime做进化树的格式,也是最为古老的格式,是一种使用圆括号和逗号代表具有节点长度理论图形树的一种方法。
library("ggtree")
library("ggplot2")
#只命名叶节点
tree1 <-"((A, B), C,D);"
x1 =read.tree(text = tree1)
ggtree(x1)+geom_tiplab()
#没有任何一个节点被命名
tree2 <-"((,),,);"
x2 =read.tree(text = tree2)
ggtree(x2)+geom_tiplab()
#所有节点均命名
tree3 <-"((A, B)E, C,D)F;"
x3 =read.tree(text = tree3)
#这里我显示叶节点并调整了叶节点标签,但是如何显示节点标签,我就不会
ggtree(x3)+geom_tiplab(offset=1,size=2,align=1)
#
#除了根节点,所有节点均命名
tree3 <-"((A, B)E, C,D);"
x3 =read.tree(text = tree3)
#这里我显示叶节点并调整了叶节点标
ggtree(x3)+geom_tiplab(offset=1,size=2,align=1)
#这里给加上了节点之间的距离
tree3<-"(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;"
x3 =read.tree(text = tree3)
#这里我显示叶节点并调整了叶节点标签,但是如何显示节点标签,我就不会
ggtree(x3)+geom_tiplab()+xlim(NA, 1) + ylim(NA, 5)#offset=1,size=2,
ggtree(x3)+geom_tiplab()+xlim(0, 1) + ylim(0, 5)
我们在上面简单了解了一下树的知识,下面我们来做一颗树
# 这次我使用番茄otu
setwd("E:/Shared_Folder/HG_usearch_HG/ggtree_hg")
# 读入分析相关文件
# 读取树文件
tree3 <-read.tree("rep_seqs_k1.tree")
# 读取树物种注释信息
tax <-read.table("rep_seqs_k1g.txt",row.names=1)
head(tax)
tax$id=rownames(tax)
head(tax)
# 物种注释等级标签,六级
colnames(tax) =c("kingdom","phylum","class","order","family","genus","id")
## 给每个OTU按门分类分组
groupInfo <-split(row.names(tax), tax$phylum)
tree3 <-groupOTU(tree3, groupInfo)
我们来出一张默认树图看看,可以看到树建立在活动窗口,当前设备,这是一颗无根树
ggtree(tree3)
#使用branch.length= "none",叶对齐,默认不是对齐的
ggtree(tree3,branch.length = "none")
#调整树的结构ladderize= F
ggtree(tree3,branch.length = "none",ladderize = F)
#改变树的形式,按照跟组上色
ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group))
#调整树的开合角度和按照一定角度旋转
p<-ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group))
open_tree(p,angle=100)%>% rotate_tree(-45)
#我们是在没有必要写这么多注释,这里我们意思一下
p<-ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group)) %<+% tax1+
geom_tiplab2(size=1.5)+ theme(legend.position= "right") +
geom_tiplab2(aes(label=phylum,color=group),size=2, offset=6)
#查看节点,为什么查看这个,因为我们要按照节点给他高亮
p+geom_label2(aes(subset=!isTip,label=node),size=2)
##加上阴影
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
今天的推送就到这这里,但是没有完,树的这些东西还有好多,我正在准备,下次发
这两天的确是发文不如之前迅猛了,不是疲软,还没有到那个年纪,只是在做一些其他的分析而已
最后还是那句话:
学习永无止境,分享永不停歇!
联系客服