基因流分析流程

基因流计算文档说明一．基因流大小估算群体结构模式是海岛结构模式，基因流在群体间是随机的、均一的而且群体也要达到漂变-迁移平衡，迁移个体来自所有其他群体中随机的一个遗传变异群体。该模型的基本思想是假设１个群体分化为无限多个亚群体，亚群体在空间呈离散分布，每个亚群体接受一小部分来自整个群体的迁移个体。迁移率与迁移基因频率在任一世代内假设为常数。岛屿模型中的Nm为每代迁入的有效个体数，即基因流的估计值。其计算公式为：Nm=1−Fst4FstWright,S.Thegeneticalstructureofpopulations.AnnEugen,1951.15(4):p.323-54.二．长期基因流流向计算采用MIGRATE-N软件中Bayesianinference的策略估算群体间基因流流向，对于群体A和群体B之间基因流流向来说，存在四种可能的模型：1）基因流在A和B之间双向流动；2）基因流从A流向B；3）基因流从B流向A；4）两者为同一群体。该软件通过分别计算前三种模型的marginallikelihood值，再转换为BayesFactors进行比较来确定各模型的可能性。（1）软件运行命令/share/nas1/zhougang/soybeans/migrate/migrate-3.6.4/bin/bin/migrate-nparmfileparmfile为参数文件，具体设置详见（3）。（2）输入文件格式MIGRATE-N软件支持多种输入格式，但此类软件运行时间均非常漫长，在利用检测到的全基因组SNP作基因流流向估计时，为减少运行时间，我们通常把所有位点SNP连成DNA序列来进行运算，在此我们仅介绍DNA序列模式的输入文件格式：第一行：分为4列数据信息，第一列留空，第二列为群体数目，第三列为位点数目，第四列为数据注释；第二行：列数与位点数目等同，每列为该位点所含碱基数；第三行：群体1信息，列出群体内个体数及群体名字；第四行起：每个位点该群体内个体信息。（3）运行参数主要为parmfile参数文件的设置，该设置可以在命令行中修改，也可以直接采用文本编辑器直接编辑。①设置输入输出文件路径上图中1和7分别为输入输出路径设置选项，也可在parmfile文件中直接对infile和outfile行进行更改。②选择需要分析的模型：图中选项8即为模型选择项，三种不同模型选择如下图示意，*表示计算，0表示忽略，因此选项值****表示模型一；*0**表示模型二，**0*表示模型三，而*00*则不进行任何计算，该选项也可在parmfile文件中直接对custom-migration行进行更改。③运行参数选择主要设置选项：10，更改迭代次数，默认5000，次数越大越好，但运行时间也相应增加，建议最少10,000次以上；12，重复计算，默认关闭，建议设置3-5次为佳；13，heating选项，默认关闭，建议打开，直接采用其打开后的heating的参数。以上参数也可在parmfile文件long-sample、replicate及heating行中修改，如：long-sample=100000，replicate=YES:3，heating=YES:0:{1.000000,1.500000,3.000000,1000000.000000}。注意：在命令行窗口中，所有参数设置完成后，必须输入W生成parmfile文件才会保存所有参数，在parmfile文件中修改则直接保存即可。（4）结果输出每个模型分别计算后，均可输出pdf结果文件，也可直接用文本编辑器打开同名文件提取结果，主要为Thermodynamicintegration值和Harmonicmean值，Thermodynamicintegration值有两个，分为1a和1b，简单来说1a,1b（Thermodynamicintegration）就是采用heating参数后的结果，而2（Harmonicmean）是没有采用heating参数的结果；1a,1b之间的区别是作者认为1b采用了贝赛尔曲线从而获得了更好地结果。要统计各个模型的可能性，要转换为BayesFactors进行比较，我们现在结果中的数值是marginallikelihood值，BF=Exp[log(P(D|thisModel)-log(P(D|otherModel)]。举例来说，我有一组1b的结果三个模型的marginallikelihood值分别为-3，-2，-1，那么我们先取常数e-3，e-2和e-1的值，然后这三个值相加每个值与和的比值就是每个模型的可能性。通常来说，因为我们输入的SNP序列会较长，计算得到的值负数太小，要是取e的负几十万次方根本无法计算，所以不好算这个P值，但是可以简单的知道互相之间的可能性大小，值越大，概率越大，即e-1e-2e-3。（5）结果解释因为我们采用的算法是基于coalescentmodel，这个模型得到的结果，跟平时的解释是相反的，因为coalescentmodel在时间上是从现在到过去，而正常情况下，我们解释都从过去到现在，也就是说，从现在到过去，是有popB流向popA，那么从过去到现在就是从popA流向popB。Beerli,P.andM.Palczewski.Unifiedframeworktoevaluatepanmixiaandmigrationdirectionamongmultiplesamplinglocations.Genetics,2010.185(1):p.313-26.三．近期基因流计算BayesAss软件被用来计算近期基因流，主要流程参数如下：（1）软件运行命令/share/nas1/zhougang/soybeans/BA3-3.0.3/BA3-v-s5339-i50000000-b5000000-n500-oresult/NDSD_nonFD.1.outexamples/NDSD_nonFD.haplotype.1-a0.55-m0.15-f0.06（2）输入文件格式第一列为个体名称，第二列为群体名称，第三列为位点名称，第四第五列为双倍体基因组基因型。（3）软件参数选择以下参数供参数：-v屏幕输出详细计算过程；-s随机种子初始数，可为任意整数，在重复运算时，可设置不同的种子数来验证结果；-i迭代次数，设置越大越好，但同样会延长计算时间；-bburning数，丢弃开始部分的迭代计算结果；-n选择迭代结果，每隔n个选择一个用于后续参数计算；-o输出结果路径-a,-m,-f混合模型参数，分别表示等位基因频率、近交系数和基因流频率，数值为0到1之间，此三个参数依照不同数据集需特异调整，如：上图为运算过程中间结果输出，%accepted后括号中第1、3、4列数值即由此三个参数分别控制，建议这三列数值为0.2-0.6之间为佳。（4）结果输出结果输出文件可用文本编辑器打开，示例如下：MigrationRates下面m[i][j]即为每一代群体i从群体j中迁移过来的个体比率。RannalaB,YangZ(2003)BayesestimationofspeciesdivergencetimesandancestralpopulationsizesusingDNAsequencesfrommultipleloci.Genetics164:1645–1656.