生信技工
0%

3s软件

Posted on Edited on In ,
4.4k 4 mins.
介绍了检测基因流的软件3s的原理和使用，3s基于三个物种的基因组检测两个近缘种的基因流的方向和强度。

1. 检测基因流的软件

  1. 基于MCMC算法的：IM,IMA
  2. 基于似然法的：除了3s外的其他软件只能分析两条序列
  3. Huristic method：ABBA(基于D统计量的计算，无数学模型，无法估计),Hyde

2. 3s软件简介

3s，3 species的简称。

软件是杨子恒和朱天琪团队开发，基于C语言的，通过计算似然率来检测两个近缘物种/群体间基因流的软件。推荐添加一个外类群，所以称为三个物种（3 species），适用于基因组数据。

  1. 3s简介
  • 3s利用似然率来推断两个物种/群体间的基因流方向和强度
  1. 3s输入
  • 输入：基因组或其他测序序列phylip文件
  • 输出：基因流方向和强度
  1. 3s优势和不足
  • 随着数据量线性增加运算时间，运算快，适合基因组数据。
  • 一次只能检测三个物种/群体，无法建立系统发育网。

3. 开发历史

  1. 3s软件主页
  1. 2012年，朱天琪和杨子恒开发的检测基因流的模型SIM3s，发布了软件3s.v2.1版本，包含6个参数。
  1. 2017年，提出了Introgression模型，发布了版本3s.v3.1版本，包含9个参数。
  1. 2020年，据作者朱天琪说开发了新的more general migration model模型，包含15个参数，2022.09还未发表。
  • 假设所有物种间都可能有基因，1-2间，2-3间，1-3间，3-5之间都可能存在基因流。
  • 这个模型参数太多，最好根据先验知识尽可能设置一些无基因流的群体间参数为0，再进行计算更准确。
  • 如果实在没有先验知识，那先打开所有可能基因流，运行一遍，再把基因流强度很低的位置关闭，再运行一遍获取更准确结果。
  • 2020年09月据说代码已有，2022.09说还未发表，可以关注着。

4. 3s软件原理

4.1. 模型基础

  1. Wright-Fisher模型：溯祖模型的逆向过程
  2. 多物种溯祖模型
  • 不区分物种间和群体间，这里的物种=群体
  • 无法跨物种边界溯祖

4.2. 3s模型的参数定义

  1. 三个物种：两个近缘种（1和2）和一个外类群（3）
  2. 序列数据由中性位点组成
  3. θ1=4N1μ，θ2=4N2μ，θ4=4N4μ，θ5=4N5μ是物种1，2，4，5的群体大小参数。这里没有物种3是因为假设物种3没有基因流。
  4. τ0, τ1是两个物种的分化时间
  5. 迁移率(migration rate)Mij=Njmij。定义是在真实世界，从群体i向群体j的迁移个体的数量，其中Nj代表群体j的群体大小，mij代表迁移比例。

Figure 1. 3个物种的基因流模型参数
图源：https://academic.oup.com/sysbio/article/66/3/379/2670069?login=false

4.3. SIM3s模型

  1. 通过假设两个参数相等，θ1=θ2=θ；M12=M21=M
  2. 结果是共有6个参数需要评估：Θ={τ0,τ1,θ4,θ5,θ,M}
  3. 这个模型是2012年发表的老版本模型SIM3s(symmetric IM model for three species)

4.4. 3s模型背后的计算

基于最大似然法估计参数：先算似然函数，对参数做最优，似然值最大的对应参数即为估计结果。

  1. 每个位点/序列计算一个似然值，似然值代表的是给定参数后基因树的概率。
  2. 似然值的计算是通过马尔科夫链来描述和实现的，三个样本的初始状态是123。
  3. 在比对好的三个样本的序列上对不同模式进行碱基计数。三个样本的碱基一致为同一种模式，12一致3不一致为一种模式，等等。
  4. 3s软件的突变模型使用JC69模型，即各种碱基的任意方向和类型的替换速率是相等的，适合近缘物种/群体。
  5. 通过似然比检验（Likelihood Ratio Test，LRT）来判断有无基因流：当按照95%显著临界值，2Δℓ>5.99时，则认为检验是显著的，即存在基因流，基因流的大小即最大似然估计值。

5. 3s的优势

  1. 3s的计算基于极大似然法，极大似然法的优势在于计算量随着样本数是线性增加的，而比如贝叶斯,MCMC算法则是指数增加。
  2. 使用的segment/locus数量越多，检验功效（指有基因流存在情况下检出率）越好，所以适合检验基因组这种大数据。
  3. 假阳性很低（<5%可接受）
  4. 基因流的方向和强度都可以被估计。

6. 3s模型假设和数据准备

6.1. 3s模型假设和数据选择

  1. 使用中性位点避免重组的影响。
  2. 3s假设segment之间是自由重组的，segment内部是不存在重组的。所以建议在基因组上取样的序列不要离得太近（建议50kb，至少大于10kb），序列长度不要太长（50bp-1kb）。
  3. 使用的segment/locus数量越多，检验功效（指有基因流存在情况下检出率）越好，所以适合检验基因组这种大数据。

6.2. 常用的数据准备方式

  1. 适用于组装好的基因组数据
  • 将基因组切分为1kb的片段
  • 每隔50kb取1kb的片段
  • 将含有基因和基因上下游10kb的片段去除（基因通常被认为是非中性的）
  1. 适用于获取的别人的partition数据
  • 去除间隔<10kb的partition
  • 删除序列长度<50bp或者>1000bp的partition
  1. 适用于微生物等基因组非常小的样品
  • 有许多微生物样品的基因组很小，能取的片段数量有限，影响分析结果。
  • 可以通过非独立重抽样增加片段数量。
  • 每个物种测序多个个体，共同作为一个数据源，以增加片段数量。比如ABC三个物种都测四个个体，A1B1C1,A2B2C2,A3B3C3,A4B4C4。

7. 3s软件使用

7.1. 3s软件安装

1
2
3
4
wget http://abacus.gene.ucl.ac.uk/software/3s.v3.1.tgz # 下载软件，只有6.74MB
tar -zxf 3s.v3.1.tgz #解压缩和解包
cd 3s.v3.1 #进入文件夹
gcc -O3 -o 3s 3s.c tools.c lfun3s.c -lm #编译

7.2. 输入文件

  1. 三个样本的碱基序列文件
  • phylip格式
  • 比对好的(aligned)
  1. 3s.ctl：参数文件基本使用默认即可，不需要变动。
1
2
3
4
5
6
7
8
9
10
11
seed = -1 #随机数种子，-1使用时间种子
outfile = path/to/out #存储运行结果的路径
seqfile = sample_align.phylip #比对好的phylip格式的序列
Imapfile = mapfile.txt #map文件，即下面准备的输入文件；如果没有map文件，默认出现的前3个物种为物种1，2，3。
ratefile = ChenLi3s.Rate.txt # 这行作者建议删除，因为使用的JC69模型，无法改变不同碱基的替换率。
nloci = 1000 # 基因位点的数量。如果小于实际位点数，只计算前1000个；如果大于实际位点数，则报错。
npoints = 16 # 做数值积分的格点数，越大精度越高，但速度慢很多。用8，16，32；作者建议16，或者尝试32。
getSE = 0 # 通常设置为0，不算方差；如果设置成1，计算方差，运算速度大幅下降，但总体上来说都很快。算不算都行。
Small_Diff = 0.5e-9 #最大似然估计的时候，似然值到多小就停止计算，建议0.5e-9，在0.5e-8到1e-9范围内调整。
simmodel = 0 #0代表不用对称模型（symmetric models），1代表使用对称模型。对称模型是指所有群体大小参数都一样，一般不用。
models = 0 2 3 #计算的模型，3代表introgression模型。
  1. mapfile
  • 用于指定三个物种分别是1，2，3的文本文件。
  • 每行两列，空格隔开。第一列物种名，第二列1，2，3。保存mapfile.txt
1
2
3
A 1
B 2
C 3

7.3. 3s软件运行

运行命令3s 3s.ctl即可

7.4. 运行结果

  1. 结果示例

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    ***Model 0 (M0)***
    lnL = -342929.817009
    MLEs
        theta4 theta5 tau0 tau1 theta1 theta2 theta3
        0.002213 0.000337 0.000866 0.000431 0.001156 0.003993 0.000502


    ***Model 2***
    lnL = -342888.249022
    2DlnL = +83.135973
    MLEs
        theta4 theta5 tau0 tau1 theta1 theta2 theta3 M12
        0.002248 0.000076 0.000853 0.000759 0.000987 0.001408 0.000500 0.934564

  2. 结果的解释

  • Model 0是没有基因流的情况下估计得到的参数，theta是θ，tau是τ。
  • Model 2是存在基因流的情况下估计得到的参数，M12代表的是从物种1到物种2的基因流，强度是0.934564，代表平均每一代有0.934564个个体从物种1迁移到物种2，是很强的基因流。

7.5. 结果有效性

7.5.1. 如何保证收敛

  1. 同一数据，运行多次保证结果一致，且2dlnL>0。当多次结果不一致时，可以考虑调小small_diff，不行再调大。
  2. 遇到收敛困难时，可以先在M0（即无基因流的模式）下分析，把M0的MLE作为其他模型的初值。
  3. 如果估计值触碰到上下界且明显不合理，需多次分析；可以尝试通过修改in.3s文件将初值设置在更合理的地方。
  4. 将不可能出现迁移的方向迁移率设置为0，减少参数个数。

7.5.2. 如何选择正确模型

  1. 如估计的树接近星状树，需调整物种树再次分析。这个调整物种树就是调整mapfile.txt文件，使得不同物种成为外类群。如果调整成任意物种树都是星状树，则可以作为辐射进化（即τ0,τ1时间很短）的证据。
  2. 如果迁移率Mij触碰下界可将该迁移率固定为0后重新分析，如两次分析-lnL相差在1.92内，则认为参数少的模型更加准确。
  3. 多次分析后找出最优模型

8. reference

这篇笔记主要是参考“朱天琪老师在中国科学院微生物研究所真菌学国家重点实验室主办的2020种群遗传学与基因组学高级培训班(第二期)和2022组学驱动的进化生态高级培训班的讲座和PPT”。

  • 欢迎关注微信公众号：生信技工
  • 公众号主要分享生信分析、生信软件、基因组学、转录组学、植物进化、生物学概念等相关内容，包括生物信息学工具的基本原理、操作步骤和学习心得。
真诚赞赏，手留余香