How to extract the GT from a given sample in a multisample VCF with cyvcf2?

Question

[LGMgeo]

Here is my VCF multisamples:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO                                           FORMAT   sample1   sample2 

chr1           801950  .       T    <DEL>   2235      PASS   SVTYPE=DEL;END=838297      GT:GQ:FT       0/0:.:. 0/0:.:.

By using cyvcf2, how can I get GT(sample1) or GT(sample2)?

I tried something like that:

from cyvcf2 import VCF

vcf = VCF(vcfFile)

for variant in vcf:

      vcf.sample

      variant.gt_types

Thank you for your help

[iraun]

vcf.samples ? (I have not tried it myself, therefore I leave this as a comment).

[LGMgeo]

it's useful to get a list of all sample names ;o)

[Arup Ghosh]

Get all the genotypes with variant.genotypes. It will return a list containing the allele and phasing information for each sample. Access the sample-specific data with the respective index.

Ref: https://brentp.github.io/cyvcf2/docstrings.html#cyvcf2.cyvcf2.Variant.genotypes

Answer 1

[Dave Carlson]

You can access a sample-specific GT based on the index of the sample. For example, to print all the GTs for "sample1":

from cyvcf2 import VCF

vcf = VCF(vcfFile)

for variant in vcf:
    print(variant.gt_types[0])

Likewise, for "sample2":

from cyvcf2 import VCF

vcf = VCF(vcfFile)

for variant in vcf:
    print(variant.gt_types[1])

[LGMgeo]

Thank you all for your responses.

In fact, I have a list of 900 samples. And it would be easier for me to access the GT from the sample name rather than the sample index. But if it doesn't exist, no problem, I'll use the index.

[Arup Ghosh]

from cyvcf2 import VCF
myvcf = VCF('1kgp.vcf.gz')

my_target_list = ['HG00124', 'HG00501', 'HG00635', 'HG02024'] #list of sample names

def get_sample_gts(vcf,target_list):
    target_index = [vcf.samples.index(x) for x in target_list]
    for v in vcf:
        target_gts = [v.genotypes[i] for i in target_index]
        print(v.REF,v.ALT,v.POS,target_gts)

get_sample_gts(myvcf,my_target_list)

[LGMgeo]

thanks Arup!