526互联

GLNexus进行joint calling时的"half-calls"（如./0, ./1）问题

发布时间 2023-09-17 23:56:47作者: 生物信息与育种

目录

关于GLNexus

GLnexus是由DNAnexus开发，用于可扩展的gVCF合并和联合变异（joint calling）要求群体测序项目，GL即genotype likelihood之意。

GATK作为变异检测金标准软件，缺点在于速度很慢。尽管在分析上游的比对、单样本HaplotypeCaller检测等环节已经有很多替代品，如Sentieon、Parabricks等，但下游joint calling这一步优化仍然有限。而这正是整个GATK流程的最限速的步骤，在GATK中只能通过分区的方法来加速，效果非常有限。

GLnexus的开发解决了这个痛点问题，在速度上不说几十上百倍的提升，至少也有十多倍。对于大规模群体/队列而言（主要针对人类基因组开发），是个非常好的工具。

Deepvariant和Clara Parabricks都推荐它来做联合变异。

详见：

https://github.com/dnanexus-rnd/GLnexus

GLnexus: joint variant calling for large cohort sequencing

由于重叠变异产生的half-calls

但是GLNexus和GATK joint calling得到的最终vcf的信息略有不同。如INFO列，GLNexus中只有AF和AQ，FORMAT列则为GT:DP:AD:SB:GQ:PL:RNC。示例如下：

      1 ##fileformat=VCFv4.2
      2 ##FILTER=<ID=MONOALLELIC,Description="Site represents one ALT allele in a region with multiple variants that could not be unified into non-overlapping multi-allelic sites">
      3 ##FILTER=<ID=PASS,Description="All filters passed">
      4 ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
      5 ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
      6 ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
      7 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
      8 ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Phred-scaled genotype Likelihoods">
      9 ##FORMAT=<ID=RNC,Number=2,Type=Character,Description="Reason for No Call in GT: . = n/a, M = Missing data, P = Partial data, I = gVCF input site is non-called, D = insufficient Dept
     10 ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fishers Exact Test to detect strand bias.">
     11 ##GLnexusConfig={unifier_config: {drop_filtered: false, min_allele_copy_number: 1, min_AQ1: 70, min_AQ2: 40, min_GQ: 40, max_alleles_per_site: 32, monoallelic_sites_for_lost_alleles
     12 ##GLnexusConfigCRC32C=1926883223
     13 ##GLnexusConfigName=gatk
     14 ##GLnexusVersion=v1.4.1-0-g68e25e5
     15 ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes">
     16 ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency estimate for each alternate allele">
     17 ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
     18 ##INFO=<ID=AQ,Number=A,Type=Integer,Description="Allele Quality score reflecting evidence for each alternate allele (Phred scale)">
     19 ##bcftools_viewCommand=view -Oz -o liangke_Zm164.B73_RefGen_v3.refine.Chr1.vcf.gz liangke_Zm164.B73_RefGen_v3.refine.Chr1.bcf; Date=Wed Sep  6 05:27:17 2023
     20 ##bcftools_viewVersion=1.13+htslib-1.13+ds
     21 ##contig=<ID=Chr1,length=305476924>
     22 ##contig=<ID=Chr2,length=159038028>
     23 #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1      Sample2     Sample3     Sample4     Sample5
     24 Chr1   3440    Chr1_3440_G_A  G       A       378     .       AF=0.036585;AQ=378      GT:DP:AD:SB:GQ:PL:RNC   0/0:51:51,0:.,.,.,.:99:.:..     ./.:118:118,0:.,.,.,.:99:.:..   0/0:8
     25 Chr1   3545    Chr1_3545_T_C  T       C       700     .       AF=0.22561;AQ=700       GT:DP:AD:SB:GQ:PL:RNC   0/0:48:48,0:.,.,.,.:99:.:..     ./.:118:95,22:58,36,13,10:99:320,0,33
     26 Chr1   3560    Chr1_3560_G_A  G       A       908     .       AF=0.268293;AQ=908      GT:DP:AD:SB:GQ:PL:RNC   0/0:51:51,0:.,.,.,.:99:.:..     ./.:112:86,25:43,43,19,7:99:596,0,311
     27 Chr1   3658    Chr1_3658_A_C  A       C       3109    .       AF=0.320122;AQ=3109     GT:DP:AD:SB:GQ:PL:RNC   0/1:47:10,36:6,4,20,17:99:1287,0,194:.. 1/1:90:6,83:2,4,34,50:39:2937

因此，它不存在类似GATK硬过滤指标，只能基于它的结果和指标做过滤。一般来说，最终得到的位点会比GATK多不少。由于很多软件都是适配GATK结果而开发，因此如果直接用GLNexus时，可能会报错。比如Plink1.9处理vcf转化为plink时报错：

$plink --vcf test.vcf --recode --out snp --double-id --allow-extra-chr

Error: Line 53 of .vcf file has a GT half-call.
Use --vcf-half-call to specify how these should be processed.

原因是vcf文件中存在half-call，如./0,./1,./2等形式的基因型。Plink默认处理为error，若想不报错，需设置其他模式（如haploid、missing、reference等），每种模式的具体含义如下所示：

为什么会产生这种基因型？需要了解下joint calling的原理，我估计我很难比官方解释得更加清楚。这里有一个图，举例说明4个人的joint calling过程，看了之后基本清楚了。

简言之，GLnexus首先将尽可能多的等位基因统一到不重叠的多等位基因位点，像往常一样调用每个样本的二倍体基因型。对于其余不能很好地统一的等位基因（如上图chr22的100位置），通常是因为它们与多个其他位点重叠，GLnexus生成单等位基因位点（MONOALLELIC）来代表每个等位基因。这些单等位基因位点在FILTER字段中标记，以便于识别，并且它们的格式和解释略有不同：GT（基因型）字段仅指示ALT等位基因的调用拷贝数（./.表示零拷贝，./1表示一个拷贝，1/1表示两个拷贝），不断言 REF 等位基因或任何其他等位基因的存在/不存在，这避免了在重叠的多等位基因位点中传达的重复信息。

vcf中对于单等位基因位点的解释如下：

MONOALLELIC,Description="Site represents one ALT allele in a region with multiple variants that could not be unified into non-overlapping multi-allelic sites"

小编推测，half-calls一般是snp和indel同时在某一个位置存在才会出现，而它的存在可能是GLNexus结果多于GATK的重要原因之一。

关于GLNexus half-calls的更多内容，请参考：

GATK joint calling对于half-calls的处理

那么，在GATK的joint calling中，GenomicsDBImport或CombineGVCFs是怎么处理的呢？

GATK的vcf中，我们可能会看到某些SNP的ALT列中出现星号* 的情况，这是由于该SNP附近有indel缺失造成，也就是它对于half-calls的处理方法。

在vcf v4.3文档中，对于ALT的*号解释是：

关于GATK的处理，请参考：

建议处理

若是GLNexus的结果，小编建议对这些位点进行处理，不然下游可能很多分析会出错。

一是直接删掉。

对于群体分析，通常拿SNP来做，因此先把SNP提取出来。其次，一般这样的位点不会太多，如果不是很核心的位点（如功能位点、芯片位点等），删除即可。

二是指定替换方式。

如plink软件的处理方式，可指定half-calls的模式，如missing（小编测试，即便在小数据情况下，也可能会出现out of memory的错误）或haploid（正常）。