mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Fix memory issues when dealing with MAF tools.
Add new bx-python eggs. Reorganize all MAF tools into a unified directory, and update tool_conf. Remember to update tool_conf.xml from one of the samples (tool_conf.xml.sample or tool_conf.xml.main).
This commit is contained in:
Daniel Blankenberg
+14
-11
@@ -21,7 +21,7 @@
|
||||
<tool file="data_source/encode_import_gencode.xml" />
|
||||
</section>
|
||||
<section name="ENCODE Tools" id="EncodeTools">
|
||||
<!-- <tool file="extract/interval2maf.xml" />
|
||||
<!-- <tool file="maf/interval2maf.xml" />
|
||||
<tool file="extract/phastOdds/phastOdds_tool.xml" />
|
||||
<tool file="stats/aggregate_binned_scores_in_intervals.xml" /> -->
|
||||
<tool file="encode/gencode_partition.xml" />
|
||||
@@ -52,8 +52,8 @@
|
||||
<tool file="stats/grouping.xml" />
|
||||
</section>
|
||||
<section name="Convert Formats" id="convert">
|
||||
<tool file="filters/maf/maf_to_fasta.xml" />
|
||||
<tool file="filters/maf/maf_to_bed.xml" />
|
||||
<tool file="maf/maf_to_fasta.xml" />
|
||||
<tool file="maf/maf_to_bed.xml" />
|
||||
<tool file="filters/gff2bed.xml" />
|
||||
<tool file="filters/bed2gff.xml" />
|
||||
<!-- <tool file="filters/axt_to_fasta.xml" />
|
||||
@@ -73,14 +73,17 @@
|
||||
<tool file="extract/twoBitToFa_wrapper.xml" />
|
||||
</section>
|
||||
<section name="Fetch Alignments" id="fetchAlign">
|
||||
<tool file="extract/interval2maf_pairwise.xml" />
|
||||
<tool file="extract/interval2maf.xml" />
|
||||
<tool file="extract/interval_maf_to_merged_fasta.xml" />
|
||||
<!-- <tool file="extract/genebed_maf_to_fasta.xml"/>
|
||||
<tool file="filters/maf/maf_stats.xml"/> -->
|
||||
<tool file="filters/maf/maf_limit_to_species.xml"/>
|
||||
<tool file="filters/maf/maf_limit_size.xml"/>
|
||||
<tool file="filters/maf/maf_by_block_number.xml"/>
|
||||
<tool file="maf/interval2maf_pairwise.xml" />
|
||||
<tool file="maf/interval2maf.xml" />
|
||||
<tool file="maf/interval_maf_to_merged_fasta.xml" />
|
||||
<tool file="maf/genebed_maf_to_fasta.xml"/>
|
||||
<tool file="maf/maf_stats.xml"/>
|
||||
<tool file="maf/maf_thread_for_species.xml"/>
|
||||
<tool file="maf/maf_limit_to_species.xml"/>
|
||||
<tool file="maf/maf_limit_size.xml"/>
|
||||
<tool file="maf/maf_by_block_number.xml"/>
|
||||
<!-- <tool file="maf/maf_reverse_complement.xml"/>
|
||||
<tool file="maf/maf_filter.xml"/> -->
|
||||
</section>
|
||||
<section name="Get Genomic Scores" id="scores">
|
||||
<tool file="stats/wiggle_to_simple.xml" />
|
||||
|
||||
+14
-14
@@ -21,7 +21,7 @@
|
||||
<tool file="data_source/encode_import_gencode.xml" />
|
||||
</section>
|
||||
<section name="ENCODE Tools" id="EncodeTools">
|
||||
<!-- <tool file="extract/interval2maf.xml" />
|
||||
<!-- <tool file="maf/interval2maf.xml" />
|
||||
<tool file="extract/phastOdds/phastOdds_tool.xml" />
|
||||
<tool file="stats/aggregate_binned_scores_in_intervals.xml" /> -->
|
||||
<tool file="encode/gencode_partition.xml" />
|
||||
@@ -52,8 +52,8 @@
|
||||
<tool file="stats/grouping.xml" />
|
||||
</section>
|
||||
<section name="Convert Formats" id="convert">
|
||||
<tool file="filters/maf/maf_to_fasta.xml" />
|
||||
<tool file="filters/maf/maf_to_bed.xml" />
|
||||
<tool file="maf/maf_to_fasta.xml" />
|
||||
<tool file="maf/maf_to_bed.xml" />
|
||||
<tool file="filters/gff2bed.xml" />
|
||||
<tool file="filters/bed2gff.xml" />
|
||||
<tool file="filters/axt_to_fasta.xml" />
|
||||
@@ -73,17 +73,17 @@
|
||||
<tool file="extract/twoBitToFa_wrapper.xml" />
|
||||
</section>
|
||||
<section name="Fetch Alignments" id="fetchAlign">
|
||||
<tool file="extract/interval2maf_pairwise.xml" />
|
||||
<tool file="extract/interval2maf.xml" />
|
||||
<tool file="extract/interval_maf_to_merged_fasta.xml" />
|
||||
<tool file="extract/genebed_maf_to_fasta.xml"/>
|
||||
<tool file="filters/maf/maf_stats.xml"/>
|
||||
<tool file="filters/maf/maf_thread_for_species.xml"/>
|
||||
<tool file="filters/maf/maf_limit_to_species.xml"/>
|
||||
<tool file="filters/maf/maf_limit_size.xml"/>
|
||||
<tool file="filters/maf/maf_by_block_number.xml"/>
|
||||
<tool file="filters/maf/maf_reverse_complement.xml"/>
|
||||
<tool file="filters/maf/maf_filter.xml"/>
|
||||
<tool file="maf/interval2maf_pairwise.xml" />
|
||||
<tool file="maf/interval2maf.xml" />
|
||||
<tool file="maf/interval_maf_to_merged_fasta.xml" />
|
||||
<tool file="maf/genebed_maf_to_fasta.xml"/>
|
||||
<tool file="maf/maf_stats.xml"/>
|
||||
<tool file="maf/maf_thread_for_species.xml"/>
|
||||
<tool file="maf/maf_limit_to_species.xml"/>
|
||||
<tool file="maf/maf_limit_size.xml"/>
|
||||
<tool file="maf/maf_by_block_number.xml"/>
|
||||
<tool file="maf/maf_reverse_complement.xml"/>
|
||||
<tool file="maf/maf_filter.xml"/>
|
||||
</section>
|
||||
<section name="Get Genomic Scores" id="scores">
|
||||
<tool file="stats/wiggle_to_simple.xml" />
|
||||
|
||||
@@ -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__()
|
||||
@@ -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__()
|
||||
@@ -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__()
|
||||
@@ -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__()
|
||||
@@ -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__()
|
||||
@@ -1,7 +1,7 @@
|
||||
<tool id="GeneBed_Maf_Fasta2" name="Stitch Gene blocks">
|
||||
<description>given a set of coding exon intervals</description>
|
||||
<command interpreter="python2.4">#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
|
||||
<command interpreter="python2.4">#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
|
||||
</command>
|
||||
<inputs>
|
||||
Executable
+125
@@ -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__()
|
||||
@@ -1,6 +1,6 @@
|
||||
<tool id="Interval2Maf_pairwise1" name="Extract Pairwise MAF blocks">
|
||||
<description>given a set of genomic intervals</description>
|
||||
<command interpreter="python2.4">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</command>
|
||||
<command interpreter="python2.4">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</command>
|
||||
<inputs>
|
||||
<param name="input1" type="data" format="interval" label="Interval File"/>
|
||||
<param name="mafType" type="select" label="Choose MAF source">
|
||||
@@ -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__()
|
||||
+2
-2
@@ -1,7 +1,7 @@
|
||||
<tool id="Interval_Maf_Merged_Fasta2" name="Stitch MAF blocks">
|
||||
<description>given a set of genomic intervals</description>
|
||||
<command interpreter="python2.4">#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
|
||||
<command interpreter="python2.4">#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
|
||||
</command>
|
||||
<inputs>
|
||||
@@ -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."
|
||||
@@ -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__()
|
||||
@@ -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()
|
||||
@@ -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
|
||||
|
||||
Executable
+41
@@ -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__()
|
||||
+9
-16
@@ -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:
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user