scout/parse/variant/transcript.py
def parse_transcripts(raw_transcripts, allele=None):
    """Parse transcript information from VCF variants

    Args:
        raw_transcripts(iterable(dict)): An iterable with raw transcript 
                                         information

    Yields:
        transcript(dict) A dictionary with transcript information
    """
    for entry in raw_transcripts:
        transcript = {}
        # There can be several functional annotations for one variant
        functional_annotations = entry.get('CONSEQUENCE', '').split('&')
        transcript['functional_annotations'] = functional_annotations
        # Get the transcript id (ensembl gene id)
        transcript_id = entry.get('FEATURE', '').split(':')[0]
        transcript['transcript_id'] = transcript_id

        # Add the hgnc gene identifiers
        # The HGNC ID is prefered and will be used if it exists
        hgnc_id = entry.get('HGNC_ID')
        if hgnc_id:
            hgnc_id = hgnc_id.split(':')[-1]
            transcript['hgnc_id'] = int(hgnc_id)
        else:
            transcript['hgnc_id'] = None

        hgnc_symbol = entry.get('SYMBOL')
        if hgnc_symbol:
            transcript['hgnc_symbol'] = hgnc_symbol
        else:
            transcript['hgnc_symbol'] = None

        ########### Fill it with the available information ###########

        ### Protein specific annotations ###

        ## Protein ID ##
        transcript['protein_id'] = entry.get('ENSP')
        
        ## Polyphen prediction ##
        polyphen_prediction = entry.get('POLYPHEN')
        # Default is 'unknown'
        prediction_term = 'unknown'
        if polyphen_prediction:
            prediction_term = polyphen_prediction.split('(')[0]
        transcript['polyphen_prediction'] = prediction_term

        ## Sift prediction ##
        # Check with other key if it does not exist
        sift_prediction = entry.get('SIFT')
        # Default is 'unknown'
        prediction_term = 'unknown'
        
        if not sift_prediction:
            sift_prediction = entry.get('SIFT_PRED')
            
        if sift_prediction:
            prediction_term = sift_prediction.split('(')[0]

        transcript['sift_prediction'] = prediction_term
        
        transcript['swiss_prot'] = entry.get('SWISSPROT') or 'unknown'

        if entry.get('DOMAINS', None):
            pfam_domains = entry['DOMAINS'].split('&')

            for annotation in pfam_domains:
                annotation = annotation.split(':')
                domain_name = annotation[0]
                domain_id = annotation[1]
                if domain_name == 'Pfam_domain':
                    transcript['pfam_domain'] = domain_id
                elif domain_name == 'PROSITE_profiles':
                    transcript['prosite_profile'] = domain_id
                elif domain_name == 'SMART_domains':
                    transcript['smart_domain'] = domain_id

        coding_sequence_entry = entry.get('HGVSC', '').split(':')
        protein_sequence_entry = entry.get('HGVSP', '').split(':')

        coding_sequence_name = None
        if len(coding_sequence_entry) > 1:
            coding_sequence_name = coding_sequence_entry[-1]

        transcript['coding_sequence_name'] = coding_sequence_name

        protein_sequence_name = None
        if len(protein_sequence_entry) > 1:
            protein_sequence_name = protein_sequence_entry[-1]

        transcript['protein_sequence_name'] = protein_sequence_name

        transcript['biotype'] = entry.get('BIOTYPE')

        transcript['exon'] = entry.get('EXON')
        transcript['intron'] = entry.get('INTRON')

        if entry.get('STRAND'):
            if entry['STRAND'] == '1':
                transcript['strand'] = '+'
            elif entry['STRAND'] == '-1':
                transcript['strand'] = '-'
        else:
            transcript['strand'] = None

        functional = []
        regional = []
        for annotation in functional_annotations:
            functional.append(annotation)
            regional.append(SO_TERMS[annotation]['region'])

        transcript['functional_annotations'] = functional
        transcript['region_annotations'] = regional

        # Check if the transcript is marked cannonical by vep
        transcript['is_canonical'] = (entry.get('CANONICAL') == 'YES')
        
        # Check if the CADD score is available on transcript level
        cadd_phred = entry.get('CADD_PHRED')
        if cadd_phred:
            transcript['cadd'] = float(cadd_phred)

        # Check frequencies
        # There are different keys for different versions of VEP
        # We only support version 90+
        thousandg_freqs = []
        gnomad_freqs = []
        try:
            # The keys for VEP v90+:
            # 'AF' or '1000GAF' - 1000G all populations combined
            # 'xxx_AF' - 1000G (or NHLBI-ESP) individual populations
            # 'gnomAD_AF' - gnomAD exomes, all populations combined
            # 'gnomAD_xxx_AF' - gnomAD exomes, individual populations
            # 'MAX_AF' - Max of all populations (1000G, gnomAD exomes, ESP)
            # https://www.ensembl.org/info/docs/tools/vep/vep_formats.html
            
            # Loop over all keys to find frequency entries
            for key in entry:
                #All frequencies endswith AF
                if not key.endswith('AF'):
                    continue
                
                value = entry[key]
                if not value:
                    continue
                
                # This is the 1000G max af information
                if (key == 'AF' or key == '1000GAF'):
                    transcript['thousand_g_maf'] = float(value)
                    continue
                
                if key == 'GNOMAD_AF':
                    transcript['gnomad_maf'] = float(value)
                    continue

                if key == 'EXAC_MAX_AF':
                    transcript['exac_max'] = float(value)
                    transcript['exac_maf'] = float(value)
                    continue
                
                if 'GNOMAD' in key:
                    gnomad_freqs.append(float(value))
                
                else:
                    thousandg_freqs.append(float(value))
                
            if thousandg_freqs:
                transcript['thousandg_max'] = max(thousandg_freqs)

            if gnomad_freqs:
                transcript['gnomad_max'] = max(gnomad_freqs)
            
        except Exception as err:
            LOG.debug("Something went wrong when parsing frequencies")
            LOG.debug("Only splitted and normalised VEP v90+ is supported")

        clinsig = entry.get('CLIN_SIG')
        if clinsig:
            transcript['clinsig'] = clinsig.split('&')

        transcript['dbsnp'] = []
        transcript['cosmic'] = []
        variant_ids = entry.get('EXISTING_VARIATION')
        if variant_ids:
            for variant_id in variant_ids.split('&'):
                if variant_id.startswith('rs'):
                    transcript['dbsnp'].append(variant_id)
                elif variant_id.startswith('COSM'):
                    transcript['cosmic'].append(int(variant_id[4:]))

        yield transcript
def parse_transcripts(raw_transcripts, allele=None):
    """Parse transcript information from VCF variants

    Args:
        raw_transcripts(iterable(dict)): An iterable with raw transcript 
                                         information

    Yields:
        transcript(dict) A dictionary with transcript information
    """
    for entry in raw_transcripts:
        transcript = {}
        # There can be several functional annotations for one variant
        functional_annotations = entry.get('CONSEQUENCE', '').split('&')
        transcript['functional_annotations'] = functional_annotations
        # Get the transcript id (ensembl gene id)
        transcript_id = entry.get('FEATURE', '').split(':')[0]
        transcript['transcript_id'] = transcript_id

        # Add the hgnc gene identifiers
        # The HGNC ID is prefered and will be used if it exists
        hgnc_id = entry.get('HGNC_ID')
        if hgnc_id:
            hgnc_id = hgnc_id.split(':')[-1]
            transcript['hgnc_id'] = int(hgnc_id)
        else:
            transcript['hgnc_id'] = None

        hgnc_symbol = entry.get('SYMBOL')
        if hgnc_symbol:
            transcript['hgnc_symbol'] = hgnc_symbol
        else:
            transcript['hgnc_symbol'] = None

        ########### Fill it with the available information ###########

        ### Protein specific annotations ###

        ## Protein ID ##
        transcript['protein_id'] = entry.get('ENSP')
        
        ## Polyphen prediction ##
        polyphen_prediction = entry.get('POLYPHEN')
        # Default is 'unknown'
        prediction_term = 'unknown'
        if polyphen_prediction:
            prediction_term = polyphen_prediction.split('(')[0]
        transcript['polyphen_prediction'] = prediction_term

        ## Sift prediction ##
        # Check with other key if it does not exist
        sift_prediction = entry.get('SIFT')
        # Default is 'unknown'
        prediction_term = 'unknown'
        
        if not sift_prediction:
            sift_prediction = entry.get('SIFT_PRED')
            
        if sift_prediction:
            prediction_term = sift_prediction.split('(')[0]

        transcript['sift_prediction'] = prediction_term
        
        transcript['swiss_prot'] = entry.get('SWISSPROT') or 'unknown'

        if entry.get('DOMAINS', None):
            pfam_domains = entry['DOMAINS'].split('&')

            for annotation in pfam_domains:
                annotation = annotation.split(':')
                domain_name = annotation[0]
                domain_id = annotation[1]
                if domain_name == 'Pfam_domain':
                    transcript['pfam_domain'] = domain_id
                elif domain_name == 'PROSITE_profiles':
                    transcript['prosite_profile'] = domain_id
                elif domain_name == 'SMART_domains':
                    transcript['smart_domain'] = domain_id

        coding_sequence_entry = entry.get('HGVSC', '').split(':')
        protein_sequence_entry = entry.get('HGVSP', '').split(':')

        coding_sequence_name = None
        if len(coding_sequence_entry) > 1:
            coding_sequence_name = coding_sequence_entry[-1]

        transcript['coding_sequence_name'] = coding_sequence_name

        protein_sequence_name = None
        if len(protein_sequence_entry) > 1:
            protein_sequence_name = protein_sequence_entry[-1]

        transcript['protein_sequence_name'] = protein_sequence_name

        transcript['biotype'] = entry.get('BIOTYPE')

        transcript['exon'] = entry.get('EXON')
        transcript['intron'] = entry.get('INTRON')

        if entry.get('STRAND'):
            if entry['STRAND'] == '1':
                transcript['strand'] = '+'
            elif entry['STRAND'] == '-1':
                transcript['strand'] = '-'
        else:
            transcript['strand'] = None

        functional = []
        regional = []
        for annotation in functional_annotations:
            functional.append(annotation)
            regional.append(SO_TERMS[annotation]['region'])

        transcript['functional_annotations'] = functional
        transcript['region_annotations'] = regional

        # Check if the transcript is marked cannonical by vep
        transcript['is_canonical'] = (entry.get('CANONICAL') == 'YES')
        
        # Check if the CADD score is available on transcript level
        cadd_phred = entry.get('CADD_PHRED')
        if cadd_phred:
            transcript['cadd'] = float(cadd_phred)

        # Check frequencies
        # There are different keys for different versions of VEP
        # We only support version 90+
        thousandg_freqs = []
        gnomad_freqs = []
        try:
            # The keys for VEP v90+:
            # 'AF' or '1000GAF' - 1000G all populations combined
            # 'xxx_AF' - 1000G (or NHLBI-ESP) individual populations
            # 'gnomAD_AF' - gnomAD exomes, all populations combined
            # 'gnomAD_xxx_AF' - gnomAD exomes, individual populations
            # 'MAX_AF' - Max of all populations (1000G, gnomAD exomes, ESP)
            # https://www.ensembl.org/info/docs/tools/vep/vep_formats.html
            
            # Loop over all keys to find frequency entries
            for key in entry:
                #All frequencies endswith AF
                if not key.endswith('AF'):
                    continue
                
                value = entry[key]
                if not value:
                    continue
                
                # This is the 1000G max af information
                if (key == 'AF' or key == '1000GAF'):
                    transcript['thousand_g_maf'] = float(value)
                    continue
                
                if key == 'GNOMAD_AF':
                    transcript['gnomad_maf'] = float(value)
                    continue

                if key == 'EXAC_MAX_AF':
                    transcript['exac_max'] = float(value)
                    transcript['exac_maf'] = float(value)
                    continue
                
                if 'GNOMAD' in key:
                    gnomad_freqs.append(float(value))
                
                else:
                    thousandg_freqs.append(float(value))
                
            if thousandg_freqs:
                transcript['thousandg_max'] = max(thousandg_freqs)

            if gnomad_freqs:
                transcript['gnomad_max'] = max(gnomad_freqs)
            
        except Exception as err:
            LOG.debug("Something went wrong when parsing frequencies")
            LOG.debug("Only splitted and normalised VEP v90+ is supported")

        clinsig = entry.get('CLIN_SIG')
        if clinsig:
            transcript['clinsig'] = clinsig.split('&')

        transcript['dbsnp'] = []
        transcript['cosmic'] = []
        variant_ids = entry.get('EXISTING_VARIATION')
        if variant_ids:
            for variant_id in variant_ids.split('&'):
                if variant_id.startswith('rs'):
                    transcript['dbsnp'].append(variant_id)
                elif variant_id.startswith('COSM'):
                    transcript['cosmic'].append(int(variant_id[4:]))

        yield transcript
def parse_transcripts ( raw_transcripts , allele = None ) : for entry in raw_transcripts : transcript = { } # There can be several functional annotations for one variant functional_annotations = entry . get ( 'CONSEQUENCE' , '' ) . split ( '&' ) transcript [ 'functional_annotations' ] = functional_annotations # Get the transcript id (ensembl gene id) transcript_id = entry . get ( 'FEATURE' , '' ) . split ( ':' ) [ 0 ] transcript [ 'transcript_id' ] = transcript_id # Add the hgnc gene identifiers # The HGNC ID is prefered and will be used if it exists hgnc_id = entry . get ( 'HGNC_ID' ) if hgnc_id : hgnc_id = hgnc_id . split ( ':' ) [ - 1 ] transcript [ 'hgnc_id' ] = int ( hgnc_id ) else : transcript [ 'hgnc_id' ] = None hgnc_symbol = entry . get ( 'SYMBOL' ) if hgnc_symbol : transcript [ 'hgnc_symbol' ] = hgnc_symbol else : transcript [ 'hgnc_symbol' ] = None ########### Fill it with the available information ########### ### Protein specific annotations ### ## Protein ID ## transcript [ 'protein_id' ] = entry . get ( 'ENSP' ) ## Polyphen prediction ## polyphen_prediction = entry . get ( 'POLYPHEN' ) # Default is 'unknown' prediction_term = 'unknown' if polyphen_prediction : prediction_term = polyphen_prediction . split ( '(' ) [ 0 ] transcript [ 'polyphen_prediction' ] = prediction_term ## Sift prediction ## # Check with other key if it does not exist sift_prediction = entry . get ( 'SIFT' ) # Default is 'unknown' prediction_term = 'unknown' if not sift_prediction : sift_prediction = entry . get ( 'SIFT_PRED' ) if sift_prediction : prediction_term = sift_prediction . split ( '(' ) [ 0 ] transcript [ 'sift_prediction' ] = prediction_term transcript [ 'swiss_prot' ] = entry . get ( 'SWISSPROT' ) or 'unknown' if entry . get ( 'DOMAINS' , None ) : pfam_domains = entry [ 'DOMAINS' ] . split ( '&' ) for annotation in pfam_domains : annotation = annotation . split ( ':' ) domain_name = annotation [ 0 ] domain_id = annotation [ 1 ] if domain_name == 'Pfam_domain' : transcript [ 'pfam_domain' ] = domain_id elif domain_name == 'PROSITE_profiles' : transcript [ 'prosite_profile' ] = domain_id elif domain_name == 'SMART_domains' : transcript [ 'smart_domain' ] = domain_id coding_sequence_entry = entry . get ( 'HGVSC' , '' ) . split ( ':' ) protein_sequence_entry = entry . get ( 'HGVSP' , '' ) . split ( ':' ) coding_sequence_name = None if len ( coding_sequence_entry ) > 1 : coding_sequence_name = coding_sequence_entry [ - 1 ] transcript [ 'coding_sequence_name' ] = coding_sequence_name protein_sequence_name = None if len ( protein_sequence_entry ) > 1 : protein_sequence_name = protein_sequence_entry [ - 1 ] transcript [ 'protein_sequence_name' ] = protein_sequence_name transcript [ 'biotype' ] = entry . get ( 'BIOTYPE' ) transcript [ 'exon' ] = entry . get ( 'EXON' ) transcript [ 'intron' ] = entry . get ( 'INTRON' ) if entry . get ( 'STRAND' ) : if entry [ 'STRAND' ] == '1' : transcript [ 'strand' ] = '+' elif entry [ 'STRAND' ] == '-1' : transcript [ 'strand' ] = '-' else : transcript [ 'strand' ] = None functional = [ ] regional = [ ] for annotation in functional_annotations : functional . append ( annotation ) regional . append ( SO_TERMS [ annotation ] [ 'region' ] ) transcript [ 'functional_annotations' ] = functional transcript [ 'region_annotations' ] = regional # Check if the transcript is marked cannonical by vep transcript [ 'is_canonical' ] = ( entry . get ( 'CANONICAL' ) == 'YES' ) # Check if the CADD score is available on transcript level cadd_phred = entry . get ( 'CADD_PHRED' ) if cadd_phred : transcript [ 'cadd' ] = float ( cadd_phred ) # Check frequencies # There are different keys for different versions of VEP # We only support version 90+ thousandg_freqs = [ ] gnomad_freqs = [ ] try : # The keys for VEP v90+: # 'AF' or '1000GAF' - 1000G all populations combined # 'xxx_AF' - 1000G (or NHLBI-ESP) individual populations # 'gnomAD_AF' - gnomAD exomes, all populations combined # 'gnomAD_xxx_AF' - gnomAD exomes, individual populations # 'MAX_AF' - Max of all populations (1000G, gnomAD exomes, ESP) # https://www.ensembl.org/info/docs/tools/vep/vep_formats.html # Loop over all keys to find frequency entries for key in entry : #All frequencies endswith AF if not key . endswith ( 'AF' ) : continue value = entry [ key ] if not value : continue # This is the 1000G max af information if ( key == 'AF' or key == '1000GAF' ) : transcript [ 'thousand_g_maf' ] = float ( value ) continue if key == 'GNOMAD_AF' : transcript [ 'gnomad_maf' ] = float ( value ) continue if key == 'EXAC_MAX_AF' : transcript [ 'exac_max' ] = float ( value ) transcript [ 'exac_maf' ] = float ( value ) continue if 'GNOMAD' in key : gnomad_freqs . append ( float ( value ) ) else : thousandg_freqs . append ( float ( value ) ) if thousandg_freqs : transcript [ 'thousandg_max' ] = max ( thousandg_freqs ) if gnomad_freqs : transcript [ 'gnomad_max' ] = max ( gnomad_freqs ) except Exception as err : LOG . debug ( "Something went wrong when parsing frequencies" ) LOG . debug ( "Only splitted and normalised VEP v90+ is supported" ) clinsig = entry . get ( 'CLIN_SIG' ) if clinsig : transcript [ 'clinsig' ] = clinsig . split ( '&' ) transcript [ 'dbsnp' ] = [ ] transcript [ 'cosmic' ] = [ ] variant_ids = entry . get ( 'EXISTING_VARIATION' ) if variant_ids : for variant_id in variant_ids . split ( '&' ) : if variant_id . startswith ( 'rs' ) : transcript [ 'dbsnp' ] . append ( variant_id ) elif variant_id . startswith ( 'COSM' ) : transcript [ 'cosmic' ] . append ( int ( variant_id [ 4 : ] ) ) yield transcript
Parse transcript information from VCF variants

    Args:
        raw_transcripts(iterable(dict)): An iterable with raw transcript 
                                         information

    Yields:
        transcript(dict) A dictionary with transcript information
Parse transcript information from VCF variants