搜索
查看: 1421|回复: 5

系统进化树构建的不简单实战

[复制链接]

18

主题

55

帖子

343

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
343
发表于 2017-5-21 21:29:04 | 显示全部楼层 |阅读模式
本帖最后由 hoptop 于 2017-5-21 21:46 编辑

在我看了好多一些进化树相关内容后,需要实战让大脑对构建进化树有一定印象,不然就会陷入知识阻塞,无法继续前行。不要企图一开始就很完美,可以说任何技能都是从不完美开始的。
此外,这里不涉及软件安装,我相信你们会的,百度和谷歌,总能解决问题。



序列查找
工具推荐:BLAST, PSI-BLAST, hmmalign/hmmsearch
由于我只是想快速上手建树,所以我直接去planttfdb下载拟南芥AP2转录因子家族30个蛋白序列
[AppleScript] 纯文本查看 复制代码
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快速但是不太精确的比对参数(也就是一切都是默认)。你可能会问,万一比对的不好怎么办?
刚开始,不要在意这些细节,等你看完说明树再去操作,你都忘了你本来的目的。(不要问我为什么知道,说出来都是泪,还有我是处女座)

[Bash shell] 纯文本查看 复制代码
mafft ap2 > ap2.aln






alignments 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 的推荐,我决定用最快速度的命令。
[Bash shell] 纯文本查看 复制代码
raxmlHPC-PTHREADS-AVX -T 4 -p 12345 -m PROTGAMMAWAG -s ap2_gb.aln -n T8

然后执行指令的时候我发现,他说有11条序列是一模一样,根据manual,这些序列可能不会用于建树。
最后我得到了如下结果:


计算自展值(bootstrp)
自展值是干什么用的?我不知道,这里写出来只是因为GGTREE需要这个内容,所以我加了这一步。当然用的还是速度最快的那个指令。不过这一次可以同时进行ML搜索和bootstraping。

[AppleScript] 纯文本查看 复制代码
raxmlHPC-PTHREADS-AVX -T 4 -f a -p 12345 -x 12345 -# 100  -m PROTGAMMAWAG -s ap2_gb.aln -n T9





可视化
这里我推荐用Y叔的GGTREE,谁用谁知道,就是需要一点R语言基础。至于其他,我又不会只更新这一篇关于进化树,毕竟以后都要做,学一点发一点了,所以麻烦大家在留言区踊跃推荐哦。
代码如下:
[AppleScript] 纯文本查看 复制代码
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))

如果raxml没有计算bootstrap,那么就是标准的newick格式,所以可以直接用read.tree导入文件。都不要前面计算bootstrap了

结果如下:
为什么那么丑,和Y叔给我们看的不一样呀。

因为不是拿来发表的呀,差不多能看就成了。还有我去看了一下PlantTFDB的进化树,发现他们是根据蛋白和domain都进行建树了。



结语
这是第一棵不适用MEGA建的树,我对这个结果肯定是存在怀疑,甚至是根本不信,但是我至少知道建树至少有以下几步
  • 寻找序列
  • 序列比对
  • 序列比对内容管理,或者叫序列比对调优
  • 根据序列进行系统发育树的极大似然法推测
  • 可视化展示
然后在这个过程中,我遇到了无数的问题,任何一个问题可能都需要我半天或者更久去寻找答案(还好有哪些愿意帮助我的人)。为了不在任何一个细节上浪费太多时间,我磕磕碰碰,迷迷糊糊的把结果做了出来。目前有以下问题要进行解决:
  • 序列到底应该如何寻找,是只需要motif么,还是全部氨基酸序列
  • 比对参数应该如何优化
  • 比对结果应该如何调整
  • 系统进化树应该选用什么方法,如何如何选择每个方法中的模型
  • 最后的可视化,如何更好看
每一个问题都是无数的细节,都需要很多经验的累积。
在之前,我提问的方式是: 我有XX序列,如何建系统发育树呢?
而现在,我相信我的提问将会更加具体。基本上就拿我上面的疑问了。

我建议新手最好去找一些比较少的序列,先用MEGA图形界面,了解基本的流程。然后如果你要处理大量数据,使用更加专业的软件,提高性能。






本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:安装.sh文件是遇到的问题解决
下一篇:python制作第一个package测试包
回复

使用道具 举报

633

主题

1172

帖子

3947

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3947
发表于 2017-5-21 23:18:23 | 显示全部楼层
你说查看序列比对结果很麻烦,所以我给你写了个网页工具:http://www.biotrainee.com:3838/jmzeng/msaR/
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

18

主题

55

帖子

343

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
343
 楼主| 发表于 2017-5-22 14:01:51 | 显示全部楼层
Jimmy 发表于 2017-5-21 23:18
你说查看序列比对结果很麻烦,所以我给你写了个网页工具:http://www.biotrainee.com:3838/jmzeng/msaR/ ...

厉害了! 我也要学
回复 支持 反对

使用道具 举报

6

主题

28

帖子

195

积分

版主

Rank: 7Rank: 7Rank: 7

积分
195
发表于 2017-8-20 19:56:25 | 显示全部楼层
1.比对结果调整:如果是少量数据的话,比对结果可以人工校对(自己测序的可以看测序峰图,网上数据一般直接看gap数量,太多gap一般直接切掉),如果是大数据的话,用gblock设置参数后选择保守区域,进行调整。
2. 进化树构建选择模型,一般是用模型检测软件,核酸:modeltest,jmodeltest等等,氨基酸:protest
3. 可视化我一般用figtree
“Nothing in Biology Makes Sense Except in the Light of Evolution”  --Theodosius Dobzhansky
回复 支持 反对

使用道具 举报

18

主题

55

帖子

343

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
343
 楼主| 发表于 2017-8-22 11:08:41 | 显示全部楼层
shuboy 发表于 2017-8-20 19:56
1.比对结果调整:如果是少量数据的话,比对结果可以人工校对(自己测序的可以看测序峰图,网上数据一般直接 ...

多谢提醒
回复 支持 反对

使用道具 举报

1

主题

8

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-12-7 17:04:39 | 显示全部楼层
有没有基因组数据做进化分析的内容哇,我手上的数据量太大了,感觉多序列比对极其耗时,后续步骤又不知道怎么进行下去了
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-12-18 18:52 , Processed in 10.115553 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.