美文网首页重点关注
二代群体遗传与重测序(中)

二代群体遗传与重测序(中)

作者: 曹草BioInfo | 来源:发表于2022-05-05 20:46 被阅读0次

    二.变异结果注释与统计

    书接上回https://www.jianshu.com/p/ab6a35502786
    如果是使用转录组重测序的话需要修改过滤条件,会发现大多数变异都发生在基因区,其他的差别都不大。比对参考基因组之后进行变异检测,同样可以用GATK的流程来做。会多一个cigar字符串的处理步骤。因为在基因组上会比对出跨区域的reads,所以要把这样一整条reads(中间是非基因区)当做两条来处理。
    三代测序由于错误率高,本身是不适合用于变异检测的,尤其是SNP和InDel。
    我们要统计的SNP变异类型包括:
    1.变异位点是否存在于基因区or上下2kbp范围内,如果在2k之内的话可能存在一些转录元件。太远的话就可以视为基因间区了
    2.如果在基因区的话,就要区分是否在intron区。在UTR区的变异影响不大,而在CDS区的话又要进一步进行分析。
    3.如果变异影响了氨基酸类型的话,就会对基因功能(蛋白)造成影响了。
    4.更加严重的变异类型有:对起始密码子的突变,突变后形成终止密码子(提前终止),干扰可变剪切位点。
    在InDel变异中,造成的非3倍数的移码突变会非常严重。即使是3倍数的移码突变也是较严重的突变。

    annovar 注释

    • 文件准备:
    基因组的位置文件gtf
    变异位点文件vcf
    gtf第八列是相位(phase)信息,意为处在CDS交界处的碱基偏移

    SNP 和 INDEL
    # 格式转换
    $ gtfToGenePred -genePredExt \
      genome.gtf \
      genome_refGene.txt#下划线后的东西不能变
    # 提取转录本(来自annovar软件)
    $ retrieve_seq_from_fasta.pl --format refGene \
      --seqfile genome.fasta \
      --out genome_refGeneMrna.fa \
      genome_refGene.txt
    

    输出得到fa文件之后进行注释。annovar可以基于基因位置进行注释,也可以基于指定区域注释。

    #格式转换
    $ convert2annovar.pl \
    -format vcf4 \ #输入文件格式
    -allsample \ #不加的话只有第一个样本的信息
    -withfreq \ #输出alt频率
    all.filtered.snp.vcf \ #输入文件,来自GATK的输出
    >all.snp.vcf.annovar.input #输出文件
    #注释
    $ annotate_variation.pl \
     -geneanno \ #基于基因进行注释
     --neargene 2000 \ #需要考虑的gene上下游范围
     -buildver genome \ #基因组数据库名称
     -dbtype refGene \ #数据库类型
     -outfile all.anno.snp \ #输出文件前缀
     -exonsort \ # 结果按exon进行排序
     all.snp.vcf.annovar.input ./#最后一项是数据库所在的目录
    
    #统计基因组各区域中snp的个数
    $ less -N all.anno.variant_function |cut -f 1|sort | uniq -c
    

    分别会得到全基因组变异文件和外显子变异文件,后者详细记载有变异类型。InDel变异只需要改变文件的名字就可以了。最后得到的输出文件中,碱基的突变会不止一个,而变异类型的种类也会更多。
    当一个变异可能属于多种变异类型的时候,annovar会按变异造成的后果严重性进行优先级排序,只输出最严重的。如果需要所有类型的输出结果可以用SnpEff软件或者加参数。

    SV与CNV

    可以继续用上一步构建好的数据库,命令几乎一样。

    $ convert2annovar.pl -format vcf4 -allsample -withfreq  all_SV.vcf >SV_annovar.input
    
    $ annotate_variation.pl -geneanno --neargene 2000 -buildver genome -dbtype refGene -outfile SV_annovar.output -exonsort SV_annovar.input ./
    

    三.群体结构分析

    包括structure,PCA,进化树分析和LD decay等

    01.snp数据过滤

    需要过滤maf次等位基因频率,missing,还有LD过滤。前二者用于去掉低质量snp位点,后者用于生成符合哈温平衡的群体。

    $ plink --vcf all.vcf \ # 输入文件
     --geno 0.1 \ # 设置缺失率阈值
     --maf 0.01 \ # 设置maf阈值
     --biallelic-only strict \ # 去掉三等位位点
     --out all.missing_maf \ # 输出文件前缀
     --recode vcf-iid \ # 输出vcf格式
     --allow-extra-chr \ # 允许其他类型的染色体
     --set-missing-var-ids @:# \ # 给没有id的snp取一个id
     --keep-allele-order #让plink不要给ref和alt换位置
    

    plink过滤maf和missing,也可以用vcftools来完成。

    $ plink --vcf all.missing_maf.vcf \
     --indep-pairwise 50 10 0.2 \ # LD过滤阈值
     --out tmp.ld \#前缀
     --allow-extra-chr \
     --set-missing-var-ids @:#
    
    $ plink --vcf all.missing_maf.vcf \
     --extract tmp.ld.prune.in \ # 保留的SNP列表
     --out all.LDfilter \
     --recode vcf-iid \#输出文件的格式
     --keep-allele-order \
     --allow-extra-chr \
     --set-missing-var-ids @:#
    

    第一步输出是可以保留下来的的位点的id信息,所以上一步过滤的时候最好可以给它加一个id。第二步再利用位点id信息来提取对应的vcf文件。
    显然,LD过滤会越来越少。如果不需要那么多snp位点去进行分析的话可以把阈值卡得更严格,从而节省计算资源。

    02.进化树构建

    • 文件准备
    vcf 文件:all.LDfilter.vcf

    $ run_pipeline.pl -Xms1G -Xmx5G \
     -SortGenotypeFilePlugin \
     -inputFile all.LDfilter.vcf \
     -outputFile all.LDfilter.sort.vcf \
     -fileType VCF
    
    $ run_pipeline.pl -Xms1G -Xmx5G \
     -importGuess all.LDfilter.sort.vcf \
     -ExportPlugin \
     -saveAs sequences.phy \
     -format Phylip_Inter
    

    phylip格式转换之前需要先排序。虽然表面上是一个pl程序,但其中有调用java的部分,所以仍然需要设置一些java参数
    dnadist可以用于nj法交互式建树,不过由于群体遗传的数据量往往非常大,脚本式利用配置文件在后台运行会更加合适。

    $ echo -e "sequences.phy\nY" > dnadist.cfg #-e打入换行符
    $ dnadist < dnadist.cfg >dnadist.log ##距离矩阵
    $ mv outfile infile.dist
    $ echo -e "infile.dist\nY" > neighbor.cfg
    $ neighbor  < neighbor.cfg >nj.log
    
    #距离矩阵
    $ less infile.dist | tr '\n' '|'| sed 's/| / /g' | tr '|' '\n' >infile.dist.table
    #去掉换行符
    $ less outtree | tr '\n' ' '|sed 's/ //g' > outtree.nwk
    

    dnadist发现outfile会报错而非覆盖,所以最好把它名字改掉。
    这样进化树文件就做好了,再改用其他软件来作图。听说ggtree包很好看。

    03.structure分析

    # vcf格式转bed格式
    $ plink --vcf ../00.filter/all.LDfilter.vcf \
    --make-bed --out all --allow-extra-chr \
    --keep-allele-order --set-missing-var-ids @:#
    # 生成k=2~4的脚本
    $ seq 2 4 | \
     awk '{print "admixture --cv -j2 all.bed "$1" 1>admix."$1".log 2>&1"}'\
     > admixture.sh
    $ cat admixture.sh 
    admixture --cv -j2 all.bed 2 1>admix.2.log 2>&1
    admixture --cv -j2 all.bed 3 1>admix.3.log 2>&1
    admixture --cv -j2 all.bed 4 1>admix.4.log 2>&1
    $ nohup sh admixture.sh &
    #绘图
    $ mkdir result
    $ cp  ./*.Q result/
    $ less outtree.nwk |sed 's/:/\n/g'|awk -F '[,(]' '{print $NF}'|awk '$1 !~")"' >sample.pop.order 
    $ Rscript draw_admixture.R result all.fam  sample.pop.order  structure
    

    在log文件里面看各K值的CV error
    Rscript第一个参数是Q矩阵;第二个参数来自格式转换时生成的fam文件,也可能用别的含样品顺序的文件;第三个参数是输出的样品顺序文件;最后是输出的前缀。可以基于树文件中给出的大概分类情况来调整顺序。
    其他的看这篇文章吧https://www.jianshu.com/p/21c2a683c6fb

    04.PCA分析

    • 文件准备:
    vcf 文件;all.LDfilter.vcf
    分组信息表:sample.pop

    plink --vcf all.LDfilter.vcf \
     --pca 10 \ # 主成分个数
     --out PCA_out \
     --allow-extra-chr --set-missing-var-ids @:#
    

    利用输出的PCA_plink.out.eigenvec和分组信息,在R中进行绘图。如果要分析三个主成分的话可以输出两张图,一张3dPCA的效果并不好。

    05.LDdecay和LDBlock

    随着两个site在基因组上的距离变远,连锁也会变弱,二者间的R²会变小。将所有SNP的LD系数按距离求平均作为纵轴,将距离作为横轴


    LD系数计算

    LD衰减距离取决于群体类型、世代间隔和染色体相对位置.
    • 文件准备
    亚群样品列表文件:
    pop.SC.table
    pop.YZR.table

    $ head -n 3 pop.SC.table
    GA0008
    GA0035
    GA0038
    
    #两个亚群
    $ PopLDdecay -InVCF all.vcf \ # 输入vcf文件
     -SubPop pop.SC.table \ # 指定要分析的亚群
     -MaxDist 500 \ # SNP对距离上限500kb
     -OutStat pop.SC.stat
    $ PopLDdecay -InVCF ../data/all.vcf  -SubPop ../data/pop.YZR.table  -MaxDist 500 -OutStat pop.YZR.stat
    #多个亚群共同绘图
    ### 准备配置文件,可以从文件名中awk出来
    $ ls pop.*.stat.gz |awk -F"." '{ print $0"\t"$2 }' > ld_stat.list
    $ Plot_MultiPop.pl -inList ld_stat.list \
     -output ld_stat.multi \
     -keepR#用于生成R脚本
    

    LDBlock用于展示染色体局部的相关性,所以需要自己先有一个关注的预选区域

    $ LDBlockShow -InVCF all.missing_maf.vcf \
    -OutPut LDBlockShow \ # 输出图片前缀
    -Region 1:300000-500000 \ # 绘图范围
    -SeleVar 2 \ # 绘制Rˆ2
    -OutPdf \ # 输出pdf图片
    -OutPng # 输出png图片
    
    LDBlock

    LDBlock展示了连锁遗传信息,一般会和曼哈顿图一起用

    四.群体选择分析

    01.π值,Fst,Tajima'D检验

    • 文件准备
    亚群样品列表文件 pop.SC.table pop.YZR.table 群体变异检测结果文件:all.vcf

    #π值计算
    $ vcftools --vcf all.vcf \
     --window-pi 100000\ #设置窗口大小
     --window-pi-step 10000\ #设置步长
     --keep pop.SC.table \ #亚群样品列表
     --out ./Pi.SC #输出文件前缀
    ###同理指定其他亚群
    $ vcftools --vcf all.vcf  --window-pi 100000 --window-pi-step 10000\  --keep pop.YZR.table  --out ./Pi.YZR
    #π值沿染色体的分布
    $ ls Pi.*pi | awk -F "." '{print $2"\t" $0}' > Pi.list
    

    接下来就可以利用Pi.list来进行CMplot绘图。染色体中间π值很低的区域往往是着丝粒区域。同样,ROD的曼哈顿图也可以用CMplot绘制。

    #TajimaD
    $ vcftools --vcf all.vcf \
    --TajimaD 100000 \
    --keep ../../data/pop.SC.table \
    --out TajimaD.SC#前缀名称
    $ vcftools --vcf all.vcf --TajimaD 100000 --keep pop.YZR.table --out TajimaD.YZR
    $ ls TajimaD.*.Tajima.D |awk -F"." '{print $2"\t"$0}' > TajimaD.list
    #Fst
    $ vcftools --vcf all.vcf \
    --fst-window-size $window \ # 窗口大小
    --fst-window-step $step \ # 设置step距离
    --weir-fst-pop pop.SC.table \
    --weir-fst-pop pop.YZR.table \
    --out ./Fst.pop1.pop2
    

    一系列的统计值计算都是用vcftools完成的,在log中有Fst的平均值和加权值的统计,可以用加权平均数来进行CMplot绘图。输出文件为.windowed.weir.fst,两个样本分化越大,这个值也会越大。
    接下来可以进行Fst和ROD的联合分析,以及Fst,D和π的联合绘图。

    02.XP-CLR

    这里使用后发表的Python版本的XPCLR软件

    $ xpclr -I ../data/all.vcf \
            -O Chr01.xpclr-python.out \
            -C Chr01 \
            -Sa  pop.SC.table \#两个亚群信息
            -Sb  pop.YZR.table  \
            --rrate  1e-8 \#单位碱基对应的遗传距离
            --ld 0.95 \
            --maxsnps 200 \
            --size 100000 \#滑窗
            --step 20000
    $ cut -f 2,3,4,12 Chr01.xpclr-python.out | awk 'NF == 4 '  > Chr01.xpclr-python.table
    #提取两个chr中xpclr的结果(第12列)
    $ xpclr -I ../data/all.vcf -O Chr02.xpclr-python.out -C Chr02 -Sa  pop.SC.table -Sb  pop.YZR.table --rrate  1e-8 --ld 0.95 --maxsnps 200 --size 100000 --step 20000
    $ cut -f 2,3,4,12 Chr02.xpclr-python.out | awk 'NF == 4 ' > Chr02.xpclr-python.table
    #如果是第一行,或者最后一列是xpclr,而且没有空值的列输出
    $ cat Chr*.xpclr-python.table | awk '{if((NR==1 || $NF != "xpclr") && NF == 4 ){print $0 }}' >  all.xpclr-python.table
    

    最后的结果也可以在R中绘图,受到选择越强的区域xpclr值会越大。

    相关文章

      网友评论

        本文标题:二代群体遗传与重测序(中)

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