解读快色木炭基因流和种群历史动态可视化

参考

大型基因组SLAF-seq
蕨类植物孢子强扩散能力下生态适应塑造的遗传分化格局
Ecological adaptation shaped the genetic structure of homoploid ferns against strong dispersal capacity

https://zhuanlan.zhihu.com/p/485185926

方法：

为了理解A. sinense复合群的进化历史，我们使用了fastsimcoal v2.6（Excoffier等，2013）中的共祖模拟方法来测试不同的假说，以确定分化的顺序和时间，每对邻近谱系之间是否存在非对称基因流，以及近期有效种群大小的变化。通过遗传聚类确定的六个谱系中（见结果部分），SK2和ALT为NEA谱系分布范围边缘的地方特有种，并且在ADMIXTURE中显示出混合祖先的模式。因此，我们将它们从模型中移除以简化模型。首先使用easySFS.py（https://github.com/isaacovercast/easySFS）从中性SNPs中构建了一个多维位点频率谱（SFS）。投射中NEA、JAP、QIN、HEN谱系的基因拷贝数分别为24、10、22和18。考虑到A. sinense在第一年无法产生孢子（张，2011），以及在人工环境中A. filix-femina的十年生存率很高（Schneller & Holderegger，1997），我们假设一代时间为五年。使用r8s（Sanderson，2003）估计了每位点每代的突变率为1.7e-9，通过设置13.5至19.3 Ma范围内分类群根的年龄限制（数据未发表），包括A. filix-femina群组和其他类群的1,056个SLAF位点的比对（包含14个个体的255,745 bp）。

我们设计了总共14个场景（见图S1）来测试不同的假设。我们对每个模型运行了50次，每次包括100,000次模拟，用于计算复合似然值，并进行了40个期望-条件最大化（ECM）周期。选择具有最小赤池信息准则（AIC）值的模型作为最优模型。最终，模型13被选为最佳模型。我们选择了具有最高最大复合似然值的运行中的点估计，并通过从点估计的模拟SFS中进行100次参数化自举重复实验，获得参数估计的置信区间（CIs）。

结果文件

可视化

解读

参考运行

fastsimcoal2 是一个复杂的软件，使用它需要编写多个文件并执行多步操作。以下是对于一个包含五个亚群的基本分析流程的详细指导，包括示例代码：

1. 准备数据文件
   首先，您需要准备遗传数据文件。这通常是序列数据、SNP数据或者其他分子标记数据。数据文件应该是 fastsimcoal2 可以接受的格式，如 .arp 文件（Arlequin项目文件）。您需要转换您的原始数据到这个格式。

2. 创建模拟参数文件 (.par)
   参数文件是一个描述您的种群模型和模拟设置的文本文件。对于五个亚群，您的 .par 文件可能看起来如下：

// Parameter file for a model with five subpopulations
5  number of population samples

// Population effective sizes (Ne)
1000  pop1_Size
2000  pop2_Size
1500  pop3_Size
1200  pop4_Size
800   pop5_Size

// Sample sizes
30  pop1_SampleSize
30  pop2_SampleSize
30  pop3_SampleSize
30  pop4_SampleSize
30  pop5_SampleSize

// Growth rates: 0 for all populations (no growth)
0   pop1_GrowthRate
0   pop2_GrowthRate
0   pop3_GrowthRate
0   pop4_GrowthRate
0   pop5_GrowthRate

// Number of migration matrices : set to 1 if migration rates are constant over time
1   numberOfMigrationMatrices

// Migration matrix: rows are source populations, columns are destination populations
0      M12    M13    M14    M15
M21    0      M23    M24    M25
M31    M32    0      M34    M35
M41    M42    M43    0      M45
M51    M52    M53    M54    0

// Historical events: format is time, source population, sink population, new size of sink population, migration rate change, new migration matrix
// Add your historical events here

// Number of independent runs
10   number of independent runs

每个 Mxy 都代表从群体 x 到群体 y 的迁移率，需要根据您的数据和假设来设定这些值。

3. 执行模拟
   在命令行界面中，您将调用 fastsimcoal2 并指定您的 .par 文件，如下：

fastsimcoal2 -t my_parameter_file.par -n 100000 -m -e

这里 -t 后面跟着参数文件名，-n 后面跟着模拟次数（例如这里是100000次），-m 表示在每次模拟之后最大似然估计的执行，-e 表示模拟结束后保留每次模拟产生的文件。

4. 评估结果
   模拟完成后，fastsimcoal2 将生成多个文件，包括估计的参数值、似然度和模拟的遗传数据。您可以查看这些文件来评估不同参数的拟合度。

5. 可视化结果
   您可以使用 R 或 Python 等工具来可视化这些结果。例如，您可以绘制估计的种群大小随时间的变化，或者迁移率的变化。

由于 fastsimcoal2 的使用涉及到多种复杂的生物统计方法，具体代码和参数需要根据您的具体数据和研究目标进行调整。上述示例仅仅提供一个基础的框架。实际使用时，还需要阅读 fastsimcoal2 的官方文档以了解所有命令和参数的具体含义，并根据自己的需求进行调整。对于初学者，建议先从简单的模型开始，逐步增加。