Update MAF tools to use species as provided by metadata. Selecting parameters for these tools should be significantly faster.

Still a bit of cleanup to do.
This commit is contained in:
Daniel Blankenberg
2007-08-30 18:06:39 +00:00
parent 42a90ea1f0
commit 477aa85f16
12 changed files with 72 additions and 83 deletions
+1 -31
View File
@@ -50,37 +50,7 @@ def get_available_species( maf_source_type ):
return available_sets
else:
try:
rval = []
species={}
input_filename = maf_source_type['input2'].file_name
try:
file_in = open(input_filename, 'r')
maf_reader = bx.align.maf.Reader( file_in )
for i, m in enumerate( maf_reader ):
l = m.components
for c in l:
spec,chrom = bx.align.maf.src_split( c.src )
if not spec or not chrom:
spec = chrom = c.src
if spec not in species:
species[spec]={"bases":0,"nongaps":0}
species[spec]["bases"] = species[spec]["bases"] + c.size + c.text.count("-")
species[spec]["nongaps"] = species[spec]["nongaps"] + c.size
file_in.close()
except Exception:
return [("There is a problem with your MAF file",'None',True)]
species_names = species.keys()
species_names.sort()
for spec in species_names:
#species_sequence[spec] = "".join(species_sequence[spec])
display = "%s: %i nongap, %i total bases" % (spec, species[spec]["nongaps"], species[spec]["bases"] )
rval.append( ( display,spec,True) )
return rval
return map(lambda spec: (spec, spec, True), maf_source_type['input2'].metadata.species)
except:
return [("<B>You must wait for the MAF file to be created before you can use this tool.</B>",'None',True)]
+1 -1
View File
@@ -24,7 +24,7 @@
</page>
</inputs>
<outputs>
<data format="maf" name="out_file1" metadata_source="input1"/>
<data format="maf" name="out_file1"/>
</outputs>
<tests>
<test>
@@ -49,37 +49,7 @@ def get_available_species( maf_source_type ):
return available_sets
else:
try:
rval = []
species={}
input_filename = maf_source_type['input2'].file_name
try:
file_in = open(input_filename, 'r')
maf_reader = bx.align.maf.Reader( file_in )
for i, m in enumerate( maf_reader ):
l = m.components
for c in l:
spec,chrom = bx.align.maf.src_split( c.src )
if not spec or not chrom:
spec = chrom = c.src
if spec not in species:
species[spec]={"bases":0,"nongaps":0}
species[spec]["bases"] = species[spec]["bases"] + c.size + c.text.count("-")
species[spec]["nongaps"] = species[spec]["nongaps"] + c.size
file_in.close()
except Exception:
return [("There is a problem with your MAF file",'None',True)]
species_names = species.keys()
species_names.sort()
for spec in species_names:
#species_sequence[spec] = "".join(species_sequence[spec])
display = "%s: %i nongap, %i total bases" % (spec, species[spec]["nongaps"], species[spec]["bases"] )
rval.append( ( display,spec,True) )
return rval
return map(lambda spec: (spec, spec, True), maf_source_type['input2'].metadata.species)
except:
return [("<B>You must wait for the MAF file to be created before you can use this tool.</B>",'None',True)]
+4 -1
View File
@@ -4,13 +4,16 @@
<inputs>
<page>
<param format="text" name="input1" type="data" label="Block Numbers"/>
<param name="block_col" type="integer" label="Column containing Block number" value="1" help="Use 1 if your file contains only block numbers"/>
<param format="maf" name="input2" label="MAF File" type="data"/>
</page>
<page>
<param name="block_col" type="select" label="Column containing Block number" dynamic_options="get_columns(input1)"/>
</page>
</inputs>
<outputs>
<data format="maf" name="out_file1" />
</outputs>
<code file="maf_by_block_number_code.py"/>
<tests>
<test>
<param name="input1" value="maf_by_block_numbers.dat"/>
@@ -0,0 +1,9 @@
#!/usr/bin/env python2.4
#Dan Blankenberg
"""
Returns valid column numbers, defaulting to 1 for text files
"""
def get_columns(dataset):
if not dataset.metadata.columns: return [("1", "1", True)]
return map(lambda col: (str(col), str(col), False), range(1, dataset.metadata.columns+1))
+2 -2
View File
@@ -16,7 +16,7 @@
</param>
</page>
<page>
<param name="species" type="select" label="Species to keep" dynamic_options="get_available_species( input1.file_name )" display="checkboxes" multiple="true"/>
<param name="species" type="select" label="Species to keep" dynamic_options="map(lambda spec: (spec, spec, True), input1.metadata.species)" display="checkboxes" multiple="true"/>
</page>
</inputs>
<outputs>
@@ -42,6 +42,6 @@ This tool allows the user to remove any undesired species from a MAF file. Colum
* **Exclude blocks with have only one species** - if this option is set to **YES** all single sequence alignment blocks WILL NOT be returned.
</help>
<code file="maf_limit_to_species_code.py"/>
<!--<code file="maf_limit_to_species_code.py"/>-->
</tool>
+26 -9
View File
@@ -22,6 +22,9 @@ def __main__():
chr_col = int(sys.argv[5].strip())-1
start_col = int(sys.argv[6].strip())-1
end_col = int(sys.argv[7].strip())-1
summary = sys.argv[8].strip()
if summary.lower() == "true": summary = True
else: summary = False
index = None
index_filename = None
@@ -79,7 +82,7 @@ def __main__():
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, 'interval2maf.py initialization error -> %s' % exc
print >>sys.stdout, 'maf_stats.py initialization error -> %s' % exc
sys.exit()
else:
print >>sys.stdout, 'Invalid source type specified: %s' % maf_source_type
@@ -88,6 +91,8 @@ def __main__():
out = open(output_filename, 'w')
num_region = 0
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)]}
@@ -96,6 +101,8 @@ def __main__():
start = region.start
end = region.end
total_length += (end - start)
blocks = index.get( src, start, end )
for maf in blocks:
#make sure all species are known
@@ -135,15 +142,25 @@ def __main__():
if i in gaps: gap_offset += 1
elif c.text[i] not in ['-']: sequences[spec][i-gap_offset+slice_start-start] = True
#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) ) )
keys = sequences.keys()
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) ) )
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)
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) ) )
keys = sequences.keys()
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) ) )
num_region += 1
print "%i regions were processed." % num_region
if summary:
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.close()
if index_filename is not None:
os.unlink(index_filename)
+19 -2
View File
@@ -3,9 +3,9 @@
<command interpreter="python2.4">
maf_stats.py
#if $maf_source_type.maf_source == "user":
$maf_source_type.maf_source $input2 $input1 $out_file1 $dbkey $input1_chromCol $input1_startCol $input1_endCol
$maf_source_type.maf_source $input2 $input1 $out_file1 $dbkey $input1_chromCol $input1_startCol $input1_endCol $summary
#else:
$maf_source_type.maf_source $maf_source_type.mafType $input1 $out_file1 $dbkey $input1_chromCol $input1_startCol $input1_endCol
$maf_source_type.maf_source $maf_source_type.mafType $input1 $out_file1 $dbkey $input1_chromCol $input1_startCol $input1_endCol $summary
#end if
</command>
<inputs>
@@ -25,6 +25,10 @@ $maf_source_type.maf_source $maf_source_type.mafType $input1 $out_file1 $dbkey $
<param name="mafType" label="MAF Type" type="select" dynamic_options="get_available_data( input1.dbkey )"/>
</when>
</conditional>
<param name="summary" type="select" label="Type of Output">
<option value="false" selected="true">Coverage by Region</option>
<option value="true">Summarize Coverage</option>
</param>
</page>
</inputs>
<outputs>
@@ -36,6 +40,7 @@ $maf_source_type.maf_source $maf_source_type.mafType $input1 $out_file1 $dbkey $
<param name="maf_source" value="cached"/>
<param name="mafType" value="8_WAY_MULTIZ_hg17"/>
<output name="out_file1" file="maf_stats_interval_out.dat"/>
<param name="summary" value="false"/>
</test>
</tests>
<help>
@@ -61,7 +66,19 @@ Consider the interval: "chrX 1000 1100 myInterval"
YYY = number of gaps
----
Alternatively, you can request only summary information for a set of intervals:
======== =========== ========
#species nucleotides coverage
======== =========== ========
hg18 30639 0.2372
rheMac2 7524 0.0582
panTro2 30390 0.2353
======== =========== ========
where **coverage** is the number of nucleotides divided by the total length of the provided intervals.
</help>
<code file="maf_stats_code.py"/>
</tool>
+3
View File
@@ -45,3 +45,6 @@ def exec_before_job(app, inp_data, out_data, param_dict, tool):
data.name = data.name + " [" + maf_sets[str(param_dict['maf_source_type']['mafType'])]['description'] + "]"
except KeyError:
data.name = data.name + " [unknown MAF source specified]"
if param_dict['summary'].lower() == "true":
for name, data in out_data.items():
data.change_datatype('tabular')
+2 -2
View File
@@ -6,7 +6,7 @@
<param format="maf" name="input1" type="data" label="MAF file"/>
</page>
<page>
<param name="species" type="select" label="Species to keep" dynamic_options="get_available_species( input1.file_name )" display="checkboxes" multiple="true"/>
<param name="species" type="select" label="Species to keep" dynamic_options="map(lambda spec: (spec, spec, True), input1.metadata.species)" display="checkboxes" multiple="true"/>
</page>
</inputs>
<outputs>
@@ -51,6 +51,6 @@ results in::
</help>
<code file="maf_thread_for_species_code.py"/>
<!--<code file="maf_thread_for_species_code.py"/>-->
</tool>
+1 -1
View File
@@ -6,7 +6,7 @@
<param format="maf" name="input1" type="data" label="MAF file to convert"/>
</page>
<page>
<param name="species" type="select" label="Select species" dynamic_options="get_available_species( input1.dbkey, input1.file_name )" display="checkboxes" multiple="true" help="a separate history item will be created for each checked species" />
<param name="species" type="select" label="Select species" dynamic_options="map(lambda spec: (spec, spec, True), input1.metadata.species)" display="checkboxes" multiple="true" help="a separate history item will be created for each checked species" />
<param name="complete_blocks" type="select">
<label>Exclude blocks which have a requested species missing</label>
<option value="partial_allowed">include blocks with missing species</option>
+3 -3
View File
@@ -15,7 +15,7 @@
<option value="concatenated">One Sequence per Species</option>
</param>
<when value="multiple">
<param name="species" type="select" label="Select species" dynamic_options="get_available_species( input1.file_name )" display="checkboxes" multiple="true" help="checked taxa will be included in the output"/>
<param name="species" type="select" label="Select species" dynamic_options="map(lambda spec: (spec, spec, True), input1.metadata.species)" display="checkboxes" multiple="true" help="checked taxa will be included in the output"/>
<param name="complete_blocks" type="select">
<label>Choose to</label>
<option value="partial_allowed">include blocks with missing species</option>
@@ -23,7 +23,7 @@
</param>
</when>
<when value="concatenated">
<param name="species" type="select" label="Species to extract" dynamic_options="get_available_species( input1.file_name )" display="checkboxes" multiple="true"/>
<param name="species" type="select" label="Species to extract" dynamic_options="map(lambda spec: (spec, spec, True), input1.metadata.species)" display="checkboxes" multiple="true"/>
</when>
</conditional>
</page>
@@ -184,5 +184,5 @@ will be converted to (**note** that the second MAF block, which does not have mm
</help>
<code file="maf_to_fasta_code.py"/>
<!--<code file="maf_to_fasta_code.py"/>-->
</tool>