系统进化树构建的简单实践
在我看了好多一些系统发育树相关内容后，是时候实战让大脑对构建进化树有一定印象，不然就会陷入知识阻塞，无法继续前行。
不要企图一开始就很完美，可以说任何技能都是从不完美开始的。
此外，这里不涉及软件安装，我相信你们会的，如果不会，先Google，后百度，总有一款适合你
序列查找
工具推荐：BLAST, PSI-BLAST, hmmalign/hmmsearch
由于我只是想快速上手建树，所以我直接去planttfdb下载拟南芥AP2转录因子家族30个蛋白序列
mkdir AP2_Ath
cd AP2_Ath
wget 'http://planttfdb.cbi.pku.edu.cn/download_seq.php?sp=Ath&fam=AP2'
mv download_seq.php\?sp\=Ath\&fam\=AP2 ap2序列比对
使用mafft快速但是不太精确的比对参数（也就是一切都是默认）。你可能会问，万一比对的不好怎么办？
刚开始，不要在意这些细节，等你看完说明树再去操作，你都忘了你本来的目的。（不要问我为什么知道，说出来都是泪，还有我是处女座）
mafft ap2 > ap2.alnalignments curation/refinement
抱歉，我不知道如何翻译， 比对结果管理？比对结果修缮？
用BioEdit查看序列比对结果。当然BioEdit还有许多好用的功能，这里就用它的序列编辑功能吧。（Y叔搞出了一个MacOS版，有兴趣的小伙伴，添加biobabble去历史记录里找找吧)
这时候出现了第一个问题，我应该如何处理’-‘部分，如下是我求助后的解答：
gap分两种情况，一种是这部分序列没测到，为了alignment，用“-”补齐；另一种情况是经过alignment后，为了同源区段的完全匹配，填补上的gap，在这种情况下，gap也是具有系统发育信号的。
构树过程中有多种选择，你可以多做尝试。如果感觉序列align之后对齐效果并不好，就用Gblock软件切一下，把 含gap多的位点和对齐不好的位点切掉，再构树。
构树不同的方法和不同的软件对gap的处理也不一样，可以选择有gap的位点不用，也可以选择部分包括或完全包括gap，看看构树的结果是否合理，支持率够不够高 —From栾云霞首先建树的理论基础是要求来自于同原序列，所以他的输入是一个alignment。因此把它删掉也有一定道理，因为这样子我们可以保留那些保守的区域来建树。但是也不可以简单粗暴的就这么做了我们需要看看gap是怎么来的：比如说有少量序列特别长，在两头引入gap，那么你可以放心把它切掉，如果说你有一些序列特别短，然后引入了很多gap的话，那么这些序列就要删掉，不能用于alignment。还有一点，比对结果不一定靠谱，它可能会实际额外引入gap，实际上没有那么多gap，可以去看我公众号里面推送的BioEdit里面的截图。 —From Y叔
对于我这个刚开始入门系统进化领域的人，我只有如下表情。
于是我下载PlantTFDB上他们的比对结果，如图：
什么他居然直接拿这些序列去建树了？想想也是，毕竟这是一个pipeline，哪有时间手工检测。
我尝试一下用Gblocks切一下吧。这是切完后的结果（重名为ap2_gb.aln，看起来似乎很不错，那这段序列去进行下一步把。
其实看了一下他生成的网页解释他选择的区域，我感觉一点都不好。但是不要在意这些细节，你又不可能一次就建好进化树。
分子进化分析
目前有如下几种方法：
距离矩阵法
最大简约法
极大似然法
贝叶斯法
目前主流比较喜欢用极大似然法和贝叶斯法。推荐工具如下：
Phyml(可信度最高，但是速度丧心病狂的慢，请在最后使用)
RAxml（次之，速度不错）
MrBayes
更多工具见http://evolution.genetics.washington.edu/phylip/software.html ， 相信我，选择太多反而等于没有选择。
所以我就拿RAxml建树吧。根据https://sco.h-its.org/exelixis/web/software/raxml/hands_on.html 的推荐，我决定用最快速度的命令。
raxmlHPC-PTHREADS-AVX -T 4 -p 12345 -m PROTGAMMAWAG -s ap2_gb.aln -n T8然后执行指令的时候我发现，他说有11条序列是一模一样，根据manual，这些序列可能不会用于建树。
最后我得到了如下结果：
计算自展值(bootstrp)
自展值是干什么用的？我不知道，这里写出来只是因为GGTREE需要这个内容，所以我加了这一步。当然用的还是速度最快的那个指令。不过这一次可以同时进行ML搜索和bootstraping。
raxmlHPC-PTHREADS-AVX -T 4 -f a -p 12345 -x 12345 -# 100 -m PROTGAMMAWAG -s ap2_gb.aln -n T9可视化
这里我推荐用Y叔的GGTREE,谁用谁知道，就是需要一点R语言基础。至于其他，我又不会只更新这一篇关于进化树，毕竟以后都要做，学一点发一点了，所以麻烦大家在留言区踊跃推荐哦。
代码如下：
library(ggtree)
raxml <- read.raxml("C:/Users/Xu/Desktop/RAxML_bipartitionsBranchLabels.T9")
lb = get.tree(raxml)$tip.label
d = data.frame(label=lb, label2 = paste("AA", substring(lb, 1, 5)))
ggtree(raxml) %<+% d + geom_tiplab(aes(label=label2))结果如下：
为什么那么丑，和Y叔给我们看的不一样呀。
因为不是拿来发表的呀，差不多能看就成了。还有我去看了一下PlantTFDB的进化树，发现他们是根据蛋白和domain都进行建树了。
结语
这是第一棵不适用MEGA建的树，我对这个结果肯定是存在怀疑，甚至是根本不信，但是我至少知道建树至少有以下几步
寻找序列
序列比对
序列比对内容管理，或者叫序列比对调优
根据序列进行系统发育树的极大似然法推测
可视化展示
然后在这个过程中，我遇到了无数的问题，任何一个问题可能都需要我半天或者更久去寻找答案（还好有哪些愿意帮助我的人）。为了不在任何一个细节上浪费太多时间，我磕磕碰碰，迷迷糊糊的把结果做了出来。目前有以下问题要进行解决：
序列到底应该如何寻找，是只需要motif么，还是全部氨基酸序列
比对参数应该如何优化
比对结果应该如何调整
系统进化树应该选用什么方法，如何如何选择每个方法中的模型
最后的可视化，如何更好看
每一个问题都是无数的细节，都需要很多经验的累积。
在之前，我提问的方式是: 我有XX序列，如何建系统发育树呢？
而现在，我相信我的提问将会更加具体。
我建议新手最好去找一些比较少的序列，先用MEGA图形界面，了解基本的流程。然后如果你要处理大量数据，使用更加专业的软件，提高性能。
这是我学习系统发育树的第一篇，根据我以前的习惯，后面肯定会出现一波相关的文章的。