说不定有用（补充）基因家族分析生信分析工具包

Phylogenomic_Tutorial || ML_Tree

2021-02-15  本文已影响0人  Dawn_WangTP

Github/mmatschiner的phylogenetic & phylogenomic学习教程记录【一】多序列比对；核算替换模型的选择；最大似然法建树的学习

[TOC]

软件准备Preparation

Basics

  1. Bash
  2. Ruby 2; 版本> 2
  3. Python 3
  4. R

安装包Libraries

软件Programs

多序列比对

Outline

多序列比对是系统发育分析的基础；其目的是为了确定序列之间同源的区段（homologous nucleotides可比较的序列）；本节为了解决以下问题：

  1. 利用MAFFT软件做多序列比对；
  2. 鉴定并排除可能错误比对的序列区段；
  3. 利用公共数据库NCBI GeneBank补全已有的数据集合；

Datasets

数据来源于Matschiner et al. 2017。此为研究cichlid(慈鲷科鱼种)在各大洋的历史演化及分布。物种包括41species，分为不同地理的groups（Non-cichlid, Neotropical cichlids, African; Indian; Malagasy）。此处利用到两个基因：线粒体mitocondrial 16S, 核基因RAG1

Softwares

利用到的软件包括：

利用MAFFT比对，AliView可视化

## 1
mafft --auto 16s.fasta >16s_aln.fasta

## 2 Gap penalty
mafft --auto --op 2 16s.fasta >16s_op2_aln.fasta

利用BMGE去除低质量比对区域

多序列比对软件对于一些区域比对质量较差，可能会出现比对错误的情况，对后续系统发育分析会产生影响，因此需要将其剔除。我们利用BMGE软件去除这些低质量比对的区域。BMGE(block mapping and gathering with entropy)由JAVA所写。

## 查看帮助
java -jar BMGE.jar -?

## 运行
java -jar BMGE.jar -i 16s_aln.fasta -t DNA -of 16s_filtered.fasta -oh 16s_filtered.html

程序before和after会告诉排除了多少问题比对的序列。后续同样的对核基因rag1也可以做相同的操作；

image

核酸替换模型的选择Substitution Model Selection

当完成多序列比对之后，进行似然法likelihood的系统发育分析之前，随后要进行的就是核酸替换模型的选择，包括有Jukes-Cantor(JC)模型，HKY模型，GTR模型等等。通常是根据核酸替换模型软件计算得出最佳的AIC值(Akaike Information Criterion)的模型，再根据此模型进行后续的系统发育分析。

Preparation

利用PAUP模型筛选

  1. 下载GUI图形化版本的PAUP软件，打开文件16s_filtered.nex

  2. PAUP计算核酸替换模型需要首先利用NJ法快速建树；

  3. 根据"Automated Model Selection"选项进行模型筛选；其中模型中包含有较多的术语：通常筛选标准包括有AIC(Akaike information criterion); AICc(Akaike information criterion corrected for small sample sizes); BIC(Bayesian information criterion); DT(decision-theoretic criterion)。通常是根据AIC值最小的Model选为最佳Model。

    Model for among-site rate variation中G(GAMMA)分布; I(Invariable sites)

  4. PAUP结果中-lnL(-log likelihood); K(the number of free parameters in the model), 代表branch length

  5. 结果显示最佳AIC模型是GTR模型；GTR模型是系统发育分析中最常用的model。

pic

已有的核酸替换模型的筛选都是认为一个模型适用于所有位点，所有的位点进化速率都是相同的。但实际上一些位点的进化速率是不同的，例如cds序列中的third-codon位点是相较于第一第二位点速率是更快的。根据"Automated Partioning"结果显示

最大似然法建树ML phylogenetic Inference

Preparation

利用热门的软件IQTREE做系统发育分析，Figtree做基本的可视化展示。Datasets利用16s_filtered.nex and rag1_filtered.nex

利用IQTREE做ML系统发育分析

iqtree --help

## iqtree 默认最简参数；
iqtree -s 16s_filtered.nex

利用bootstrap法获取ML的各个node支持率

根据默认参数的-s所得到的系统发育树，我们尚无法确定各个支系的支持率程度，可靠度(reliability)，因此我们需要借用bootstrap法获取ml树的支持率；

## -B 为ultrafast的bootstrap值，--prefix为预设输出名字
iqtree -s 16s_filtered.nex -B 1000 --prefix 16s_filtered.bs.nex

对cds核基因分段(Partitioned)的ML系统推断

由于前面做模型选择时我们发现，对于cds的codon(1,2,3)序列来说，不同位置会用不同的替换模型，在IQTREE中也可以同样根据此进行Partitioned inference；

#

#NEXUS
  BEGIN SETS;
    CHARSET codon1 = 1-1368\3;
    CHARSET codon2 = 2-1368\3;
    CHARSET codon3 = 3-1368\3;
  END;
  
  ### 将以上输出到partitions.txt文件中
  
## run iqtree
iqtree -s rag1_filtered.nex -p partitions.txt -B 1000 --prefix rag1_filtered.bs.nex

比较不同ML树之间的差异

通过不同的基因我们会获得不同的系统发育树。而为检查两个树之间的差异以及评估整体的差别，我们可以利用Robinson-Foulds distance来鉴定树之间拓扑结构差异。在IQTREE中的-t和--tree-dist2

iqtree -t 16s_filtered.bs.nex.treefile --tree-dist2 rag1_filtered.bs.nex.treefile

同样也可以计算不同树文件中的平均bootstrap值判断支持率。

串联比对(concatenated alignments)下的系统推断

一个基本前提是，序列位点越多，最后获得各个支系的支持率也就越大，所以我们可以将多个基因串联起来统一进行系统发育推断。当然另一个前提是不同的基因会随着物种进行同样的进化历程，这个前提在对一些近缘物种推断时是会存在问题的。


#NEXUS
  BEGIN SETS;
        CHARSET 16S = 16s_filtered.nex: *;
        CHARSET rag1_codon1 = rag1_filtered.nex: 1-1368\3;
        CHARSET rag1_codon2 = rag1_filtered.nex: 2-1368\3;
        CHARSET rag1_codon3 = rag1_filtered.nex: 3-1368\3;
  END;

## 根据组合方式进行统一iqtree计算 -p
iqtree -p partitions.txt -B 1000 --prefix concatenated.bs.nex

结果确实会比单基因建树的支持率效果更好。多基因串联建树也是后期由基因树推断物种树的一个重要的方法。

参考资料

Substitution Model Selection
Maximum-Likelihood Phylogenetic Inference

上一篇下一篇

猜你喜欢

热点阅读