R语言编程入门进阶：谈谈可视化（二）

R语言编程入门（进化树的可视化）

写在前面的

        上篇推文以RNA-seq的样本间相关性分析为例，以样本间相关性系数矩阵为基础，展示了3种不同的可视化方法。从根本上理解可视化的数学含义！可视化并不单单只是作图，可视化是通过数学图形或者颜色，以一定的映射关系呈现并反映数据的变化趋势。同一个数据有多种不同的方式可视化，不同的数据也可以使用同一种方式进行可视化。所以不要去盲目学习可视化软件，也不要轻信各种讲R语言可视化的培训班，因为没有良好的数学理解能力，学再多的可视化技巧都是等于0

进化树的可视化

        由于可视化的类型实在很多，上次在我发布推文后不久就有人私聊我说看见我在推文里用了一颗聚类树来展示样本间相关性距离。这个聚类树横着倒过来很像进化树。她问我能不能做进化树。其实进化树也是一种聚类形式，进化树是将物种的遗传距离进行聚类。由于涉及到一些基本的生物学假设，进化树又和普通的聚类树稍有不同。今天我们来学习学习如何做进化树。今天的内容不只是单纯画个树，我们一起来学习如何从NCBI的数据->可视化树。
我记得在我上学期的生物信息学课上老师说过，同样是进化树，有的人能画个10+20+的paper，有的人只能画个1,2分的paper。这句话是真的把进化树的内涵说的一清二楚。话不多说，我们直接看门见山讲流程。要画进化树，首先要有目的，而本次内容仅仅只是教学，我们以一部分Ebola病毒的全基因组为基础，研究研究这些埃博拉病毒的系统发生关系。首先在NCBI的Nucleotide搜索Ebola virus, complete genome。

然后弹出如下窗口，在你想要的基因组前画勾，点击send，选择fasta格式。

        拿到序列以后，我们可以进行下一步分析了。首先是做多序列比对，多序列比对的软件有很多，如mafft，clustalW，phylip等等。在这里我们选择mafft软件，因为我们要做Ebola病毒的全基因组比对，这个软件做多序列比对非常快。先看看mafft的教程

        主要参数是输出的格式还有线程，以及比对的模型。localpair是局部比对，利用的 L-INS-i 算法，是相对最准确的比对模式。genafpair利用E-INS-i算法选择序列中有较多不可比对区域的序列, 以及globalpair,利用G-INS-i算法，适合于序列长度相似的序列。但是mafft在default模式下会自动根据序列的特征选择合适的比对模型。所以我们进行多序列比对使用如下命令

mafft --auto --thread 12 Ebola.fasta > mafftEbola.fasta

        多序列比对完成后就可以构建进化树了，为了快速。我们选择FastTree软件，这个软件使用最大似然法进行快速系统发生分析。而且准确率也很高。由于是核酸序列，我们选择gtr模型。

FastTree -gtr -fastest -n 1000 mafftEbola.fasta > Ebola_tree

        FastTree会生成一个树文件，该文件的内容如(((A,B),(C,D)),E)类型。这是一种用符号语言来描述聚类关系的形式。我们要使用R语言将其可视化。R语言有一个函数包叫ggtree，这是Y叔写的函数包，这个函数包是十分强大的进化树可视化工具，能与ggplot2兼容。我个人认为用R语言画树有三个难点，1是读取文件，2是定根（reroot），3是画bootstrap值。目前网上对于这三个难点的综合运用不全，我在这里post出我解决这三个问题的代码。ggtree设置了很多函数用于读取不同进化树软件输出的树文件，例如read.raxml函数是读取raxml软件输出的树文件，也有read.tree函数读取树文件。由于可读的软件类型非常多，在这里我用一个列表画出可读的文件和对应函数

        但是，我们这里恰好需要用另外一个函数读取，我们是FastTree输出的带有Bootstrap值的树文件，这个是newick格式的文件。要使用treeio函数包的read.newick函数读取。用read.newick函数读取FastTree树文件后会产生一个数据格式为treedata数据格式的变量，这个变量中包含的内容可以用@符号跟踪查看，可以使用is函数查看数据格式。例如下面的代码中，tree是储存树文件的初始变量，tree@phylo则是tree变量中的phylo变量.
        Reroot：ggtree在默认模式下输出的可视化树不一定是我们想要的，因为不一定是我们想要的外群定根，所以这里我们需要手动选择外群定根。由于我只下载了Ebola的序列，我们这里假设一个外群，即编号为KR063671.1的外群。（实际情况下外群一定要是跟你的目标序列相似但又不能太相似的序列，例如选择同一属下不同科的物种序列）。定根使用的函数是root函数，这个函数是ape函数包里的函数，当然ggtree函数包里也有reroot函数，phytools里也有reroot函数。网上（Biostar）我看见有教程用ggtree的reroot有报错，但是我相信应该是使用者没有理解好函数。由于我在这里使用的是read.newick函数读取的树，已经是phylo格式了，我需要使用能够处理phylo格式的函数，所以我选择的是ape函数包的root函数。这个函数首先要输入你的树变量，其次有一个outgroup参数输入你想定外群的节点（node），具体的节点你可以在tip.label中查看每个序列对应的节点。tip.label是储存在phylo中，而phylo又储存在tree里，使用tree@phylo号是什么意思，如果你们想找答案，我建议你们is(tree),is(tree@phylo)就知道答案了，因为他们的数据类型不同
        画bootstrap值，要将bootstrap值画出来就要从读取文件那一步开始做准备。我在第一步读文件时，使用了node.label标记了bootstrap值，用support表示。（如果换成bootstrap就报错了，我争取读读源码以后解释这是为什么）事先标记了bootstrap的label后就可以在后面使用ggplot2的语法geom_label将bootstrap值可视化了。接下来我贴上我的全部代码

library(ggtree)
library(phytools)
library(treeio)
library(ggplot2)
#使用treeio包的read.newick函数读取树文件。
tree<-treeio::read.newick('Ebola_tree',node.label='support')
#定根，选择KR063671.1作为外群
tree@phylo<-root(
  tree@phylo,
  outgroup = tree@phylo$tip.label[which(tree@phylo$tip.label=='KR063671.1')],
  resolve.root = T
)
#画树，树枝颜色调橙色，调大小以及线条，geom_treescale()则是画出树的标尺，ggtitle是写入树的名称，xlim则是调整显示的长度。
p<-ggtree(tree,color='orange',size=1,linetype=1,ladderize = T)+
  geom_treescale()+ggtitle(label = 'Ebola Phylogenetic Tree')+xlim(0,0.02)
#画bootstrap值
p<-p+geom_label(aes(label=support,fill=support),size=2)
#调整树末端字体大小，以及将bootstrap值按照高值为红色，低值为暗绿色调整
p<-p+geom_tiplab(size=3.0)+
  scale_fill_continuous(low='darkgreen',high = 'red')
#调整树的主题，将图例显示在右端
p<-p+theme_tree2(legend.position='right')

树的结果如下

重点

     R语言画树有三个难点，1是读取文件，2是定根（reroot），3是画bootstrap值对应的函数要根据实际的系统发生分析软件输出的数据选择合适的函数，如果使用的是我上面列表里写的软件输出的结果，读取函数就用那里面对应的函数，定根便使用ggtree的reroot函数，如果使用的FastTree输出的结果则可以使用我的代码可视化。但是另外强调一点，像使用RAxML软件做完系统发生分析后会输出不同的树文件，有的文件里面没有bootstrap值，一定要选择有bootstrap值的文件才能可视化出bootstrap值。