view translate_bed_sequences.py @ 4:aa93f7910259

Fix test case

#!/usr/bin/env python
"""
#
#------------------------------------------------------------------------------
#                         University of Minnesota
#         Copyright 2014, Regents of the University of Minnesota
#------------------------------------------------------------------------------
# Author:
#
#  James E Johnson
#
#------------------------------------------------------------------------------
"""

"""
Input:  BED file (12 column) + 13th sequence column appended by extract_genomic_dna
Output: Fasta of 3-frame translations of the spliced sequence
  
"""

import sys,re,os.path
import optparse
from optparse import OptionParser
from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate

class BedEntry( object ):
  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]
      self.chrom = chrom
      self.chromStart = int(chromStart)
      self.chromEnd = int(chromEnd)
      self.name = name
      self.score = int(score)
      self.strand = strand
      self.thickStart = int(thickStart)
      self.thickEnd = int(thickEnd)
      self.itemRgb = itemRgb
      self.blockCount = int(blockCount)
      self.blockSizes = [int(x) for x in blockSizes.split(',')]
      self.blockStarts = [int(x) for x in blockStarts.split(',')]
      self.seq = seq
    except Exception, e:
      print >> sys.stderr, "Unable to read Bed entry" % e
      exit(1)
  def get_splice_junctions(self): 
    splice_juncs = []
    for i in range(self.blockCount  - 1):
      splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1])
      splice_juncs.append(splice_junc)
    return splice_juncs
  def get_exon_seqs(self):
    exons = []
    for i in range(self.blockCount):
      # splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1])
      exons.append(self.seq[self.blockStarts[i]:self.blockStarts[i] + self.blockSizes[i]])
    if self.strand == '-':  #reverse complement
      exons.reverse()
      for i,s in enumerate(exons):
        exons[i] = reverse_complement(s)
    return exons
  def get_spliced_seq(self):
    seq = ''.join(self.get_exon_seqs())
    return seq
  def get_translation(self,sequence=None):
    translation = None
    seq = sequence if sequence else self.get_spliced_seq()
    if seq:
      seqlen = len(seq) / 3 * 3;
      if seqlen >= 3:
        translation = translate(seq[:seqlen])
    return translation
  def get_translations(self):
    translations = []
    seq = self.get_spliced_seq()
    if seq:
      for i in range(3):
        translation = self.get_translation(sequence=seq[i:])
        if translation:
          translations.append(translation)
    return translations
  ## [[start,end,seq],[start,end,seq],[start,end,seq]]
  ## 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]
    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
    if self.strand == '-':
      exon_sizes.reverse()
    splice_sites = [sum(exon_sizes[:x]) / 3 for x in range(1,len(exon_sizes))]
    junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0]
    if seq:
      for i in range(3):
        translation = self.get_translation(sequence=seq[i:])
        if translation:
          tstart = 0
          tstop = len(translation)
          if not untrimmed:
            tstart = translation.rfind('*',0,junc) + 1
            stop = translation.find('*',junc)
            tstop = stop if stop >= 0 else len(translation)
          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]
    return translations
  def get_seq_id(self,seqtype='unk:unk',reference='',frame=None):
    ## Ensembl fasta ID format
    # >ID SEQTYPE:STATUS LOCATION GENE TRANSCRIPT
    # >ENSP00000328693 pep:splice chromosome:NCBI35:1:904515:910768:1 gene:ENSG00000158815:transcript:ENST00000328693 gene_biotype:protein_coding transcript_biotype:protein_coding
    frame_name = ''
    chromStart = self.chromStart
    chromEnd = self.chromEnd
    strand = 1 if self.strand == '+' else -1
    if frame != None:
      block_sum = sum(self.blockSizes)
      offset = (block_sum - frame) % 3
      frame_name = '_' + str(frame + 1)
      if self.strand == '+':
        chromStart += frame
        chromEnd -= offset
      else:
        chromStart += offset
        chromEnd -= frame
    location = "chromosome:%s:%s:%s:%s:%s" % (reference,self.chrom,chromStart,chromEnd,strand)
    seq_id = "%s%s %s %s" % (self.name,frame_name,seqtype,location)
    return seq_id
  def get_line(self, start_offset = 0, end_offset = 0):
    if start_offset or end_offset:
      s_offset = start_offset if start_offset else 0
      e_offset = end_offset if end_offset else 0
      if s_offset > self.chromStart:
        s_offset = self.chromStart
      chrStart = self.chromStart - s_offset
      chrEnd = self.chromEnd + e_offset
      blkSizes = self.blockSizes
      blkSizes[0] += s_offset
      blkSizes[-1] += e_offset
      blkStarts = self.blockStarts
      for i in range(1,self.blockCount):
        blkStarts[i] += s_offset
      items = [str(x) for x in [self.chrom,chrStart,chrEnd,self.name,self.score,self.strand,self.thickStart,self.thickEnd,self.itemRgb,self.blockCount,','.join([str(x) for x in blkSizes]),','.join([str(x) for x in blkStarts])]]
      return '\t'.join(items) + '\n'
    return self.line

def __main__():
  #Parse Command Line
  parser = optparse.OptionParser()
  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( '-S', '--seqtype', dest='seqtype', default='pep:splice', help='SEQTYPE:STATUS for fasta ID line'  )
  parser.add_option( '-R', '--reference', dest='reference', default=None, help='Genome Reference Name for fasta ID location '  )
  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' )
  parser.add_option( '-U', '--unfiltered', dest='filtering', action='store_false', default=True, help='Do NOT filterout translation with stop codon in the first exon'  )
  parser.add_option( '-u', '--untrimmed', dest='untrimmed', action='store_true', default=False, help='Do NOT trim from splice site to stop codon'  )
  parser.add_option( '-L', '--min_length', dest='min_length', type='int', default=None, help='Minimun length (to first stop codon)'  )
  parser.add_option( '-M', '--max_stop_codons', dest='max_stop_codons', type='int', default=None, help='Filter out translations with more than max_stop_codons'  )
  parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stdout'  )
  (options, args) = parser.parse_args()
  # Input files
  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
  # Output files
  outFile = None
  if options.output == None:
    #write to stdout
    outFile = sys.stdout
  else:
    try:
      outPath = os.path.abspath(options.output)
      outFile = open(outPath, 'w')
    except Exception, e:
      print >> sys.stderr, "failed: %s" % e
      exit(3)
  leading_bp = 0
  trailing_bp = 0
  if options.leading_bp:
    if options.leading_bp >= 0:
      leading_bp = options.leading_bp
    else:
      print >> sys.stderr, "failed: leading_bp must be positive"
      exit(5)
  if options.trailing_bp:
    if  options.trailing_bp >= 0:
      trailing_bp = options.trailing_bp
    else:
      print >> sys.stderr, "failed: trailing_bp must be positive"
      exit(5)
  # Scan bed file 
  try:
    for i, line in enumerate( inputFile ):
      if line.startswith('track'):
        if outFile and options.bed_format:
          outFile.write(line)
        continue
      entry = BedEntry(line)
      strand = 1 if entry.strand == '+' else -1
      translations = entry.get_translations()
      if options.debug:
        exon_seqs = entry.get_exon_seqs()
        exon_sizes = [len(seq) for seq in exon_seqs]
        splice_sites = [sum(exon_sizes[:x]) / 3 for x in range(1,len(exon_sizes))]
        print >> sys.stderr, entry.name
        print >> sys.stderr, line.rstrip('\r\n')
        print >> sys.stderr, "exons:  %s" % exon_seqs
        print >> sys.stderr, "%s" % splice_sites
        for i,translation in enumerate(translations):
          print >> sys.stderr, "frame %d:  %s" % (i+1,translation)
          print >> sys.stderr, "splice:   %s" % (''.join(['^' if (((j*3)+i)/3) in splice_sites else '-' for j in range(len(translation))]))
        print >> sys.stderr, ""
      if options.bed_format:
        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)
        for i,tx in enumerate(translations):
          if tx:
            (chromStart,chromEnd,translation) = 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)
            location = "chromosome:%s:%s:%s:%s:%s" % (options.reference,entry.chrom,chromStart,chromEnd,strand)
            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')
  except Exception, e:
    print >> sys.stderr, "failed: Error reading %s - %s" % (options.input if options.input else 'stdin',e)

if __name__ == "__main__" : __main__()