cyvcf2

To add new information to the INFO field of a VCF, you must first update the VCF header using add_info_to_header. The header update requires a dictionary containing the keys 'ID', 'Description', 'Type', and 'Number'. After updating the header, use a Writer initialized with the original VCF as a template to write the modified records to a new file.

from cyvcf2 import VCF, Writer
vcf = VCF(VCF_PATH)

# adjust the header to contain the new field
# the keys 'ID', 'Description', 'Type', and 'Number' are required.
vcf.add_info_to_header({'ID': 'gene', 'Description': 'overlapping gene',
    'Type':'Character', 'Number': '1'})

# create a new vcf Writer using the input vcf as a template.
fname = "out.vcf"
w = Writer(fname, vcf)

for v in vcf:
    # Perform user-defined logic to find new info
    genes = get_gene_intersections(v)
    if genes is not None:
        v.INFO["gene"] = ",".join(genes)
    w.write_record(v)

w.close(); vcf.close()

To add new data to the FORMAT field or modify genotypes, follow these steps:

1. Update Header: Use add_format_to_header with the required keys: 'ID', 'Description', 'Type', and 'Number'.
2. Initialize Writer: Create a Writer using the input VCF as a template.
3. Set Format Data: Use v.set_format('FIELD_NAME', data_array) to add data to the FORMAT field. If the field type is Integer and Number=1, the data must be an $n imes 1$ numpy array of integers (where $n$ is the number of samples).
4. Modify Genotypes: You can overwrite genotypes by assigning values to v.genotypes[index]. For example, to set a sample to 'no calls', use [-1]*v.ploidy + [False].
5. Refresh Genotypes: After modifying individual genotype indices, you must reassign the field to itself (v.genotypes = v.genotypes) to ensure the Variant object reprocesses the changes and subsequent calls (like v.genotype.array()) reflect the updates.

Note: Modifying v.genotypes[index] directly leaves the Variant object in an inconsistent state until the reassignment step occurs.

from cyvcf2 import VCF, Writer
vcf = VCF(VCF_PATH)

# adjust the header to contain the new field
vcf.add_format_to_header({
    'ID': 'FILTER_CODE',
    'Description': 'Numeric code for filtering reason',
    'Type': 'Integer',
    'Number': '1'
})

# create a new vcf Writer using the input vcf as a template.
fname = "out.vcf"
w = Writer(fname, vcf)

for v in vcf:
    # Example logic to find samples to filter
    indicies, reasons = filter_samples(v)
    if indicies:
        # add the reasons array to the format dictionary at this locus
        v.set_format('FILTER_CODE', reasons)
        for index in indicies:
            # overwrite the genotypes of each filtered locus to be nocalls
            v.genotypes[index] = [-1]*v.ploidy + [False]
        
        # it is necessary to reassign the genotypes field
        # so that the v Variant object reprocess it
        v.genotypes = v.genotypes

    w.write_record(v)

w.close(); vcf.close()

The VCF class is the primary entry point for parsing VCF and BCF files. You can iterate over variants in a file using a standard loop.

Important Note on Numpy Arrays: Attributes like variant.gt_ref_depths return numpy arrays that are backed by underlying C data. Once the variant object goes out of scope, these arrays will contain invalid data. To persist the data, you must create a copy using np.array().

from cyvcf2 import VCF
import numpy as np

for variant in VCF('some.vcf.gz'):
    # Basic attributes
    print(variant.CHROM, variant.start, variant.end, variant.REF, variant.ALT)
    
    # Accessing sample fields as numpy arrays
    # Note: Use np.array() to copy if you need the data after the loop iteration
    gt_types = variant.gt_types
    gt_depths = np.array(variant.gt_ref_depths) 
    
    # Accessing INFO fields
    dp = variant.INFO.get('DP') # returns int
    fs = variant.INFO.get('FS') # returns float

    # Accessing FORMAT fields per sample
    dp_per_sample = variant.format('DP')

from cyvcf2 import VCF

for variant in VCF('some.vcf.gz'):
    print(variant.CHROM, variant.start, variant.end)