Get PL values according to GT in VCFs
I'm trying to parse out the PL information from a VCF according to genotype. I tried applying some of the pseduocode from the VCF format pages, but I think something is wrong.
Can someone please check my code? I'm uncertain if it's functioning properly. I tested it and thought it was working correctly. But after some analysis of the output, I think the code is wrong.
def genotype_pl_index(gt=None):
"""
In presence of the GT field the same
ploidy is expected and the canonical order is used;
without GT field, diploidy is assumed. If A is the allele in
REF and B,C,... are the alleles as ordered in ALT, the ordering
of genotypes for the likelihoods is given by:
F(j/k) = (k*(k+1)/2)+j
"""
_a = gt.replace('|','/').split('/')
if len(_a)==2:
a1,a2 = int(_a[0]),int(_a[1])
t = a1
if a2 < a1: a1,a2=a2,t
return int((a2*(a2+1)/2.)+a1)
else:
return -1
def phred_quals(entry=None,format=None):
"""
returns PL values
entry is a sample entry for a VCF (like GT:DP:PL)
format is a dict of FORMAT strings to index numbers
format['GT']=0
"""
_entry = str(entry).split(':')
if len(_entry)<format['PL']: return -1
gt = _entry[format['GT']]
else:
_pl = _entry[format['PL']].split(',')
idx = genotype_pl_index(gt)
if idx==-1: return -1
else:
pl = try_index(_pl,idx)
if pl==None: return -1
else: return float(pl)
• 124 views
•
link
0 answers
No answers yet.
Log in to answer this question.
Hello,
could you please give an example how your input looks like and what output you expected?
fin swimmer