Mercurial > repos > jjohnson > ensembl_variant_report

comparison ensembl_variant_report.py @ 7:d5cb252c68da draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/jj-umn/galaxytools/tree/master/ensembl_variant_report commit b62e927469f96a857e880c77530c50c885ad92fc-dirty
author jjohnson
date Thu, 14 Jun 2018 18:07:10 -0400
parents d9fad18ffdb1
children fd612f8119a2
comparison
equal deleted inserted replaced
6:d9fad18ffdb1 7:d5cb252c68da
150 print >> sys.stderr, "Parsing gene model failed: %s" % e 150 print >> sys.stderr, "Parsing gene model failed: %s" % e
151 exit(2) 151 exit(2)
152 try: 152 try:
153 parse_input = parse_tabular if options.format == 'tabular' else parse_snpeff_vcf 153 parse_input = parse_tabular if options.format == 'tabular' else parse_snpeff_vcf
154 for tid,pos1,ref,alt,dp,freq in parse_input(): 154 for tid,pos1,ref,alt,dp,freq in parse_input():
155 if options.debug: print >> sys.stderr, "%s\t%d\t%s\t%s\t%d\t%f" % (tid,pos1,ref,alt,dp,freq)
155 if not tid: 156 if not tid:
156 continue 157 continue
157 if options.min_depth and dp is not None and dp < options.min_depth: 158 if options.min_depth and dp is not None and dp < options.min_depth:
158 continue 159 continue
159 if options.min_freq and freq is not None and freq < options.min_freq: 160 if options.min_freq and freq is not None and freq < options.min_freq:
160 continue 161 continue
161 ## transcript_id, pos, ref, alt, dp, dpr 162 try:
162 tx = ens_ref.get_gtf_transcript(tid) 163 ## transcript_id, pos, ref, alt, dp, dpr
163 if not tx: 164 tx = ens_ref.get_gtf_transcript(tid)
164 continue 165 if not tx and tid.find('.') > 0:
165 coding = ens_ref.transcript_is_coding(tid) 166 tid = tid.split('.')[0]
166 if not coding: 167 tx = ens_ref.get_gtf_transcript(tid)
167 continue 168 if not tx:
168 frame_shift = len(ref) != len(alt) 169 continue
169 cds = ens_ref.get_cds(tid) 170 if options.debug: print >> sys.stderr, "%s" % (tx)
170 pos0 = pos1 - 1 # zero based position 171 coding = ens_ref.transcript_is_coding(tid)
171 spos = pos0 if tx.gene.strand else pos0 + len(ref) - 1 172 if not coding:
172 alt_seq = alt if tx.gene.strand else reverse_complement(alt) 173 continue
173 ref_seq = ref if tx.gene.strand else reverse_complement(ref) 174 frame_shift = len(ref) != len(alt)
174 cds_pos = ens_ref.genome_to_cds_pos(tid, spos) 175 cds = ens_ref.get_cds(tid)
175 alt_cds = cds[:cds_pos] + alt_seq + cds[cds_pos+len(ref):] if cds_pos+len(ref) < len(cds) else '' 176 if options.debug or not cds: print >> sys.stderr, "cds:\n%s" % (cds)
176 offset = 0 177 pos0 = pos1 - 1 # zero based position
177 if tx.gene.strand: 178 spos = pos0 if tx.gene.strand else pos0 + len(ref) - 1
178 for i in range(min(len(ref),len(alt))): 179 alt_seq = alt if tx.gene.strand else reverse_complement(alt)
179 if ref[i] == alt[i]: 180 ref_seq = ref if tx.gene.strand else reverse_complement(ref)
180 offset = i 181 cds_pos = ens_ref.genome_to_cds_pos(tid, spos)
181 else: 182 if options.debug: print >> sys.stderr, "cds_pos: %s" % (str(cds_pos))
182 break 183 alt_cds = cds[:cds_pos] + alt_seq + cds[cds_pos+len(ref):] if cds_pos+len(ref) < len(cds) else ''
183 else: 184 offset = 0
184 for i in range(-1,-min(len(ref),len(alt)) -1,-1): 185 if tx.gene.strand:
185 if ref[i] == alt[i]: 186 for i in range(min(len(ref),len(alt))):
186 offset = i 187 if ref[i] == alt[i]:
187 else: 188 offset = i
188 break 189 else:
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 190 break
212 else: 191 else:
213 reported_stop_codon = "%s-%s" % (alt_cds[peplen*3-4],alt_cds[peplen*3-3:peplen*3]) 192 for i in range(-1,-min(len(ref),len(alt)) -1,-1):
214 reported_peptide = "%s_%s_%s" % (pep[max(aa_pos-options.leading_aa,0):aa_pos], 193 if ref[i] == alt[i]:
215 pep[aa_pos], 194 offset = i
216 pep[aa_pos+1:min(aa_pos+1+options.trailing_aa,len(pep))]) 195 else:
217 cs_pos = aa_pos * 3 196 break
218 aa_pos = (cds_pos + offset) / 3 197 refpep = translate(cds[:len(cds)/3*3])
219 ref_codon = cds[cs_pos:cs_pos+3] 198 pep = translate(alt_cds[:len(alt_cds)/3*3])
220 ref_aa = translate(ref_codon) 199 peplen = len(pep)
221 alt_codon = alt_cds[cs_pos:cs_pos+3] 200 aa_pos = (cds_pos + offset) / 3
222 alt_aa = translate(alt_codon) 201 reported_stop_codon = ''
223 gene = tx.gene.names[0] 202 if aa_pos >= len(pep):
224 report_fields = [tx.gene.names[0], 203 print >> sys.stderr, "aa_pos %d >= peptide length %d : %s %d %s %s" % (aa_pos,len(pep),tid,pos1,ref,alt)
225 '%s:%d %s' % (tx.gene.contig,pos1,'+' if tx.gene.strand else '-'), 204 continue
226 ref_seq, 205 if frame_shift:
227 alt_seq, 206 #find stop_codons
228 "%1.2f" % freq if freq is not None else '', 207 nstops = 0
229 str(dp), 208 stop_codons = []
230 "%s|%s" % (tx.gene.gene_id,tx.cdna_id), 209 for i in range(aa_pos,peplen):
231 "%d" % (aa_pos + 1), 210 if refpep[i] != pep[i]:
232 "%s%d%s" % (ref_aa,aa_pos + 1,alt_aa), 211 aa_pos = i
233 "%d" % len(pep), 212 break
234 reported_stop_codon, 213 reported_peptide = pep[aa_pos:]
235 reported_peptide 214 for i in range(aa_pos,peplen):
236 ] 215 if pep[i] == '*':
237 if options.debug: 216 nstops += 1
238 report_fields.append("%d %d %d %d %s %s" % (cds_pos, offset, cs_pos,aa_pos,ref_codon,alt_codon)) 217 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 ''))
239 outputFile.write('\t'.join(report_fields)) 218 if nstops > options.readthrough:
240 if options.debug: 219 reported_peptide = pep[aa_pos:i+1]
241 print >> sys.stderr, "%s %s\n%s\n%s\n%s %s" % ( 220 reported_stop_codon = ','.join(stop_codons)
242 cds[cs_pos-6:cs_pos], cds[cs_pos:cs_pos+15], 221 break
243 translate(cds[cs_pos-6:cs_pos+15]), 222 else:
244 translate(alt_cds[cs_pos-6:cs_pos+15]), 223 reported_stop_codon = "%s-%s" % (alt_cds[peplen*3-4],alt_cds[peplen*3-3:peplen*3])
245 alt_cds[cs_pos-6:cs_pos], alt_cds[cs_pos:cs_pos+15]) 224 reported_peptide = "%s_%s_%s" % (pep[max(aa_pos-options.leading_aa,0):aa_pos],
246 outputFile.write('\n') 225 pep[aa_pos],
226 pep[aa_pos+1:min(aa_pos+1+options.trailing_aa,len(pep))])
227 cs_pos = aa_pos * 3
228 aa_pos = (cds_pos + offset) / 3
229 ref_codon = cds[cs_pos:cs_pos+3]
230 ref_aa = translate(ref_codon)
231 alt_codon = alt_cds[cs_pos:cs_pos+3]
232 alt_aa = translate(alt_codon)
233 gene = tx.gene.names[0]
234 report_fields = [tx.gene.names[0],
235 '%s:%d %s' % (tx.gene.contig,pos1,'+' if tx.gene.strand else '-'),
236 ref_seq,
237 alt_seq,
238 "%1.2f" % freq if freq is not None else '',
239 str(dp),
240 "%s|%s" % (tx.gene.gene_id,tx.cdna_id),
241 "%d" % (aa_pos + 1),
242 "%s%d%s" % (ref_aa,aa_pos + 1,alt_aa),
243 "%d" % len(pep),
244 reported_stop_codon,
245 reported_peptide,
246 tx.get_transcript_type_name()
247 ]
248 if options.debug:
249 report_fields.append("%d %d %d %d %s %s" % (cds_pos, offset, cs_pos,aa_pos,ref_codon,alt_codon))
250 outputFile.write('\t'.join(report_fields))
251 if options.debug:
252 print >> sys.stderr, "%s %s\n%s\n%s\n%s %s" % (
253 cds[cs_pos-6:cs_pos], cds[cs_pos:cs_pos+15],
254 translate(cds[cs_pos-6:cs_pos+15]),
255 translate(alt_cds[cs_pos-6:cs_pos+15]),
256 alt_cds[cs_pos-6:cs_pos], alt_cds[cs_pos:cs_pos+15])
257 outputFile.write('\n')
258 except Exception, e:
259 print >> sys.stderr, "failed: %s" % e
260 print >> sys.stderr, "%s\t%d\t%s\t%s\t%d\t%f\t%s" % (tid,pos1,ref,alt,dp,freq,e)
247 except Exception, e: 261 except Exception, e:
248 print >> sys.stderr, "failed: %s" % e 262 print >> sys.stderr, "failed: %s" % e
249 exit(1) 263 exit(1)
250 264
251 265