美文网首页
PSMC软件分析群体历史有效群体大小步骤

PSMC软件分析群体历史有效群体大小步骤

作者: iBioinformatics | 来源:发表于2022-09-29 06:34 被阅读0次

    https://blog.csdn.net/zaprily/article/details/108764219

    现在出现了更好的算法,能用多个样本估计一个群体的历史群体大小:MSMC详情可见此链接

    @PSMC软件分析群体历史有效群体大小流程
    首先是软件下载,需要用到PSMC和samtool
    PSMC下载地址

    1、文件转换
    基因组文件格式为.bam,callsnp之后,再转换为.fq.gz格式,关于callSNP博客写的比较明确,我们在此只展示bwt之后已经生成bam的代码
    在软件readme中有代码:

    samtools mpileup -C50 -uf ref.fa aln.bam | bcftools view -c
    | vcfutils.pl vcf2fq -d 10 -D 100 | gzip > diploid.fq.gz
    1
    2
    这段代码重点需要注意两个问题,第一个是参考基因组ref.fa文件一定要用生成bam的参考文件,否则很可能在psmc.fa中跑出空文档
    第二个可能是由于samtool版本的改变,bcftools view -c需要改成bcftools call -c才能运行成功。
    运行代码

    samtools mpileup -C50 -uf pig.fa sample.recal.bam | bcftools call -c | vcfutils.pl vcf2fq -d 10 -D 100 | gzip > sample.fq.gz
    1
    2、生成psmc.fa和psmc文件
    这一步在README中也有

    utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa
    psmc -N25 -t15 -r5 -p "4+252+4+6" -o diploid.psmc diploid.psmcfa
    1
    2
    -t 的意思在文章里面应该对应着Tmax,和时间间隔有关的一个参数
    t i = 0.1 ∗ e i / n ∗ l o g ( 1 + 10 ∗ T m a x ) − 0.1 t_i=0.1
    e^{i/nlog(1+10T_{max})}-0.1
    t
    i

    =0.1∗e
    i/n∗log(1+10∗T
    max

    )
    −0.1

    i = 0 , … , n i = 0,…, ni=0,…,n n是设定的时间间隔(64),和下一步中的-p有关系
    最重要的参数-p,文中给出的解释为

    option specifies that there are 64 atomic time intervals and 28 (=1+25+1+1) free interval parameters.The first parameter spans the first 4 atomic time intervals, each of the next 25 parameters spans 2 intervals, the 27th spans 4 intervals and the last parameter spans the last 6 time intervals.

    我的理解就是64个时间间隔,在计算时需要估计28个参数,将64个时间间隔划分一下,第一个划分了4个时间间隔,在这四个时间间隔估计1个Ne的参数;25*2划分了25组2个时间间隔,每两个时间间隔估计一个和Ne有关的参数,后面的都是这个意思。大家如果有兴趣可以数一下:

    两边的参数分别是psmc -N25 -t15 -r5 -p "4+25 * 2+4+6"和psmc -N30 -t5 -r5 -p “4+30 * 2+4+6+10”
    其他的两个参数“-N25 -r5”,这个在参数里面我找了文档和原文章都找不到,但是应该不会有很大的影响,如果不放心的话也可设置几个参数组测试一下。

    3、画图
    生成psmc文件之后就是画图的步骤,

    perl utils/psmc_plot.pl -u 2e-09 -g 1 -p sample_plot sample.psmc
    1
    -u为突变率,-g为世代间隔,用perl编写的脚本需要补充一些库,可以根据错误提示补充更新库。

    我当时就遇到库没有更新的问题,(上图)应该是libtif库

    可以使用他给出的perl直接绘制多个个样本的历史群体大小,这个回答中 davidaray 的解决方法很细致。
    说实话,这个问题我看了很多回答,我觉得github里面的这段代码应该是可以的,虽然我用这个跑不出来。

    seq 100 | xargs -i echo psmc -N25 -t15 -r5 -b -p "4+252+4+6" -o round-{}.psmc split.fa | sh
    cat diploid.psmc round-
    .psmc > combined.psmc
    perl utils/psmc_plot.pl -u 2e-09 -g 1 -p sample_plot sample.psmc
    1
    2
    3
    不同物种绘制到同一张图中
    这个回答要感谢@幻影刺客的私信回复,具体大家可以验证一下,有问题欢迎一起探讨。
    下面是如何使用软件自带脚本拼接的过程

    psmc_plot.pl -u mutation_rate -g generation -Y years -M 'sample1,sample2,sample3' combine_plot combine.psmc
    1
    只需将需要出图的 *.psmc合并在一个psmc下,然后在-M参数后面按照样本顺序添加样本名(要和你合并的数量样本数一致哦) 就可以实现合并出图了

    同一物种多个样本推算PSMC
    这个解答 回答了用多个样本可以怎么算

    有一些坑
    psmc_plot.pl -u -g 参数:-u 的说明是 absolute mutation rate per nucleotide [2.5e-08]
    注意,这里填的是per generation per site的mutation,不是per year per site 的mutation,否则会导致PSMC曲线平移一个数量级。

    这个来自于简书
    做的时候也参考了一些别人的回答
    种群密度怎么算&PSMC推测种群大小历史动态

    原文链接:https://blog.csdn.net/zaprily/article/details/108764219

    相关文章

      网友评论

          本文标题:PSMC软件分析群体历史有效群体大小步骤

          本文链接:https://www.haomeiwen.com/subject/ctuvortx.html