scout/parse/variant/genotype.py
def parse_genotype(variant, ind, pos):
    """Get the genotype information in the proper format

    Sv specific format fields:

    ##FORMAT=<ID=DV,Number=1,Type=Integer,
    Description="Number of paired-ends that support the event">

    ##FORMAT=<ID=PE,Number=1,Type=Integer,
    Description="Number of paired-ends that support the event">

    ##FORMAT=<ID=PR,Number=.,Type=Integer,
    Description="Spanning paired-read support for the ref and alt alleles
                in the order listed">

    ##FORMAT=<ID=RC,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the SV">

    ##FORMAT=<ID=RCL,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the left control region">

    ##FORMAT=<ID=RCR,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the right control region">

    ##FORMAT=<ID=RR,Number=1,Type=Integer,
    Description="# high-quality reference junction reads">

    ##FORMAT=<ID=RV,Number=1,Type=Integer,
    Description="# high-quality variant junction reads">

    ##FORMAT=<ID=SR,Number=1,Type=Integer,
    Description="Number of split reads that support the event">

    Args:
        variant(cyvcf2.Variant)
        ind_id(dict): A dictionary with individual information
        pos(int): What position the ind has in vcf

    Returns:
        gt_call(dict)

    """
    gt_call = {}
    ind_id = ind['individual_id']

    gt_call['individual_id'] = ind_id
    gt_call['display_name'] = ind['display_name']

    # Fill the object with the relevant information:
    genotype = variant.genotypes[pos]

    ref_call = genotype[0]
    alt_call = genotype[1]

    gt_call['genotype_call'] = '/'.join([GENOTYPE_MAP[ref_call],
                                         GENOTYPE_MAP[alt_call]])

    paired_end_alt = None
    paired_end_ref = None
    split_read_alt = None
    split_read_ref = None

    # Check if PE is annotated
    # This is the number of paired end reads that supports the variant
    if 'PE' in variant.FORMAT:
        try:
            value = int(variant.format('PE')[pos])
            if not value < 0:
                paired_end_alt = value
        except ValueError as e:
            pass

    # Check if PR is annotated
    # Number of paired end reads that supports ref and alt
    if 'PR' in variant.FORMAT:
        values = variant.format('PR')[pos]
        try:
            alt_value = int(values[1])
            ref_value = int(values[0])
            if not alt_value < 0:
                paired_end_alt = alt_value
            if not ref_value < 0:
                paired_end_ref = ref_value
        except ValueError as r:
            pass

    # Check if 'SR' is annotated
    if 'SR' in variant.FORMAT:
        values = variant.format('SR')[pos]
        alt_value = 0
        ref_value = 0
        if len(values) == 1:
            alt_value = int(values[0])
        elif len(values) == 2:
            alt_value = int(values[1])
            ref_value = int(values[0])
        if not alt_value < 0:
            split_read_alt = alt_value
        if not ref_value < 0:
            split_read_ref = ref_value

    # Number of paired ends that supports the event
    if 'DV' in variant.FORMAT:
        values = variant.format('DV')[pos]
        alt_value = int(values[0])
        if not alt_value < 0:
            paired_end_alt = alt_value

    # Number of paired ends that supports the reference
    if 'DR' in variant.FORMAT:
        values = variant.format('DR')[pos]
        ref_value = int(values[0])
        if not alt_value < 0:
            paired_end_ref = ref_value

    # Number of split reads that supports the event
    if 'RV' in variant.FORMAT:
        values = variant.format('RV')[pos]
        alt_value = int(values[0])
        if not alt_value < 0:
            split_read_alt = alt_value

    # Number of split reads that supports the reference
    if 'RR' in variant.FORMAT:
        values = variant.format('RR')[pos]
        ref_value = int(values[0])
        if not ref_value < 0:
            split_read_ref = ref_value


    alt_depth = int(variant.gt_alt_depths[pos])
    if alt_depth == -1:
        if 'VD' in variant.FORMAT:
            alt_depth = int(variant.format('VD')[pos][0])

        if (paired_end_alt != None or split_read_alt != None):
            alt_depth = 0
            if paired_end_alt:
                alt_depth += paired_end_alt
            if split_read_alt:
                alt_depth += split_read_alt

    gt_call['alt_depth'] = alt_depth

    ref_depth = int(variant.gt_ref_depths[pos])
    if ref_depth == -1:
        if (paired_end_ref != None or split_read_ref != None):
            ref_depth = 0
            if paired_end_ref:
                ref_depth += paired_end_ref
            if split_read_ref:
                ref_depth += split_read_ref

    gt_call['ref_depth'] = ref_depth

    alt_frequency = float(variant.gt_alt_freqs[pos])
    if alt_frequency == -1:
        if 'AF' in variant.FORMAT:
            alt_frequency = float(variant.format('AF')[pos][0])

    read_depth = int(variant.gt_depths[pos])
    if read_depth == -1:
        # If read depth could not be parsed by cyvcf2, try to get it manually
        if 'DP' in variant.FORMAT:
            read_depth = int(variant.format('DP')[pos][0])
        elif (alt_depth != -1 or ref_depth != -1):
            read_depth = 0
            if alt_depth != -1:
                read_depth += alt_depth
            if ref_depth != -1:
                read_depth += alt_depth

    gt_call['read_depth'] = read_depth

    gt_call['alt_frequency'] = alt_frequency

    gt_call['genotype_quality'] = int(variant.gt_quals[pos])

    return gt_call
def parse_genotype(variant, ind, pos):
    """Get the genotype information in the proper format

    Sv specific format fields:

    ##FORMAT=<ID=DV,Number=1,Type=Integer,
    Description="Number of paired-ends that support the event">

    ##FORMAT=<ID=PE,Number=1,Type=Integer,
    Description="Number of paired-ends that support the event">

    ##FORMAT=<ID=PR,Number=.,Type=Integer,
    Description="Spanning paired-read support for the ref and alt alleles
                in the order listed">

    ##FORMAT=<ID=RC,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the SV">

    ##FORMAT=<ID=RCL,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the left control region">

    ##FORMAT=<ID=RCR,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the right control region">

    ##FORMAT=<ID=RR,Number=1,Type=Integer,
    Description="# high-quality reference junction reads">

    ##FORMAT=<ID=RV,Number=1,Type=Integer,
    Description="# high-quality variant junction reads">

    ##FORMAT=<ID=SR,Number=1,Type=Integer,
    Description="Number of split reads that support the event">

    Args:
        variant(cyvcf2.Variant)
        ind_id(dict): A dictionary with individual information
        pos(int): What position the ind has in vcf

    Returns:
        gt_call(dict)

    """
    gt_call = {}
    ind_id = ind['individual_id']

    gt_call['individual_id'] = ind_id
    gt_call['display_name'] = ind['display_name']

    # Fill the object with the relevant information:
    genotype = variant.genotypes[pos]

    ref_call = genotype[0]
    alt_call = genotype[1]

    gt_call['genotype_call'] = '/'.join([GENOTYPE_MAP[ref_call],
                                         GENOTYPE_MAP[alt_call]])

    paired_end_alt = None
    paired_end_ref = None
    split_read_alt = None
    split_read_ref = None

    # Check if PE is annotated
    # This is the number of paired end reads that supports the variant
    if 'PE' in variant.FORMAT:
        try:
            value = int(variant.format('PE')[pos])
            if not value < 0:
                paired_end_alt = value
        except ValueError as e:
            pass

    # Check if PR is annotated
    # Number of paired end reads that supports ref and alt
    if 'PR' in variant.FORMAT:
        values = variant.format('PR')[pos]
        try:
            alt_value = int(values[1])
            ref_value = int(values[0])
            if not alt_value < 0:
                paired_end_alt = alt_value
            if not ref_value < 0:
                paired_end_ref = ref_value
        except ValueError as r:
            pass

    # Check if 'SR' is annotated
    if 'SR' in variant.FORMAT:
        values = variant.format('SR')[pos]
        alt_value = 0
        ref_value = 0
        if len(values) == 1:
            alt_value = int(values[0])
        elif len(values) == 2:
            alt_value = int(values[1])
            ref_value = int(values[0])
        if not alt_value < 0:
            split_read_alt = alt_value
        if not ref_value < 0:
            split_read_ref = ref_value

    # Number of paired ends that supports the event
    if 'DV' in variant.FORMAT:
        values = variant.format('DV')[pos]
        alt_value = int(values[0])
        if not alt_value < 0:
            paired_end_alt = alt_value

    # Number of paired ends that supports the reference
    if 'DR' in variant.FORMAT:
        values = variant.format('DR')[pos]
        ref_value = int(values[0])
        if not alt_value < 0:
            paired_end_ref = ref_value

    # Number of split reads that supports the event
    if 'RV' in variant.FORMAT:
        values = variant.format('RV')[pos]
        alt_value = int(values[0])
        if not alt_value < 0:
            split_read_alt = alt_value

    # Number of split reads that supports the reference
    if 'RR' in variant.FORMAT:
        values = variant.format('RR')[pos]
        ref_value = int(values[0])
        if not ref_value < 0:
            split_read_ref = ref_value


    alt_depth = int(variant.gt_alt_depths[pos])
    if alt_depth == -1:
        if 'VD' in variant.FORMAT:
            alt_depth = int(variant.format('VD')[pos][0])

        if (paired_end_alt != None or split_read_alt != None):
            alt_depth = 0
            if paired_end_alt:
                alt_depth += paired_end_alt
            if split_read_alt:
                alt_depth += split_read_alt

    gt_call['alt_depth'] = alt_depth

    ref_depth = int(variant.gt_ref_depths[pos])
    if ref_depth == -1:
        if (paired_end_ref != None or split_read_ref != None):
            ref_depth = 0
            if paired_end_ref:
                ref_depth += paired_end_ref
            if split_read_ref:
                ref_depth += split_read_ref

    gt_call['ref_depth'] = ref_depth

    alt_frequency = float(variant.gt_alt_freqs[pos])
    if alt_frequency == -1:
        if 'AF' in variant.FORMAT:
            alt_frequency = float(variant.format('AF')[pos][0])

    read_depth = int(variant.gt_depths[pos])
    if read_depth == -1:
        # If read depth could not be parsed by cyvcf2, try to get it manually
        if 'DP' in variant.FORMAT:
            read_depth = int(variant.format('DP')[pos][0])
        elif (alt_depth != -1 or ref_depth != -1):
            read_depth = 0
            if alt_depth != -1:
                read_depth += alt_depth
            if ref_depth != -1:
                read_depth += alt_depth

    gt_call['read_depth'] = read_depth

    gt_call['alt_frequency'] = alt_frequency

    gt_call['genotype_quality'] = int(variant.gt_quals[pos])

    return gt_call
def parse_genotype ( variant , ind , pos ) : gt_call = { } ind_id = ind [ 'individual_id' ] gt_call [ 'individual_id' ] = ind_id gt_call [ 'display_name' ] = ind [ 'display_name' ] # Fill the object with the relevant information: genotype = variant . genotypes [ pos ] ref_call = genotype [ 0 ] alt_call = genotype [ 1 ] gt_call [ 'genotype_call' ] = '/' . join ( [ GENOTYPE_MAP [ ref_call ] , GENOTYPE_MAP [ alt_call ] ] ) paired_end_alt = None paired_end_ref = None split_read_alt = None split_read_ref = None # Check if PE is annotated # This is the number of paired end reads that supports the variant if 'PE' in variant . FORMAT : try : value = int ( variant . format ( 'PE' ) [ pos ] ) if not value < 0 : paired_end_alt = value except ValueError as e : pass # Check if PR is annotated # Number of paired end reads that supports ref and alt if 'PR' in variant . FORMAT : values = variant . format ( 'PR' ) [ pos ] try : alt_value = int ( values [ 1 ] ) ref_value = int ( values [ 0 ] ) if not alt_value < 0 : paired_end_alt = alt_value if not ref_value < 0 : paired_end_ref = ref_value except ValueError as r : pass # Check if 'SR' is annotated if 'SR' in variant . FORMAT : values = variant . format ( 'SR' ) [ pos ] alt_value = 0 ref_value = 0 if len ( values ) == 1 : alt_value = int ( values [ 0 ] ) elif len ( values ) == 2 : alt_value = int ( values [ 1 ] ) ref_value = int ( values [ 0 ] ) if not alt_value < 0 : split_read_alt = alt_value if not ref_value < 0 : split_read_ref = ref_value # Number of paired ends that supports the event if 'DV' in variant . FORMAT : values = variant . format ( 'DV' ) [ pos ] alt_value = int ( values [ 0 ] ) if not alt_value < 0 : paired_end_alt = alt_value # Number of paired ends that supports the reference if 'DR' in variant . FORMAT : values = variant . format ( 'DR' ) [ pos ] ref_value = int ( values [ 0 ] ) if not alt_value < 0 : paired_end_ref = ref_value # Number of split reads that supports the event if 'RV' in variant . FORMAT : values = variant . format ( 'RV' ) [ pos ] alt_value = int ( values [ 0 ] ) if not alt_value < 0 : split_read_alt = alt_value # Number of split reads that supports the reference if 'RR' in variant . FORMAT : values = variant . format ( 'RR' ) [ pos ] ref_value = int ( values [ 0 ] ) if not ref_value < 0 : split_read_ref = ref_value alt_depth = int ( variant . gt_alt_depths [ pos ] ) if alt_depth == - 1 : if 'VD' in variant . FORMAT : alt_depth = int ( variant . format ( 'VD' ) [ pos ] [ 0 ] ) if ( paired_end_alt != None or split_read_alt != None ) : alt_depth = 0 if paired_end_alt : alt_depth += paired_end_alt if split_read_alt : alt_depth += split_read_alt gt_call [ 'alt_depth' ] = alt_depth ref_depth = int ( variant . gt_ref_depths [ pos ] ) if ref_depth == - 1 : if ( paired_end_ref != None or split_read_ref != None ) : ref_depth = 0 if paired_end_ref : ref_depth += paired_end_ref if split_read_ref : ref_depth += split_read_ref gt_call [ 'ref_depth' ] = ref_depth alt_frequency = float ( variant . gt_alt_freqs [ pos ] ) if alt_frequency == - 1 : if 'AF' in variant . FORMAT : alt_frequency = float ( variant . format ( 'AF' ) [ pos ] [ 0 ] ) read_depth = int ( variant . gt_depths [ pos ] ) if read_depth == - 1 : # If read depth could not be parsed by cyvcf2, try to get it manually if 'DP' in variant . FORMAT : read_depth = int ( variant . format ( 'DP' ) [ pos ] [ 0 ] ) elif ( alt_depth != - 1 or ref_depth != - 1 ) : read_depth = 0 if alt_depth != - 1 : read_depth += alt_depth if ref_depth != - 1 : read_depth += alt_depth gt_call [ 'read_depth' ] = read_depth gt_call [ 'alt_frequency' ] = alt_frequency gt_call [ 'genotype_quality' ] = int ( variant . gt_quals [ pos ] ) return gt_call
Get the genotype information in the proper format

    Sv specific format fields:

    ##FORMAT=<ID=DV,Number=1,Type=Integer,
    Description="Number of paired-ends that support the event">

    ##FORMAT=<ID=PE,Number=1,Type=Integer,
    Description="Number of paired-ends that support the event">

    ##FORMAT=<ID=PR,Number=.,Type=Integer,
    Description="Spanning paired-read support for the ref and alt alleles
                in the order listed">

    ##FORMAT=<ID=RC,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the SV">

    ##FORMAT=<ID=RCL,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the left control region">

    ##FORMAT=<ID=RCR,Number=1,Type=Integer,
    Description="Raw high-quality read counts for the right control region">

    ##FORMAT=<ID=RR,Number=1,Type=Integer,
    Description="# high-quality reference junction reads">

    ##FORMAT=<ID=RV,Number=1,Type=Integer,
    Description="# high-quality variant junction reads">

    ##FORMAT=<ID=SR,Number=1,Type=Integer,
    Description="Number of split reads that support the event">

    Args:
        variant(cyvcf2.Variant)
        ind_id(dict): A dictionary with individual information
        pos(int): What position the ind has in vcf

    Returns:
        gt_call(dict)
Get the genotype information in the proper format