Skip to content

VCF preprocessor converts an explicit missing genotype (./.) to 0/0 during catalog ALT completion #240

Description

@Luker121

Hi, I am running the PharmCAT VCF Preprocessor 3.4.0. Here, it converts an explicitly missing genotype
(./.) to homozygous reference (0/0) while completing the ALT alleles
from pharmcat_positions.vcf.

This is related to Issue #218, but the important difference is that the input
genotype in this example is missing, not already 0/0. The direct
preprocessor output therefore changes the genotype call.

My environment has:

  • PharmCAT VCF Preprocessor: 3.4.0
  • bcftools: 1.22
  • htslib: 1.22.1
  • Python: 3.10.6
  • Input: ordinary single-sample VCF, not a block gVCF

None of the following options were used:

  • --missing-to-ref
  • --absent-to-ref
  • --unspecified-to-ref

Here I have some example data:

The official PharmCAT 3.4.0 position catalog contains:

chr6  18130687  rs1142345  T  C,G  .  PASS  PX=TPMT  GT  0/0

Example input vcf:

##fileformat=VCFv4.2
##contig=<ID=chr6>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS      ID REF ALT QUAL FILTER INFO FORMAT SAMPLE
chr6   18130687 .  T   C   .    PASS   .    GT     ./.

I used this command:

pharmcat_vcf_preprocessor \
  -vcf missing_gt.vcf \
  -o output \
  -refFna GRCh38_reference.fa \
  -refVcf pharmcat_positions.vcf.bgz

And my output is:

#CHROM POS      ID          REF ALT QUAL FILTER INFO    FORMAT SAMPLE
chr6   18130687 rs1142345   T   C,G .    PASS   PX=TPMT GT     0/0

The missing genotype ./. has become 0/0.

Completing the ALT representation from C to C,G is understandable, but it should not promote an missing genotype to a confident homozygous-reference genotype, am I correct?

The output genotype should remain missing:
GT=./.
This is also consistent with the PharmCAT VCF requirements, which state that reference (0/0) and missing (./.) positions should be specified explicitly and that missing positions should not be interpreted as
reference.
The behavior appears to originate from the catalog ALT completion block:
https://github.com/PharmGKB/PharmCAT/blob/v3.4.0/preprocessor/pcat/utilities.py#L1206-L1215

If upstream QC marks a site as missing, this can turn that unknown position into reference. The Named Allele Matcher can then use information that was not present in the input.

Could you therefore please confirm whether the ./. -> 0/0 conversion is intended?

Best,
Lukas

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions