mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-10-25 00:12:08 +08:00
Mostly cleaning up MAF tool code. maf_to_fasta_concat.py no longer stores sequences in memory.
This commit is contained in:
Daniel Blankenberg
@@ -14,32 +14,27 @@ def __main__():
|
||||
input_block_filename = sys.argv[1].strip()
|
||||
input_maf_filename = sys.argv[2].strip()
|
||||
output_filename1 = sys.argv[3].strip()
|
||||
block_col = int(sys.argv[4].strip())-1
|
||||
block_col = int( sys.argv[4].strip() ) - 1
|
||||
if block_col < 0:
|
||||
print >> sys.stderr, "Invalid column specified"
|
||||
sys.exit(0)
|
||||
|
||||
maf_writer = bx.align.maf.Writer( open(output_filename1, 'w') )
|
||||
maf_writer = bx.align.maf.Writer( open( output_filename1, 'w' ) )
|
||||
#we want to maintain order of block file and write blocks as many times as they are listed
|
||||
ctr = 0
|
||||
failed_lines = []
|
||||
for line in open(input_block_filename, 'r'):
|
||||
ctr += 1
|
||||
for ctr, line in enumerate( open( input_block_filename, 'r' ) ):
|
||||
try:
|
||||
block_wanted = int(line.split("\t")[block_col].strip())
|
||||
block_wanted = int( line.split( "\t" )[block_col].strip() )
|
||||
except:
|
||||
failed_lines.append(str(ctr))
|
||||
failed_lines.append( str( ctr ) )
|
||||
continue
|
||||
try:
|
||||
maf_reader = bx.align.maf.Reader( open(input_maf_filename, 'r') )
|
||||
count = 0
|
||||
for m in maf_reader:
|
||||
for count, m in enumerate( bx.align.maf.Reader( open( input_maf_filename, 'r' ) ) ):
|
||||
if count == block_wanted:
|
||||
maf_writer.write(m)
|
||||
maf_writer.write( m )
|
||||
break
|
||||
count += 1
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
if len(failed_lines) > 0: print "Failed to extract from %i lines (%s)." % (len(failed_lines), ",".join(failed_lines))
|
||||
if len( failed_lines ) > 0: print "Failed to extract from %i lines (%s)." % ( len( failed_lines ), ",".join( failed_lines ) )
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
@@ -7,41 +7,41 @@ import bx.align.maf
|
||||
def main():
|
||||
#Read command line arguments
|
||||
try:
|
||||
script_file = sys.argv.pop(1)
|
||||
maf_file = sys.argv.pop(1)
|
||||
out_file = sys.argv.pop(1)
|
||||
additional_files_path = sys.argv.pop(1)
|
||||
script_file = sys.argv.pop( 1 )
|
||||
maf_file = sys.argv.pop( 1 )
|
||||
out_file = sys.argv.pop( 1 )
|
||||
additional_files_path = sys.argv.pop( 1 )
|
||||
except:
|
||||
print >>sys.stderr, "One or more arguments is missing.\nUsage: maf_filter.py maf_filter_file input_maf output_maf path_to_save_debug"
|
||||
sys.exit()
|
||||
|
||||
#Open input and output MAF files
|
||||
try:
|
||||
maf_reader = bx.align.maf.Reader( open(maf_file,'r') )
|
||||
maf_writer = bx.align.maf.Writer( open(out_file,'w') )
|
||||
maf_reader = bx.align.maf.Reader( open( maf_file,'r' ) )
|
||||
maf_writer = bx.align.maf.Writer( open( out_file,'w' ) )
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
|
||||
#Save script file for debuging/verification info later
|
||||
os.mkdir(additional_files_path)
|
||||
shutil.copy(script_file, os.path.join(additional_files_path,'debug.txt'))
|
||||
os.mkdir( additional_files_path )
|
||||
shutil.copy( script_file, os.path.join( additional_files_path, 'debug.txt' ) )
|
||||
|
||||
#Loop through blocks, running filter on each
|
||||
#'maf_block' and 'ret_val' are used/shared in the provided code file
|
||||
#'ret_val' should be set to True if the block is to be kept
|
||||
i = 0
|
||||
blocks_kept = 0
|
||||
for i, maf_block in enumerate(maf_reader):
|
||||
for i, maf_block in enumerate( maf_reader ):
|
||||
local = {'maf_block':maf_block, 'ret_val':False}
|
||||
execfile(script_file, {}, local)
|
||||
execfile( script_file, {}, local )
|
||||
if local['ret_val']:
|
||||
maf_writer.write(maf_block)
|
||||
maf_writer.write( maf_block )
|
||||
blocks_kept += 1
|
||||
maf_writer.close()
|
||||
maf_reader.close()
|
||||
if i == 0: print "Your file contains no valid maf_blocks."
|
||||
else: print 'Kept %s of %s blocks (%.2f%%).' % (blocks_kept, i+1, float(blocks_kept) / float(i+1) * 100.0)
|
||||
else: print 'Kept %s of %s blocks (%.2f%%).' % ( blocks_kept, i + 1, float( blocks_kept ) / float( i + 1 ) * 100.0 )
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
|
||||
@@ -12,20 +12,22 @@ def __main__():
|
||||
|
||||
input_maf_filename = sys.argv[1].strip()
|
||||
output_filename1 = sys.argv[2].strip()
|
||||
min_size = int(sys.argv[3].strip())
|
||||
max_size = int(sys.argv[4].strip())
|
||||
min_size = int( sys.argv[3].strip() )
|
||||
max_size = int( sys.argv[4].strip() )
|
||||
if max_size < 1: max_size = sys.maxint
|
||||
maf_writer = bx.align.maf.Writer( open(output_filename1, 'w') )
|
||||
maf_writer = bx.align.maf.Writer( open( output_filename1, 'w' ) )
|
||||
try:
|
||||
maf_reader = bx.align.maf.Reader( open(input_maf_filename, 'r') )
|
||||
maf_reader = bx.align.maf.Reader( open( input_maf_filename, 'r' ) )
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
maf_count = 0
|
||||
for m in maf_reader:
|
||||
|
||||
blocks_kept = 0
|
||||
i = 0
|
||||
for i, m in enumerate( maf_reader ):
|
||||
if min_size <= m.components[0].size <= max_size:
|
||||
maf_writer.write(m)
|
||||
maf_count += 1
|
||||
print "%i MAF blocks kept." % maf_count
|
||||
maf_writer.write( m )
|
||||
blocks_kept += 1
|
||||
print 'Kept %s of %s blocks (%.2f%%).' % ( blocks_kept, i + 1, float( blocks_kept ) / float( i + 1 ) * 100.0 )
|
||||
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
@@ -7,43 +7,38 @@ columns containing only gaps.
|
||||
|
||||
usage: %prog species,species2,... input_maf output_maf allow_partial min_species_per_block
|
||||
"""
|
||||
|
||||
#Dan Blankenberg
|
||||
import pkg_resources; pkg_resources.require( "bx-python" )
|
||||
import bx.align.maf
|
||||
import sys
|
||||
|
||||
|
||||
def main():
|
||||
|
||||
species = sys.argv[1].split( ',' )
|
||||
try:
|
||||
maf_reader = bx.align.maf.Reader( open(sys.argv[2],'r') )
|
||||
maf_writer = bx.align.maf.Writer( open(sys.argv[3],'w') )
|
||||
maf_reader = bx.align.maf.Reader( open( sys.argv[2],'r' ) )
|
||||
maf_writer = bx.align.maf.Writer( open( sys.argv[3],'w' ) )
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
allow_partial = False
|
||||
if int(sys.argv[4]): allow_partial = True
|
||||
min_species_per_block = int(sys.argv[5])
|
||||
if int( sys.argv[4] ): allow_partial = True
|
||||
min_species_per_block = int( sys.argv[5] )
|
||||
|
||||
maf_blocks_kept = 0
|
||||
|
||||
for m in maf_reader:
|
||||
new_components = []
|
||||
for comp in m.components:
|
||||
if comp.src.split( '.' )[0] in species:
|
||||
new_components.append( comp )
|
||||
m.components = new_components
|
||||
for m in maf_reader:
|
||||
m = m.limit_to_species( species )
|
||||
m.remove_all_gap_columns()
|
||||
output_block = False
|
||||
if allow_partial: output_block = True
|
||||
elif len( m.components ) == len(species): output_block = True
|
||||
elif len( m.components ) == len( species ): output_block = True
|
||||
if output_block and len( m.components ) > min_species_per_block:
|
||||
maf_writer.write( m )
|
||||
maf_blocks_kept += 1
|
||||
|
||||
|
||||
maf_reader.close()
|
||||
maf_writer.close()
|
||||
|
||||
print "%i MAF blocks have been kept." % maf_blocks_kept
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
@@ -11,26 +11,25 @@ import pkg_resources; pkg_resources.require( "bx-python" )
|
||||
import bx.align.maf
|
||||
import sys
|
||||
|
||||
|
||||
def __main__():
|
||||
|
||||
#Parse Command Line
|
||||
input_file = sys.argv.pop(1)
|
||||
output_file = sys.argv.pop(1)
|
||||
input_file = sys.argv.pop( 1 )
|
||||
output_file = sys.argv.pop( 1 )
|
||||
|
||||
try:
|
||||
maf_writer = bx.align.maf.Writer( open(output_file, 'w') )
|
||||
maf_writer = bx.align.maf.Writer( open( output_file, 'w' ) )
|
||||
except:
|
||||
print sys.stderr, "Unable to open output file"
|
||||
sys.exit()
|
||||
try:
|
||||
count = 0
|
||||
for maf in bx.align.maf.Reader( open( input_file ) ):
|
||||
maf_writer.write(maf.reverse_complement())
|
||||
count += 1
|
||||
for count, maf in enumerate( bx.align.maf.Reader( open( input_file ) ) ):
|
||||
maf_writer.write( maf.reverse_complement() )
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
print "%i regions were reverse complemented." % count
|
||||
maf_writer.close()
|
||||
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
@@ -12,16 +12,57 @@ import bx.intervals.io
|
||||
import bx.interval_index_file
|
||||
import psyco_full
|
||||
|
||||
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__():
|
||||
maf_source_type = sys.argv.pop(1)
|
||||
maf_source_type = sys.argv.pop( 1 )
|
||||
input_maf_filename = sys.argv[1].strip()
|
||||
input_interval_filename = sys.argv[2].strip()
|
||||
output_filename = sys.argv[3].strip()
|
||||
dbkey = sys.argv[4].strip()
|
||||
try:
|
||||
chr_col = int(sys.argv[5].strip())-1
|
||||
start_col = int(sys.argv[6].strip())-1
|
||||
end_col = int(sys.argv[7].strip())-1
|
||||
chr_col = int( sys.argv[5].strip() ) - 1
|
||||
start_col = int( sys.argv[6].strip() ) - 1
|
||||
end_col = int( sys.argv[7].strip() ) - 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()
|
||||
@@ -29,63 +70,18 @@ def __main__():
|
||||
if summary.lower() == "true": summary = True
|
||||
else: summary = False
|
||||
|
||||
index = None
|
||||
index_filename = None
|
||||
index = index_filename = None
|
||||
if maf_source_type == "user":
|
||||
#index maf for use here
|
||||
indexes = bx.interval_index_file.Indexes()
|
||||
try:
|
||||
maf_reader = bx.align.maf.Reader( open( input_maf_filename ) )
|
||||
# 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:
|
||||
index, index_filename = build_maf_index( input_maf_filename, species = [dbkey] )
|
||||
if index is None:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
index = bx.align.maf.Indexed(input_maf_filename, index_filename = index_filename, keep_open=True)
|
||||
elif maf_source_type == "cached":
|
||||
#access existing indexes
|
||||
try:
|
||||
maf_sets = {}
|
||||
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_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
|
||||
index = bx.align.maf.MultiIndexed( maf_sets[input_maf_filename]['paths'], keep_open=True, parse_e_rows=True )
|
||||
except Exception, exc:
|
||||
print >>sys.stdout, 'maf_stats.py initialization error -> %s' % exc
|
||||
index = 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()
|
||||
else:
|
||||
print >>sys.stdout, 'Invalid source type specified: %s' % maf_source_type
|
||||
@@ -97,8 +93,8 @@ def __main__():
|
||||
species_summary = {}
|
||||
total_length = 0
|
||||
#loop through interval file
|
||||
for region in bx.intervals.io.NiceReaderWrapper( open(input_interval_filename, 'r' ), chrom_col=chr_col, start_col=start_col, end_col=end_col, fix_strand=True, return_header=False, return_comments=False):
|
||||
sequences = {dbkey: [ False for i in range( region.end - region.start)]}
|
||||
for region in bx.intervals.io.NiceReaderWrapper( open( input_interval_filename, 'r' ), chrom_col = chr_col, start_col = start_col, end_col = end_col, fix_strand = True, return_header = False, return_comments = False ):
|
||||
sequences = {dbkey: [ False for i in range( region.end - region.start ) ]}
|
||||
|
||||
src = dbkey + "." + region.chrom
|
||||
start = region.start
|
||||
@@ -110,14 +106,14 @@ def __main__():
|
||||
for maf in blocks:
|
||||
#make sure all species are known
|
||||
for c in maf.components:
|
||||
spec = bx.align.src_split(c.src)[0]
|
||||
if spec not in sequences: sequences[spec] = [ False for i in range( region.end - region.start)]
|
||||
spec = c.src.split( '.' )[0]
|
||||
if spec not in sequences: sequences[spec] = [ False for i in range( region.end - region.start ) ]
|
||||
|
||||
#save old score here for later use, since slice results score==0
|
||||
old_score = maf.score
|
||||
|
||||
#slice maf by start and end
|
||||
ref = maf.get_component_by_src( bx.align.src_merge(dbkey, region.chrom) )
|
||||
ref = maf.get_component_by_src( bx.align.src_merge( dbkey, region.chrom ) )
|
||||
# If the reference component is on the '-' strand we should complement the interval
|
||||
if ref.strand == '-':
|
||||
slice_start = max( ref.src_size - end, ref.start )
|
||||
@@ -133,38 +129,38 @@ def __main__():
|
||||
ref = sliced.get_component_by_src( ref.src )
|
||||
#look for gaps (indels) in primary sequence, we do not include these columns in our stats
|
||||
gaps = []
|
||||
for i in range(len(ref.text)):
|
||||
for i in range( len( ref.text ) ):
|
||||
if ref.text[i] in ['-']:
|
||||
gaps.append(i)
|
||||
gaps.append( i )
|
||||
|
||||
#Set nucleotide containing columns
|
||||
for c in sliced.components:
|
||||
spec = bx.align.src_split(c.src)[0]
|
||||
spec = c.src.split( '.' )[0]
|
||||
gap_offset = 0
|
||||
for i in range(len(c.text)):
|
||||
for i in range( len( c.text ) ):
|
||||
if i in gaps: gap_offset += 1
|
||||
elif c.text[i] not in ['-']: sequences[spec][i-gap_offset+slice_start-start] = True
|
||||
elif c.text[i] not in ['-']: sequences[spec][i - gap_offset + slice_start - start] = True
|
||||
|
||||
if summary:
|
||||
#record summary
|
||||
for key in sequences.keys():
|
||||
if key not in species_summary: species_summary[key] = 0
|
||||
species_summary[key] = species_summary[key] + sequences[key].count(True)
|
||||
species_summary[key] = species_summary[key] + sequences[key].count( True )
|
||||
else:
|
||||
#print sequences
|
||||
out.write("%s\t%s\t%s\t%s\n" % ( "\t".join(region.fields), dbkey, sequences[dbkey].count(True), sequences[dbkey].count(False) ) )
|
||||
out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), dbkey, sequences[dbkey].count( True ), sequences[dbkey].count( False ) ) )
|
||||
keys = sequences.keys()
|
||||
keys.remove(dbkey)
|
||||
keys.remove( dbkey )
|
||||
keys.sort()
|
||||
for key in keys:
|
||||
out.write("%s\t%s\t%s\t%s\n" % ( "\t".join(region.fields), key, sequences[key].count(True), sequences[key].count(False) ) )
|
||||
out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), key, sequences[key].count( True ), sequences[key].count( False ) ) )
|
||||
num_region += 1
|
||||
if summary:
|
||||
out.write("#species\tnucleotides\tcoverage\n")
|
||||
out.write( "#species\tnucleotides\tcoverage\n" )
|
||||
for spec in species_summary:
|
||||
out.write("%s\t%s\t%.4f\n" % ( spec, species_summary[spec], float(species_summary[spec]) / total_length ))
|
||||
print "%i regions were processed with a total length of %i." % (num_region, total_length)
|
||||
out.write( "%s\t%s\t%.4f\n" % ( spec, species_summary[spec], float( species_summary[spec] ) / total_length ) )
|
||||
print "%i regions were processed with a total length of %i." % ( num_region, total_length )
|
||||
out.close()
|
||||
if index_filename is not None:
|
||||
os.unlink(index_filename)
|
||||
os.unlink( index_filename )
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
@@ -8,7 +8,7 @@ which are adjacent after the unwanted species have been dropped.
|
||||
|
||||
usage: %prog input_maf output_maf species1,species2
|
||||
"""
|
||||
|
||||
#Dan Blankenberg
|
||||
import sys
|
||||
import pkg_resources; pkg_resources.require( "bx-python" )
|
||||
import bx.align.maf
|
||||
@@ -17,18 +17,17 @@ from bx.align.tools.thread import *
|
||||
from bx.align.tools.fuse import *
|
||||
|
||||
def main():
|
||||
|
||||
input_file = sys.argv.pop(1)
|
||||
output_file = sys.argv.pop(1)
|
||||
species = sys.argv.pop(1).split(',')
|
||||
input_file = sys.argv.pop( 1 )
|
||||
output_file = sys.argv.pop( 1 )
|
||||
species = sys.argv.pop( 1 ).split( ',' )
|
||||
|
||||
try:
|
||||
maf_reader = bx.align.maf.Reader( open(input_file) )
|
||||
maf_reader = bx.align.maf.Reader( open( input_file ) )
|
||||
except:
|
||||
print >> sys.stderr, "Unable to open source MAF file"
|
||||
sys.exit()
|
||||
try:
|
||||
maf_writer = FusingAlignmentWriter( bx.align.maf.Writer( open(output_file, 'w') ) )
|
||||
maf_writer = FusingAlignmentWriter( bx.align.maf.Writer( open( output_file, 'w' ) ) )
|
||||
except:
|
||||
print >> sys.stderr, "Unable to open output file"
|
||||
sys.exit()
|
||||
|
||||
@@ -27,12 +27,12 @@ def __main__():
|
||||
if "None" in species:
|
||||
species = {}
|
||||
try:
|
||||
for i, m in enumerate( maf.Reader( open(input_filename, 'r') ) ):
|
||||
for i, m in enumerate( maf.Reader( open( input_filename, 'r' ) ) ):
|
||||
for c in m.components:
|
||||
spec,chrom = maf.src_split( c.src )
|
||||
if not spec or not chrom:
|
||||
spec = chrom = c.src
|
||||
species[spec]=""
|
||||
species[spec] = ""
|
||||
species = species.keys()
|
||||
except:
|
||||
print >>sys.stderr, "Invalid MAF file specified"
|
||||
@@ -43,24 +43,24 @@ def __main__():
|
||||
return
|
||||
|
||||
|
||||
for i in range(0,len(species)):
|
||||
for i in range( 0, len( species ) ):
|
||||
spec = species[i]
|
||||
if i==0:
|
||||
out_files[spec] = open(output_filename, 'w')
|
||||
if i == 0:
|
||||
out_files[spec] = open( output_filename, 'w' )
|
||||
primary_spec = spec
|
||||
else:
|
||||
out_files[spec] = tempfile.NamedTemporaryFile(mode='w',dir=database_tmp_dir, suffix='.maf_to_bed')
|
||||
out_files[spec] = tempfile.NamedTemporaryFile( mode = 'w', dir = database_tmp_dir, suffix = '.maf_to_bed' )
|
||||
filename = out_files[spec].name
|
||||
out_files[spec].close()
|
||||
out_files[spec] = open(filename, 'w')
|
||||
num_species = len(species)
|
||||
out_files[spec] = open( filename, 'w' )
|
||||
num_species = len( species )
|
||||
|
||||
print "Restricted to species:",",".join(species)
|
||||
print "Restricted to species:", ",".join( species )
|
||||
|
||||
file_in = open(input_filename, 'r')
|
||||
file_in = open( input_filename, 'r' )
|
||||
maf_reader = maf.Reader( file_in )
|
||||
|
||||
block_num=-1
|
||||
block_num = -1
|
||||
|
||||
for i, m in enumerate( maf_reader ):
|
||||
block_num += 1
|
||||
@@ -73,15 +73,15 @@ def __main__():
|
||||
if not spec or not chrom:
|
||||
spec = chrom = c.src
|
||||
if spec not in out_files.keys():
|
||||
out_files[spec] = tempfile.NamedTemporaryFile(mode='w',dir=database_tmp_dir, suffix='.maf_to_bed')
|
||||
out_files[spec] = tempfile.NamedTemporaryFile( mode='w', dir = database_tmp_dir, suffix = '.maf_to_bed' )
|
||||
filename = out_files[spec].name
|
||||
out_files[spec].close()
|
||||
out_files[spec] = open(filename, 'w')
|
||||
out_files[spec] = open( filename, 'w' )
|
||||
|
||||
if c.strand == "-":
|
||||
out_files[spec].write(chrom+"\t"+str(c.src_size - c.end)+"\t"+str(c.src_size - c.start)+"\t"+spec+"_"+str(block_num)+"\t"+"0\t"+c.strand+"\n")
|
||||
out_files[spec].write( chrom + "\t" + str( c.src_size - c.end ) + "\t" + str( c.src_size - c.start ) + "\t" + spec + "_" + str( block_num ) + "\t" + "0\t" + c.strand + "\n" )
|
||||
else:
|
||||
out_files[spec].write(chrom+"\t"+str(c.start)+"\t"+str(c.end)+"\t"+spec+"_"+str(block_num)+"\t"+"0\t"+c.strand+"\n")
|
||||
out_files[spec].write( chrom + "\t" + str( c.start ) + "\t" + str( c.end ) + "\t" + spec + "_" + str( block_num ) + "\t" + "0\t" + c.strand + "\n" )
|
||||
|
||||
file_in.close()
|
||||
for file_out in out_files.keys():
|
||||
@@ -89,8 +89,8 @@ def __main__():
|
||||
|
||||
for spec in out_files.keys():
|
||||
if spec != primary_spec:
|
||||
print "#FILE\t"+spec+"\t"+os.path.join(database_tmp_dir,os.path.split(out_files[spec].name)[1])
|
||||
print "#FILE\t" + spec + "\t" + os.path.join( database_tmp_dir, os.path.split( out_files[spec].name )[1] )
|
||||
else:
|
||||
print "#FILE1\t"+spec+"\t"+out_files[spec].name
|
||||
print "#FILE1\t" + spec + "\t" + out_files[spec].name
|
||||
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
@@ -5,7 +5,7 @@ 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
|
||||
@@ -14,58 +14,45 @@ import pkg_resources; pkg_resources.require( "bx-python" )
|
||||
from bx.align import maf
|
||||
|
||||
def __main__():
|
||||
print "Restricted to species:",sys.argv[1]
|
||||
print "Restricted to species:", sys.argv[1]
|
||||
|
||||
texts = {}
|
||||
|
||||
input_filename = sys.argv[2]
|
||||
output_filename = sys.argv[3]
|
||||
species = sys.argv[1].split(',')
|
||||
species = sys.argv[1].split( ',' )
|
||||
|
||||
if "None" in species:
|
||||
species = get_species(input_filename)
|
||||
species = get_species( input_filename )
|
||||
|
||||
file_in = open(input_filename, 'r')
|
||||
file_out = open(output_filename, 'w')
|
||||
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()
|
||||
|
||||
for s in species: texts[s] = []
|
||||
try:
|
||||
maf_reader = maf.Reader( file_in )
|
||||
for m in maf_reader:
|
||||
for s in species:
|
||||
c = m.get_component_by_src_start( s )
|
||||
if c: texts[s].append( c.text )
|
||||
else: texts[s].append( "-" * m.text_size )
|
||||
for s in species:
|
||||
file_out.write(">"+s+"\n")
|
||||
file_out.write("".join(texts[s])+"\n")
|
||||
|
||||
file_in.close()
|
||||
file_out.close()
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
|
||||
def print_n( s, n, f = sys.stdout ):
|
||||
p = 0
|
||||
while p < len( s ):
|
||||
print >> f, s[p:min(p+n,len(s))]
|
||||
p += n
|
||||
|
||||
def get_species(maf_filename):
|
||||
def get_species( maf_filename ):
|
||||
try:
|
||||
species={}
|
||||
|
||||
file_in = open(maf_filename, 'r')
|
||||
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 )
|
||||
spec, chrom = maf.src_split( c.src )
|
||||
if not spec or not chrom:
|
||||
spec = chrom = c.src
|
||||
species[spec]=spec
|
||||
species[spec] = spec
|
||||
|
||||
file_in.close()
|
||||
|
||||
@@ -75,5 +62,4 @@ def get_species(maf_filename):
|
||||
except:
|
||||
return []
|
||||
|
||||
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
@@ -3,7 +3,7 @@
|
||||
"""
|
||||
Read a maf and print the text as a fasta file.
|
||||
"""
|
||||
|
||||
#Dan Blankenberg
|
||||
from __future__ import division
|
||||
|
||||
import textwrap
|
||||
@@ -12,21 +12,21 @@ import pkg_resources; pkg_resources.require( "bx-python" )
|
||||
from bx.align import maf
|
||||
|
||||
def __main__():
|
||||
print "Restricted to species:",sys.argv[3]
|
||||
print "Restricted to species:", sys.argv[3]
|
||||
|
||||
input_filename = sys.argv[1]
|
||||
output_filename = sys.argv[2]
|
||||
species = sys.argv[3].split(',')
|
||||
species = sys.argv[3].split( ',' )
|
||||
partial = sys.argv[4]
|
||||
num_species = len(species)
|
||||
num_species = len( species )
|
||||
|
||||
file_in = open(input_filename, 'r')
|
||||
file_in = open( input_filename, 'r' )
|
||||
try:
|
||||
maf_reader = maf.Reader( file_in )
|
||||
|
||||
file_out = open(output_filename, 'w')
|
||||
file_out = open( output_filename, 'w' )
|
||||
|
||||
block_num=-1
|
||||
block_num = -1
|
||||
|
||||
for i, m in enumerate( maf_reader ):
|
||||
block_num += 1
|
||||
@@ -35,22 +35,16 @@ def __main__():
|
||||
l = m.components
|
||||
if len(l) < num_species and partial == "partial_disallowed": continue
|
||||
for c in l:
|
||||
spec,chrom = maf.src_split( c.src )
|
||||
spec, chrom = maf.src_split( c.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")
|
||||
file_out.write("\n")
|
||||
file_out.write( ">" + c.src + "(" + c.strand + "):" + str( c.start ) + "-" + str( c.end ) + "|" + spec + "_" + str( block_num ) + "\n" )
|
||||
file_out.write( c.text + "\n" )
|
||||
file_out.write( "\n" )
|
||||
file_in.close()
|
||||
except:
|
||||
print >>sys.stderr, "Your MAF file appears to be malformed."
|
||||
sys.exit()
|
||||
file_out.close()
|
||||
|
||||
def print_n( s, n, f = sys.stdout ):
|
||||
p = 0
|
||||
while p < len( s ):
|
||||
print >> f, s[p:min(p+n,len(s))]
|
||||
p += n
|
||||
|
||||
if __name__ == "__main__": __main__()
|
||||
|
||||
Reference in New Issue
Block a user