Mercurial > repos > galaxyp > map_peptides_to_bed
diff map_peptides_to_bed.py @ 3:704ea6303c4c draft default tip
"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 11:41:15 -0400 |
| parents | db90662d26f9 |
| children |
line wrap: on
line diff
--- a/map_peptides_to_bed.py Thu Dec 15 18:36:55 2016 -0500 +++ b/map_peptides_to_bed.py Tue Apr 07 11:41:15 2020 -0400 @@ -10,324 +10,342 @@ # James E Johnson # #------------------------------------------------------------------------------ -""" - -""" Input: list of protein_accessions, peptide_sequence - GFF3 with fasta + GFF3 with fasta Output: GFF3 of peptides Filter: Must cross splice boundary - """ -import sys,re,os.path -import tempfile import optparse -from optparse import OptionParser -from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate +import os.path +import sys + +from Bio.Seq import ( + reverse_complement, + translate +) + + +class BedEntry(object): + def __init__(self, line): + self.line = line + try: + 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) + 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 as e: + sys.stderr.write("Unable to read Bed entry %s \n" % 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): + 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 -class BedEntry( object ): - def __init__(self, line): - self.line = line - try: - 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) - 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 __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): - 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) - 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] - ## print >> sys.stderr, "%d chromStart: %d chromEnd: %s bStart: %s bEnd: %d" % (x,chromStart,chromEnd,bStart,bEnd) - if bStart <= tstart < bEnd: - if self.strand == '+': - chromStart = self.chromStart + self.blockStarts[x] + (tstart - bStart) - else: - chromEnd = self.chromStart + self.blockStarts[x] + self.blockSizes[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,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))] - 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): - 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 - print >> sys.stderr, "tstart: %d tstop: %d offset: %d" % (tstart,tstop,offset) - 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) - #get the blocks for this translation - tblockCount = 0 - tblockSizes = [] - tblockStarts = [] - for x in range(self.blockCount): + 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 = int(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 + + def get_subrange(self, tstart, tstop): + """ + (start, end) + """ + chromStart = self.chromStart + chromEnd = self.chromEnd + r = range(self.blockCount) + if self.strand == '-': + r = list(r) + r.reverse() + bStart = 0 + for x in r: + bEnd = bStart + self.blockSizes[x] + # print >> sys.stderr, "%d chromStart: %d chromEnd: %s bStart: %s bEnd: %d\n" % (x, chromStart, chromEnd, bStart, bEnd) + if bStart <= tstart < bEnd: + if self.strand == '+': + chromStart = self.chromStart + self.blockStarts[x] + (tstart - bStart) + else: + chromEnd = self.chromStart + self.blockStarts[x] + self.blockSizes[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) + + def get_blocks(self, chromStart, chromEnd): + """ + get the blocks for sub range + """ + 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 + break if bEnd < chromStart: - continue - cStart = max(chromStart,bStart) + continue + cStart = max(chromStart, bStart) tblockStarts.append(cStart - chromStart) - tblockSizes.append(min(chromEnd,bEnd) - cStart) + tblockSizes.append(min(chromEnd, bEnd) - cStart) tblockCount += 1 - 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 - # >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 + sys.stderr.write("tblockCount: %d tblockStarts: %s tblockSizes: %s\n" % (tblockCount, tblockStarts, tblockSizes)) + return (tblockCount, tblockSizes, tblockStarts) + + def get_filterd_translations(self, untrimmed=False, filtering=True, ignore_left_bp=0, ignore_right_bp=0): + """ + [[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 + """ + translations = [None, None, None, None, None, None] + seq = self.get_spliced_seq() + ignore = int((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 = [int(sum(exon_sizes[:x]) / 3) for x in range(1, len(exon_sizes))] + sys.stderr.write("splice_sites: %s\n" % splice_sites) + 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 + sys.stderr.write("tstart: %d tstop: %d offset: %d\n" % (tstart, tstop, offset)) + 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) + # get the blocks for this translation + 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 + sys.stderr.write("tblockCount: %d tblockStarts: %s tblockSizes: %s\n" % (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 + >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 is not 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( '-t', '--translated_bed', dest='translated_bed', default=None, help='A bed file with added 13th column having a translation' ) - parser.add_option( '-i', '--input', dest='input', default=None, help='Tabular file with peptide_sequence column' ) - parser.add_option( '-p', '--peptide_column', type='int', dest='peptide_column', default=1, help='column ordinal with peptide sequence' ) - parser.add_option( '-n', '--name_column', type='int', dest='name_column', default=2, help='column ordinal with protein name' ) - parser.add_option( '-s', '--start_column', type='int', dest='start_column', default=None, help='column with peptide start position in protein' ) - 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( '-f', '--fasta', dest='fasta', default=None, help='Protein fasta' ) - parser.add_option( '-T', '--gffTags', dest='gffTags', action='store_true', default=False, help='Add #gffTags to bed output for IGV' ) - parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stderr' ) - (options, args) = parser.parse_args() - # Input files - if options.input != None: + # Parse Command Line + parser = optparse.OptionParser() + parser.add_option('-t', '--translated_bed', dest='translated_bed', default=None, help='A bed file with added 13th column having a translation') + parser.add_option('-i', '--input', dest='input', default=None, help='Tabular file with peptide_sequence column') + parser.add_option('-p', '--peptide_column', type='int', dest='peptide_column', default=1, help='column ordinal with peptide sequence') + parser.add_option('-n', '--name_column', type='int', dest='name_column', default=2, help='column ordinal with protein name') + parser.add_option('-s', '--start_column', type='int', dest='start_column', default=None, help='column with peptide start position in protein') + 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('-f', '--fasta', dest='fasta', default=None, help='Protein fasta') + parser.add_option('-T', '--gffTags', dest='gffTags', action='store_true', default=False, help='Add #gffTags to bed output for IGV') + parser.add_option('-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stderr') + (options, args) = parser.parse_args() + # Input files + if options.input is not None: + try: + inputPath = os.path.abspath(options.input) + inputFile = open(inputPath, 'r') + except Exception as e: + sys.stderr("failed: %s\n" % e) + exit(2) + else: + inputFile = sys.stdin + inputBed = None + if options.translated_bed is not None: + inputBed = open(os.path.abspath(options.translated_bed), 'r') + peptide_column = options.peptide_column - 1 + name_column = options.name_column - 1 if options.name_column else None + start_column = options.start_column - 1 if options.start_column else None + # Read in peptides + # peps[prot_name] = [seq] + prot_peps = dict() + unassigned_peps = set() 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 - inputBed = None - if options.translated_bed != None: - inputBed = open(os.path.abspath(options.translated_bed),'r') - peptide_column = options.peptide_column - 1 - name_column = options.name_column - 1 if options.name_column else None - start_column = options.start_column - 1 if options.start_column else None - # Read in peptides - # peps[prot_name] = [seq] - prot_peps = dict() - unassigned_peps = set() - try: - for i, line in enumerate( inputFile ): - ## print >> sys.stderr, "%3d\t%s" % (i,line) - if line.startswith('#'): - continue - fields = line.rstrip('\r\n').split('\t') - ## print >> sys.stderr, "%3d\t%s" % (i,fields) - if peptide_column < len(fields): - peptide = fields[peptide_column] - prot_name = fields[name_column] if name_column is not None and name_column < len(fields) else None - if prot_name: - offset = fields[start_column] if start_column is not None and start_column < len(fields) else -1 - if prot_name not in prot_peps: - prot_peps[prot_name] = dict() - prot_peps[prot_name][peptide] = offset - else: - unassigned_peps.add(peptide) - if options.debug: - print >> sys.stderr, "prot_peps: %s" % prot_peps - print >> sys.stderr, "unassigned_peps: %s" % unassigned_peps - except Exception, e: - print >> sys.stderr, "failed: Error reading %s - %s" % (options.input if options.input else 'stdin',e) - exit(1) - # Output files - bed_fh = None - ## gff_fh = None - ## gff_fa_file = None - gff_fa = None - outFile = None - if options.bed: - bed_fh = open(options.bed,'w') - bed_fh.write('track name="%s" type=bedDetail description="%s" \n' % ('novel_junction_peptides','test')) - if options.gffTags: - bed_fh.write('#gffTags\n') - ## 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)) - try: - for i, line in enumerate( inputBed ): - ## print >> sys.stderr, "%3d:\t%s" % (i,line) - if line.startswith('track'): - continue - entry = BedEntry(line) - if entry.name in prot_peps: - for (peptide,offset) in prot_peps[entry.name].iteritems(): - if offset < 0: - offset = entry.seq.find(peptide) - if options.debug: - print >> sys.stderr, "%s\t%s\t%d\t%s\n" % (entry.name, peptide,offset,entry.seq) - if offset >= 0: - tstart = offset * 3 - tstop = tstart + len(peptide) * 3 - if options.debug: - print >> sys.stderr, "%d\t%d\t%d" % (offset,tstart,tstop) - (pepStart,pepEnd) = entry.get_subrange(tstart,tstop) - if options.debug: - print >> sys.stderr, "%d\t%d\t%d" % (offset,pepStart,pepEnd) - if bed_fh: - entry.thickStart = pepStart - entry.thickEnd = pepEnd - bedfields = str(entry).split('\t') - if options.gffTags: - bedfields[3] = "ID=%s;Name=%s" % (entry.name,peptide) - bed_fh.write("%s\t%s\t%s\n" % ('\t'.join(bedfields[:12]),peptide,entry.seq)) - except Exception, e: - print >> sys.stderr, "failed: Error reading %s - %s" % (options.input if options.input else 'stdin',e) + for i, line in enumerate(inputFile): + # print >> sys.stderr, "%3d\t%s\n" % (i, line) + if line.startswith('#'): + continue + fields = line.rstrip('\r\n').split('\t') + # print >> sys.stderr, "%3d\t%s\n" % (i, fields) + if peptide_column < len(fields): + peptide = fields[peptide_column] + prot_name = fields[name_column] if name_column is not None and name_column < len(fields) else None + if prot_name: + offset = fields[start_column] if start_column is not None and start_column < len(fields) else -1 + if prot_name not in prot_peps: + prot_peps[prot_name] = dict() + prot_peps[prot_name][peptide] = offset + else: + unassigned_peps.add(peptide) + if options.debug: + sys.stderr.write("prot_peps: %s\n" % prot_peps) + sys.stderr.write("unassigned_peps: %s\n" % unassigned_peps) + except Exception as e: + sys.stderr.write("failed: Error reading %s - %s\n" % (options.input if options.input else 'stdin', e)) + exit(1) + # Output files + bed_fh = None + if options.bed: + bed_fh = open(options.bed, 'w') + bed_fh.write('track name="%s" type=bedDetail description="%s" \n' % ('novel_junction_peptides', 'test')) + if options.gffTags: + bed_fh.write('#gffTags\n') + # 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)) + try: + for i, line in enumerate(inputBed): + # print >> sys.stderr, "%3d:\t%s\n" % (i, line) + if line.startswith('track'): + continue + entry = BedEntry(line) + if entry.name in prot_peps: + for (peptide, offset) in prot_peps[entry.name].items(): + if offset < 0: + offset = entry.seq.find(peptide) + if options.debug: + sys.stderr.write("%s\t%s\t%d\t%s\n" % (entry.name, peptide, offset, entry.seq)) + if offset >= 0: + tstart = offset * 3 + tstop = tstart + len(peptide) * 3 + if options.debug: + sys.stderr.write("%d\t%d\t%d\n" % (offset, tstart, tstop)) + (pepStart, pepEnd) = entry.get_subrange(tstart, tstop) + if options.debug: + sys.stderr.write("%d\t%d\t%d\n" % (offset, pepStart, pepEnd)) + if bed_fh: + entry.thickStart = pepStart + entry.thickEnd = pepEnd + bedfields = str(entry).split('\t') + if options.gffTags: + bedfields[3] = "ID=%s;Name=%s" % (entry.name, peptide) + bed_fh.write("%s\t%s\t%s\n" % ('\t'.join(bedfields[:12]), peptide, entry.seq)) + except Exception as e: + sys.stderr.write("failed: Error reading %s - %s\n" % (options.input if options.input else 'stdin', e)) + raise -if __name__ == "__main__" : __main__() +if __name__ == "__main__": + __main__()