群体模拟软件ms使用说明

ms是中性进化条件下,群体数据模拟最常用的软件。群体数据模拟在群体遗传学研究中非常重要,可以用来评估各统计量的可靠性,研究某种群体历史下各统计量的估计值,帮助我们理解真实数据所经历的群体历史。

基本用法:

1
ms nsam nreps -t θ

nsam为模拟的样本数量(单倍型),nreps为模拟的DNA片段数量,各片段之间不存在连锁。这两个参数为必需的。θ为突变参数,θ=4N0μL,详见下面参数解释。

输出文件如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
ms 4 2 -t 5.0
27473 36154 10290
//
segsites: 4
positions: 0.0110 0.0765 0.6557 0.7571
0010
0100
0000
1001
//
segsites: 5
positions: 0.0491 0.2443 0.2923 0.5984 0.8312
00001
00000
00010
11110

第一行为命令行。
第二行为随机数种子,用该随机数种子可以重复模拟的结果。
之后的数据以‘//’分隔各DNA片段。每个片段第一部分为分离位点数。第二部分为这些分离位点的相对位置,0~1之间,默认保留4位小数,可以用-p n 来修改精度,n为保留的小数位数。第三部分为基因型,每一行为一个单倍型,祖先型用0表示,突变用1表示。

各参数含义

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
-f filename # 参数文件,ms可以从文件读取除nsam, nreps之外的其它参数
-seeds x1 x2 x3 # 随机数种子,3个整数,一样的随机数种子可以重复结果,缺省后ms使用默认随机数种子
-t θ # 突变率,θ=4NμL, N为二倍体群体大小, μ为突变率,L为模拟的DNA片段长度。
-s int # 每个DNA片段生成固定的分离位点数,如果我们要模拟不连锁的SNPs,可以设置-s 1。
-T # 输出系统发育树,newick格式。
-L # 输出最近共同祖先的时间及总枝长,单位是4N generation。
-p int # 保留小数位数,默认是4
-r ρ nsites # ρ=4Nr,r为模拟DNA片段两末端之间的重组率,nsites为DNA片段长度。
-c f λ # 基因转换, f=g/r, g为每碱基转换率,如果r=0,则f=4Ng,λ为平均基因转换长度,即使ρ=0,nsites也要通过-c指定。
-G α # 将所有群体的增长率设为α,群体呈指数增长,α为自然数e的幂。时间的单位为4N generation。
-I npop n1 n2 ...[4Nm] # 设定一共npop个群体,每个群体的样本数分别为n1,n2 ...,群体的总基因流为4Nm,默认为0。注意,每个群体之间的基因流为4Nm/(npop-1)。
-n i x # 设置第i个群体大小为x*N。
-g i a_i # 设置第i个群体的增长率为a_i。
-m i j Mij # 设置群体j到群体i的基因流为Mij,单位为4N generation。
-ma M11, M12...M21, M22... # 设置基因流矩阵。i=j的时候,是没有任何意义的,ms会忽略,可以用x代替。
-eG t a # 设置时间t的时候,所有群体的增长率,注意ms模拟的逻辑是backward,所有的时间都是距离当前的时间,单位为4N generation。
-eg t i ai # 设置时间t的时候,i群体的增长率为ai。
-eN t x # 设置时间t的时候,所有群体的大小为x*N,注意,此设置默认效果是所有群体的增长率变为0。因此,要使同一时间的增长率生效,增长率设置要在群体大小设置之后。
-en t i x # 设置i亚群的增长率,效果同上。
-eM t x # 设置t时间所有群体之间的基因流为x/(npop-1)。
-em t i j x # 设置时间t,j 到i的基因流4Nmij=x。
-ema t npop M11, M12... # 设置t时间的基因流矩阵。
-es t i p # t时间i群体分为i和npop+1两个亚群,i中每个个体分到新i的概率为p,分到npop+1的概率为1-p。新群体的基因流和增长率都为0,群体大小设为N,i群体的基因流增长率不变。从正向历史来看,这是一个admixture事件。
-ej t i j # t时间,i群体所有个体合并到j群体,i群体到其它群体的基因流都变为0,增长率不变。从正向历史来看,这是一个群体分化事件。

注: 以上N皆指 N0因为在markdown的代码块里面,html标记不生效

例子

1. 群体大小瞬间收缩然后指数扩增



图中N1=10,000,N2=5,000,N3=20,000,如果突变率为每代10-8,如果模拟的DNA片段长度为8,000 bp,将N0设置为20,000,则θ=4*20000*10-8*8000=6.40,如果T1是16000代,则转为4N0代单位为16000/(4*20000)=0.2,T2为24000代,即0.3以4N0代为单位。5000=20000*exp-a0.2,算出a为6.93。如果该群体取15个样本,模拟1000个DNA片段,ms命令如下:

1
ms 15 1000 -t 6.4 -G 6.93 -eG 0.2 0.0 -eN 0.3 0.5

2. 两个群体分化后经历不同的群体大小变化



各个参数的计算参考1,不再赘述,只给每个模型的ms命令。

1
2
ms 15 100 -t 11.2 -I 2 3 12 -g 1 44.36 -n 2 0.125 -eg 0.03125 1
0.0 -en 0.0625 2 0.05 -ej 0.09375 2 1

3. 近期隔离的踏脚石模型



1
2
3
ms 15 100 -t 3.0 -I 6 0 7 0 0 8 0 -m 1 2 2.5 -m 2 1 2.5 -m 2 3 2.5
-m 3 2 2.5 -m 4 5 2.5 -m 5 4 2.5 -m 5 6 2.5 -m 6 5 2.5 -em 2.0 3 4
2.5 -em 2.0 4 3 2.5

4. 近期Admixture事件



1
2
3
4
5
6
7
8
9
10
11
ms 76 50000 -s 1 -I 17 24 22 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 \
-n 1 1 -n 2 0.3333 -n 3 0.25 -n 4 0.25 -n 5 0.2 -n 6 0.2 -n 7 0.2 -n 8 0.2 \
-n 9 0.2 -n 10 0.2 -n 11 0.2 -n 12 0.2 -n 13 0.2 -n 14 0.2 -n 15 0.2 \
-n 16 0.2 -n 17 0.2 \
-es 0.0003125 5 0.5 -en 0.0003125 5 0.2 -ej 0.0003125 18 2 \
-es 0.0003125 6 0.5 -en 0.0003125 6 0.2 -ej 0.0003125 19 2 \
-ej 0.01041667 6 5 -ej 0.01041667 7 5 -ej 0.01041667 8 5 -ej 0.01041667 9 5 \
-ej 0.01041667 10 5 -ej 0.01041667 11 5 -ej 0.01041667 12 5 \
-ej 0.01041667 13 5 -ej 0.01041667 14 5 -ej 0.01041667 15 5 \
-ej 0.01041667 16 5 -ej 0.01041667 17 5 -ej 0.01354167 5 4 \
-ej 0.01458333 4 3 -ej 0.015625 3 2 -ej 0.034375 2 1
如果我的博客能帮你节约时间,我会非常高兴你请我喝杯咖啡。