easySFS计算sfs
安装easySFS适用的conda环境
1
2
3
4
5
6
7
8
9
10
11
| conda create -n easySFS
conda activate easySFS
#Install dependencies:
conda install -c bioconda dadi pandas
#Clone this repo:
git clone https://github.com/isaacovercast/easySFS.git
cd easySFS
chmod 777 easySFS.py
./easySFS.py
|
试运行，并判断是否需要down sample
这一步是为missing data很多的情况准备的。easySFS对missing data的考虑并不如angsd那么周到，所以一般需要进行这一步，根据输出的文件判断是否要缩小后续采样数目：
1
2
3
4
| /path/to/easySFS.py \
-i /path/to/file.vcf \
-p pops_file.txt \
--preview
|
根据官网的解释，这里生成的文件里，每个“（）”里面前面的数字代表采样投影数目（即projection），后面的代表这个数目下的frequency，要在考虑到样本数目的情况下选择frequency最大的projection数目。
之前查到一个说法，说因为研究的对象是二倍体，所以project数目要全部乘二，但是【后续发现其实不需要乘二，因为是“down” sample，乘二会超出最高样本量限制，导致算出的sfs里所有数值都是0！】
再次运行并生成sfs文件
上一步得到的projection数目会作为这一步里的一个参数：
1
2
3
4
| /path/to/easySFS.py \
-i /path/to/file.vcf \
-p pops_file.txt \
--proj 27,22,13,13 #上一步得出的最佳projection数目，按顺序排列
|
这一步会自动生成output文件夹，里面按照具体需求分成了dadi和fastsimcoal会用到的数据。fastsimcoal2需要的文件是obs和sfs后缀的那些。
fastsimcoal2种群历史模拟
基本原理和思路 ▲
fastsimcoal2的原理其实就是从现在到过去的溯祖过程，所以一切都需要反过来思考，比如瓶颈效应按照进化的时间顺序应该是从较大的种群规模变得更小，Ne是【减小】的，但在溯祖模型里就是从正在经历的瓶颈，即更小的种群规模回到了更大的规模，所以Ne是【增大】的。
同理，种群的分化/split在这里就应该理解为“两个原先分开的种群合并为了一个”，再精确一些，应用在实际模型运行中，就是“两个原本分开的种群，其中一个种群中的全部个体都迁移到了另一个当中”，因此可以用migration参数为1的迁移事件表示一次历史上的种群split。
这里的思路转换非常重要，因为所有模型构建都必须按照这个思路来进行。简单来讲，就是【从现在回到过去】模拟种群的历史。
输入文件
obs文件
上一步运行easySFS得到的结果文件。如果使用其他软件比如angsd计算sfs，需要自己手动写很多obs后缀的文件，也就是写一行类似这样的开头：
1
2
| 1 observation
d0_0 d0_1 d0_2 d0_3 d0_4 d0_5 d0_6 d0_7 d0_8 d0_9 d0_10 d0_11 d0_12 d0_13 d0_14 d0_15 d0_16 d0_17 d0_18 d0_19 d0_20 d0_21 d0_22 d0_23 d0_24 d0_25 d0_26 d0_27
|
然后把folded sfs的值（包含后半部分全是0的那些）复制到下一行就行了。你有几个种群，就需要写多少个两两比对的obs文件。为了便于这一步的计算，还是直接用easySFS算sfs比较好。
tpl文件
也就是template模板的意思。这个文件标明了一些很基础的参数，比如是否有基因流/瓶颈效应/种群扩张，有多少次，各发生在哪个时间段里。参数解释：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
| //Number of population samples (demes) #这里的demes数目是指从现在往回溯祖之前，【现今】总共有多少个种群
4
//Population effective sizes (number of genes) #因为这里想要计算的未知参数之一就是现今的Ne，所以在这个文件里都会用字母表示这类需要软件计算的参数，相当于解方程时设置x和y的意思
Nwfp
Nckp
Nhyb
Ntw
//Sample sizes #因为研究物种是二倍体，所以这里是每个种群个体数×2
54
44
26
26
//Growth rates : negative growth implies population expansion #我没有用到这个参数
0
0
0
0
//Number of migration matrices : 0 implies no migration between demes #这里没有用到基因流矩阵，所以没设置
0
//historical event: time, source, sink, migrants, new size, new growth rate, migr. matrix
#这七个参数分别的意思是：历史事件发生的时间，起源于哪个种群，终止于哪个种群，有多少比例的迁移/基因流，该事件结束后终止种群的新规模是之前的多少倍，该事件结束后终止种群的新增长率，该事件发生时使用哪个基因流矩阵
7 historical event
T2BOT 2 2 0 RE2BOT 0 0
T3BOT 3 3 0 RE3BOT 0 0
T0BOT 0 0 0 RE0BOT 0 0
T1BOT 1 1 0 RE1BOT 0 0
Thyb 2 0 1 REhyb 0 0
Ttw 3 1 1 REtw 0 0
TANC 0 1 1 REANC 0 0
//Number of independent loci [chromosome] #这里我没改
1 0
//Per chromosome: Number of linkage blocks #这里也没改
1
//per Block: data type, num loci, rec. rate and mut rate + optional parameters #如果输入文件用的是empirical的sfs数据就用FREQ这个类型表示，突变率我设置了和stairway plot 2一样的数值，最后那个参数也没改
FREQ 1 0 8.1e-10 OUTEXP
|
*关于这里的migration matrix：基因流矩阵，就是定义基因流方向的一个参数。我的个人理解如下表（以一个四种群模型为例）：
| sink\source |
A |
B |
C |
D |
| A |
/ |
B→A |
C→A |
D→A |
| B |
A→B |
/ |
C→B |
D→B |
| C |
A→C |
B→C |
/ |
D→C |
| D |
A→D |
B→D |
C→D |
/ |
按照这个矩阵，如果是一个只包含从A到C基因流的矩阵，写在tpl文件里就是：
1
2
3
4
| 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000
migAtoC 0.000 0.000 0.000
0.000 0.000 0.000 0.000
|
如果同时包含A和C之间的双向基因流，又包含从B到D的单向基因流，那就是：
1
2
3
4
| 0.000 0.000 migCtoA 0.000
0.000 0.000 0.000 0.000
migAtoC 0.000 0.000 0.000
0.000 migBtoD 0.000 0.000
|
不包含任何基因流的矩阵就是：
1
2
3
4
| 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000
|
est文件
向软件解释每个参数的估算方式和范围的文件。从fsc2.7开始往后的版本里【不会再需要“RULES”这一栏】了，所以这个文件只会有两栏，“PARAMETERS”和“COMPLEX PARAMETERS”，其中前者表示简单的、没有彼此包含关系的参数，后者表示具有依赖和包含关系的参数，可以理解为后者是前者的因变量。示例：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
| // Priors and rules file
// *********************
[PARAMETERS]
#这里的七个参数分别代表：是否为整数（1代表是，0代表否，我一般估算种群大小和时间的时候就设为1，如果估算的是增长率/resize比例这种就设为0），
#参数的名字，这个unif具体含义我没查到，最低值，最高值，是否输出结果，是否对该参数限制范围
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 Nwfp unif 100 100000 output bounded
1 Nckp unif 1000 100000 output bounded
1 Nhyb unif 100 100000 output bounded
1 Ntw unif 100 100000 output bounded
1 N2preBOT unif 100 150000 output bounded
1 N3preBOT unif 100 150000 output bounded
1 N0preBOT unif 100 150000 output bounded
1 N1preBOT unif 100 600000 output bounded
1 NAwfp unif 100 100000 output bounded
1 NAkp unif 1000 100000 output bounded
1 NANC unif 1000 100000 output bounded
1 Ttw unif 10 600000 output bounded
1 Thyb unif 10 600000 output bounded
1 TANC unif 10 6000000 output bounded
1 T2BOT unif 10 600000 output bounded
1 T3BOT unif 10 600000 output bounded
1 T0BOT unif 10 600000 output bounded
1 T1BOT unif 10 600000 output bounded
[RULES] #fsc2.7之后的版本都不再需要这个了
[COMPLEX PARAMETERS] #这里定义的都是需要依赖以上简单参数才能计算的那些参数，你需要在这里给定几个参数间的相关关系
#有一点非常值得注意的是，这些参数的名字不能互为子集关系，比如两个参数的名字不能是“Nhyb”和“NhybBOT”，因为前者是后者的子集，所以软件会报错
0 REhyb = NAwfp/Nwfp hide
0 REtw = NAkp/Nckp hide
0 REANC = NANC/NAkp hide
0 RE2BOT = N2preBOT/Nhyb hide
0 RE3BOT = N3preBOT/Ntw hide
0 RE0BOT = N0preBOT/Nwfp hide
0 RE1BOT = N1preBOT/Nckp hide
|