Mercurial > repos > galaxyp > map_peptides_to_bed
view map_peptides_to_bed.py @ 0:51f8f9041724 draft
planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/map_peptides_to_bed commit e04ed4b4960d6109a85c1cc68a2bf4931c8751ef-dirty
author | galaxyp |
---|---|
date | Mon, 25 Jan 2016 15:32:49 -0500 |
parents | |
children | db90662d26f9 |
line wrap: on
line source
#!/usr/bin/env python """ # #------------------------------------------------------------------------------ # University of Minnesota # Copyright 2014, Regents of the University of Minnesota #------------------------------------------------------------------------------ # Author: # # James E Johnson # #------------------------------------------------------------------------------ """ """ Input: list of protein_accessions, peptide_sequence 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 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): 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) 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 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=None, 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: 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) if __name__ == "__main__" : __main__()