16s分析之进化树＋差异分析（一）

之前我便说过，16s出图主要是用R，自这两年来Y叔的ggtree出来，好评不断，引用不断攀升

这里Y树曾经详细介绍过ggtree，里面还顺带介绍了ggplot，有兴趣的亲查看（反正我收获颇丰）：

这里我参考了宏基因组公众号的代码：（参考别人的嘛，也不难为情，工具这种事情，总是不多看别人怎么做，自己学着搞）

现在我们开始：

这里我的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

x1 =read.tree(text = tree1)

ggtree(x1)+geom_tiplab()

#没有任何一个节点被命名

tree2

x2 =read.tree(text = tree2)

ggtree(x2)+geom_tiplab()

#所有节点均命名

tree3

x3 =read.tree(text = tree3)

#这里我显示叶节点并调整了叶节点标签,但是如何显示节点标签，我就不会

ggtree(x3)+geom_tiplab(offset=1,size=2,align=1)

#

#除了根节点，所有节点均命名

tree3

x3 =read.tree(text = tree3)

#这里我显示叶节点并调整了叶节点标

ggtree(x3)+geom_tiplab(offset=1,size=2,align=1)

#这里给加上了节点之间的距离

tree3

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

#读取树物种注释信息

tax

head(tax)

tax$id=rownames(tax)

head(tax)

#物种注释等级标签，六级

colnames(tax) =c("kingdom","phylum","class","order","family","genus","id")

##给每个OTU按门分类分组

groupInfo

tree3

我们来出一张默认树图看看，可以看到树建立在活动窗口，当前设备，这是一颗无根树

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

open_tree(p,angle=100)%>% rotate_tree(-45)

#我们是在没有必要写这么多注释，这里我们意思一下

p

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)

##加上阴影

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

今天的推送就到这这里，但是没有完，树的这些东西还有好多，我正在准备，下次发

这两天的确是发文不如之前迅猛了，不是疲软，还没有到那个年纪，只是在做一些其他的分析而已

最后还是那句话：

学习永无止境，分享永不停歇！