This is a test version of Biostars. For the public version, visit https://www.biostars.org.
update genotype calls with PyVCF

Hi everyone, I'm parsing my vcf files with PyVCF and I would like to update my genotype calls: for every position in which the DP field is smaller than 10, genotype should be "./.". For the moment I've generated a list with the updated calls as I wish, like this:

vcf_reader = vcf.Reader(filename='myfile.vcf.gz')
new_gt = []
for record in vcf_reader:
    if record.samples[0].data[2] > 10:
        new_gt.append(record.samples[0].data[0])
    elif record.samples[0].data[2] < 10:
        new_gt.append('./.')

However, I've also tried to directly modifying the genotypes in the vcf_reader object, but I get an error message saying 'CallData' object does not support item assignment. In this link they have a similar problem in the FORMAT field. I was wondering: has someone found a way to directly modify genotype calls in PyVCF or to insert a list with the desired output? Thanks a lot in advanced!

pyvcf genotypes
• 2,178 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by biosol biosol

1 answer

a one liner using http://lindenb.github.io/jvarkit/VcfFilterJdk.html

 java -jar dist/vcffilterjdk.jar -e 'return new VariantContextBuilder(variant).genotypes(variant.getGenotypes().stream().map(G->G.isCalled() && G.hasDP() && G.getDP()<10?GenotypeBuilder.createMissing(G.getSampleName(),G.getPloidy()):G).collect(Collectors.toList())).make();' --recalc  input.vcf
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Yep! Works super fast, thanks a lot!

• 0 views • link
written by biosol biosol

Log in to answer this question.

More posts like this