diff --git a/tool_conf.xml.main b/tool_conf.xml.main index 9ea4e1dbdcc..e294e57cc82 100644 --- a/tool_conf.xml.main +++ b/tool_conf.xml.main @@ -21,7 +21,7 @@
- @@ -52,8 +52,8 @@
- - + + - - - + + + + + + + + + +
diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 2d0abaf6fd2..1436b60bdb5 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -21,7 +21,7 @@
- @@ -52,8 +52,8 @@
- - + + @@ -73,17 +73,17 @@
- - - - - - - - - - - + + + + + + + + + + +
diff --git a/tools/extract/genebed_maf_to_fasta.py b/tools/extract/genebed_maf_to_fasta.py deleted file mode 100644 index 9809c5b55ed..00000000000 --- a/tools/extract/genebed_maf_to_fasta.py +++ /dev/null @@ -1,371 +0,0 @@ -#!/usr/bin/env python2.4 - -""" -Reads a gene BED and an indexed MAF. Produces a FASTA file containing -the aligned gene sequences, based upon the provided coordinates - -Alignment blocks are layered ontop of each other based upon score. - -usage: %prog dbkey_of_BED comma_separated_list_of_additional_dbkeys_to_extract comma_separated_list_of_indexed_maf_files input_gene_bed_file output_fasta_file cached|user -""" - -#Dan Blankenberg -import pkg_resources; pkg_resources.require( "bx-python" ) -import bx.align.maf -import bx.intervals.io -import bx.interval_index_file -import sys, os, tempfile, string - -MAF_LOCATION_FILE = "/depot/data2/galaxy/maf_index.loc" - -#an object corresponding to a reference layered alignment -class RegionAlignment( object ): - - DNA_COMPLEMENT = string.maketrans( "ACGTacgt", "TGCAtgca" ) - - def __init__( self, size, species = [] ): - self.size = size - self.sequences = {} - if not isinstance( species, list ): - species = [species] - for spec in species: - self.add_species( spec ) - - #add a species to the alignment - def add_species( self, species ): - #make temporary sequence files, need to delete when done - fd, file_path = tempfile.mkstemp() - self.sequences[species] = { 'file':os.fdopen( fd, 'w+' ), 'path':file_path } - self.sequences[species]['file'].write( "-" * self.size ) - - #returns the names for species found in alignment, skipping names as requested - def get_species_names( self, skip = [] ): - if not isinstance( skip, list ): skip = [skip] - names = self.sequences.keys() - for name in skip: - try: names.remove( name ) - except: pass - return names - - #returns the sequence for a species - def get_sequence( self, species ): - self.sequences[species]['file'].seek( 0 ) - return( self.sequences[species]['file'].read() ) - - #returns the reverse complement of the sequence for a species - def get_sequence_reverse_complement( self, species ): - complement = [base for base in self.get_sequence( species ).translate( self.DNA_COMPLEMENT )] - complement.reverse() - return "".join( complement ) - - #sets a position for a species - def set_position( self, index, species, base ): - if index >= self.size or index < 0: raise "Your index (%i) is out of range (0 - %i)." % ( index, self.size - 1 ) - if len(base) != 1: raise "A genomic position can only have a length of 1." - if species not in self.sequences.keys(): self.add_species( species ) - self.sequences[species]['file'].seek( index ) - self.sequences[species]['file'].write( base ) - - #Flush temp file of specified species, or all species - def flush( self, species = None ): - if species is None: - species = self.sequences.keys() - elif not isinstance( species, list ): - species = [species] - for spec in species: - self.sequences[spec]['file'].flush() - - #object cleanup, delete temporary files - def __del__( self ): - for species, sequence in self.sequences.items(): - sequence['file'].close() - os.unlink( sequence['path'] ) - - -class GenomicRegionAlignment( RegionAlignment ): - - def __init__( self, start, end, species = [] ): - RegionAlignment.__init__( self, end - start, species ) - self.start = start - self.end = end - -class SplicedAlignment( object ): - - DNA_COMPLEMENT = string.maketrans( "ACGTacgt", "TGCAtgca" ) - - def __init__( self, exon_starts, exon_ends, species = [] ): - if not isinstance( exon_starts, list ): - exon_starts = [exon_starts] - if not isinstance( exon_ends, list ): - exon_ends = [exon_ends] - assert len( exon_starts ) == len( exon_ends ), "The number of starts does not match the number of sizes." - self.exons = [] - for i in range( len( exon_starts ) ): - self.exons.append( GenomicRegionAlignment( exon_starts[i], exon_ends[i], species ) ) - - #returns the names for species found in alignment, skipping names as requested - def get_species_names( self, skip = [] ): - if not isinstance( skip, list ): skip = [skip] - names = [] - for exon in self.exons: - for name in exon.get_species_names( skip = skip ): - if name not in names: - names.append( name ) - return names - - #returns the sequence for a species - def get_sequence( self, species ): - sequence = tempfile.TemporaryFile() - for exon in self.exons: - if species in exon.get_species_names(): - sequence.write( exon.get_sequence( species ) ) - else: - sequence.write( "-" * exon.size ) - sequence.seek( 0 ) - return sequence.read() - - #returns the reverse complement of the sequence for a species - def get_sequence_reverse_complement( self, species ): - complement = [base for base in self.get_sequence( species ).translate( self.DNA_COMPLEMENT )] - complement.reverse() - return "".join( complement ) - - #Start and end of coding region - @property - def start( self ): - return self.exons[0].start - @property - def end( self ): - return self.exons[-1].end - -def maf_index_by_uid( maf_uid ): - for line in open( MAF_LOCATION_FILE ): - try: - #read each line, if not enough fields, go to next line - if line[0:1] == "#" : continue - fields = line.split('\t') - if maf_uid == fields[1]: - try: - maf_files = fields[3].replace( "\n", "" ).replace( "\r", "" ).split( "," ) - return bx.align.maf.MultiIndexed( maf_files, keep_open = True, parse_e_rows = True ) - except Exception, e: - raise 'MAF UID (%s) found, but configuration appears to be malformed: %s' % ( maf_uid, e ) - except: - pass - return None - -#builds and returns (index, index_filename) for specified maf_file -def build_maf_index( maf_file, species = None ): - indexes = bx.interval_index_file.Indexes() - try: - maf_reader = bx.align.maf.Reader( open( maf_file ) ) - # Need to be a bit tricky in our iteration here to get the 'tells' right - while True: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - if species is not None and c.src.split( "." )[0] not in species: - continue - indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) - fd, index_filename = tempfile.mkstemp() - out = os.fdopen( fd, 'w' ) - indexes.write( out ) - out.close() - return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = True ), index_filename ) - except: - return ( None, None ) - -def __main__(): - #Parse Command Line - primary_species = sys.argv.pop( 1 ) - secondary_species = sys.argv.pop( 1 ).split( "," ) - include_primary = True - if secondary_species == ["None"]: - secondary_species = None - else: - try: secondary_species.remove( primary_species ) - except: include_primary = False - maf_identifier = sys.argv.pop( 1 ) - interval_file = sys.argv.pop( 1 ) - output_file = sys.argv.pop( 1 ) - maf_source_type = sys.argv.pop( 1 ) - - #ensure primary_species is set - if primary_species == "?": - print >> sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." - sys.exit() - - #get index for mafs based on type - index = index_filename = None - #using specified uid for locally cached - if maf_source_type.lower() in ["cached"]: - index = maf_index_by_uid( maf_identifier ) - if index is None: - print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( maf_uid ) - sys.exit() - elif maf_source_type.lower() in ["user"]: - #index maf for use here, need to remove index_file when finished - index, index_filename = build_maf_index( maf_identifier, species = [primary_species] ) - if index is None: - print >> sys.stderr, "Your MAF file appears to be malformed." - sys.exit() - else: - print >> sys.stderr, "Invalid MAF source type specified." - sys.exit() - - #open output file - output = open( output_file, "w" ) - - #Step through gene bed - genes_extracted = 0 - line_count = 0 - for line_count, line in enumerate( open( interval_file, "r" ).readlines() ): - try: - if line[0:1]=="#": - continue - - #load gene bed fields - try: - #Starts and ends for exons - starts = [] - ends = [] - - fields = line.split() - #Requires atleast 12 BED columns - if len(fields) < 12: - continue - chrom = fields[0] - tx_start = int( fields[1] ) - tx_end = int( fields[2] ) - name = fields[3] - strand = fields[5] - if strand != '-': strand='+' #Default strand is + - cds_start = int( fields[6] ) - cds_end = int( fields[7] ) - - #Calculate and store starts and ends of coding exons - region_start, region_end = cds_start, cds_end - exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) ) - exon_starts = map( ( lambda x: x + tx_start ), exon_starts ) - exon_ends = map( int, fields[10].rstrip( ',' ).split( ',' ) ) - exon_ends = map( ( lambda x, y: x + y ), exon_starts, exon_ends ); - for start, end in zip( exon_starts, exon_ends ): - start = max( start, region_start ) - end = min( end, region_end ) - if start < end: - starts.append( start ) - ends.append( end ) - #create spliced alignment object - if secondary_species is not None: alignment = SplicedAlignment( starts, ends, [primary_species] + secondary_species ) - else: alignment = SplicedAlignment( starts, ends, [primary_species] ) - except Exception, e: - print "Error loading exon positions from input line %i: %s" % ( line_count, e ) - continue - - for exon in alignment.exons: - try: - primary_src = "%s.%s" % ( primary_species, chrom ) - start = exon.start - end = exon.end - - #Get blocks overlaping this position - blocks = index.get( primary_src, start, end ) - #Order the blocks by score, lowest first - blocks_order = [] - for i, block in enumerate( blocks ): - for j in range( 0, len( blocks_order ) ): - if float( block.score ) < float( blocks[blocks_order[j]].score ): - blocks_order.insert( j, i ) - break - else: - blocks_order.append( i ) - - #Loop through ordered block indexes and layer them - for block_index in blocks_order: - #Get maf block - maf = blocks[block_index] - #Limit maf block to desired species - if secondary_species is not None: - maf = maf.limit_to_species( [primary_species] + secondary_species ) - #Colapse extraneous gap columns - maf.remove_all_gap_columns() - #Positions and strand are in reference to ref - ref = maf.get_component_by_src( primary_src ) - #We want our block coordinates to be from positive strand, if region is on negative strand, we will reverse compliment it at the end - if ref.strand == "-": - maf = maf.reverse_complement() - ref = maf.get_component_by_src( primary_src ) - - #slice maf by start and end - slice_start = max( start, ref.start ) - slice_end = min( end, ref.end ) - try: - maf = maf.slice_by_component( ref, slice_start, slice_end ) - ref = maf.get_component_by_src( primary_src ) - except: - continue - - #skip gap locations due to insertions in secondary species relative to primary species - start_offset = slice_start - start - num_gaps = 0 - for i in range( len( ref.text.rstrip().rstrip("-") ) ): - if ref.text[i] in ["-"]: - num_gaps += 1 - continue - #Set base for all species - for spec in [ c.src.split( '.' )[0] for c in maf.components ]: - try: - #NB: If a gap appears in higher scoring secondary species block, - #it will overwrite any bases that have been set by lower scoring blocks - #this seems more proper than allowing, e.g. a single base from lower scoring alignment to exist outside of its genomic context - exon.set_position( start_offset + i - num_gaps, spec, maf.get_component_by_src_start( spec ).text[i] ) - except: - #species/sequence for species does not exist - pass - except Exception, e: - print "Error filling exons with MAFs from input line %i: %s" % ( line_count, e ) - continue - - #Write alignment to output file - #Output primary species first, if requested - if include_primary: - output.write( ">%s.%s\n" %( primary_species, name ) ) - if strand == "-": - output.write( alignment.get_sequence_reverse_complement( primary_species ) ) - else: - output.write( alignment.get_sequence( primary_species ) ) - output.write( "\n" ) - #Output all remainging species - for spec in secondary_species or alignment.get_species_names( skip = primary_species ): - output.write( ">%s.%s\n" % ( spec, name ) ) - if strand == "-": - output.write( alignment.get_sequence_reverse_complement( spec ) ) - else: - output.write( alignment.get_sequence( spec ) ) - output.write( "\n" ) - - output.write( "\n" ) - - genes_extracted += 1 - - except Exception, e: - print "Unexpected error from input line %i: %s" % ( line_count, e ) - continue - - #close output file - output.close() - - #remove index file if created during run - if index_filename is not None: - os.unlink( index_filename ) - - #Print message about success for user - if genes_extracted > 0: - print "%i genes were extracted successfully." % ( genes_extracted ) - else: - print "No genes were extracted." - if line_count > 0: - print "This tool requires your input file to conform to the 12 column BED standard." - -if __name__ == "__main__": __main__() diff --git a/tools/extract/interval2maf.py b/tools/extract/interval2maf.py deleted file mode 100755 index 7f91ba224e9..00000000000 --- a/tools/extract/interval2maf.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/env python2.4 - -""" -Reads a list of intervals and a maf. Produces a new maf containing the -blocks or parts of blocks in the original that overlapped the intervals. - -If a MAF file, not UID, is provided the MAF file is indexed before being processed. - -NOTE: If two intervals overlap the same block it will be written twice. - -usage: %prog maf_file [options] - -d, --dbkey=d: Database key, ie hg17 - -c, --chromCol=c: Column of Chr - -s, --startCol=s: Column of Start - -e, --endCol=e: Column of End - -S, --strandCol=S: Column of Strand - -t, --mafType=t: Type of MAF source to use - -m, --mafFile=m: Path of source MAF file, if not using cached version - -i, --interval_file=i: Input interval file - -o, --output_file=o: Output MAF file - -p, --species=p: Species to include in output -""" - -#Dan Blankenberg -import pkg_resources; pkg_resources.require( "bx-python" ) -from bx.cookbook import doc_optparse -import bx.align.maf -import bx.intervals.io -import bx.interval_index_file -import sys, os, tempfile - -MAF_LOCATION_FILE = "/depot/data2/galaxy/maf_index.loc" - -def maf_index_by_uid( maf_uid ): - for line in open( MAF_LOCATION_FILE ): - try: - #read each line, if not enough fields, go to next line - if line[0:1] == "#" : continue - fields = line.split('\t') - if maf_uid == fields[1]: - try: - maf_files = fields[3].replace( "\n", "" ).replace( "\r", "" ).split( "," ) - return bx.align.maf.MultiIndexed( maf_files, keep_open = True, parse_e_rows = True ) - except Exception, e: - raise 'MAF UID (%s) found, but configuration appears to be malformed: %s' % ( maf_uid, e ) - except: - pass - return None - -#builds and returns (index, index_filename) for specified maf_file -def build_maf_index( maf_file, species = None ): - indexes = bx.interval_index_file.Indexes() - try: - maf_reader = bx.align.maf.Reader( open( maf_file ) ) - # Need to be a bit tricky in our iteration here to get the 'tells' right - while True: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - if species is not None and c.src.split( "." )[0] not in species: - continue - indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) - fd, index_filename = tempfile.mkstemp() - out = os.fdopen( fd, 'w' ) - indexes.write( out ) - out.close() - return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = True ), index_filename ) - except: - return ( None, None ) - -def __main__(): - index = index_filename = None - mincols = 0 - - # Parse Command Line - options, args = doc_optparse.parse( __doc__ ) - - if options.dbkey: dbkey = options.dbkey - else: dbkey = None - if dbkey in [None, "?"]: - print >>sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." - sys.exit() - - species = None - if options.species: - species = options.species.split( ',' ) - if "None" in species: species = None - - if options.chromCol: chromCol = int( options.chromCol ) - 1 - else: - print >>sys.stderr, "Chromosome column has not been specified." - sys.exit() - - if options.startCol: startCol = int( options.startCol ) - 1 - else: - print >>sys.stderr, "Start column has not been specified." - sys.exit() - - if options.endCol: endCol = int( options.endCol ) - 1 - else: - print >>sys.stderr, "End column has not been specified." - sys.exit() - - if options.strandCol: strandCol = int( options.strandCol ) - 1 - else: - print >>sys.stderr, "Strand column has not been specified." - sys.exit() - - if options.interval_file: interval_file = options.interval_file - else: - print >>sys.stderr, "Input interval file has not been specified." - sys.exit() - - if options.output_file: output_file = options.output_file - else: - print >>sys.stderr, "Output file has not been specified." - sys.exit() - - #Open indexed access to MAFs - if options.mafType: - index = maf_index_by_uid( options.mafType ) - if index is None: - print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( options.mafType ) - sys.exit() - elif options.mafFile: - index, index_filename = build_maf_index( options.mafFile, species = [dbkey] ) - if index is None: - print >> sys.stderr, "Your MAF file appears to be malformed." - sys.exit() - else: - print >>sys.stderr, "Desired source MAF type has not been specified." - sys.exit() - - out = bx.align.maf.Writer( open(output_file, "w") ) - - # Iterate over input regions - num_blocks = 0 - num_lines = 0 - for num_lines, region in enumerate( bx.intervals.io.NiceReaderWrapper( open( interval_file, 'r' ), chrom_col = chromCol, start_col = startCol, end_col = endCol, strand_col = strandCol, fix_strand = True, return_header = False, return_comments = False ) ): - try: - src = "%s.%s" % ( dbkey, region.chrom ) - - blocks = index.get( src, region.start, region.end ) - - for block in blocks: - ref = block.get_component_by_src( src ) - #We want our block coordinates to be from positive strand - if ref.strand == "-": - block = block.reverse_complement() - ref = block.get_component_by_src( src ) - - #save old score here for later use - old_score = block.score - slice_start = max( region.start, ref.start ) - slice_end = min( region.end, ref.end ) - - #when interval is out-of-range (not in maf index), fail silently: else could create tons of scroll - try: - block = block.slice_by_component( ref, slice_start, slice_end ) - except: - continue - - if block.text_size > mincols: - if region.strand != ref.strand: block = block.reverse_complement() - # restore old score, may not be accurate, but it is better than 0 for everything - block.score = old_score - if species is not None: - block = block.limit_to_species( species ) - block.remove_all_gap_columns() - out.write( block ) - num_blocks += 1 - except Exception, e: - print "Error found on input line %s: %s." % ( num_lines, e ) - continue - - # Close output MAF - out.close() - - #remove index file if created during run - if index_filename is not None: - os.unlink( index_filename ) - - print "%s MAF blocks extracted." % num_blocks - -if __name__ == "__main__": __main__() diff --git a/tools/extract/interval2maf_pairwise.py b/tools/extract/interval2maf_pairwise.py deleted file mode 100644 index 8c407d6efce..00000000000 --- a/tools/extract/interval2maf_pairwise.py +++ /dev/null @@ -1,200 +0,0 @@ -#!/usr/bin/env python2.4 - -""" -Reads a list of intervals and a maf. Produces a new maf containing the -blocks or parts of blocks in the original that overlapped the intervals. - -If index_file is not provided maf_file.index is used. - -NOTE: If two intervals overlap the same block it will be written twice. With - non-overlapping intervals and --chop this is never a problem. - Chop has been made always on. - -usage: %prog maf_file index_file [options] < interval_file - -d, --dbkey=d: Database key, ie hg17 - -c, --chromCol=c: Column of Chr - -s, --startCol=s: Column of Start - -e, --endCol=e: Column of End - -S, --strandCol=S: Column of Strand - -t, --mafType=t: Type of MAF source to use - -i, --interval_file=i: Input interval file - -o, --output_file=o: Output MAF file -""" - -#import psyco_full - -import pkg_resources; pkg_resources.require( "bx-python" ) -from bx.cookbook import doc_optparse - -import bx.align.maf -from bx import interval_index_file -import bx.intervals.io -from bx import misc -import os -import sys - -def __main__(): - - # Parse Command Line - - options, args = doc_optparse.parse( __doc__ ) - - #dictionary of available maf files - maf_sets = {} - try: - for line in open( "/depot/data2/galaxy/maf_pairwise.loc" ): - if line[0:1] == "#" : continue - fields = line.split('\t') - #read each line, if not enough fields, go to next line - try: - maf_desc = fields[0] - maf_uid = fields[1] - builds = fields[2] - build_to_common_list = {} - common_to_build_list = {} - split_builds = builds.split(",") - for build in split_builds: - this_build = build.split("=")[0] - try: - this_common = build.split("=")[1] - except: - this_common = this_build - build_to_common_list[this_build]=this_common - common_to_build_list[this_common]=this_build - - paths = fields[3].replace("\n","").replace("\r","") - maf_sets[maf_uid]={} - maf_sets[maf_uid]['description']=maf_desc - maf_sets[maf_uid]['builds']=build_to_common_list - maf_sets[maf_uid]['common']=common_to_build_list - maf_sets[maf_uid]['paths']=paths.split(",") - except: - continue - - except Exception, exc: - print >>sys.stdout, 'interval2maf_pairwise.py initialization error -> %s' % exc - - - try: - mincols=0 - - if options.dbkey: dbkey = options.dbkey - else: dbkey="?" - - - if options.chromCol: chromCol= int(options.chromCol) - 1 - else: - print >>sys.stderr, "Chromosome column has not been specified." - sys.exit() - - if options.startCol: startCol= int(options.startCol) - 1 - else: - print >>sys.stderr, "Start column has not been specified." - sys.exit() - - if options.endCol: endCol= int(options.endCol) - 1 - else: - print >>sys.stderr, "End column has not been specified." - sys.exit() - - if options.strandCol: strandCol= int(options.strandCol) - 1 - else: - print >>sys.stderr, "Strand column has not been specified." - sys.exit() - - if options.mafType: mafType= options.mafType - else: - print >>sys.stderr, "Desired source MAF type has not been specified." - sys.exit() - - if options.interval_file: interval_file= options.interval_file - else: - print >>sys.stderr, "Input interval file has not been specified." - sys.exit() - - if options.output_file: output_file= options.output_file - else: - print >>sys.stderr, "Output file has not been specified." - sys.exit() - except: - sys.exit() - - if dbkey == "?": - print >>sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." - sys.exit() - - - #Open MAF Files, with indexes - try: - maf_files = maf_sets[mafType]['paths'] - except: - print >>sys.stderr, "The MAF source specified appears to be invalid." - sys.exit() - - try: - # Open indexed access to mafs - index = bx.align.maf.MultiIndexed( maf_files, keep_open=True, parse_e_rows=True ) - except: - print >>sys.stderr, "The MAF source specified [", mafType ,"] appears to be missing." - sys.exit() - - #convert dbkey to name in maf file, if no db->name entry, use build - try: - dbkey = maf_sets[mafType]['builds'][dbkey] - except: - print >>sys.stderr, "This MAF set is not available for this build." - sys.exit() - - out = bx.align.maf.Writer( open(output_file, "w") ) - - # Iterate over input ranges - num_blocks=0 - num_lines = 0 - for region in bx.intervals.io.NiceReaderWrapper( open(interval_file, 'r' ), chrom_col=chromCol, start_col=startCol, end_col=endCol, strand_col=strandCol, fix_strand=True, return_header=False, return_comments=False): - try: - num_lines += 1 - src = "%s.%s" % (dbkey,region.chrom) - start = region.start - end = region.end - strand = region.strand - - # Find overlap with reference component - blocks = index.get( src, start, end ) - for block in blocks: - ref = block.get_component_by_src( src ) - #We want our block coordinates to be from positive strand - if ref.strand == "-": - block = block.reverse_complement() - ref = block.get_component_by_src( src ) - #save old score here for later use - old_score = block.score - #slice maf by start and end - slice_start = max( start, ref.start ) - slice_end = min( end, ref.end ) - - #when interval is out-of-range (not in maf index), fail silently: else could create tons of scroll - try: - sliced = block.slice_by_component( ref, slice_start, slice_end ) - except: - continue - - if sliced.text_size > mincols: - if strand != ref.strand: sliced = sliced.reverse_complement() - # restore old score, may not be accurate, but it is better than 0 for everything - sliced.score = old_score - for c in sliced.components: - spec,chrom = bx.align.src_split( c.src ) - if not spec or not chrom: - spec = chrom = c.src - if spec in maf_sets[mafType]['common']: - c.src = bx.align.src_merge(maf_sets[mafType]['common'][spec],chrom) - out.write( sliced ) - num_blocks+=1 - except Exception, e: - print "Error found on input line:",num_lines - print e - continue - - # Close output MAF - out.close() - print num_blocks, "MAF blocks extracted." -if __name__ == "__main__": __main__() diff --git a/tools/extract/interval_maf_to_merged_fasta.py b/tools/extract/interval_maf_to_merged_fasta.py deleted file mode 100644 index ba41dda0d15..00000000000 --- a/tools/extract/interval_maf_to_merged_fasta.py +++ /dev/null @@ -1,271 +0,0 @@ -#!/usr/bin/env python2.4 - -""" -Reads an interval file and an indexed MAF. Produces a FASTA file containing -the aligned sequences, based upon the provided coordinates - -Alignment blocks are layered ontop of each other based upon score. - -usage: %prog dbkey_of_interval_file comma_separated_list_of_additional_dbkeys_to_extract maf_file|cached_maf_uid input_interval_file output_fasta_file chromCol startCol endCol strandCol user|cached -""" - -#Dan Blankenberg -import pkg_resources; pkg_resources.require( "bx-python" ) -import bx.align.maf -import bx.intervals.io -import bx.interval_index_file -import sys, os, tempfile, string - -MAF_LOCATION_FILE = "/depot/data2/galaxy/maf_index.loc" - -#an object corresponding to a reference layered alignment -class RegionAlignment( object ): - - DNA_COMPLEMENT = string.maketrans( "ACGTacgt", "TGCAtgca" ) - - def __init__( self, size, species = [] ): - self.size = size - self.sequences = {} - if not isinstance( species, list ): - species = [species] - for spec in species: - self.add_species( spec ) - - #add a species to the alignment - def add_species( self, species ): - #make temporary sequence files, need to delete when done - fd, file_path = tempfile.mkstemp() - self.sequences[species] = { 'file':os.fdopen( fd, 'w+' ), 'path':file_path } - self.sequences[species]['file'].write( "-" * self.size ) - - #returns the names for species found in alignment, skipping names as requested - def get_species_names( self, skip = [] ): - if not isinstance( skip, list ): skip = [skip] - names = self.sequences.keys() - for name in skip: - try: names.remove( name ) - except: pass - return names - - #returns the sequence for a species - def get_sequence( self, species ): - self.sequences[species]['file'].seek( 0 ) - return self.sequences[species]['file'].read() - - #returns the reverse complement of the sequence for a species - def get_sequence_reverse_complement( self, species ): - complement = [base for base in self.get_sequence( species ).translate( self.DNA_COMPLEMENT )] - complement.reverse() - return "".join( complement ) - - #sets a position for a species - def set_position( self, index, species, base ): - if index >= self.size or index < 0: raise "Your index (%i) is out of range (0 - %i)." % ( index, self.size - 1 ) - if len(base) != 1: raise "A genomic position can only have a length of 1." - if species not in self.sequences.keys(): self.add_species( species ) - self.sequences[species]['file'].seek( index ) - self.sequences[species]['file'].write( base ) - - #Flush temp file of specified species, or all species - def flush( self, species = None ): - if species is None: - species = self.sequences.keys() - elif not isinstance( species, list ): - species = [species] - for spec in species: - self.sequences[spec]['file'].flush() - - #object cleanup, delete temporary files - def __del__( self ): - for species, sequence in self.sequences.items(): - sequence['file'].close() - os.unlink( sequence['path'] ) - -def maf_index_by_uid( maf_uid ): - for line in open( MAF_LOCATION_FILE ): - try: - #read each line, if not enough fields, go to next line - if line[0:1] == "#" : continue - fields = line.split('\t') - if maf_uid == fields[1]: - try: - maf_files = fields[3].replace( "\n", "" ).replace( "\r", "" ).split( "," ) - return bx.align.maf.MultiIndexed( maf_files, keep_open = True, parse_e_rows = True ) - except Exception, e: - raise 'MAF UID (%s) found, but configuration appears to be malformed: %s' % ( maf_uid, e ) - except: - pass - return None - -#builds and returns (index, index_filename) for specified maf_file -def build_maf_index( maf_file, species = None ): - indexes = bx.interval_index_file.Indexes() - try: - maf_reader = bx.align.maf.Reader( open( maf_file ) ) - # Need to be a bit tricky in our iteration here to get the 'tells' right - while True: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - if species is not None and c.src.split( "." )[0] not in species: - continue - indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) - fd, index_filename = tempfile.mkstemp() - out = os.fdopen( fd, 'w' ) - indexes.write( out ) - out.close() - return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = True ), index_filename ) - except: - return ( None, None ) - -def __main__(): - #Parse Command Line - primary_species = sys.argv.pop( 1 ) - secondary_species = sys.argv.pop( 1 ).split( "," ) - include_primary = True - if secondary_species == ["None"]: - secondary_species = None - else: - try: secondary_species.remove( primary_species ) - except: include_primary = False - maf_identifier = sys.argv.pop( 1 ) - interval_file = sys.argv.pop( 1 ) - output_file = sys.argv.pop( 1 ) - try: - chr_col = int( sys.argv.pop( 1 ).strip() ) - 1 - start_col = int( sys.argv.pop( 1 ).strip() ) - 1 - end_col = int( sys.argv.pop( 1 ).strip() ) - 1 - strand_col = int( sys.argv.pop( 1 ).strip() ) - 1 - maf_source_type = sys.argv.pop( 1 ) - except: - print >> sys.stderr, "You appear to be missing metadata. You can specify your metadata by clicking on the pencil icon associated with your interval file." - sys.exit() - - #ensure primary_species is set - if primary_species == "?": - print >> sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." - sys.exit() - - #get index for mafs based on type - index = index_filename = None - #using specified uid for locally cached - if maf_source_type.lower() in ["cached"]: - index = maf_index_by_uid( maf_identifier ) - if index is None: - print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( maf_identifier ) - sys.exit() - elif maf_source_type.lower() in ["user"]: - #index maf for use here, need to remove index_file when finished - index, index_filename = build_maf_index( maf_identifier, species = [primary_species] ) - if index is None: - print >> sys.stderr, "Your MAF file appears to be malformed." - sys.exit() - else: - print >> sys.stderr, "Invalid MAF source type specified." - sys.exit() - - #open output file - output = open( output_file, "w" ) - - #Step through interval file - intervals_extracted = 0 - line_count = 0 - for line_count, region in enumerate( bx.intervals.io.NiceReaderWrapper( open(interval_file, 'r' ), chrom_col = chr_col, start_col = start_col, end_col = end_col, strand_col = strand_col, fix_strand = True, return_header = False, return_comments = False ) ): - #create alignment object - if secondary_species is not None: alignment = RegionAlignment( region.end - region.start, [primary_species] + secondary_species ) - else: alignment = RegionAlignment( region.end - region.start, primary_species ) - primary_src = "%s.%s" % ( primary_species, region.chrom ) - - #Get blocks overlaping this position - blocks = index.get( primary_src, region.start, region.end ) - #Order the blocks by score, lowest first - blocks_order = [] - for i, block in enumerate( blocks ): - for j in range( 0, len( blocks_order ) ): - if float( block.score ) < float( blocks[blocks_order[j]].score ): - blocks_order.insert( j, i ) - break - else: - blocks_order.append( i ) - - #Loop through ordered block indexes and layer them - for block_index in blocks_order: - #Get maf block - maf = blocks[block_index] - #Limit maf block to desired species - if secondary_species is not None: - maf = maf.limit_to_species( [primary_species] + secondary_species ) - #Colapse extraneous gap columns - maf.remove_all_gap_columns() - #Positions and strand are in reference to ref - ref = maf.get_component_by_src( primary_src ) - #We want our block coordinates to be from positive strand, if region is on negative strand, we will reverse compliment it at the end - if ref.strand == "-": - maf = maf.reverse_complement() - ref = maf.get_component_by_src( primary_src ) - - #slice maf by start and end - slice_start = max( region.start, ref.start ) - slice_end = min( region.end, ref.end ) - try: - maf = maf.slice_by_component( ref, slice_start, slice_end ) - ref = maf.get_component_by_src( primary_src ) - except: - continue - - #skip gap locations due to insertions in secondary species relative to primary species - start_offset = slice_start - region.start - num_gaps = 0 - for i in range( len( ref.text.rstrip().rstrip( "-" ) ) ): - if ref.text[i] in ["-"]: - num_gaps += 1 - continue - #Set base for all species - for spec in [ c.src.split( '.' )[0] for c in maf.components ]: - try: - #NB: If a gap appears in higher scoring secondary species block, - #it will overwrite any bases that have been set by lower scoring blocks - #this seems more proper than allowing, e.g. a single base from lower scoring alignment to exist outside of its genomic context - alignment.set_position( start_offset + i - num_gaps, spec, maf.get_component_by_src_start( spec ).text[i] ) - except: - #species/sequence for species does not exist - pass - - #Write alignment to output file - #Output primary species first, if requested - if include_primary: - output.write( ">%s.%s(%s):%s-%s\n" %( primary_species, region.chrom, region.strand, region.start, region.end ) ) - if region.strand == "-": - output.write( alignment.get_sequence_reverse_complement( primary_species ) ) - else: - output.write( alignment.get_sequence( primary_species ) ) - output.write( "\n" ) - #Output all remainging species - for spec in secondary_species or alignment.get_species_names( skip = primary_species ): - output.write( ">%s\n" % ( spec ) ) - if region.strand == "-": - output.write( alignment.get_sequence_reverse_complement( spec ) ) - else: - output.write( alignment.get_sequence( spec ) ) - output.write( "\n" ) - - output.write( "\n" ) - intervals_extracted += 1 - - output.close() - - #remove index file if created during run - if index_filename is not None: - os.unlink( index_filename ) - - #Print message about success for user - if intervals_extracted > 0: - print "%i regions were extracted successfully." % ( intervals_extracted ) - else: - print "No regions were extracted." - if line_count > 0: - print "Make sure your metadata is properly set by clicking the pencil icon associated with your interval file." - - -if __name__ == "__main__": __main__() diff --git a/tools/filters/maf/maf_to_fasta_concat.py b/tools/filters/maf/maf_to_fasta_concat.py deleted file mode 100755 index 4d2aeb1cb25..00000000000 --- a/tools/filters/maf/maf_to_fasta_concat.py +++ /dev/null @@ -1,65 +0,0 @@ -#!/usr/bin/env python2.3 - -""" -Read a maf and print the text as a fasta file, concatenating blocks - -usage %prog species1,species2 maf_file out_file -""" -#Dan Blankenberg -from __future__ import division - -import textwrap -import sys -import pkg_resources; pkg_resources.require( "bx-python" ) -from bx.align import maf - -def __main__(): - print "Restricted to species:", sys.argv[1] - - texts = {} - - input_filename = sys.argv[2] - output_filename = sys.argv[3] - species = sys.argv[1].split( ',' ) - - if "None" in species: - species = get_species( input_filename ) - - file_out = open( output_filename, 'w' ) - for spec in species: - file_out.write( ">" + spec + "\n" ) - try: - for m in maf.Reader( open( input_filename, 'r' ) ): - c = m.get_component_by_src_start( spec ) - if c: file_out.write( c.text ) - else: file_out.write( "-" * m.text_size ) - except: - print >>sys.stderr, "Your MAF file appears to be malformed." - sys.exit() - file_out.write( "\n" ) - file_out.close() - -def get_species( maf_filename ): - try: - species={} - - file_in = open( maf_filename, 'r' ) - maf_reader = maf.Reader( file_in ) - - for i, m in enumerate( maf_reader ): - l = m.components - for c in l: - spec, chrom = maf.src_split( c.src ) - if not spec or not chrom: - spec = chrom = c.src - species[spec] = spec - - file_in.close() - - species = species.keys() - species.sort() - return species - except: - return [] - -if __name__ == "__main__": __main__() diff --git a/tools/extract/genebed_maf_to_fasta.xml b/tools/maf/genebed_maf_to_fasta.xml similarity index 86% rename from tools/extract/genebed_maf_to_fasta.xml rename to tools/maf/genebed_maf_to_fasta.xml index 80afe71c7b7..5d951a522b8 100644 --- a/tools/extract/genebed_maf_to_fasta.xml +++ b/tools/maf/genebed_maf_to_fasta.xml @@ -1,7 +1,7 @@ given a set of coding exon intervals - #if $maf_source_type.maf_source == "user":#genebed_maf_to_fasta.py $dbkey $maf_source_type.species $maf_source_type.maf_file $input1 $out_file1 $maf_source_type.maf_source -#else:#genebed_maf_to_fasta.py $dbkey $maf_source_type.species $maf_source_type.maf_identifier $input1 $out_file1 $maf_source_type.maf_source + #if $maf_source_type.maf_source == "user":#interval_maf_to_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_file --interval_file=$input1 --output_file=$out_file1 --mafSourceType=$maf_source_type.maf_source --geneBED +#else:#interval_maf_to_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_identifier --interval_file=$input1 --output_file=$out_file1 --mafSourceType=$maf_source_type.maf_source --geneBED #end if diff --git a/tools/maf/interval2maf.py b/tools/maf/interval2maf.py new file mode 100755 index 00000000000..a82b24b315f --- /dev/null +++ b/tools/maf/interval2maf.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python2.4 + +""" +Reads a list of intervals and a maf. Produces a new maf containing the +blocks or parts of blocks in the original that overlapped the intervals. + +If a MAF file, not UID, is provided the MAF file is indexed before being processed. + +NOTE: If two intervals overlap the same block it will be written twice. + +usage: %prog maf_file [options] + -d, --dbkey=d: Database key, ie hg17 + -c, --chromCol=c: Column of Chr + -s, --startCol=s: Column of Start + -e, --endCol=e: Column of End + -S, --strandCol=S: Column of Strand + -t, --mafType=t: Type of MAF source to use + -m, --mafFile=m: Path of source MAF file, if not using cached version + -i, --interval_file=i: Input interval file + -o, --output_file=o: Output MAF file + -p, --species=p: Species to include in output + -l, --indexLocation=l: Override default maf_index.loc file +""" + +#Dan Blankenberg +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.cookbook import doc_optparse +import bx.align.maf +import bx.intervals.io +import maf_utilities +import sys + +def __main__(): + index = index_filename = None + mincols = 0 + + #Parse Command Line + options, args = doc_optparse.parse( __doc__ ) + + if options.dbkey: dbkey = options.dbkey + else: dbkey = None + if dbkey in [None, "?"]: + print >>sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." + sys.exit() + + species = None + if options.species: + species = options.species.split( ',' ) + if "None" in species: species = None + + if options.chromCol: chromCol = int( options.chromCol ) - 1 + else: + print >>sys.stderr, "Chromosome column has not been specified." + sys.exit() + + if options.startCol: startCol = int( options.startCol ) - 1 + else: + print >>sys.stderr, "Start column has not been specified." + sys.exit() + + if options.endCol: endCol = int( options.endCol ) - 1 + else: + print >>sys.stderr, "End column has not been specified." + sys.exit() + + if options.strandCol: strandCol = int( options.strandCol ) - 1 + else: + print >>sys.stderr, "Strand column has not been specified." + sys.exit() + + if options.interval_file: interval_file = options.interval_file + else: + print >>sys.stderr, "Input interval file has not been specified." + sys.exit() + + if options.output_file: output_file = options.output_file + else: + print >>sys.stderr, "Output file has not been specified." + sys.exit() + #Finish parsing command line + + #Open indexed access to MAFs + if options.mafType: + if options.indexLocation: + index = maf_utilities.maf_index_by_uid( options.mafType, options.indexLocation ) + else: + index = maf_utilities.maf_index_by_uid( options.mafType ) + if index is None: + print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( options.mafType ) + sys.exit() + elif options.mafFile: + index, index_filename = maf_utilities.build_maf_index( options.mafFile, species = [dbkey] ) + if index is None: + print >> sys.stderr, "Your MAF file appears to be malformed." + sys.exit() + else: + print >>sys.stderr, "Desired source MAF type has not been specified." + sys.exit() + + #Create MAF writter + out = bx.align.maf.Writer( open(output_file, "w") ) + + #Iterate over input regions + num_blocks = 0 + num_regions = None + for num_regions, region in enumerate( bx.intervals.io.NiceReaderWrapper( open( interval_file, 'r' ), chrom_col = chromCol, start_col = startCol, end_col = endCol, strand_col = strandCol, fix_strand = True, return_header = False, return_comments = False ) ): + src = "%s.%s" % ( dbkey, region.chrom ) + for block in maf_utilities.get_chopped_blocks_for_region( index, src, region, species, mincols ): + out.write( block ) + num_blocks += 1 + + #Close output MAF + out.close() + + #remove index file if created during run + maf_utilities.remove_temp_index_file( index_filename ) + + if num_blocks: + print "%i MAF blocks extracted for %i regions." % ( num_blocks, ( num_regions + 1 ) ) + elif num_regions is not None: + print "No MAF blocks could be extracted for %i regions." % ( num_regions + 1 ) + else: + print "No valid regions have been provided." + +if __name__ == "__main__": __main__() diff --git a/tools/extract/interval2maf.xml b/tools/maf/interval2maf.xml similarity index 100% rename from tools/extract/interval2maf.xml rename to tools/maf/interval2maf.xml diff --git a/tools/extract/interval2maf_pairwise.xml b/tools/maf/interval2maf_pairwise.xml similarity index 81% rename from tools/extract/interval2maf_pairwise.xml rename to tools/maf/interval2maf_pairwise.xml index 0f0fb52dd9e..9ce61776274 100644 --- a/tools/extract/interval2maf_pairwise.xml +++ b/tools/maf/interval2maf_pairwise.xml @@ -1,6 +1,6 @@ given a set of genomic intervals - interval2maf_pairwise.py --dbkey=$dbkey --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafType=$mafType --interval_file=$input1 --output_file=$out_file1 + interval2maf.py --dbkey=$input1_dbkey --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafType=$mafType --interval_file=$input1 --output_file=$out_file1 --indexLocation=/depot/data2/galaxy/maf_pairwise.loc diff --git a/tools/maf/interval_maf_to_merged_fasta.py b/tools/maf/interval_maf_to_merged_fasta.py new file mode 100644 index 00000000000..41f370f64e5 --- /dev/null +++ b/tools/maf/interval_maf_to_merged_fasta.py @@ -0,0 +1,183 @@ +#!/usr/bin/env python2.4 + +""" +Reads an interval or gene BED and a MAF Source. +Produces a FASTA file containing the aligned intervals/gene sequences, based upon the provided coordinates + +Alignment blocks are layered ontop of each other based upon score. + +usage: %prog maf_file [options] + -d, --dbkey=d: Database key, ie hg17 + -c, --chromCol=c: Column of Chr + -s, --startCol=s: Column of Start + -e, --endCol=e: Column of End + -S, --strandCol=S: Column of Strand + -G, --geneBED: Input is a Gene BED file, process and join exons as one region + -t, --mafSourceType=t: Type of MAF source to use + -m, --mafSource=m: Path of source MAF file, if not using cached version + -i, --interval_file=i: Input interval file + -o, --output_file=o: Output MAF file + -p, --species=p: Species to include in output + +usage: %prog dbkey_of_BED comma_separated_list_of_additional_dbkeys_to_extract comma_separated_list_of_indexed_maf_files input_gene_bed_file output_fasta_file cached|user +""" + +#Dan Blankenberg +import maf_utilities +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.cookbook import doc_optparse +import bx.intervals.io +import sys + +def __main__(): + + #Parse Command Line + options, args = doc_optparse.parse( __doc__ ) + mincols = 0 + + if options.dbkey: primary_species = options.dbkey + else: primary_species = None + if primary_species in [None, "?", "None"]: + print >>sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." + sys.exit() + + include_primary = True + if options.species: + secondary_species = options.species.split( ',' ) + if "None" in secondary_species: + secondary_species = None + species = None + else: + try: secondary_species.remove( primary_species ) + except: include_primary = False + species = [primary_species] + secondary_species + + if options.interval_file: interval_file = options.interval_file + else: + print >>sys.stderr, "Input interval file has not been specified." + sys.exit() + + if options.output_file: output_file = options.output_file + else: + print >>sys.stderr, "Output file has not been specified." + sys.exit() + + if not options.geneBED: + if options.chromCol: chr_col = int( options.chromCol ) - 1 + else: + print >>sys.stderr, "Chromosome column has not been specified." + sys.exit() + + if options.startCol: start_col = int( options.startCol ) - 1 + else: + print >>sys.stderr, "Start column has not been specified." + sys.exit() + + if options.endCol: end_col = int( options.endCol ) - 1 + else: + print >>sys.stderr, "End column has not been specified." + sys.exit() + + if options.strandCol: strand_col = int( options.strandCol ) - 1 + else: + print >>sys.stderr, "Strand column has not been specified." + sys.exit() + #Finish parsing command line + + #get index for mafs based on type + index = index_filename = None + #using specified uid for locally cached + if options.mafSourceType.lower() in ["cached"]: + index = maf_utilities.maf_index_by_uid( options.mafSource ) + if index is None: + print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( maf_identifier ) + sys.exit() + elif options.mafSourceType.lower() in ["user"]: + #index maf for use here, need to remove index_file when finished + index, index_filename = maf_utilities.build_maf_index( options.mafSource, species = [primary_species] ) + if index is None: + print >> sys.stderr, "Your MAF file appears to be malformed." + sys.exit() + else: + print >> sys.stderr, "Invalid MAF source type specified." + sys.exit() + + #open output file + output = open( output_file, "w" ) + + + if options.geneBED: + region_enumerator = maf_utilities.line_enumerator( open( interval_file, "r" ).readlines() ) + else: + region_enumerator = enumerate( bx.intervals.io.NiceReaderWrapper( open( interval_file, 'r' ), chrom_col = chr_col, start_col = start_col, end_col = end_col, strand_col = strand_col, fix_strand = True, return_header = False, return_comments = False ) ) + + #Step through intervals + regions_extracted = 0 + line_count = 0 + for line_count, line in region_enumerator: + try: + if options.geneBED: #Process as Gene BED + try: + starts, ends, fields = maf_utilities.get_starts_ends_fields_from_gene_bed( line ) + #create spliced alignment object + alignment = maf_utilities.get_spliced_region_alignment( index, primary_species, fields[0], starts, ends, strand = '+', species = species, mincols = mincols ) + primary_name = secondary_name = fields[3] + alignment_strand = fields[5] + except Exception, e: + print "Error loading exon positions from input line %i: %s" % ( line_count, e ) + continue + else: #Process as standard intervals + try: + #create spliced alignment object + alignment = maf_utilities.get_region_alignment( index, primary_species, line.chrom, line.start, line.end, strand = '+', species = species, mincols = mincols ) + primary_name = "%s(%s):%s-%s" % ( line.chrom, line.strand, line.start, line.end ) + secondary_name = "" + alignment_strand = line.strand + except Exception, e: + print "Error loading region positions from input line %i: %s" % ( line_count, e ) + continue + + #Write alignment to output file + #Output primary species first, if requested + if include_primary: + output.write( ">%s.%s\n" %( primary_species, primary_name ) ) + if alignment_strand == "-": + output.write( alignment.get_sequence_reverse_complement( primary_species ) ) + else: + output.write( alignment.get_sequence( primary_species ) ) + output.write( "\n" ) + #Output all remainging species + for spec in secondary_species or alignment.get_species_names( skip = primary_species ): + if secondary_name: + output.write( ">%s.%s\n" % ( spec, secondary_name ) ) + else: + output.write( ">%s\n" % ( spec ) ) + if alignment_strand == "-": + output.write( alignment.get_sequence_reverse_complement( spec ) ) + else: + output.write( alignment.get_sequence( spec ) ) + output.write( "\n" ) + + output.write( "\n" ) + + regions_extracted += 1 + + except Exception, e: + print "Unexpected error from input line %i: %s" % ( line_count, e ) + continue + + #close output file + output.close() + + #remove index file if created during run + maf_utilities.remove_temp_index_file( index_filename ) + + #Print message about success for user + if regions_extracted > 0: + print "%i regions were processed successfully." % ( regions_extracted ) + else: + print "No regions were processed successfully." + if line_count > 0 and options.geneBED: + print "This tool requires your input file to conform to the 12 column BED standard." + +if __name__ == "__main__": __main__() diff --git a/tools/extract/interval_maf_to_merged_fasta.xml b/tools/maf/interval_maf_to_merged_fasta.xml similarity index 82% rename from tools/extract/interval_maf_to_merged_fasta.xml rename to tools/maf/interval_maf_to_merged_fasta.xml index 0cbe851089e..d035bf98820 100644 --- a/tools/extract/interval_maf_to_merged_fasta.xml +++ b/tools/maf/interval_maf_to_merged_fasta.xml @@ -1,7 +1,7 @@ given a set of genomic intervals - #if $maf_source_type.maf_source == "user":#interval_maf_to_merged_fasta.py $dbkey $maf_source_type.species $maf_source_type.maf_file $input1 $out_file1 $input1_chromCol $input1_startCol $input1_endCol $input1_strandCol $maf_source_type.maf_source -#else:#interval_maf_to_merged_fasta.py $dbkey $maf_source_type.species $maf_source_type.maf_identifier $input1 $out_file1 $input1_chromCol $input1_startCol $input1_endCol $input1_strandCol $maf_source_type.maf_source + #if $maf_source_type.maf_source == "user":#interval_maf_to_merged_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_file --interval_file=$input1 --output_file=$out_file1 --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafSourceType=$maf_source_type.maf_source +#else:#interval_maf_to_merged_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_identifier --interval_file=$input1 --output_file=$out_file1 --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafSourceType=$maf_source_type.maf_source #end if diff --git a/tools/filters/maf/maf_by_block_number.py b/tools/maf/maf_by_block_number.py similarity index 90% rename from tools/filters/maf/maf_by_block_number.py rename to tools/maf/maf_by_block_number.py index abf266b9136..a5a386b8208 100644 --- a/tools/filters/maf/maf_by_block_number.py +++ b/tools/maf/maf_by_block_number.py @@ -29,9 +29,9 @@ def __main__(): failed_lines.append( str( ctr ) ) continue try: - for count, m in enumerate( bx.align.maf.Reader( open( input_maf_filename, 'r' ) ) ): + for count, block in enumerate( bx.align.maf.Reader( open( input_maf_filename, 'r' ) ) ): if count == block_wanted: - maf_writer.write( m ) + maf_writer.write( block ) break except: print >>sys.stderr, "Your MAF file appears to be malformed." diff --git a/tools/filters/maf/maf_by_block_number.xml b/tools/maf/maf_by_block_number.xml similarity index 100% rename from tools/filters/maf/maf_by_block_number.xml rename to tools/maf/maf_by_block_number.xml diff --git a/tools/filters/maf/maf_filter.py b/tools/maf/maf_filter.py similarity index 100% rename from tools/filters/maf/maf_filter.py rename to tools/maf/maf_filter.py diff --git a/tools/filters/maf/maf_filter.xml b/tools/maf/maf_filter.xml similarity index 100% rename from tools/filters/maf/maf_filter.xml rename to tools/maf/maf_filter.xml diff --git a/tools/filters/maf/maf_limit_size.py b/tools/maf/maf_limit_size.py similarity index 100% rename from tools/filters/maf/maf_limit_size.py rename to tools/maf/maf_limit_size.py diff --git a/tools/filters/maf/maf_limit_size.xml b/tools/maf/maf_limit_size.xml similarity index 100% rename from tools/filters/maf/maf_limit_size.xml rename to tools/maf/maf_limit_size.xml diff --git a/tools/filters/maf/maf_limit_to_species.py b/tools/maf/maf_limit_to_species.py similarity index 100% rename from tools/filters/maf/maf_limit_to_species.py rename to tools/maf/maf_limit_to_species.py diff --git a/tools/filters/maf/maf_limit_to_species.xml b/tools/maf/maf_limit_to_species.xml similarity index 100% rename from tools/filters/maf/maf_limit_to_species.xml rename to tools/maf/maf_limit_to_species.xml diff --git a/tools/filters/maf/maf_reverse_complement.py b/tools/maf/maf_reverse_complement.py similarity index 100% rename from tools/filters/maf/maf_reverse_complement.py rename to tools/maf/maf_reverse_complement.py diff --git a/tools/filters/maf/maf_reverse_complement.xml b/tools/maf/maf_reverse_complement.xml similarity index 100% rename from tools/filters/maf/maf_reverse_complement.xml rename to tools/maf/maf_reverse_complement.xml diff --git a/tools/filters/maf/maf_stats.py b/tools/maf/maf_stats.py similarity index 57% rename from tools/filters/maf/maf_stats.py rename to tools/maf/maf_stats.py index 0fe5a9b2a70..31640ce2e59 100644 --- a/tools/filters/maf/maf_stats.py +++ b/tools/maf/maf_stats.py @@ -4,54 +4,11 @@ Reads a list of intervals and a maf. Outputs a new set of intervals with statistics appended. """ -import sys, tempfile, os +import sys import pkg_resources; pkg_resources.require( "bx-python" ) -import bx.align.maf import bx.intervals.io -import bx.interval_index_file -import psyco_full from numpy import zeros - -MAF_LOCATION_FILE = "/depot/data2/galaxy/maf_index.loc" - -def maf_index_by_uid( maf_uid ): - for line in open( MAF_LOCATION_FILE ): - try: - #read each line, if not enough fields, go to next line - if line[0:1] == "#" : continue - fields = line.split('\t') - if maf_uid == fields[1]: - try: - maf_files = fields[3].replace( "\n", "" ).replace( "\r", "" ).split( "," ) - return bx.align.maf.MultiIndexed( maf_files, keep_open = True, parse_e_rows = False ) - except Exception, e: - raise 'MAF UID (%s) found, but configuration appears to be malformed: %s' % ( maf_uid, e ) - except: - pass - return None - -#builds and returns (index, index_filename) for specified maf_file -def build_maf_index( maf_file, species = None ): - indexes = bx.interval_index_file.Indexes() - try: - maf_reader = bx.align.maf.Reader( open( maf_file ) ) - # Need to be a bit tricky in our iteration here to get the 'tells' right - while True: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - if species is not None and c.src.split( "." )[0] not in species: - continue - indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) - fd, index_filename = tempfile.mkstemp() - out = os.fdopen( fd, 'w' ) - indexes.write( out ) - out.close() - return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = False ), index_filename ) - except: - return ( None, None ) - +import maf_utilities def __main__(): maf_source_type = sys.argv.pop( 1 ) @@ -73,13 +30,13 @@ def __main__(): index = index_filename = None if maf_source_type == "user": #index maf for use here - index, index_filename = build_maf_index( input_maf_filename, species = [dbkey] ) + index, index_filename = maf_utilities.build_maf_index( input_maf_filename, species = [dbkey] ) if index is None: print >>sys.stderr, "Your MAF file appears to be malformed." sys.exit() elif maf_source_type == "cached": #access existing indexes - index = maf_index_by_uid( input_maf_filename ) + index = maf_utilities.maf_index_by_uid( input_maf_filename ) if index is None: print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( input_maf_filename ) sys.exit() @@ -98,35 +55,21 @@ def __main__(): total_length += ( region.end - region.start ) coverage = { dbkey: zeros( region.end - region.start, dtype = bool ) } - blocks = index.get( src, region.start, region.end ) - for maf in blocks: + for block in maf_utilities.get_chopped_blocks_for_region( index, src, region, force_strand='+' ): #make sure all species are known - for c in maf.components: + for c in block.components: spec = c.src.split( '.' )[0] if spec not in coverage: coverage[spec] = zeros( region.end - region.start, dtype = bool ) - #slice maf by start and end - ref = maf.get_component_by_src( src ) - # If the reference component is on the '-' strand we should complement the interval - if ref.strand == '-': - maf = maf.reverse_complement() - ref = maf.get_component_by_src( src ) - slice_start = max( region.start, ref.start ) - slice_end = min( region.end, ref.end ) - try: - maf = maf.slice_by_component( ref, slice_start, slice_end ) - except: - continue - ref = maf.get_component_by_src( ref.src ) - + ref = block.get_component_by_src( src ) #skip gap locations due to insertions in secondary species relative to primary species - start_offset = slice_start - region.start + start_offset = ref.start - region.start num_gaps = 0 for i in range( len( ref.text.rstrip().rstrip( "-" ) ) ): if ref.text[i] in ["-"]: num_gaps += 1 continue #Toggle base if covered - for comp in maf.components: + for comp in block.components: spec = comp.src.split( '.' )[0] if comp.text and comp.text[i] not in ['-']: coverage[spec][start_offset + i - num_gaps] = True @@ -149,6 +92,6 @@ def __main__(): out.write( "%s\t%s\t%.4f\n" % ( spec, species_summary[spec], float( species_summary[spec] ) / total_length ) ) out.close() print "%i regions were processed with a total length of %i." % ( num_region, total_length ) - if index_filename is not None: - os.unlink( index_filename ) + maf_utilities.remove_temp_index_file( index_filename ) + if __name__ == "__main__": __main__() diff --git a/tools/filters/maf/maf_stats.xml b/tools/maf/maf_stats.xml similarity index 100% rename from tools/filters/maf/maf_stats.xml rename to tools/maf/maf_stats.xml diff --git a/tools/filters/maf/maf_stats_code.py b/tools/maf/maf_stats_code.py similarity index 100% rename from tools/filters/maf/maf_stats_code.py rename to tools/maf/maf_stats_code.py diff --git a/tools/filters/maf/maf_thread_for_species.py b/tools/maf/maf_thread_for_species.py similarity index 94% rename from tools/filters/maf/maf_thread_for_species.py rename to tools/maf/maf_thread_for_species.py index 44348468f4f..d29bdbe2f37 100644 --- a/tools/filters/maf/maf_thread_for_species.py +++ b/tools/maf/maf_thread_for_species.py @@ -41,8 +41,8 @@ def main(): m.components = new_components m.score = 0.0 maf_writer.write( m ) - except: - print >> sys.stderr, "Error steping through MAF File" + except Exception, e: + print >> sys.stderr, "Error steping through MAF File: %s" % e sys.exit() maf_reader.close() maf_writer.close() diff --git a/tools/filters/maf/maf_thread_for_species.xml b/tools/maf/maf_thread_for_species.xml similarity index 100% rename from tools/filters/maf/maf_thread_for_species.xml rename to tools/maf/maf_thread_for_species.xml diff --git a/tools/filters/maf/maf_to_bed.py b/tools/maf/maf_to_bed.py similarity index 97% rename from tools/filters/maf/maf_to_bed.py rename to tools/maf/maf_to_bed.py index d1a4e0398d3..fc921ccf99a 100644 --- a/tools/filters/maf/maf_to_bed.py +++ b/tools/maf/maf_to_bed.py @@ -3,11 +3,7 @@ """ Read a maf and output intervals for specified list of species. """ - -from __future__ import division - -import textwrap -import sys, tempfile, os +import sys, os import pkg_resources; pkg_resources.require( "bx-python" ) from bx.align import maf diff --git a/tools/filters/maf/maf_to_bed.xml b/tools/maf/maf_to_bed.xml similarity index 100% rename from tools/filters/maf/maf_to_bed.xml rename to tools/maf/maf_to_bed.xml diff --git a/tools/filters/maf/maf_to_bed_code.py b/tools/maf/maf_to_bed_code.py similarity index 100% rename from tools/filters/maf/maf_to_bed_code.py rename to tools/maf/maf_to_bed_code.py diff --git a/tools/filters/maf/maf_to_fasta.xml b/tools/maf/maf_to_fasta.xml similarity index 100% rename from tools/filters/maf/maf_to_fasta.xml rename to tools/maf/maf_to_fasta.xml diff --git a/tools/maf/maf_to_fasta_concat.py b/tools/maf/maf_to_fasta_concat.py new file mode 100755 index 00000000000..2437afccd2f --- /dev/null +++ b/tools/maf/maf_to_fasta_concat.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python2.3 + +""" +Read a maf and print the text as a fasta file, concatenating blocks + +usage %prog species1,species2 maf_file out_file +""" +#Dan Blankenberg +import sys +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.align import maf +import maf_utilities + +def __main__(): + print "Restricted to species:", sys.argv[1] + + texts = {} + + input_filename = sys.argv[2] + output_filename = sys.argv[3] + species = sys.argv[1].split( ',' ) + + if "None" in species: + species = maf_utilities.get_species_in_maf( input_filename ) + + file_out = open( output_filename, 'w' ) + for spec in species: + file_out.write( ">" + spec + "\n" ) + try: + for block in maf.Reader( open( input_filename, 'r' ) ): + component = block.get_component_by_src_start( spec ) + if component: file_out.write( component.text ) + else: file_out.write( "-" * m.text_size ) + except: + print >>sys.stderr, "Your MAF file appears to be malformed." + sys.exit() + file_out.write( "\n" ) + file_out.close() + + +if __name__ == "__main__": __main__() diff --git a/tools/filters/maf/maf_to_fasta_multiple_sets.py b/tools/maf/maf_to_fasta_multiple_sets.py similarity index 57% rename from tools/filters/maf/maf_to_fasta_multiple_sets.py rename to tools/maf/maf_to_fasta_multiple_sets.py index 2c496d2181e..1eb5abf0a59 100755 --- a/tools/filters/maf/maf_to_fasta_multiple_sets.py +++ b/tools/maf/maf_to_fasta_multiple_sets.py @@ -4,16 +4,13 @@ Read a maf and print the text as a fasta file. """ #Dan Blankenberg -from __future__ import division - -import textwrap import sys import pkg_resources; pkg_resources.require( "bx-python" ) from bx.align import maf def __main__(): print "Restricted to species:", sys.argv[3] - + input_filename = sys.argv[1] output_filename = sys.argv[2] species = sys.argv[3].split( ',' ) @@ -26,20 +23,16 @@ def __main__(): file_out = open( output_filename, 'w' ) - block_num = -1 - - for i, m in enumerate( maf_reader ): - block_num += 1 + for block_num, block in enumerate( maf_reader ): if "None" not in species: - m = m.limit_to_species( species ) - l = m.components - if len(l) < num_species and partial == "partial_disallowed": continue - for c in l: - spec, chrom = maf.src_split( c.src ) + block = block.limit_to_species( species ) + if len( block.components ) < num_species and partial == "partial_disallowed": continue + for component in block.components: + spec, chrom = maf.src_split( component.src ) if not spec or not chrom: - spec = chrom = c.src - file_out.write( ">" + c.src + "(" + c.strand + "):" + str( c.start ) + "-" + str( c.end ) + "|" + spec + "_" + str( block_num ) + "\n" ) - file_out.write( c.text + "\n" ) + spec = chrom = component.src + file_out.write( ">" + component.src + "(" + component.strand + "):" + str( component.start ) + "-" + str( component.end ) + "|" + spec + "_" + str( block_num ) + "\n" ) + file_out.write( component.text + "\n" ) file_out.write( "\n" ) file_in.close() except: diff --git a/tools/maf/maf_utilities.py b/tools/maf/maf_utilities.py new file mode 100644 index 00000000000..cca583ae553 --- /dev/null +++ b/tools/maf/maf_utilities.py @@ -0,0 +1,321 @@ +#!/usr/bin/env python2.4 +""" +Provides wrappers and utilities for working with MAF files and alignments. +""" +#Dan Blankenberg +import pkg_resources; pkg_resources.require( "bx-python" ) +import bx.align.maf +import bx.intervals +import bx.interval_index_file +import sys, os, string, tempfile + +MAF_LOCATION_FILE = "/depot/data2/galaxy/maf_index.loc" + +#an object corresponding to a reference layered alignment +class RegionAlignment( object ): + + DNA_COMPLEMENT = string.maketrans( "ACGTacgt", "TGCAtgca" ) + MAX_SEQUENCE_SIZE = sys.maxint #Maximum length of sequence allowed + + def __init__( self, size, species = [] ): + assert size <= self.MAX_SEQUENCE_SIZE, "Maximum length allowed for an individual sequence has been exceeded (%i > %i)." % ( size, self.MAX_SEQUENCE_SIZE ) + self.size = size + self.sequences = {} + if not isinstance( species, list ): + species = [species] + for spec in species: + self.add_species( spec ) + + #add a species to the alignment + def add_species( self, species ): + #make temporary sequence files + self.sequences[species] = tempfile.TemporaryFile() + self.sequences[species].write( "-" * self.size ) + + #returns the names for species found in alignment, skipping names as requested + def get_species_names( self, skip = [] ): + if not isinstance( skip, list ): skip = [skip] + names = self.sequences.keys() + for name in skip: + try: names.remove( name ) + except: pass + return names + + #returns the sequence for a species + def get_sequence( self, species ): + self.sequences[species].seek( 0 ) + return self.sequences[species].read() + + #returns the reverse complement of the sequence for a species + def get_sequence_reverse_complement( self, species ): + complement = [base for base in self.get_sequence( species ).translate( self.DNA_COMPLEMENT )] + complement.reverse() + return "".join( complement ) + + #sets a position for a species + def set_position( self, index, species, base ): + if index >= self.size or index < 0: raise "Your index (%i) is out of range (0 - %i)." % ( index, self.size - 1 ) + if len(base) != 1: raise "A genomic position can only have a length of 1." + if species not in self.sequences.keys(): self.add_species( species ) + self.sequences[species].seek( index ) + self.sequences[species].write( base ) + + #Flush temp file of specified species, or all species + def flush( self, species = None ): + if species is None: + species = self.sequences.keys() + elif not isinstance( species, list ): + species = [species] + for spec in species: + self.sequences[spec].flush() + +class GenomicRegionAlignment( RegionAlignment ): + + def __init__( self, start, end, species = [] ): + RegionAlignment.__init__( self, end - start, species ) + self.start = start + self.end = end + +class SplicedAlignment( object ): + + DNA_COMPLEMENT = string.maketrans( "ACGTacgt", "TGCAtgca" ) + + def __init__( self, exon_starts, exon_ends, species = [] ): + if not isinstance( exon_starts, list ): + exon_starts = [exon_starts] + if not isinstance( exon_ends, list ): + exon_ends = [exon_ends] + assert len( exon_starts ) == len( exon_ends ), "The number of starts does not match the number of sizes." + self.exons = [] + for i in range( len( exon_starts ) ): + self.exons.append( GenomicRegionAlignment( exon_starts[i], exon_ends[i], species ) ) + + #returns the names for species found in alignment, skipping names as requested + def get_species_names( self, skip = [] ): + if not isinstance( skip, list ): skip = [skip] + names = [] + for exon in self.exons: + for name in exon.get_species_names( skip = skip ): + if name not in names: + names.append( name ) + return names + + #returns the sequence for a species + def get_sequence( self, species ): + sequence = tempfile.TemporaryFile() + for exon in self.exons: + if species in exon.get_species_names(): + sequence.write( exon.get_sequence( species ) ) + else: + sequence.write( "-" * exon.size ) + sequence.seek( 0 ) + return sequence.read() + + #returns the reverse complement of the sequence for a species + def get_sequence_reverse_complement( self, species ): + complement = [base for base in self.get_sequence( species ).translate( self.DNA_COMPLEMENT )] + complement.reverse() + return "".join( complement ) + + #Start and end of coding region + @property + def start( self ): + return self.exons[0].start + @property + def end( self ): + return self.exons[-1].end + +#Open a MAF index using a UID +def maf_index_by_uid( maf_uid, index_location_file = None ): + for line in open( index_location_file or MAF_LOCATION_FILE ): + try: + #read each line, if not enough fields, go to next line + if line[0:1] == "#" : continue + fields = line.split('\t') + if maf_uid == fields[1]: + try: + maf_files = fields[3].replace( "\n", "" ).replace( "\r", "" ).split( "," ) + return bx.align.maf.MultiIndexed( maf_files, keep_open = True, parse_e_rows = True ) + except Exception, e: + raise 'MAF UID (%s) found, but configuration appears to be malformed: %s' % ( maf_uid, e ) + except: + pass + return None + +#builds and returns (index, index_filename) for specified maf_file +def build_maf_index( maf_file, species = None ): + indexes = bx.interval_index_file.Indexes() + try: + maf_reader = bx.align.maf.Reader( open( maf_file ) ) + # Need to be a bit tricky in our iteration here to get the 'tells' right + while True: + pos = maf_reader.file.tell() + block = maf_reader.next() + if block is None: break + for c in block.components: + if species is not None and c.src.split( "." )[0] not in species: + continue + indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) + fd, index_filename = tempfile.mkstemp() + out = os.fdopen( fd, 'w' ) + indexes.write( out ) + out.close() + return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = True ), index_filename ) + except: + return ( None, None ) + +#generator yielding only chopped and valid blocks for a specified region +def get_chopped_blocks_for_region( index, src, region, species = None, mincols = 0, force_strand = None ): + for block in index.get_as_iterator( src, region.start, region.end ): + ref = block.get_component_by_src( src ) + #We want our block coordinates to be from positive strand + if ref.strand == "-": + block = block.reverse_complement() + ref = block.get_component_by_src( src ) + + #save old score here for later use + old_score = block.score + slice_start = max( region.start, ref.start ) + slice_end = min( region.end, ref.end ) + + #slice block by reference species at determined limits + block = block.slice_by_component( ref, slice_start, slice_end ) + + if block.text_size > mincols: + if ( force_strand is None and region.strand != ref.strand ) or ( force_strand is not None and force_strand != ref.strand ): + block = block.reverse_complement() + # restore old score, may not be accurate, but it is better than 0 for everything + block.score = old_score + if species is not None: + block = block.limit_to_species( species ) + block.remove_all_gap_columns() + yield block + +#returns a filled region alignment for specified regions +def get_region_alignment( index, primary_species, chrom, start, end, strand = '+', species = None, mincols = 0 ): + if species is not None: alignment = RegionAlignment( end - start, species ) + else: alignment = RegionAlignment( end - start, primary_species ) + return fill_region_alignment( alignment, index, primary_species, chrom, start, end, strand, species, mincols ) + +#fills a region alignment +def fill_region_alignment( alignment, index, primary_species, chrom, start, end, strand = '+', species = None, mincols = 0 ): + #first step through blocks, save index and score in array, then order by score (array will start as 0=index0,scoreX) + #step through ordered list, step through maf blocks, stopping at index, store, then break inner loop + region = bx.intervals.Interval( start, end ) + region.chrom = chrom + region.strand = strand + primary_src = "%s.%s" % ( primary_species, chrom ) + + #Order blocks overlaping this position by score, lowest first + blocks_order = [] + for i, block in enumerate( get_chopped_blocks_for_region( index, primary_src, region, species, mincols ) ): + for j in range( 0, len( blocks_order ) ): + if float( block.score ) < float( blocks_order[j]['score'] ): + blocks_order.insert( j, {'index':i, 'score':block.score} ) + break + else: + blocks_order.append( {'index':i, 'score':block.score} ) + + #Loop through ordered block indexes and layer blocks by score + for block_dict in blocks_order: + for block_index, block in enumerate( get_chopped_blocks_for_region( index, primary_src, region, species, mincols ) ): + if block_index == block_dict['index']: + ref = block.get_component_by_src( primary_src ) + #skip gap locations due to insertions in secondary species relative to primary species + start_offset = ref.start - start + num_gaps = 0 + for i in range( len( ref.text.rstrip().rstrip("-") ) ): + if ref.text[i] in ["-"]: + num_gaps += 1 + continue + #Set base for all species + for spec in [ c.src.split( '.' )[0] for c in block.components ]: + try: + #NB: If a gap appears in higher scoring secondary species block, + #it will overwrite any bases that have been set by lower scoring blocks + #this seems more proper than allowing, e.g. a single base from lower scoring alignment to exist outside of its genomic context + alignment.set_position( start_offset + i - num_gaps, spec, block.get_component_by_src_start( spec ).text[i] ) + except: + #species/sequence for species does not exist + pass + break + return alignment + +#returns a filled spliced region alignment for specified region with start and end lists +def get_spliced_region_alignment( index, primary_species, chrom, starts, ends, strand = '+', species = None, mincols = 0 ): + #create spliced alignment object + if species is not None: alignment = SplicedAlignment( starts, ends, species ) + else: alignment = SplicedAlignment( starts, ends, [primary_species] ) + for exon in alignment.exons: + fill_region_alignment( exon, index, primary_species, chrom, exon.start, exon.end, strand, species, mincols) + return alignment + +#loop through string array, only return non-commented lines +def line_enumerator( lines, comment_start = '#' ): + i = 0 + for line in lines: + if not line.startswith( comment_start ): + i += 1 + yield ( i, line ) + +#read a GeneBed file, return list of starts, ends, raw fields +def get_starts_ends_fields_from_gene_bed( line ): + #Starts and ends for exons + starts = [] + ends = [] + + fields = line.split() + #Requires atleast 12 BED columns + if len(fields) < 12: + raise Exception, "Not a proper 12 column BED line (%s)." % line + chrom = fields[0] + tx_start = int( fields[1] ) + tx_end = int( fields[2] ) + name = fields[3] + strand = fields[5] + if strand != '-': strand='+' #Default strand is + + cds_start = int( fields[6] ) + cds_end = int( fields[7] ) + + #Calculate and store starts and ends of coding exons + region_start, region_end = cds_start, cds_end + exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) ) + exon_starts = map( ( lambda x: x + tx_start ), exon_starts ) + exon_ends = map( int, fields[10].rstrip( ',' ).split( ',' ) ) + exon_ends = map( ( lambda x, y: x + y ), exon_starts, exon_ends ); + for start, end in zip( exon_starts, exon_ends ): + start = max( start, region_start ) + end = min( end, region_end ) + if start < end: + starts.append( start ) + ends.append( end ) + return ( starts, ends, fields ) + +def get_species_in_maf( maf_filename ): + try: + species={} + + file_in = open( maf_filename, 'r' ) + maf_reader = maf.Reader( file_in ) + + for i, m in enumerate( maf_reader ): + l = m.components + for c in l: + spec, chrom = maf.src_split( c.src ) + if not spec or not chrom: + spec = chrom = c.src + species[spec] = spec + + file_in.close() + + species = species.keys() + species.sort() + return species + except: + return [] + + + +def remove_temp_index_file( index_filename ): + try: os.unlink( index_filename ) + except: pass