Mercurial > repos > jjohnson > translate_bed_sequences
diff translate_bed_sequences.py @ 5:c626a939eef7 draft default tip
Uploaded
| author | jjohnson |
|---|---|
| date | Tue, 12 Jan 2016 14:38:03 -0500 |
| parents | 3b526a780849 |
| children |
line wrap: on
line diff
--- a/translate_bed_sequences.py Thu Jan 30 13:26:58 2014 -0600 +++ b/translate_bed_sequences.py Tue Jan 12 14:38:03 2016 -0500 @@ -19,6 +19,7 @@ """ import sys,re,os.path +import tempfile import optparse from optparse import OptionParser from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate @@ -27,7 +28,9 @@ def __init__(self, line): self.line = line try: - (chrom,chromStart,chromEnd,name,score,strand,thickStart,thickEnd,itemRgb,blockCount,blockSizes,blockStarts,seq) = line.split('\t')[0:13] + fields = line.rstrip('\r\n').split('\t') + (chrom,chromStart,chromEnd,name,score,strand,thickStart,thickEnd,itemRgb,blockCount,blockSizes,blockStarts) = fields[0:12] + seq = fields[12] if len(fields) > 12 else None self.chrom = chrom self.chromStart = int(chromStart) self.chromEnd = int(chromEnd) @@ -44,6 +47,12 @@ except Exception, e: print >> sys.stderr, "Unable to read Bed entry" % e exit(1) + def __str__(self): + return '%s\t%d\t%d\t%s\t%d\t%s\t%d\t%d\t%s\t%d\t%s\t%s%s' % ( + self.chrom, self.chromStart, self.chromEnd, self.name, self.score, self.strand, self.thickStart, self.thickEnd, self.itemRgb, self.blockCount, + ','.join([str(x) for x in self.blockSizes]), + ','.join([str(x) for x in self.blockStarts]), + '\t%s' % self.seq if self.seq else '') def get_splice_junctions(self): splice_juncs = [] for i in range(self.blockCount - 1): @@ -80,17 +89,59 @@ if translation: translations.append(translation) return translations - ## [[start,end,seq],[start,end,seq],[start,end,seq]] + ## (start,end) + def get_subrange(self,tstart,tstop): + chromStart = self.chromStart + chromEnd = self.chromEnd + r = range(self.blockCount) + if self.strand == '-': + r.reverse() + bStart = 0 + for x in r: + bEnd = bStart + self.blockSizes[x] + if bStart <= tstart < bEnd: + if self.strand == '+': + chromStart = self.chromStart + self.blockStarts[x] + (tstart - bStart) + else: + chromEnd = self.chromStart + self.blockStarts[x] + (tstart - bStart) + if bStart <= tstop < bEnd: + if self.strand == '+': + chromEnd = self.chromStart + self.blockStarts[x] + (tstop - bStart) + else: + chromStart = self.chromStart + self.blockStarts[x] + self.blockSizes[x] - (tstop - bStart) + bStart += self.blockSizes[x] + return(chromStart,chromEnd) + #get the blocks for sub range + def get_blocks(self,chromStart,chromEnd): + tblockCount = 0 + tblockSizes = [] + tblockStarts = [] + for x in range(self.blockCount): + bStart = self.chromStart + self.blockStarts[x] + bEnd = bStart + self.blockSizes[x] + if bStart > chromEnd: + break + if bEnd < chromStart: + continue + cStart = max(chromStart,bStart) + tblockStarts.append(cStart - chromStart) + tblockSizes.append(min(chromEnd,bEnd) - cStart) + tblockCount += 1 + ## print >> sys.stderr, "tblockCount: %d tblockStarts: %s tblockSizes: %s" % (tblockCount,tblockStarts,tblockSizes) + return (tblockCount,tblockSizes,tblockStarts) + ## [(start,end,seq,blockCount,blockSizes,blockStarts),(start,end,seq,blockCount,blockSizes,blockStarts),(start,end,seq,blockCount,blockSizes,blockStarts)] ## filter: ignore translation if stop codon in first exon after ignore_left_bp - def get_filterd_translations(self,untrimmed=False,filtering=True,ignore_left_bp=0,ignore_right_bp=0): - translations = [None,None,None] + def get_filterd_translations(self,untrimmed=False,filtering=True,ignore_left_bp=0,ignore_right_bp=0,debug=False): + translations = [None,None,None,None,None,None] seq = self.get_spliced_seq() ignore = (ignore_left_bp if self.strand == '+' else ignore_right_bp) / 3 block_sum = sum(self.blockSizes) - exon_sizes = self.blockSizes + exon_sizes = [x for x in self.blockSizes] if self.strand == '-': exon_sizes.reverse() splice_sites = [sum(exon_sizes[:x]) / 3 for x in range(1,len(exon_sizes))] + if debug: + print >> sys.stderr, "splice_sites: %s" % splice_sites junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0] if seq: for i in range(3): @@ -98,22 +149,32 @@ if translation: tstart = 0 tstop = len(translation) + offset = (block_sum - i) % 3 + if debug: + print >> sys.stderr, "frame: %d\ttstart: %d tstop: %d offset: %d\t%s" % (i,tstart,tstop,offset,translation) if not untrimmed: tstart = translation.rfind('*',0,junc) + 1 stop = translation.find('*',junc) tstop = stop if stop >= 0 else len(translation) + offset = (block_sum - i) % 3 + trimmed = translation[tstart:tstop] + if debug: + print >> sys.stderr, "frame: %d\ttstart: %d tstop: %d offset: %d\t%s" % (i,tstart,tstop,offset,trimmed) if filtering and tstart > ignore: continue - trimmed = translation[tstart:tstop] #get genomic locations for start and end - offset = (block_sum - i) % 3 if self.strand == '+': chromStart = self.chromStart + i + (tstart * 3) chromEnd = self.chromEnd - offset - (len(translation) - tstop) * 3 else: chromStart = self.chromStart + offset + (len(translation) - tstop) * 3 chromEnd = self.chromEnd - i - (tstart * 3) - translations[i] = [chromStart,chromEnd,trimmed] + #get the blocks for this translation + (tblockCount,tblockSizes,tblockStarts) = self.get_blocks(chromStart,chromEnd) + translations[i] = (chromStart,chromEnd,trimmed,tblockCount,tblockSizes,tblockStarts) + if debug: + print >> sys.stderr, "tblockCount: %d tblockStarts: %s tblockSizes: %s" % (tblockCount,tblockStarts,tblockSizes) + # translations[i] = (chromStart,chromEnd,trimmed,tblockCount,tblockSizes,tblockStarts) return translations def get_seq_id(self,seqtype='unk:unk',reference='',frame=None): ## Ensembl fasta ID format @@ -160,8 +221,14 @@ parser.add_option( '-i', '--input', dest='input', help='BED file (tophat junctions.bed) with sequence column added' ) parser.add_option( '-o', '--output', dest='output', help='Translations of spliced sequence') parser.add_option( '-b', '--bed_format', dest='bed_format', action='store_true', default=False, help='Append translations to bed file instead of fasta' ) + parser.add_option( '-D', '--fa_db', dest='fa_db', default=None, help='Prefix DB identifier for fasta ID line, e.g. generic' ) + parser.add_option( '-s', '--fa_sep', dest='fa_sep', default='|', help='fasta ID separator defaults to pipe char, e.g. generic|ProtID|description' ) + parser.add_option( '-B', '--bed', dest='bed', default=None, help='Output a bed file with added 13th column having translation' ) + parser.add_option( '-G', '--gff3', dest='gff', default=None, help='Output translations to a GFF3 file' ) parser.add_option( '-S', '--seqtype', dest='seqtype', default='pep:splice', help='SEQTYPE:STATUS for fasta ID line' ) + parser.add_option( '-P', '--id_prefix', dest='id_prefix', default='', help='prefix for the sequence ID' ) parser.add_option( '-R', '--reference', dest='reference', default=None, help='Genome Reference Name for fasta ID location ' ) + parser.add_option( '-r', '--refsource', dest='refsource', default=None, help='Source for Genome Reference, e.g. Ensembl, UCSC, or NCBI' ) parser.add_option( '-Q', '--score_name', dest='score_name', default=None, help='include in the fasta ID line score_name:score ' ) parser.add_option( '-l', '--leading_bp', dest='leading_bp', type='int', default=None, help='leading number of base pairs to ignore when filtering' ) parser.add_option( '-t', '--trailing_bp', dest='trailing_bp', type='int', default=None, help='trailing number of base pairs to ignore when filtering' ) @@ -182,10 +249,17 @@ else: inputFile = sys.stdin # Output files + bed_fh = None + gff_fh = None + gff_fa_file = None + gff_fa = None outFile = None if options.output == None: #write to stdout outFile = sys.stdout + if options.gff: + gff_fa_file = tempfile.NamedTemporaryFile(prefix='gff_fasta_',suffix=".fa",dir=os.getcwd()).name + gff_fa = open(gff_fa_file,'w') else: try: outPath = os.path.abspath(options.output) @@ -193,6 +267,16 @@ except Exception, e: print >> sys.stderr, "failed: %s" % e exit(3) + if options.gff: + gff_fa_file = outPath + if options.bed: + bed_fh = open(options.bed,'w') + bed_fh.write('track name="%s" description="%s" \n' % ('novel_junctioni_translations','test')) + if options.gff: + gff_fh = open(options.gff,'w') + gff_fh.write("##gff-version 3.2.1\n") + if options.reference: + gff_fh.write("##genome-build %s %s\n" % (options.refsource if options.refsource else 'unknown', options.reference)) leading_bp = 0 trailing_bp = 0 if options.leading_bp: @@ -233,21 +317,59 @@ tx_entry = "%s\t%s\n" % (line.rstrip('\r\n'),'\t'.join(translations)) outFile.write(tx_entry) else: - translations = entry.get_filterd_translations(untrimmed=options.untrimmed,filtering=options.filtering,ignore_left_bp=leading_bp,ignore_right_bp=trailing_bp) + translations = entry.get_filterd_translations(untrimmed=options.untrimmed,filtering=options.filtering,ignore_left_bp=leading_bp,ignore_right_bp=trailing_bp,debug=options.debug) for i,tx in enumerate(translations): if tx: - (chromStart,chromEnd,translation) = tx + (chromStart,chromEnd,translation,blockCount,blockSizes,blockStarts) = tx if options.min_length != None and len(translation) < options.min_length: continue if options.max_stop_codons != None and translation.count('*') > options.max_stop_codons: continue frame_name = '_%s' % (i + 1) + pep_id = "%s%s%s" % (options.id_prefix,entry.name,frame_name) + if bed_fh: + bed_fh.write('%s\t%d\t%d\t%s\t%d\t%s\t%d\t%d\t%s\t%d\t%s\t%s\t%s\n' % (str(entry.chrom),chromStart,chromEnd,pep_id,entry.score,entry.strand,chromStart,chromEnd,entry.itemRgb,blockCount,','.join([str(x) for x in blockSizes]),','.join([str(x) for x in blockStarts]),translation)) location = "chromosome:%s:%s:%s:%s:%s" % (options.reference,entry.chrom,chromStart,chromEnd,strand) + if blockCount: + location += " blockCount:%d blockSizes:%s blockStarts:%s" % (blockCount,','.join([str(x) for x in blockSizes]),','.join([str(x) for x in blockStarts])) score = " %s:%s" % (options.score_name,entry.score) if options.score_name else '' - seq_id = "%s%s %s %s%s" % (entry.name,frame_name,options.seqtype,location, score) - outFile.write(">%s\n" % seq_id) - outFile.write(translation) - outFile.write('\n') + seq_description = "%s %s%s" % (options.seqtype, location, score) + seq_id = "%s " % pep_id + if options.fa_db: + seq_id = "%s%s%s%s" % (options.fa_db,options.fa_sep,pep_id,options.fa_sep) + fa_id = "%s%s" % (seq_id,seq_description) + fa_entry = ">%s\n%s\n" % (fa_id,translation) + outFile.write(fa_entry) + if gff_fh: + if gff_fa: + gff_fa.write(fa_entry) + gff_fh.write("##sequence-region %s %d %d\n" % (entry.chrom,chromStart + 1,chromEnd - 1)) + gff_fh.write("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%d\tID=%s\n" % (entry.chrom,'splice_junc','gene',chromStart + 1,chromEnd - 1,entry.score,entry.strand,0,pep_id)) + for x in range(blockCount): + start = chromStart+blockStarts[x] + 1 + end = start + blockSizes[x] - 1 + phase = (3 - sum(blockSizes[:x]) % 3) % 3 + gff_fh.write("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%d\tParent=%s;ID=%s_%d\n" % (entry.chrom,'splice_junc','CDS',start,end,entry.score,entry.strand,phase,pep_id,pep_id,x)) + """ + ##gff-version 3 + ##sequence-region 19 1 287484 + 19 MassSpec peptide 282299 287484 10.0 - 0 ID=TEARLSFYSGHSSFGMYCMVFLALYVQ + 19 MassSpec CDS 287474 287484 . - 0 Parent=TEARLSFYSGHSSFGMYCMVFLALYVQ;transcript_id=ENST00000269812 + 19 MassSpec CDS 282752 282809 . - 1 Parent=TEARLSFYSGHSSFGMYCMVFLALYVQ;transcript_id=ENST00000269812 + 19 MassSpec CDS 282299 282310 . - 0 Parent=TEARLSFYSGHSSFGMYCMVFLALYVQ;transcript_id=ENST00000269812 + """ + if bed_fh: + bed_fh.close() + if gff_fh: + if gff_fa: + gff_fa.close() + else: + outFile.close() + gff_fa = open(gff_fa_file,'r') + gff_fh.write("##FASTA\n") + for i, line in enumerate(gff_fa): + gff_fh.write(line) + gff_fh.close() except Exception, e: print >> sys.stderr, "failed: Error reading %s - %s" % (options.input if options.input else 'stdin',e)