comparison ensembl_variant_report.py @ 0:c3a9e63e8c51 draft

planemo upload for repository https://github.com/jj-umn/galaxytools/tree/master/ensembl_variant_report commit 239c1ee096e5fc3e2e929f7bf2d4afba5c677d4b-dirty

0:c3a9e63e8c51

#!/usr/bin/env python
"""
Report variants to Ensembl Transcripts

FrameShift report line
# Gene	Variant Position	Reference	Variation	Prevalence	Sequencing Depth	Transcript	AA Position	AA change	AA Length	Stop Codon	Stop Region	AA Variation
TIGD6	chr5:149374879	AAG	AG	1.00	13	ENSG00000164296|ENST00000296736	345	Q344	522	A-TGA-T	KRWTSSRPST*

MissSense report line
# Gene	Variant Position	Reference	Variation	Prevalence	Sequencing Depth	Transcript	AA Position	AA change	AA Length	Stop Codon	Stop Region	AA Variation
FN1	chr2:216235089	G	A	1.00	7394	ENSG00000115414|ENST00000354785	2261	V2261I	2478	G-TAA	TGLTRGATYN_I_IVEALKDQQR
"""
import sys
import os.path
import re
import time
import optparse
from ensemblref import EnsemblRef
from Bio.Seq import reverse_complement, translate


def __main__():
    #Parse Command Line
    parser = optparse.OptionParser()
    #I/O
    parser.add_option( '-i', '--input', dest='input', default=None, help='Tabular file with peptide_sequence column' )
    parser.add_option( '-s', '--format', dest='format', default='tabular', choices=['tabular','snpeff'], help='Tabular file with peptide_sequence column' )
    #Columns for tabular input
    parser.add_option( '-C', '--chrom_column', type='int', dest='chrom_column', default=1, help='column ordinal with Ensembl transctip ID' )
    parser.add_option( '-P', '--pos_column', type='int', dest='pos_column', default=2, help='column ordinal with Ensembl transctip ID' )
    parser.add_option( '-R', '--ref_column', type='int', dest='ref_column', default=3, help='column ordinal with Ensembl transctip ID' )
    parser.add_option( '-A', '--alt_column', type='int', dest='alt_column', default=4, help='column ordinal with Ensembl transctip ID' )
    parser.add_option( '-T', '--transcript_column', type='int', dest='transcript_column', default=1, help='column ordinal with Ensembl transctip ID' )
    parser.add_option( '-F', '--dpr_column', type='int', dest='dpr_column', default=1, help='column with VCF: DPR or AD' )
    parser.add_option( '-D', '--dp_column', type='int', dest='dp_column', default=1, help='column with VCF: DP' )
    parser.add_option( '-g', '--gene_model', dest='gene_model', default=None, help='GTF gene model file. Used to annotate NSJ peptide entries.')
    parser.add_option( '-2', '--twobit', dest='twobit', default=None, help='Reference genome in UCSC twobit format')
    #Output file
    parser.add_option( '-o', '--output', dest='output', default=None, help='The output report (else write to stdout)' )
    #filters
    parser.add_option( '-d', '--min_depth', type='int', dest='min_depth', default=None, help='Minimum read depth to report' )
    parser.add_option( '-f', '--min_freq', type='float', dest='min_freq', default=None, help='Minimum variant frequency to report' )
    #peptide options
    parser.add_option( '-l', '--leading_aa', type='int', dest='leading_aa', default=10, help='Number AAs before missense variant' )
    parser.add_option( '-t', '--trailing_aa', type='int', dest='trailing_aa', default=10, help='Number AAs after missense variant' )
    parser.add_option( '-r', '--readthrough', type='int', dest='readthrough', default=0, help='' )
    # 
    parser.add_option('--debug', dest='debug', action='store_true', default=False, help='Print debugging messages')
    (options, args) = parser.parse_args()
    
    ##INPUTS##
    if options.input != None:
        try:
            inputPath = os.path.abspath(options.input)
            inputFile = open(inputPath, 'r')
        except Exception, e:
            print >> sys.stderr, "failed: %s" % e
            exit(2)
    else:
        inputFile = sys.stdin

    if options.output != None:
        try:
            outputPath = os.path.abspath(options.output)
            outputFile = open(outputPath, 'w')
        except Exception, e:
            print >> sys.stderr, "failed: %s" % e
            exit(3)
    else:
        outputFile = sys.stdout

    def parse_tabular():
        ci = options.chrom_column - 1
        pi = options.pos_column - 1
        ri = options.ref_column - 1
        ai = options.alt_column - 1
        ti = options.transcript_column - 1
        di = options.dp_column - 1
        fi = options.dpr_column - 1
        for linenum,line in enumerate(inputFile):
            if options.debug:
                print >> sys.stderr, "%d: %s\n" % (linenum,line)
            if line.startswith('#'):
                continue
            if line.strip() == '':
                continue
            fields = line.rstrip('\r\n').split('\t')
            transcript = fields[ti]
            if not transcript:
                print >> sys.stderr, "%d: %s\n" % (linenum,line)
                continue
            chrom = fields[ci]
            pos = int(fields[pi])
            ref = fields[ri]
            alt = fields[ai]
            dp = int(fields[di])
            dpr = [int(x) for x in fields[fi].split(',')]
            yield (transcript,pos,ref,alt,dp,dpr)

    def parse_snpeff_vcf():
        for linenum,line in enumerate(inputFile):
            if line.startswith('##'): 
                if line.find('SnpEffVersion=') > 0:
                    SnpEffVersion = re.search('SnpEffVersion="?(\d+\.\d+)',line).groups()[0]
            elif line.startswith('#CHROM'): 
                pass
            else:
                fields = line.strip('\r\n').split('\t')
                if options.debug: print >> sys.stderr, "\n%s" % (fields)
                (chrom, pos, id, ref, alts, qual, filter, info) = fields[0:8]
                pos = int(pos)
                qual = float(qual)
                for info_item in info.split(';'):
                    if info_item.find('=') < 0: continue
                    (key, val) = info_item.split('=', 1)
                    if key == 'DP':
                        dp = int(val)
                    if key == 'DPR':
                        dpr = dpr = [int(x) for x in val.split(',')]
                    if key in ['EFF','ANN']:
                        for effect in val.split(','):
                            if options.debug: print >> sys.stderr, "\n%s" % (effect.split('|'))
                            if key == 'ANN':
                                (alt,eff,impact,gene_name,gene_id,feature_type,transcript,biotype,exon,c_hgvs,p_hgvs,cdna,cds,aa,distance,info) = effect.split('|')
                            elif key == 'EFF':
                                (eff, effs) = effect.rstrip(')').split('(')
                                (impact, functional_class, codon_change, aa_change, aa_len, gene_name, biotype, coding, transcript, exon, alt) = effs.split('|')[0:11]
                            yield (transcript,pos,ref,alt,dp,dpr)


    #Process gene model
    ens_ref = None
    if options.gene_model != None:
        try:
            geneModelFile = os.path.abspath(options.gene_model)
            twoBitFile = os.path.abspath(options.twobit)
            print >> sys.stderr, "Parsing ensembl ref: %s %s" % (options.gene_model,options.twobit)
            time1 = time.time()
            ens_ref = EnsemblRef(geneModelFile,twoBitFile)
            time2 = time.time()
            print >> sys.stderr, "Parsing ensembl ref: %d seconds" % (int(time2-time1))
        except Exception, e:
            print >> sys.stderr, "Parsing gene model failed: %s" % e
            exit(2)
    try:
        parse_input = parse_tabular if options.format == 'tabular' else parse_snpeff_vcf
        for tid,pos1,ref,alt,dp,dpr in parse_input():
            freq = float(dpr[1])/float(sum(dpr)) if dpr else None
            if options.min_depth and dp is not None and dp < options.min_depth:
                continue
            if options.min_freq and freq is not None and freq < options.min_freq:
                continue
            ## transcript_id, pos, ref, alt, dp, dpr
            tx = ens_ref.get_gtf_transcript(tid)
            coding = ens_ref.transcript_is_coding(tid)
            if not coding:
                continue
            frame_shift = len(ref) != len(alt)
            cds = ens_ref.get_cds(tid)
            pos0 = pos1 - 1 # zero based position
            spos = pos0 if tx.gene.strand else pos0 + len(ref) - 1
            alt_seq = alt if tx.gene.strand else reverse_complement(alt)
            cds_pos = ens_ref.genome_to_cds_pos(tid, spos)
            alt_cds = cds[:cds_pos] + alt_seq + cds[cds_pos+len(ref):]
            offset = 0
            if tx.gene.strand:
                for i in range(min(len(ref),len(alt))):
                    if ref[i] == alt[i]:
                        offset = i
                    else:
                        break
            else:
                for i in range(-1,-min(len(ref),len(alt)) -1,-1):
                    if ref[i] == alt[i]:
                        offset = i
                    else:
                        break
            refpep = translate(cds[:len(cds)/3*3])
            pep = translate(alt_cds[:len(alt_cds)/3*3])
            peplen = len(pep)
            aa_pos = (cds_pos + offset) / 3
            if aa_pos >= len(pep):
                print >> sys.stderr, "%d: %s\n" % (linenum,line)
                continue
            if frame_shift:
                #find stop_codons
                nstops = 0
                stop_codons = []
                for i in range(aa_pos,peplen):
                    if refpep[i] != pep[i]:
                        aa_pos = i
                        break
                for i in range(aa_pos,peplen):
                    if pep[i] == '*':
                        nstops += 1
                        stop_codons.append("%s-%s%s" % (alt_cds[i*3-1],alt_cds[i*3:i*3+3],"-%s" % alt_cds[i*3+4] if len(alt_cds) > i*3 else ''))
                        if nstops > options.readthrough:
                            reported_peptide = pep[aa_pos:i+1]
                            reported_stop_codon = ','.join(stop_codons)
                            break
            else:
                reported_stop_codon = "%s-%s" % (alt_cds[peplen*3-4],alt_cds[peplen*3-3:peplen*3])
                reported_peptide = "%s_%s_%s" % (pep[max(aa_pos-options.leading_aa,0):aa_pos],
                                                 pep[aa_pos],
                                                 pep[aa_pos+1:min(aa_pos+1+options.trailing_aa,len(pep))])
            cs_pos = aa_pos * 3
            aa_pos = (cds_pos + offset) / 3
            ref_codon = cds[cs_pos:cs_pos+3]
            ref_aa = translate(ref_codon)
            alt_codon = alt_cds[cs_pos:cs_pos+3]
            alt_aa = translate(alt_codon)
            gene = tx.gene.names[0]
            report_fields = [tx.gene.names[0],
             '%s:%d %s' % (tx.gene.contig,pos1,'+' if tx.gene.strand else '-'),
             ref,
             alt,
             "%1.2f" % freq if freq is not None else '',
             str(dp),
             "%s|%s" % (tx.gene.gene_id,tx.cdna_id),
             "%d" % (aa_pos + 1),
             "%s%d%s" % (ref_aa,aa_pos + 1,alt_aa),
             "%d" % len(pep),
             reported_stop_codon,
             reported_peptide
             ]
            if options.debug:
                report_fields.append("%d %d %d %d  %s  %s" % (cds_pos, offset, cs_pos,aa_pos,ref_codon,alt_codon))
            outputFile.write('\t'.join(report_fields))
            if options.debug:
                print >> sys.stderr, "%s %s\n%s\n%s\n%s %s" % (
                         cds[cs_pos-6:cs_pos], cds[cs_pos:cs_pos+15],
                         translate(cds[cs_pos-6:cs_pos+15]),
                         translate(alt_cds[cs_pos-6:cs_pos+15]),
                         alt_cds[cs_pos-6:cs_pos], alt_cds[cs_pos:cs_pos+15])
            outputFile.write('\n')
    except Exception, e:
        print >> sys.stderr, "failed: %s" % e
        exit(1)


if __name__ == "__main__" : __main__()