From 6eaa6d511cc44eefaec2b24e728b8db4e269c9d6 Mon Sep 17 00:00:00 2001 From: Daniel Blankenberg Date: Mon, 12 Nov 2007 21:06:24 +0000 Subject: [PATCH] Update interval2maf tools to use one script file. --- tools/extract/interval2maf.py | 211 +++++++++++++---------------- tools/extract/interval2maf.xml | 4 +- tools/extract/interval2maf_code.py | 47 ------- tools/extract/user_interval2maf.py | 160 ---------------------- 4 files changed, 98 insertions(+), 324 deletions(-) delete mode 100644 tools/extract/interval2maf_code.py delete mode 100644 tools/extract/user_interval2maf.py diff --git a/tools/extract/interval2maf.py b/tools/extract/interval2maf.py index 8f9b4afc2b9..436ebfc13c0 100755 --- a/tools/extract/interval2maf.py +++ b/tools/extract/interval2maf.py @@ -4,83 +4,87 @@ 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. +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. With - non-overlapping intervals and --chop this is never a problem. - Chop has been made always on. +NOTE: If two intervals overlap the same block it will be written twice. -usage: %prog maf_file index_file [options] < interval_file +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 """ -#import psyco_full - +#Dan Blankenberg 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 +import bx.interval_index_file +import sys, os -def __main__(): +MAF_LOCATION_FILE = "/depot/data2/galaxy/maf_index.loc" - # 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_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') - #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 + 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 - except Exception, exc: - print >>sys.stdout, 'interval2maf.py initialization error -> %s' % exc +#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 + options, args = doc_optparse.parse( __doc__ ) + index = index_filename = None + mincols = 0 try: - mincols=0 - if options.dbkey: dbkey = options.dbkey - else: 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() + if options.species: species = options.species.split( ',' ) + else: species = None if options.chromCol: chromCol= int(options.chromCol) - 1 else: @@ -101,12 +105,7 @@ def __main__(): 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." @@ -117,53 +116,34 @@ def __main__(): 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() + #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() + except Exception, exc: + print >>sys.stdout, 'interval2maf.py initialization error -> %s' % exc out = bx.align.maf.Writer( open(output_file, "w") ) - # Iterate over input ranges - num_blocks=0 + # Iterate over input regions + 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): + 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: - num_lines += 1 - src = "%s.%s" % (dbkey,region.chrom) - start = region.start - end = region.end - strand = region.strand - - #do a fix on src, chromosome for mlagan alignments (they lack chr) - if mafType == "ENCODE_MLAGAN": - src = src.replace(".chr",".") + src = "%s.%s" % ( dbkey, region.chrom ) + + blocks = index.get( src, region.start, region.end ) - 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 @@ -173,34 +153,35 @@ def __main__(): #save old score here for later use old_score = block.score - slice_start = max( start, ref.start ) - slice_end = min( end, ref.end ) + 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: - sliced = block.slice_by_component( ref, slice_start, slice_end ) + block = block.slice_by_component( ref, slice_start, slice_end ) except: continue if sliced.text_size > mincols: - if strand != ref.strand: sliced = sliced.reverse_complement() + if region.strand != ref.strand: block = block.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 - c.src = bx.align.src_merge(spec,chrom) - 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 + 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:",num_lines - print e + print "Error found on input line %s: %s." % ( num_lines, e ) continue # Close output MAF out.close() - print num_blocks, "MAF blocks extracted." + + #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.xml b/tools/extract/interval2maf.xml index 462307879ee..236b022e936 100644 --- a/tools/extract/interval2maf.xml +++ b/tools/extract/interval2maf.xml @@ -1,6 +1,6 @@ given a set of genomic intervals - #if $maf_source_type.maf_source == "user":#user_interval2maf.py --dbkey=$input1_dbkey --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafFile=$maf_source_type.mafFile --interval_file=$input1 --output_file=$out_file1 + #if $maf_source_type.maf_source == "user":#interval2maf.py --dbkey=$input1_dbkey --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafFile=$maf_source_type.mafFile --interval_file=$input1 --output_file=$out_file1 #else:#interval2maf.py --dbkey=$input1_dbkey --chromCol=$input1_chromCol --startCol=$input1_startCol --endCol=$input1_endCol --strandCol=$input1_strandCol --mafType=$maf_source_type.mafType --interval_file=$input1 --output_file=$out_file1 #end if @@ -15,7 +15,7 @@ - + diff --git a/tools/extract/interval2maf_code.py b/tools/extract/interval2maf_code.py deleted file mode 100644 index 5a5b3dbad24..00000000000 --- a/tools/extract/interval2maf_code.py +++ /dev/null @@ -1,47 +0,0 @@ -#build list of available data -import os, sys -maf_sets = {} - -try: - for line in open( "/depot/data2/galaxy/maf_index.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_list =[] - split_builds = builds.split(",") - for build in split_builds: - this_build = build.split("=")[0] - build_list.append(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_list - except: - continue - -except Exception, exc: - print >>sys.stdout, 'interval2maf_code.py initialization error -> %s' % exc - -#return available datasets for group and build, set None option as selected for hg16 -def get_available_data( build ): - available_sets = [] - for key in maf_sets: - if build in maf_sets[key]['builds']: - available_sets.append((maf_sets[key]['description'],key,False)) - if len(available_sets) < 1: - available_sets.append(('No data available for this build','None',True)) - return available_sets - - -def exec_before_job(app, inp_data, out_data, param_dict, tool): - if param_dict['maf_source_type']['maf_source'] == "cached": - for name, data in out_data.items(): - try: - data.name = data.name + " [" + maf_sets[str(param_dict['maf_source_type']['mafType'])]['description'] + "]" - except KeyError: - data.name = data.name + " [unknown MAF source specified]" diff --git a/tools/extract/user_interval2maf.py b/tools/extract/user_interval2maf.py deleted file mode 100644 index 97c2daa96e2..00000000000 --- a/tools/extract/user_interval2maf.py +++ /dev/null @@ -1,160 +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 - -m, --mafFile=m: MAF file 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 -import tempfile -import os, sys - -def __main__(): - - # Parse Command Line - - options, args = doc_optparse.parse( __doc__ ) - - 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.mafFile: mafFile= options.mafFile - 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() - - #index maf for use here - indexes = interval_index_file.Indexes() - - try: - maf_reader = bx.align.maf.Reader( open( mafFile ) ) - # Need to be a bit tricky in our iteration here to get the 'tells' right - while 1: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) - index_filename = tempfile.NamedTemporaryFile().name - out = open(index_filename,'w') - indexes.write(out) - out.close() - except: - print >>sys.stderr, "Your MAF file appears to be malformed." - sys.exit() - - index = bx.align.maf.Indexed(mafFile, index_filename = index_filename, keep_open=True, parse_e_rows=True) - 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 - 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_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 - c.src = bx.align.src_merge(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() - #Delete index file - os.unlink(index_filename) - print num_blocks, "MAF blocks extracted." -if __name__ == "__main__": __main__() \ No newline at end of file