Mercurial > repos > jjohnson > ensembl_variant_report
comparison ensembl_variant_report.py @ 0:9f4ea174ce3d draft
planemo upload for repository https://github.com/jj-umn/galaxytools/tree/master/ensembl_variant_report commit e6aa05bbbee3cc7d98f16354fc41c674f439ff1b-dirty
| author | jjohnson |
|---|---|
| date | Thu, 14 Jun 2018 17:51:39 -0400 |
| parents | |
| children | a67b4de184c2 |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:9f4ea174ce3d |
|---|---|
| 1 #!/usr/bin/env python | |
| 2 """ | |
| 3 Report variants to Ensembl Transcripts | |
| 4 | |
| 5 FrameShift report line | |
| 6 # Gene Variant Position Reference Variation Prevalence Sequencing Depth Transcript AA Position AA change AA Length Stop Codon Stop Region AA Variation | |
| 7 TIGD6 chr5:149374879 AAG AG 1.00 13 ENSG00000164296|ENST00000296736 345 Q344 522 A-TGA-T KRWTSSRPST* | |
| 8 | |
| 9 MissSense report line | |
| 10 # Gene Variant Position Reference Variation Prevalence Sequencing Depth Transcript AA Position AA change AA Length Stop Codon Stop Region AA Variation | |
| 11 FN1 chr2:216235089 G A 1.00 7394 ENSG00000115414|ENST00000354785 2261 V2261I 2478 G-TAA TGLTRGATYN_I_IVEALKDQQR | |
| 12 """ | |
| 13 import sys | |
| 14 import os.path | |
| 15 import re | |
| 16 import time | |
| 17 import optparse | |
| 18 from ensemblref import EnsemblRef | |
| 19 from Bio.Seq import reverse_complement, translate | |
| 20 | |
| 21 | |
| 22 def __main__(): | |
| 23 #Parse Command Line | |
| 24 parser = optparse.OptionParser() | |
| 25 #I/O | |
| 26 parser.add_option( '-i', '--input', dest='input', default=None, help='Tabular file with peptide_sequence column' ) | |
| 27 parser.add_option( '-s', '--format', dest='format', default='tabular', choices=['tabular','snpeff'], help='Tabular file with peptide_sequence column' ) | |
| 28 #Columns for tabular input | |
| 29 parser.add_option( '-C', '--chrom_column', type='int', dest='chrom_column', default=1, help='column ordinal with Ensembl transctip ID' ) | |
| 30 parser.add_option( '-P', '--pos_column', type='int', dest='pos_column', default=2, help='column ordinal with Ensembl transctip ID' ) | |
| 31 parser.add_option( '-R', '--ref_column', type='int', dest='ref_column', default=3, help='column ordinal with Ensembl transctip ID' ) | |
| 32 parser.add_option( '-A', '--alt_column', type='int', dest='alt_column', default=4, help='column ordinal with Ensembl transctip ID' ) | |
| 33 parser.add_option( '-T', '--transcript_column', type='int', dest='transcript_column', default=1, help='column ordinal with Ensembl transctip ID' ) | |
| 34 parser.add_option( '-F', '--dpr_column', type='int', dest='dpr_column', default=1, help='column with VCF: DPR or AD' ) | |
| 35 parser.add_option( '-D', '--dp_column', type='int', dest='dp_column', default=1, help='column with VCF: DP' ) | |
| 36 parser.add_option( '-g', '--gene_model', dest='gene_model', default=None, help='GTF gene model file. Used to annotate NSJ peptide entries.') | |
| 37 parser.add_option( '-2', '--twobit', dest='twobit', default=None, help='Reference genome in UCSC twobit format') | |
| 38 #Output file | |
| 39 parser.add_option( '-o', '--output', dest='output', default=None, help='The output report (else write to stdout)' ) | |
| 40 #filters | |
| 41 parser.add_option( '-d', '--min_depth', type='int', dest='min_depth', default=None, help='Minimum read depth to report' ) | |
| 42 parser.add_option( '-f', '--min_freq', type='float', dest='min_freq', default=None, help='Minimum variant frequency to report' ) | |
| 43 #peptide options | |
| 44 parser.add_option( '-l', '--leading_aa', type='int', dest='leading_aa', default=10, help='Number AAs before missense variant' ) | |
| 45 parser.add_option( '-t', '--trailing_aa', type='int', dest='trailing_aa', default=10, help='Number AAs after missense variant' ) | |
| 46 parser.add_option( '-r', '--readthrough', type='int', dest='readthrough', default=0, help='' ) | |
| 47 # | |
| 48 parser.add_option('--debug', dest='debug', action='store_true', default=False, help='Print debugging messages') | |
| 49 (options, args) = parser.parse_args() | |
| 50 | |
| 51 ##INPUTS## | |
| 52 if options.input != None: | |
| 53 try: | |
| 54 inputPath = os.path.abspath(options.input) | |
| 55 inputFile = open(inputPath, 'r') | |
| 56 except Exception, e: | |
| 57 print >> sys.stderr, "failed: %s" % e | |
| 58 exit(2) | |
| 59 else: | |
| 60 inputFile = sys.stdin | |
| 61 | |
| 62 if options.output != None: | |
| 63 try: | |
| 64 outputPath = os.path.abspath(options.output) | |
| 65 outputFile = open(outputPath, 'w') | |
| 66 except Exception, e: | |
| 67 print >> sys.stderr, "failed: %s" % e | |
| 68 exit(3) | |
| 69 else: | |
| 70 outputFile = sys.stdout | |
| 71 | |
| 72 def parse_tabular(): | |
| 73 ci = options.chrom_column - 1 | |
| 74 pi = options.pos_column - 1 | |
| 75 ri = options.ref_column - 1 | |
| 76 ai = options.alt_column - 1 | |
| 77 ti = options.transcript_column - 1 | |
| 78 di = options.dp_column - 1 | |
| 79 fi = options.dpr_column - 1 | |
| 80 for linenum,line in enumerate(inputFile): | |
| 81 if options.debug: | |
| 82 print >> sys.stderr, "%d: %s\n" % (linenum,line) | |
| 83 if line.startswith('#'): | |
| 84 continue | |
| 85 if line.strip() == '': | |
| 86 continue | |
| 87 fields = line.rstrip('\r\n').split('\t') | |
| 88 transcript = fields[ti] | |
| 89 if not transcript: | |
| 90 print >> sys.stderr, "%d: %s\n" % (linenum,line) | |
| 91 continue | |
| 92 chrom = fields[ci] | |
| 93 pos = int(fields[pi]) | |
| 94 ref = fields[ri] | |
| 95 alts = fields[ai] | |
| 96 dp = int(fields[di]) | |
| 97 dpr = [int(x) for x in fields[fi].split(',')] | |
| 98 for i,alt in enumerate(alts.split(',')): | |
| 99 freq = float(dpr[i+1])/float(sum(dpr)) if dpr else None | |
| 100 yield (transcript,pos,ref,alt,dp,freq) | |
| 101 | |
| 102 def parse_snpeff_vcf(): | |
| 103 for linenum,line in enumerate(inputFile): | |
| 104 if line.startswith('##'): | |
| 105 if line.find('SnpEffVersion=') > 0: | |
| 106 SnpEffVersion = re.search('SnpEffVersion="?(\d+\.\d+)',line).groups()[0] | |
| 107 elif line.startswith('#CHROM'): | |
| 108 pass | |
| 109 else: | |
| 110 fields = line.strip('\r\n').split('\t') | |
| 111 if options.debug: print >> sys.stderr, "\n%s" % (fields) | |
| 112 (chrom, pos, id, ref, alts, qual, filter, info) = fields[0:8] | |
| 113 alt_list = alts.split(',') | |
| 114 pos = int(pos) | |
| 115 qual = float(qual) | |
| 116 dp = None | |
| 117 dpr = None | |
| 118 for info_item in info.split(';'): | |
| 119 if info_item.find('=') < 0: continue | |
| 120 (key, val) = info_item.split('=', 1) | |
| 121 if key == 'DP': | |
| 122 dp = int(val) | |
| 123 if key == 'DPR': | |
| 124 dpr = [int(x) for x in val.split(',')] | |
| 125 if key in ['EFF','ANN']: | |
| 126 for effect in val.split(','): | |
| 127 if options.debug: print >> sys.stderr, "\n%s" % (effect.split('|')) | |
| 128 if key == 'ANN': | |
| 129 (alt,eff,impact,gene_name,gene_id,feature_type,transcript,biotype,exon,c_hgvs,p_hgvs,cdna,cds,aa,distance,info) = effect.split('|') | |
| 130 elif key == 'EFF': | |
| 131 (eff, effs) = effect.rstrip(')').split('(') | |
| 132 (impact, functional_class, codon_change, aa_change, aa_len, gene_name, biotype, coding, transcript, exon, alt) = effs.split('|')[0:11] | |
| 133 i = alt_list.index(alt) if alt in alt_list else 0 | |
| 134 freq = float(dpr[i+1])/float(sum(dpr)) if dpr else None | |
| 135 yield (transcript,pos,ref,alt,dp,freq) | |
| 136 | |
| 137 | |
| 138 #Process gene model | |
| 139 ens_ref = None | |
| 140 if options.gene_model != None: | |
| 141 try: | |
| 142 geneModelFile = os.path.abspath(options.gene_model) | |
| 143 twoBitFile = os.path.abspath(options.twobit) | |
| 144 print >> sys.stderr, "Parsing ensembl ref: %s %s" % (options.gene_model,options.twobit) | |
| 145 time1 = time.time() | |
| 146 ens_ref = EnsemblRef(geneModelFile,twoBitFile) | |
| 147 time2 = time.time() | |
| 148 print >> sys.stderr, "Parsing ensembl ref: %d seconds" % (int(time2-time1)) | |
| 149 except Exception, e: | |
| 150 print >> sys.stderr, "Parsing gene model failed: %s" % e | |
| 151 exit(2) | |
| 152 try: | |
| 153 parse_input = parse_tabular if options.format == 'tabular' else parse_snpeff_vcf | |
| 154 for tid,pos1,ref,alt,dp,freq in parse_input(): | |
| 155 if not tid: | |
| 156 continue | |
| 157 if options.min_depth and dp is not None and dp < options.min_depth: | |
| 158 continue | |
| 159 if options.min_freq and freq is not None and freq < options.min_freq: | |
| 160 continue | |
| 161 ## transcript_id, pos, ref, alt, dp, dpr | |
| 162 tx = ens_ref.get_gtf_transcript(tid) | |
| 163 if not tx: | |
| 164 continue | |
| 165 coding = ens_ref.transcript_is_coding(tid) | |
| 166 if not coding: | |
| 167 continue | |
| 168 frame_shift = len(ref) != len(alt) | |
| 169 cds = ens_ref.get_cds(tid) | |
| 170 pos0 = pos1 - 1 # zero based position | |
| 171 spos = pos0 if tx.gene.strand else pos0 + len(ref) - 1 | |
| 172 alt_seq = alt if tx.gene.strand else reverse_complement(alt) | |
| 173 ref_seq = ref if tx.gene.strand else reverse_complement(ref) | |
| 174 cds_pos = ens_ref.genome_to_cds_pos(tid, spos) | |
| 175 alt_cds = cds[:cds_pos] + alt_seq + cds[cds_pos+len(ref):] if cds_pos+len(ref) < len(cds) else '' | |
| 176 offset = 0 | |
| 177 if tx.gene.strand: | |
| 178 for i in range(min(len(ref),len(alt))): | |
| 179 if ref[i] == alt[i]: | |
| 180 offset = i | |
| 181 else: | |
| 182 break | |
| 183 else: | |
| 184 for i in range(-1,-min(len(ref),len(alt)) -1,-1): | |
| 185 if ref[i] == alt[i]: | |
| 186 offset = i | |
| 187 else: | |
| 188 break | |
| 189 refpep = translate(cds[:len(cds)/3*3]) | |
| 190 pep = translate(alt_cds[:len(alt_cds)/3*3]) | |
| 191 peplen = len(pep) | |
| 192 aa_pos = (cds_pos + offset) / 3 | |
| 193 if aa_pos >= len(pep): | |
| 194 print >> sys.stderr, "aa_pos %d >= peptide length %d : %s %d %s %s\n" % (aa_pos,len(pep),tid,pos1,ref,alt) | |
| 195 continue | |
| 196 if frame_shift: | |
| 197 #find stop_codons | |
| 198 nstops = 0 | |
| 199 stop_codons = [] | |
| 200 for i in range(aa_pos,peplen): | |
| 201 if refpep[i] != pep[i]: | |
| 202 aa_pos = i | |
| 203 break | |
| 204 for i in range(aa_pos,peplen): | |
| 205 if pep[i] == '*': | |
| 206 nstops += 1 | |
| 207 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 '')) | |
| 208 if nstops > options.readthrough: | |
| 209 reported_peptide = pep[aa_pos:i+1] | |
| 210 reported_stop_codon = ','.join(stop_codons) | |
| 211 break | |
| 212 else: | |
| 213 reported_stop_codon = "%s-%s" % (alt_cds[peplen*3-4],alt_cds[peplen*3-3:peplen*3]) | |
| 214 reported_peptide = "%s_%s_%s" % (pep[max(aa_pos-options.leading_aa,0):aa_pos], | |
| 215 pep[aa_pos], | |
| 216 pep[aa_pos+1:min(aa_pos+1+options.trailing_aa,len(pep))]) | |
| 217 cs_pos = aa_pos * 3 | |
| 218 aa_pos = (cds_pos + offset) / 3 | |
| 219 ref_codon = cds[cs_pos:cs_pos+3] | |
| 220 ref_aa = translate(ref_codon) | |
| 221 alt_codon = alt_cds[cs_pos:cs_pos+3] | |
| 222 alt_aa = translate(alt_codon) | |
| 223 gene = tx.gene.names[0] | |
| 224 report_fields = [tx.gene.names[0], | |
| 225 '%s:%d %s' % (tx.gene.contig,pos1,'+' if tx.gene.strand else '-'), | |
| 226 ref_seq, | |
| 227 alt_seq, | |
| 228 "%1.2f" % freq if freq is not None else '', | |
| 229 str(dp), | |
| 230 "%s|%s" % (tx.gene.gene_id,tx.cdna_id), | |
| 231 "%d" % (aa_pos + 1), | |
| 232 "%s%d%s" % (ref_aa,aa_pos + 1,alt_aa), | |
| 233 "%d" % len(pep), | |
| 234 reported_stop_codon, | |
| 235 reported_peptide | |
| 236 ] | |
| 237 if options.debug: | |
| 238 report_fields.append("%d %d %d %d %s %s" % (cds_pos, offset, cs_pos,aa_pos,ref_codon,alt_codon)) | |
| 239 outputFile.write('\t'.join(report_fields)) | |
| 240 if options.debug: | |
| 241 print >> sys.stderr, "%s %s\n%s\n%s\n%s %s" % ( | |
| 242 cds[cs_pos-6:cs_pos], cds[cs_pos:cs_pos+15], | |
| 243 translate(cds[cs_pos-6:cs_pos+15]), | |
| 244 translate(alt_cds[cs_pos-6:cs_pos+15]), | |
| 245 alt_cds[cs_pos-6:cs_pos], alt_cds[cs_pos:cs_pos+15]) | |
| 246 outputFile.write('\n') | |
| 247 except Exception, e: | |
| 248 print >> sys.stderr, "failed: %s" % e | |
| 249 exit(1) | |
| 250 | |
| 251 | |
| 252 if __name__ == "__main__" : __main__() | |
| 253 |