cyvcf2 icon indicating copy to clipboard operation
cyvcf2 copied to clipboard

Published 20 hours ago • verified publisher icon brentp

raw_header is a bit cooked

Open davmlaw opened this issue 11 months ago • 4 comments

raw_header is modified - loading a VCF with bad contigs in header, it has fake contigs with no length added to it.

This always happens on 1 machine, not on another 2 I tried. I suspect it's due to htslib doing some magic and adding fake contigs to the header. But for that to happen, raw_header must have been built from data structures or something rather than just being raw.

This VCF file has contigs given as "chr1" in the header, but the records are "1" - so the used contigs are not declared

The header looks like:

In [1]: import gzip
In [2]: from cyvcf2 import Reader
In [3]: filename = "/data/incoming/clinvar_20240624_GCA_009914755.4.vcf.gz"
In [4]: for line in gzip.open(filename, "rt"):
   ...:     if line.startswith("#"):
   ...:         print(line.strip())
   ...:     else:
   ...:         break
   ...: 
##fileformat=VCFv4.1
##INFO=<ID=AF_ESP,Number=1,Type=Float,Description="allele frequencies from GO-ESP">
##INFO=<ID=AF_EXAC,Number=1,Type=Float,Description="allele frequencies from ExAC">
##INFO=<ID=AF_TGP,Number=1,Type=Float,Description="allele frequencies from TGP">
##INFO=<ID=ALLELEID,Number=1,Type=Integer,Description="the ClinVar Allele ID">
##INFO=<ID=CLNDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
##INFO=<ID=CLNDNINCL,Number=.,Type=String,Description="For included Variant : ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
##INFO=<ID=CLNDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier submitted for germline classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=CLNDISDBINCL,Number=.,Type=String,Description="For included Variant: Tag-value pairs of disease database name and identifier for germline classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=CLNHGVS,Number=.,Type=String,Description="Top-level (primary assembly, alt, or patch) HGVS expression.">
##INFO=<ID=CLNREVSTAT,Number=.,Type=String,Description="ClinVar review status of germline classification for the Variation ID">
##INFO=<ID=CLNSIG,Number=.,Type=String,Description="Aggregate germline classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=CLNSIGCONF,Number=.,Type=String,Description="Conflicting germline classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=CLNSIGINCL,Number=.,Type=String,Description="Germline classification for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:classification; multiple values are separated by a vertical bar">
##INFO=<ID=CLNVC,Number=1,Type=String,Description="Variant type">
##INFO=<ID=CLNVCSO,Number=1,Type=String,Description="Sequence Ontology id for variant type">
##INFO=<ID=CLNVI,Number=.,Type=String,Description="the variant's clinical sources reported as tag-value pairs of database and variant identifier">
##INFO=<ID=DBVARID,Number=.,Type=String,Description="nsv accessions from dbVar for the variant">
##INFO=<ID=GENEINFO,Number=1,Type=String,Description="Gene(s) for the variant reported as gene symbol:gene id. The gene symbol and id are delimited by a colon (:) and each pair is delimited by a vertical bar (|)">
##INFO=<ID=MC,Number=.,Type=String,Description="comma separated list of molecular consequence in the form of Sequence Ontology ID|molecular_consequence">
##INFO=<ID=ONCDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in ONCDISDB">
##INFO=<ID=ONCDNINCL,Number=.,Type=String,Description="For included variant: ClinVar's preferred disease name for the concept specified by disease identifiers in ONCDISDBINCL">
##INFO=<ID=ONCDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier submitted for oncogenicity classifications, e.g. MedGen:NNNNNN">
##INFO=<ID=ONCDISDBINCL,Number=.,Type=String,Description="For included variant: Tag-value pairs of disease database name and identifier for oncogenicity classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=ONC,Number=.,Type=String,Description="Aggregate oncogenicity classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=ONCINCL,Number=.,Type=String,Description="Oncogenicity classification for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:classification; multiple values are separated by a vertical bar">
##INFO=<ID=ONCREVSTAT,Number=.,Type=String,Description="ClinVar review status of oncogenicity classification for the Variation ID">
##INFO=<ID=ONCCONF,Number=.,Type=String,Description="Conflicting oncogenicity classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=ORIGIN,Number=.,Type=String,Description="Allele origin. One or more of the following values may be added: 0 - unknown; 1 - germline; 2 - somatic; 4 - inherited; 8 - paternal; 16 - maternal; 32 - de-novo; 64 - biparental; 128 - uniparental; 256 - not-tested; 512 - tested-inconclusive; 1073741824 - other">
##INFO=<ID=RS,Number=.,Type=String,Description="dbSNP ID (i.e. rs number)">
##INFO=<ID=SCIDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in SCIDISDB">
##INFO=<ID=SCIDNINCL,Number=.,Type=String,Description="For included variant: ClinVar's preferred disease name for the concept specified by disease identifiers in SCIDISDBINCL">
##INFO=<ID=SCIDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier submitted for somatic clinial impact classifications, e.g. MedGen:NNNNNN">
##INFO=<ID=SCIDISDBINCL,Number=.,Type=String,Description="For included variant: Tag-value pairs of disease database name and identifier for somatic clinical impact classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=SCIREVSTAT,Number=.,Type=String,Description="ClinVar review status of somatic clinical impact for the Variation ID">
##INFO=<ID=SCI,Number=.,Type=String,Description="Aggregate somatic clinical impact for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=SCIINCL,Number=.,Type=String,Description="Somatic clinical impact classification for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:classification; multiple values are separated by a vertical bar">
##contig=<ID=chr1,length=248387328,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr10,length=134758134,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr11,length=135127769,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr12,length=133324548,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr13,length=113566686,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr14,length=101161492,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr15,length=99753195,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr16,length=96330374,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr17,length=84276897,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr18,length=80542538,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr19,length=61707364,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr2,length=242696752,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr20,length=66210255,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr21,length=45090682,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr22,length=51324926,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr3,length=201105948,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr4,length=193574945,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr5,length=182045439,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr6,length=172126628,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr7,length=160567428,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr8,length=146259331,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr9,length=150617247,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chrX,length=154259566,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chrY,length=62460029,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##liftOverProgram=CrossMap,version=0.6.6,website=https://crossmap.readthedocs.io/en/latest/
##liftOverChainFile=chain/grch38-T2T_CHM13_v2.chain
##originalFile=clinvar_20240624.vcf.gz
##targetRefGenome=fasta/Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz
##liftOverDate=July26,2024
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO

while the htslib/cyvcf2.raw_header has extra contigs at the bottom:

In [5]: print(Reader(filename).raw_header)
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=AF_ESP,Number=1,Type=Float,Description="allele frequencies from GO-ESP">
##INFO=<ID=AF_EXAC,Number=1,Type=Float,Description="allele frequencies from ExAC">
##INFO=<ID=AF_TGP,Number=1,Type=Float,Description="allele frequencies from TGP">
##INFO=<ID=ALLELEID,Number=1,Type=Integer,Description="the ClinVar Allele ID">
##INFO=<ID=CLNDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
##INFO=<ID=CLNDNINCL,Number=.,Type=String,Description="For included Variant : ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
##INFO=<ID=CLNDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier submitted for germline classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=CLNDISDBINCL,Number=.,Type=String,Description="For included Variant: Tag-value pairs of disease database name and identifier for germline classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=CLNHGVS,Number=.,Type=String,Description="Top-level (primary assembly, alt, or patch) HGVS expression.">
##INFO=<ID=CLNREVSTAT,Number=.,Type=String,Description="ClinVar review status of germline classification for the Variation ID">
##INFO=<ID=CLNSIG,Number=.,Type=String,Description="Aggregate germline classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=CLNSIGCONF,Number=.,Type=String,Description="Conflicting germline classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=CLNSIGINCL,Number=.,Type=String,Description="Germline classification for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:classification; multiple values are separated by a vertical bar">
##INFO=<ID=CLNVC,Number=1,Type=String,Description="Variant type">
##INFO=<ID=CLNVCSO,Number=1,Type=String,Description="Sequence Ontology id for variant type">
##INFO=<ID=CLNVI,Number=.,Type=String,Description="the variant's clinical sources reported as tag-value pairs of database and variant identifier">
##INFO=<ID=DBVARID,Number=.,Type=String,Description="nsv accessions from dbVar for the variant">
##INFO=<ID=GENEINFO,Number=1,Type=String,Description="Gene(s) for the variant reported as gene symbol:gene id. The gene symbol and id are delimited by a colon (:) and each pair is delimited by a vertical bar (|)">
##INFO=<ID=MC,Number=.,Type=String,Description="comma separated list of molecular consequence in the form of Sequence Ontology ID|molecular_consequence">
##INFO=<ID=ONCDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in ONCDISDB">
##INFO=<ID=ONCDNINCL,Number=.,Type=String,Description="For included variant: ClinVar's preferred disease name for the concept specified by disease identifiers in ONCDISDBINCL">
##INFO=<ID=ONCDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier submitted for oncogenicity classifications, e.g. MedGen:NNNNNN">
##INFO=<ID=ONCDISDBINCL,Number=.,Type=String,Description="For included variant: Tag-value pairs of disease database name and identifier for oncogenicity classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=ONC,Number=.,Type=String,Description="Aggregate oncogenicity classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=ONCINCL,Number=.,Type=String,Description="Oncogenicity classification for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:classification; multiple values are separated by a vertical bar">
##INFO=<ID=ONCREVSTAT,Number=.,Type=String,Description="ClinVar review status of oncogenicity classification for the Variation ID">
##INFO=<ID=ONCCONF,Number=.,Type=String,Description="Conflicting oncogenicity classification for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=ORIGIN,Number=.,Type=String,Description="Allele origin. One or more of the following values may be added: 0 - unknown; 1 - germline; 2 - somatic; 4 - inherited; 8 - paternal; 16 - maternal; 32 - de-novo; 64 - biparental; 128 - uniparental; 256 - not-tested; 512 - tested-inconclusive; 1073741824 - other">
##INFO=<ID=RS,Number=.,Type=String,Description="dbSNP ID (i.e. rs number)">
##INFO=<ID=SCIDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in SCIDISDB">
##INFO=<ID=SCIDNINCL,Number=.,Type=String,Description="For included variant: ClinVar's preferred disease name for the concept specified by disease identifiers in SCIDISDBINCL">
##INFO=<ID=SCIDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier submitted for somatic clinial impact classifications, e.g. MedGen:NNNNNN">
##INFO=<ID=SCIDISDBINCL,Number=.,Type=String,Description="For included variant: Tag-value pairs of disease database name and identifier for somatic clinical impact classifications, e.g. OMIM:NNNNNN">
##INFO=<ID=SCIREVSTAT,Number=.,Type=String,Description="ClinVar review status of somatic clinical impact for the Variation ID">
##INFO=<ID=SCI,Number=.,Type=String,Description="Aggregate somatic clinical impact for this single variant; multiple values are separated by a vertical bar">
##INFO=<ID=SCIINCL,Number=.,Type=String,Description="Somatic clinical impact classification for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:classification; multiple values are separated by a vertical bar">
##contig=<ID=chr1,length=248387328,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr10,length=134758134,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr11,length=135127769,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr12,length=133324548,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr13,length=113566686,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr14,length=101161492,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr15,length=99753195,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr16,length=96330374,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr17,length=84276897,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr18,length=80542538,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr19,length=61707364,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr2,length=242696752,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr20,length=66210255,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr21,length=45090682,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr22,length=51324926,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr3,length=201105948,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr4,length=193574945,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr5,length=182045439,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr6,length=172126628,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr7,length=160567428,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr8,length=146259331,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chr9,length=150617247,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chrX,length=154259566,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##contig=<ID=chrY,length=62460029,assembly=Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz>
##liftOverProgram=CrossMap,version=0.6.6,website=https://crossmap.readthedocs.io/en/latest/
##liftOverChainFile=chain/grch38-T2T_CHM13_v2.chain
##originalFile=clinvar_20240624.vcf.gz
##targetRefGenome=fasta/Homo_sapiens_gca009914755v4.T2T_CHM13_v2.dna.primary_assembly.fa.gz
##liftOverDate=July26,2024
##contig=<ID=1>
##contig=<ID=10>
##contig=<ID=11>
##contig=<ID=12>
##contig=<ID=13>
##contig=<ID=14>
##contig=<ID=15>
##contig=<ID=16>
##contig=<ID=17>
##contig=<ID=18>
##contig=<ID=19>
##contig=<ID=2>
##contig=<ID=20>
##contig=<ID=21>
##contig=<ID=22>
##contig=<ID=3>
##contig=<ID=4>
##contig=<ID=5>
##contig=<ID=6>
##contig=<ID=7>
##contig=<ID=8>
##contig=<ID=9>
##contig=<ID=X>
##contig=<ID=Y>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO

Lib versions

In [12]: import cyvcf2

In [13]: cyvcf2.__version__
Out[13]: '0.31.1'

If you can tell me what to run to find out what htslib version cyvcf2 is using e.g.:

In [14]: print(pysam.version.__htslib_version__)
    ...: 
1.18

davmlaw avatar Dec 19 '24 13:12 davmlaw

hmm. that is odd. here is the code for raw header: https://github.com/brentp/cyvcf2/blob/541ab16a255a5287c331843d8180ed6b9ef10e00/cyvcf2/cyvcf2.pyx#L641-L650

brentp avatar Dec 19 '24 16:12 brentp

is __set__ somehow getting called? https://github.com/brentp/cyvcf2/blob/541ab16a255a5287c331843d8180ed6b9ef10e00/cyvcf2/cyvcf2.pyx#L1987-L1995

brentp avatar Dec 19 '24 16:12 brentp

Hi, sorry for slow reply it's summer here and I was at the beach

I haven't been able to recreate it using:

print(Reader(filename).raw_header)

Though I am pretty sure I did just copy/paste at the time. This was on a VM I have since deleted, sorry.

I can recreate it on any machine by doing:

reader = Reader(filename)
next(reader)
print(reader.raw_header)

In which case because it did call set when it created the VCF record using a chrom that isn't in the header, so the header is now:

...
##liftOverDate=July26,2024
##contig=<ID=1>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO

Maybe I didn't notice I looped on reader before, though?

So I guess I am pointing out the behavior here - I feel that "raw" implies that it is straight out of the file, and not silently modified, even if trying to help. I can understand if you need to keep backwards compatibility though

davmlaw avatar Jan 28 '25 06:01 davmlaw

yeah, the contigs only get added when they are not found on a record. I think the solution is to document that in order to get the unaltered header, you must access raw_header before iteration (or have a VCF that already has the ##contig rows).

brentp avatar Jan 28 '25 17:01 brentp