This is a test version of Biostars. For the public version, visit https://www.biostars.org.
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)
vcf pl

Hello,

could you please give an example how your input looks like and what output you expected?

fin swimmer

0 answers

No answers yet.

Log in to answer this question.