Mercurial > repos > galaxyp > map_peptides_to_bed

comparison map_peptides_to_bed.py @ 6:aac354f8d5e8 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/map_peptides_to_bed commit 2a470e2c775a7427aa530e058510e4dc7b6d8e80"
author galaxyp
date Tue, 07 Apr 2020 15:31:21 +0000
parents 88d2c4c2cb0a
children
comparison
equal deleted inserted replaced
5:8175366a591f 6:aac354f8d5e8
8 # Author: 8 # Author:
9 # 9 #
10 # James E Johnson 10 # James E Johnson
11 # 11 #
12 #------------------------------------------------------------------------------ 12 #------------------------------------------------------------------------------
13 Input: list of protein_accessions, peptide_sequence
14 GFF3 with fasta
15 Output: GFF3 of peptides
16
17 Filter: Must cross splice boundary
13 """ 18 """
14 19
15 """
16 Input: list of protein_accessions, peptide_sequence
17 GFF3 with fasta
18 Output: GFF3 of peptides
19
20 Filter: Must cross splice boundary
21
22 """
23
24 import sys,re,os.path
25 import tempfile
26 import optparse 20 import optparse
27 from optparse import OptionParser 21 import os.path
28 from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate 22 import sys
29 23
30 class BedEntry( object ): 24 from Bio.Seq import (
31 def __init__(self, line): 25 reverse_complement,
32 self.line = line 26 translate
33 try: 27 )
34 fields = line.rstrip('\r\n').split('\t') 28
35 (chrom,chromStart,chromEnd,name,score,strand,thickStart,thickEnd,itemRgb,blockCount,blockSizes,blockStarts) = fields[0:12] 29
36 seq = fields[12] if len(fields) > 12 else None 30 class BedEntry(object):
37 self.chrom = chrom 31 def __init__(self, line):
38 self.chromStart = int(chromStart) 32 self.line = line
39 self.chromEnd = int(chromEnd) 33 try:
40 self.name = name 34 fields = line.rstrip('\r\n').split('\t')
41 self.score = int(score) 35 (chrom, chromStart, chromEnd, name, score, strand, thickStart, thickEnd, itemRgb, blockCount, blockSizes, blockStarts) = fields[0:12]
42 self.strand = strand 36 seq = fields[12] if len(fields) > 12 else None
43 self.thickStart = int(thickStart) 37 self.chrom = chrom
44 self.thickEnd = int(thickEnd) 38 self.chromStart = int(chromStart)
45 self.itemRgb = itemRgb 39 self.chromEnd = int(chromEnd)
46 self.blockCount = int(blockCount) 40 self.name = name
47 self.blockSizes = [int(x) for x in blockSizes.split(',')] 41 self.score = int(score)
48 self.blockStarts = [int(x) for x in blockStarts.split(',')] 42 self.strand = strand
49 self.seq = seq 43 self.thickStart = int(thickStart)
50 except Exception, e: 44 self.thickEnd = int(thickEnd)
51 print >> sys.stderr, "Unable to read Bed entry" % e 45 self.itemRgb = itemRgb
52 exit(1) 46 self.blockCount = int(blockCount)
53 def __str__(self): 47 self.blockSizes = [int(x) for x in blockSizes.split(',')]
54 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' % ( 48 self.blockStarts = [int(x) for x in blockStarts.split(',')]
55 self.chrom, self.chromStart, self.chromEnd, self.name, self.score, self.strand, self.thickStart, self.thickEnd, self.itemRgb, self.blockCount, 49 self.seq = seq
56 ','.join([str(x) for x in self.blockSizes]), 50 except Exception as e:
57 ','.join([str(x) for x in self.blockStarts]), 51 sys.stderr.write("Unable to read Bed entry %s \n" % e)
58 '\t%s' % self.seq if self.seq else '') 52 exit(1)
59 def get_splice_junctions(self): 53
60 splice_juncs = [] 54 def __str__(self):
61 for i in range(self.blockCount - 1): 55 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' % (
62 splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1]) 56 self.chrom, self.chromStart, self.chromEnd, self.name, self.score, self.strand, self.thickStart, self.thickEnd, self.itemRgb, self.blockCount,
63 splice_juncs.append(splice_junc) 57 ','.join([str(x) for x in self.blockSizes]),
64 return splice_juncs 58 ','.join([str(x) for x in self.blockStarts]),
65 def get_exon_seqs(self): 59 '\t%s' % self.seq if self.seq else '')
66 exons = [] 60
67 for i in range(self.blockCount): 61 def get_splice_junctions(self):
68 # splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1]) 62 splice_juncs = []
69 exons.append(self.seq[self.blockStarts[i]:self.blockStarts[i] + self.blockSizes[i]]) 63 for i in range(self.blockCount - 1):
70 if self.strand == '-': #reverse complement 64 splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i + 1])
71 exons.reverse() 65 splice_juncs.append(splice_junc)
72 for i,s in enumerate(exons): 66 return splice_juncs
73 exons[i] = reverse_complement(s) 67
74 return exons 68 def get_exon_seqs(self):
75 def get_spliced_seq(self): 69 exons = []
76 seq = ''.join(self.get_exon_seqs()) 70 for i in range(self.blockCount):
77 return seq 71 # splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1])
78 def get_translation(self,sequence=None): 72 exons.append(self.seq[self.blockStarts[i]:self.blockStarts[i] + self.blockSizes[i]])
79 translation = None 73 if self.strand == '-': # reverse complement
80 seq = sequence if sequence else self.get_spliced_seq() 74 exons.reverse()
81 if seq: 75 for i, s in enumerate(exons):
82 seqlen = len(seq) / 3 * 3; 76 exons[i] = reverse_complement(s)
83 if seqlen >= 3: 77 return exons
84 translation = translate(seq[:seqlen]) 78
85 return translation 79 def get_spliced_seq(self):
86 def get_translations(self): 80 seq = ''.join(self.get_exon_seqs())
87 translations = [] 81 return seq
88 seq = self.get_spliced_seq() 82
89 if seq: 83 def get_translation(self, sequence=None):
90 for i in range(3): 84 translation = None
91 translation = self.get_translation(sequence=seq[i:]) 85 seq = sequence if sequence else self.get_spliced_seq()
92 if translation: 86 if seq:
93 translations.append(translation) 87 seqlen = int(len(seq) / 3) * 3
94 return translations 88 if seqlen >= 3:
95 ## (start,end) 89 translation = translate(seq[:seqlen])
96 def get_subrange(self,tstart,tstop): 90 return translation
97 chromStart = self.chromStart 91
98 chromEnd = self.chromEnd 92 def get_translations(self):
99 r = range(self.blockCount) 93 translations = []
100 if self.strand == '-': 94 seq = self.get_spliced_seq()
101 r.reverse() 95 if seq:
102 bStart = 0 96 for i in range(3):
103 for x in r: 97 translation = self.get_translation(sequence=seq[i:])
104 bEnd = bStart + self.blockSizes[x] 98 if translation:
105 ## print >> sys.stderr, "%d chromStart: %d chromEnd: %s bStart: %s bEnd: %d" % (x,chromStart,chromEnd,bStart,bEnd) 99 translations.append(translation)
106 if bStart <= tstart < bEnd: 100 return translations
107 if self.strand == '+': 101
108 chromStart = self.chromStart + self.blockStarts[x] + (tstart - bStart) 102 def get_subrange(self, tstart, tstop):
109 else: 103 """
110 chromEnd = self.chromStart + self.blockStarts[x] + self.blockSizes[x] - (tstart - bStart) 104 (start, end)
111 if bStart <= tstop < bEnd: 105 """
112 if self.strand == '+': 106 chromStart = self.chromStart
113 chromEnd = self.chromStart + self.blockStarts[x] + (tstop - bStart) 107 chromEnd = self.chromEnd
114 else: 108 r = range(self.blockCount)
115 chromStart = self.chromStart + self.blockStarts[x] + self.blockSizes[x] - (tstop - bStart) 109 if self.strand == '-':
116 bStart += self.blockSizes[x] 110 r = list(r)
117 return(chromStart,chromEnd) 111 r.reverse()
118 #get the blocks for sub range 112 bStart = 0
119 def get_blocks(self,chromStart,chromEnd): 113 for x in r:
120 tblockCount = 0 114 bEnd = bStart + self.blockSizes[x]
121 tblockSizes = [] 115 # print >> sys.stderr, "%d chromStart: %d chromEnd: %s bStart: %s bEnd: %d\n" % (x, chromStart, chromEnd, bStart, bEnd)
122 tblockStarts = [] 116 if bStart <= tstart < bEnd:
123 for x in range(self.blockCount): 117 if self.strand == '+':
124 bStart = self.chromStart + self.blockStarts[x] 118 chromStart = self.chromStart + self.blockStarts[x] + (tstart - bStart)
125 bEnd = bStart + self.blockSizes[x] 119 else:
126 if bStart > chromEnd: 120 chromEnd = self.chromStart + self.blockStarts[x] + self.blockSizes[x] - (tstart - bStart)
127 break 121 if bStart <= tstop < bEnd:
128 if bEnd < chromStart: 122 if self.strand == '+':
129 continue 123 chromEnd = self.chromStart + self.blockStarts[x] + (tstop - bStart)
130 cStart = max(chromStart,bStart) 124 else:
131 tblockStarts.append(cStart - chromStart) 125 chromStart = self.chromStart + self.blockStarts[x] + self.blockSizes[x] - (tstop - bStart)
132 tblockSizes.append(min(chromEnd,bEnd) - cStart) 126 bStart += self.blockSizes[x]
133 tblockCount += 1 127 return(chromStart, chromEnd)
134 print >> sys.stderr, "tblockCount: %d tblockStarts: %s tblockSizes: %s" % (tblockCount,tblockStarts,tblockSizes) 128
135 return (tblockCount,tblockSizes,tblockStarts) 129 def get_blocks(self, chromStart, chromEnd):
136 130 """
137 ## [[start,end,seq,blockCount,blockSizes,blockStarts],[start,end,seq,blockCount,blockSizes,blockStarts],[start,end,seq,blockCount,blockSizes,blockStarts]] 131 get the blocks for sub range
138 ## filter: ignore translation if stop codon in first exon after ignore_left_bp 132 """
139 def get_filterd_translations(self,untrimmed=False,filtering=True,ignore_left_bp=0,ignore_right_bp=0): 133 tblockCount = 0
140 translations = [None,None,None,None,None,None] 134 tblockSizes = []
141 seq = self.get_spliced_seq() 135 tblockStarts = []
142 ignore = (ignore_left_bp if self.strand == '+' else ignore_right_bp) / 3 136 for x in range(self.blockCount):
143 block_sum = sum(self.blockSizes)
144 exon_sizes = self.blockSizes
145 if self.strand == '-':
146 exon_sizes.reverse()
147 splice_sites = [sum(exon_sizes[:x]) / 3 for x in range(1,len(exon_sizes))]
148 print >> sys.stderr, "splice_sites: %s" % splice_sites
149 junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0]
150 if seq:
151 for i in range(3):
152 translation = self.get_translation(sequence=seq[i:])
153 if translation:
154 tstart = 0
155 tstop = len(translation)
156 if not untrimmed:
157 tstart = translation.rfind('*',0,junc) + 1
158 stop = translation.find('*',junc)
159 tstop = stop if stop >= 0 else len(translation)
160 if filtering and tstart > ignore:
161 continue
162 trimmed = translation[tstart:tstop]
163 #get genomic locations for start and end
164 offset = (block_sum - i) % 3
165 print >> sys.stderr, "tstart: %d tstop: %d offset: %d" % (tstart,tstop,offset)
166 if self.strand == '+':
167 chromStart = self.chromStart + i + (tstart * 3)
168 chromEnd = self.chromEnd - offset - (len(translation) - tstop) * 3
169 else:
170 chromStart = self.chromStart + offset + (len(translation) - tstop) * 3
171 chromEnd = self.chromEnd - i - (tstart * 3)
172 #get the blocks for this translation
173 tblockCount = 0
174 tblockSizes = []
175 tblockStarts = []
176 for x in range(self.blockCount):
177 bStart = self.chromStart + self.blockStarts[x] 137 bStart = self.chromStart + self.blockStarts[x]
178 bEnd = bStart + self.blockSizes[x] 138 bEnd = bStart + self.blockSizes[x]
179 if bStart > chromEnd: 139 if bStart > chromEnd:
180 break 140 break
181 if bEnd < chromStart: 141 if bEnd < chromStart:
182 continue 142 continue
183 cStart = max(chromStart,bStart) 143 cStart = max(chromStart, bStart)
184 tblockStarts.append(cStart - chromStart) 144 tblockStarts.append(cStart - chromStart)
185 tblockSizes.append(min(chromEnd,bEnd) - cStart) 145 tblockSizes.append(min(chromEnd, bEnd) - cStart)
186 tblockCount += 1 146 tblockCount += 1
187 print >> sys.stderr, "tblockCount: %d tblockStarts: %s tblockSizes: %s" % (tblockCount,tblockStarts,tblockSizes) 147 sys.stderr.write("tblockCount: %d tblockStarts: %s tblockSizes: %s\n" % (tblockCount, tblockStarts, tblockSizes))
188 translations[i] = [chromStart,chromEnd,trimmed,tblockCount,tblockSizes,tblockStarts] 148 return (tblockCount, tblockSizes, tblockStarts)
189 return translations 149
190 def get_seq_id(self,seqtype='unk:unk',reference='',frame=None): 150 def get_filterd_translations(self, untrimmed=False, filtering=True, ignore_left_bp=0, ignore_right_bp=0):
191 ## Ensembl fasta ID format 151 """
192 # >ID SEQTYPE:STATUS LOCATION GENE TRANSCRIPT 152 [[start, end, seq, blockCount, blockSizes, blockStarts], [start, end, seq, blockCount, blockSizes, blockStarts], [start, end, seq, blockCount, blockSizes, blockStarts]]
193 # >ENSP00000328693 pep:splice chromosome:NCBI35:1:904515:910768:1 gene:ENSG00000158815:transcript:ENST00000328693 gene_biotype:protein_coding transcript_biotype:protein_coding 153 filter: ignore translation if stop codon in first exon after ignore_left_bp
194 frame_name = '' 154 """
195 chromStart = self.chromStart 155 translations = [None, None, None, None, None, None]
196 chromEnd = self.chromEnd 156 seq = self.get_spliced_seq()
197 strand = 1 if self.strand == '+' else -1 157 ignore = int((ignore_left_bp if self.strand == '+' else ignore_right_bp) / 3)
198 if frame != None: 158 block_sum = sum(self.blockSizes)
199 block_sum = sum(self.blockSizes) 159 exon_sizes = self.blockSizes
200 offset = (block_sum - frame) % 3 160 if self.strand == '-':
201 frame_name = '_' + str(frame + 1) 161 exon_sizes.reverse()
202 if self.strand == '+': 162 splice_sites = [int(sum(exon_sizes[:x]) / 3) for x in range(1, len(exon_sizes))]
203 chromStart += frame 163 sys.stderr.write("splice_sites: %s\n" % splice_sites)
204 chromEnd -= offset 164 junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0]
205 else: 165 if seq:
206 chromStart += offset 166 for i in range(3):
207 chromEnd -= frame 167 translation = self.get_translation(sequence=seq[i:])
208 location = "chromosome:%s:%s:%s:%s:%s" % (reference,self.chrom,chromStart,chromEnd,strand) 168 if translation:
209 seq_id = "%s%s %s %s" % (self.name,frame_name,seqtype,location) 169 tstart = 0
210 return seq_id 170 tstop = len(translation)
211 def get_line(self, start_offset = 0, end_offset = 0): 171 if not untrimmed:
212 if start_offset or end_offset: 172 tstart = translation.rfind('*', 0, junc) + 1
213 s_offset = start_offset if start_offset else 0 173 stop = translation.find('*', junc)
214 e_offset = end_offset if end_offset else 0 174 tstop = stop if stop >= 0 else len(translation)
215 if s_offset > self.chromStart: 175 if filtering and tstart > ignore:
216 s_offset = self.chromStart 176 continue
217 chrStart = self.chromStart - s_offset 177 trimmed = translation[tstart:tstop]
218 chrEnd = self.chromEnd + e_offset 178 # get genomic locations for start and end
219 blkSizes = self.blockSizes 179 offset = (block_sum - i) % 3
220 blkSizes[0] += s_offset 180 sys.stderr.write("tstart: %d tstop: %d offset: %d\n" % (tstart, tstop, offset))
221 blkSizes[-1] += e_offset 181 if self.strand == '+':
222 blkStarts = self.blockStarts 182 chromStart = self.chromStart + i + (tstart * 3)
223 for i in range(1,self.blockCount): 183 chromEnd = self.chromEnd - offset - (len(translation) - tstop) * 3
224 blkStarts[i] += s_offset 184 else:
225 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])]] 185 chromStart = self.chromStart + offset + (len(translation) - tstop) * 3
226 return '\t'.join(items) + '\n' 186 chromEnd = self.chromEnd - i - (tstart * 3)
227 return self.line 187 # get the blocks for this translation
188 tblockCount = 0
189 tblockSizes = []
190 tblockStarts = []
191 for x in range(self.blockCount):
192 bStart = self.chromStart + self.blockStarts[x]
193 bEnd = bStart + self.blockSizes[x]
194 if bStart > chromEnd:
195 break
196 if bEnd < chromStart:
197 continue
198 cStart = max(chromStart, bStart)
199 tblockStarts.append(cStart - chromStart)
200 tblockSizes.append(min(chromEnd, bEnd) - cStart)
201 tblockCount += 1
202 sys.stderr.write("tblockCount: %d tblockStarts: %s tblockSizes: %s\n" % (tblockCount, tblockStarts, tblockSizes))
203 translations[i] = [chromStart, chromEnd, trimmed, tblockCount, tblockSizes, tblockStarts]
204 return translations
205
206 def get_seq_id(self, seqtype='unk:unk', reference='', frame=None):
207 """
208 # Ensembl fasta ID format
209 >ID SEQTYPE:STATUS LOCATION GENE TRANSCRIPT
210 >ENSP00000328693 pep:splice chromosome:NCBI35:1:904515:910768:1 gene:ENSG00000158815:transcript:ENST00000328693 gene_biotype:protein_coding transcript_biotype:protein_coding
211 """
212 frame_name = ''
213 chromStart = self.chromStart
214 chromEnd = self.chromEnd
215 strand = 1 if self.strand == '+' else -1
216 if frame is not None:
217 block_sum = sum(self.blockSizes)
218 offset = (block_sum - frame) % 3
219 frame_name = '_' + str(frame + 1)
220 if self.strand == '+':
221 chromStart += frame
222 chromEnd -= offset
223 else:
224 chromStart += offset
225 chromEnd -= frame
226 location = "chromosome:%s:%s:%s:%s:%s" % (reference, self.chrom, chromStart, chromEnd, strand)
227 seq_id = "%s%s %s %s" % (self.name, frame_name, seqtype, location)
228 return seq_id
229
230 def get_line(self, start_offset=0, end_offset=0):
231 if start_offset or end_offset:
232 s_offset = start_offset if start_offset else 0
233 e_offset = end_offset if end_offset else 0
234 if s_offset > self.chromStart:
235 s_offset = self.chromStart
236 chrStart = self.chromStart - s_offset
237 chrEnd = self.chromEnd + e_offset
238 blkSizes = self.blockSizes
239 blkSizes[0] += s_offset
240 blkSizes[-1] += e_offset
241 blkStarts = self.blockStarts
242 for i in range(1, self.blockCount):
243 blkStarts[i] += s_offset
244 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])]]
245 return '\t'.join(items) + '\n'
246 return self.line
247
228 248
229 def __main__(): 249 def __main__():
230 #Parse Command Line 250 # Parse Command Line
231 parser = optparse.OptionParser() 251 parser = optparse.OptionParser()
232 parser.add_option( '-t', '--translated_bed', dest='translated_bed', default=None, help='A bed file with added 13th column having a translation' ) 252 parser.add_option('-t', '--translated_bed', dest='translated_bed', default=None, help='A bed file with added 13th column having a translation')
233 parser.add_option( '-i', '--input', dest='input', default=None, help='Tabular file with peptide_sequence column' ) 253 parser.add_option('-i', '--input', dest='input', default=None, help='Tabular file with peptide_sequence column')
234 parser.add_option( '-p', '--peptide_column', type='int', dest='peptide_column', default=1, help='column ordinal with peptide sequence' ) 254 parser.add_option('-p', '--peptide_column', type='int', dest='peptide_column', default=1, help='column ordinal with peptide sequence')
235 parser.add_option( '-n', '--name_column', type='int', dest='name_column', default=2, help='column ordinal with protein name' ) 255 parser.add_option('-n', '--name_column', type='int', dest='name_column', default=2, help='column ordinal with protein name')
236 parser.add_option( '-s', '--start_column', type='int', dest='start_column', default=None, help='column with peptide start position in protein' ) 256 parser.add_option('-s', '--start_column', type='int', dest='start_column', default=None, help='column with peptide start position in protein')
237 parser.add_option( '-B', '--bed', dest='bed', default=None, help='Output a bed file with added 13th column having translation' ) 257 parser.add_option('-B', '--bed', dest='bed', default=None, help='Output a bed file with added 13th column having translation')
238 ## parser.add_option( '-G', '--gff3', dest='gff', default=None, help='Output translations to a GFF3 file' ) 258 # parser.add_option('-G', '--gff3', dest='gff', default=None, help='Output translations to a GFF3 file')
239 ## parser.add_option( '-f', '--fasta', dest='fasta', default=None, help='Protein fasta' ) 259 # parser.add_option('-f', '--fasta', dest='fasta', default=None, help='Protein fasta')
240 parser.add_option( '-T', '--gffTags', dest='gffTags', action='store_true', default=False, help='Add #gffTags to bed output for IGV' ) 260 parser.add_option('-T', '--gffTags', dest='gffTags', action='store_true', default=False, help='Add #gffTags to bed output for IGV')
241 parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stderr' ) 261 parser.add_option('-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stderr')
242 (options, args) = parser.parse_args() 262 (options, args) = parser.parse_args()
243 # Input files 263 # Input files
244 if options.input != None: 264 if options.input is not None:
265 try:
266 inputPath = os.path.abspath(options.input)
267 inputFile = open(inputPath, 'r')
268 except Exception as e:
269 sys.stderr("failed: %s\n" % e)
270 exit(2)
271 else:
272 inputFile = sys.stdin
273 inputBed = None
274 if options.translated_bed is not None:
275 inputBed = open(os.path.abspath(options.translated_bed), 'r')
276 peptide_column = options.peptide_column - 1
277 name_column = options.name_column - 1 if options.name_column else None
278 start_column = options.start_column - 1 if options.start_column else None
279 # Read in peptides
280 # peps[prot_name] = [seq]
281 prot_peps = dict()
282 unassigned_peps = set()
245 try: 283 try:
246 inputPath = os.path.abspath(options.input) 284 for i, line in enumerate(inputFile):
247 inputFile = open(inputPath, 'r') 285 # print >> sys.stderr, "%3d\t%s\n" % (i, line)
248 except Exception, e: 286 if line.startswith('#'):
249 print >> sys.stderr, "failed: %s" % e 287 continue
250 exit(2) 288 fields = line.rstrip('\r\n').split('\t')
251 else: 289 # print >> sys.stderr, "%3d\t%s\n" % (i, fields)
252 inputFile = sys.stdin 290 if peptide_column < len(fields):
253 inputBed = None 291 peptide = fields[peptide_column]
254 if options.translated_bed != None: 292 prot_name = fields[name_column] if name_column is not None and name_column < len(fields) else None
255 inputBed = open(os.path.abspath(options.translated_bed),'r') 293 if prot_name:
256 peptide_column = options.peptide_column - 1 294 offset = fields[start_column] if start_column is not None and start_column < len(fields) else -1
257 name_column = options.name_column - 1 if options.name_column else None 295 if prot_name not in prot_peps:
258 start_column = options.start_column - 1 if options.start_column else None 296 prot_peps[prot_name] = dict()
259 # Read in peptides 297 prot_peps[prot_name][peptide] = offset
260 # peps[prot_name] = [seq] 298 else:
261 prot_peps = dict() 299 unassigned_peps.add(peptide)
262 unassigned_peps = set() 300 if options.debug:
263 try: 301 sys.stderr.write("prot_peps: %s\n" % prot_peps)
264 for i, line in enumerate( inputFile ): 302 sys.stderr.write("unassigned_peps: %s\n" % unassigned_peps)
265 ## print >> sys.stderr, "%3d\t%s" % (i,line) 303 except Exception as e:
266 if line.startswith('#'): 304 sys.stderr.write("failed: Error reading %s - %s\n" % (options.input if options.input else 'stdin', e))
267 continue 305 exit(1)
268 fields = line.rstrip('\r\n').split('\t') 306 # Output files
269 ## print >> sys.stderr, "%3d\t%s" % (i,fields) 307 bed_fh = None
270 if peptide_column < len(fields): 308 if options.bed:
271 peptide = fields[peptide_column] 309 bed_fh = open(options.bed, 'w')
272 prot_name = fields[name_column] if name_column is not None and name_column < len(fields) else None 310 bed_fh.write('track name="%s" type=bedDetail description="%s" \n' % ('novel_junction_peptides', 'test'))
273 if prot_name: 311 if options.gffTags:
274 offset = fields[start_column] if start_column is not None and start_column < len(fields) else -1 312 bed_fh.write('#gffTags\n')
275 if prot_name not in prot_peps: 313 # if options.gff:
276 prot_peps[prot_name] = dict() 314 # gff_fh = open(options.gff, 'w')
277 prot_peps[prot_name][peptide] = offset 315 # gff_fh.write("##gff-version 3.2.1\n")
278 else: 316 # if options.reference:
279 unassigned_peps.add(peptide) 317 # gff_fh.write("##genome-build %s %s\n" % (options.refsource if options.refsource else 'unknown', options.reference))
280 if options.debug: 318 try:
281 print >> sys.stderr, "prot_peps: %s" % prot_peps 319 for i, line in enumerate(inputBed):
282 print >> sys.stderr, "unassigned_peps: %s" % unassigned_peps 320 # print >> sys.stderr, "%3d:\t%s\n" % (i, line)
283 except Exception, e: 321 if line.startswith('track'):
284 print >> sys.stderr, "failed: Error reading %s - %s" % (options.input if options.input else 'stdin',e) 322 continue
285 exit(1) 323 entry = BedEntry(line)
286 # Output files 324 if entry.name in prot_peps:
287 bed_fh = None 325 for (peptide, offset) in prot_peps[entry.name].items():
288 ## gff_fh = None 326 if offset < 0:
289 ## gff_fa_file = None 327 offset = entry.seq.find(peptide)
290 gff_fa = None 328 if options.debug:
291 outFile = None 329 sys.stderr.write("%s\t%s\t%d\t%s\n" % (entry.name, peptide, offset, entry.seq))
292 if options.bed: 330 if offset >= 0:
293 bed_fh = open(options.bed,'w') 331 tstart = offset * 3
294 bed_fh.write('track name="%s" type=bedDetail description="%s" \n' % ('novel_junction_peptides','test')) 332 tstop = tstart + len(peptide) * 3
295 if options.gffTags: 333 if options.debug:
296 bed_fh.write('#gffTags\n') 334 sys.stderr.write("%d\t%d\t%d\n" % (offset, tstart, tstop))
297 ## if options.gff: 335 (pepStart, pepEnd) = entry.get_subrange(tstart, tstop)
298 ## gff_fh = open(options.gff,'w') 336 if options.debug:
299 ## gff_fh.write("##gff-version 3.2.1\n") 337 sys.stderr.write("%d\t%d\t%d\n" % (offset, pepStart, pepEnd))
300 ## if options.reference: 338 if bed_fh:
301 ## gff_fh.write("##genome-build %s %s\n" % (options.refsource if options.refsource else 'unknown', options.reference)) 339 entry.thickStart = pepStart
302 try: 340 entry.thickEnd = pepEnd
303 for i, line in enumerate( inputBed ): 341 bedfields = str(entry).split('\t')
304 ## print >> sys.stderr, "%3d:\t%s" % (i,line) 342 if options.gffTags:
305 if line.startswith('track'): 343 bedfields[3] = "ID=%s;Name=%s" % (entry.name, peptide)
306 continue 344 bed_fh.write("%s\t%s\t%s\n" % ('\t'.join(bedfields[:12]), peptide, entry.seq))
307 entry = BedEntry(line) 345 except Exception as e:
308 if entry.name in prot_peps: 346 sys.stderr.write("failed: Error reading %s - %s\n" % (options.input if options.input else 'stdin', e))
309 for (peptide,offset) in prot_peps[entry.name].iteritems(): 347 raise
310 if offset < 0: 348
311 offset = entry.seq.find(peptide) 349
312 if options.debug: 350 if __name__ == "__main__":
313 print >> sys.stderr, "%s\t%s\t%d\t%s\n" % (entry.name, peptide,offset,entry.seq) 351 __main__()
314 if offset >= 0:
315 tstart = offset * 3
316 tstop = tstart + len(peptide) * 3
317 if options.debug:
318 print >> sys.stderr, "%d\t%d\t%d" % (offset,tstart,tstop)
319 (pepStart,pepEnd) = entry.get_subrange(tstart,tstop)
320 if options.debug:
321 print >> sys.stderr, "%d\t%d\t%d" % (offset,pepStart,pepEnd)
322 if bed_fh:
323 entry.thickStart = pepStart
324 entry.thickEnd = pepEnd
325 bedfields = str(entry).split('\t')
326 if options.gffTags:
327 bedfields[3] = "ID=%s;Name=%s" % (entry.name,peptide)
328 bed_fh.write("%s\t%s\t%s\n" % ('\t'.join(bedfields[:12]),peptide,entry.seq))
329 except Exception, e:
330 print >> sys.stderr, "failed: Error reading %s - %s" % (options.input if options.input else 'stdin',e)
331
332 if __name__ == "__main__" : __main__()
333