# HG changeset patch
# User jjohnson
# Date 1452627483 18000
# Node ID c626a939eef731fcb451ce9967003f00c4c1fab4
# Parent  aa93f79102597fdf87004d8870ebf9bef0156cc9
Uploaded

diff -r aa93f7910259 -r c626a939eef7 test-data/._translated_bed_sequences.fa
Binary file test-data/._translated_bed_sequences.fa has changed
diff -r aa93f7910259 -r c626a939eef7 tool_dependencies.xml
--- a/tool_dependencies.xml	Thu Jan 30 13:26:58 2014 -0600
+++ b/tool_dependencies.xml	Tue Jan 12 14:38:03 2016 -0500
@@ -1,6 +1,6 @@
 <?xml version="1.0"?>
 <tool_dependency>
     <package name="biopython" version="1.62">
-        <repository changeset_revision="ac9cc2992b69" name="package_biopython_1_62" owner="biopython" toolshed="http://testtoolshed.g2.bx.psu.edu" />
+        <repository changeset_revision="6b6a336db91a" name="package_biopython_1_62" owner="biopython" toolshed="https://testtoolshed.g2.bx.psu.edu" />
     </package>
 </tool_dependency>
diff -r aa93f7910259 -r c626a939eef7 translate_bed_sequences.py
--- 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)
 
diff -r aa93f7910259 -r c626a939eef7 translate_bed_sequences.xml
--- a/translate_bed_sequences.xml	Thu Jan 30 13:26:58 2014 -0600
+++ b/translate_bed_sequences.xml	Tue Jan 12 14:38:03 2016 -0500
@@ -1,16 +1,29 @@
 <?xml version="1.0"?>
-<tool id="translate_bed_sequences" name="Translate BED Sequences" version="0.0.1">
+<tool id="translate_bed_sequences" name="Translate BED Sequences" version="0.1.0">
   <description>3 frame translation of BED augmented with a sequence column</description>
   <requirements>
     <requirement type="package" version="1.62">biopython</requirement>
     <requirement type="python-module">Bio</requirement>
   </requirements>
-  <command interpreter="python">translate_bed_sequences.py  --input "$input" 
+  <command interpreter="python">
+  translate_bed_sequences.py  --input "$input" 
+  #if $fa_db:
+   --fa_db='$fa_db'
+  #end if
+  #if $fa_sep:
+   --fa_sep='$fa_sep'
+  #end if
+  #if $id_prefix:
+   --id_prefix='$id_prefix'
+  #end if
   #if $reference:
    --reference $reference
   #else:
    --reference ${input.metadata.dbkey}
   #end if
+  #if $refsource:
+   --refsource $refsource
+  #end if
   #if $seqtype:
     --seqtype $seqtype
   #end if
@@ -29,18 +42,31 @@
   #end if
   #if $trim.trimseqs == 'no':
     --untrimmed
-    #if $trim.max_stop_codons.__str__ != '':
+    #if str($trim.max_stop_codons) != '':
       --max_stop_codons $trim.max_stop_codons
     #end if
   #end if
-  #if $min_length:
+  #if str($min_length) != '':
    --min_length $min_length 
   #end if
+  --bed $translated_bed
   --output "$output"
   </command>
   <inputs>
     <param name="input" type="data" format="bed" label="BED file with added sequence column" 
            help="Output from 'Extract Genomic DNA' run on tophat junctions.bed "/> 
+    <param name="fa_db" type="text" value="" optional="true" label="fasta ID source, e.g. generic"
+           help="Any Compomics application such as PeptideShaker, requires a source">
+    </param>
+    <param name="fa_sep" type="text" value="" optional="true" label="fasta ID source, e.g. generic"
+           help="Only used when a fasta ID source is given, default to the pipe character">
+    </param>
+    <param name="id_prefix" type="text" value="" optional="true" label="ID prefix for generated IDs"
+           help="Can be used to distinguish samples">
+        <validator type="regex" message="Allowed chars:a-z A-Z 0-9 _ - |">^[a-zA-Z0-9_-|]*$</validator>
+    </param>
+    <param name="refsource" type="text" value="Ensembl" optional="true" label="Genome reference source"
+           help=""/>
     <param name="reference" type="text" value="" optional="true" label="Genome reference name"
            help="By default, the database metadata will be used."/>
     <param name="seqtype" type="text" value="" optional="true" label="The SEQTYPE:STATUS to include in the fasta ID lines"
@@ -75,8 +101,9 @@
     <exit_code range="1:" level="fatal" description="Error" />
   </stdio>
   <outputs>
-    <data name="output" metadata_source="input" format="fasta" label="${tool.name} on ${on_string}">
-      <filter>'found' in str(outputs)</filter>
+    <data name="translated_bed" metadata_source="input" format="bed" label="${tool.name} on ${on_string} bed">
+    </data>
+    <data name="output" metadata_source="input" format="fasta" label="${tool.name} on ${on_string} fasta">
     </data>
   </outputs>
   <tests>
@@ -97,6 +124,5 @@
 It generates a peptide fasta file with the 3-frame translations of the spliced sequence 
 defined by each entry in the input BED file.
 
-
   </help>
 </tool>