diff --git a/tools/extract/genebed_maf_to_fasta_code.py b/tools/extract/genebed_maf_to_fasta_code.py index 211fef66f09..308204aa883 100644 --- a/tools/extract/genebed_maf_to_fasta_code.py +++ b/tools/extract/genebed_maf_to_fasta_code.py @@ -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 [("You must wait for the MAF file to be created before you can use this tool.",'None',True)] diff --git a/tools/extract/interval2maf.xml b/tools/extract/interval2maf.xml index b91ae98f81a..2db4fcece5f 100644 --- a/tools/extract/interval2maf.xml +++ b/tools/extract/interval2maf.xml @@ -24,7 +24,7 @@ - + diff --git a/tools/extract/interval_maf_to_merged_fasta_code.py b/tools/extract/interval_maf_to_merged_fasta_code.py index 00483c7e0d8..e28c7acebc3 100644 --- a/tools/extract/interval_maf_to_merged_fasta_code.py +++ b/tools/extract/interval_maf_to_merged_fasta_code.py @@ -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 [("You must wait for the MAF file to be created before you can use this tool.",'None',True)] diff --git a/tools/filters/maf/maf_by_block_number.xml b/tools/filters/maf/maf_by_block_number.xml index b702f16731b..a1d6e06cc39 100644 --- a/tools/filters/maf/maf_by_block_number.xml +++ b/tools/filters/maf/maf_by_block_number.xml @@ -4,13 +4,16 @@ - + + + + diff --git a/tools/filters/maf/maf_by_block_number_code.py b/tools/filters/maf/maf_by_block_number_code.py new file mode 100644 index 00000000000..38cc90c7277 --- /dev/null +++ b/tools/filters/maf/maf_by_block_number_code.py @@ -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)) diff --git a/tools/filters/maf/maf_limit_to_species.xml b/tools/filters/maf/maf_limit_to_species.xml index 786123cfdb0..e8e308a0062 100644 --- a/tools/filters/maf/maf_limit_to_species.xml +++ b/tools/filters/maf/maf_limit_to_species.xml @@ -16,7 +16,7 @@ - + @@ -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. - + diff --git a/tools/filters/maf/maf_stats.py b/tools/filters/maf/maf_stats.py index 89d8e309ebe..4c4f6eca7d3 100644 --- a/tools/filters/maf/maf_stats.py +++ b/tools/filters/maf/maf_stats.py @@ -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) diff --git a/tools/filters/maf/maf_stats.xml b/tools/filters/maf/maf_stats.xml index fe2af56f9ae..16bebf18856 100644 --- a/tools/filters/maf/maf_stats.xml +++ b/tools/filters/maf/maf_stats.xml @@ -3,9 +3,9 @@ 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 @@ -25,6 +25,10 @@ $maf_source_type.maf_source $maf_source_type.mafType $input1 $out_file1 $dbkey $ + + + + @@ -36,6 +40,7 @@ $maf_source_type.maf_source $maf_source_type.mafType $input1 $out_file1 $dbkey $ + @@ -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. diff --git a/tools/filters/maf/maf_stats_code.py b/tools/filters/maf/maf_stats_code.py index 43501f76a76..a991d10b7a6 100644 --- a/tools/filters/maf/maf_stats_code.py +++ b/tools/filters/maf/maf_stats_code.py @@ -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') diff --git a/tools/filters/maf/maf_thread_for_species.xml b/tools/filters/maf/maf_thread_for_species.xml index bb43d2eed33..77bd998a4f6 100644 --- a/tools/filters/maf/maf_thread_for_species.xml +++ b/tools/filters/maf/maf_thread_for_species.xml @@ -6,7 +6,7 @@ - + @@ -51,6 +51,6 @@ results in:: - + diff --git a/tools/filters/maf/maf_to_bed.xml b/tools/filters/maf/maf_to_bed.xml index 669c2851d3a..4456395db63 100644 --- a/tools/filters/maf/maf_to_bed.xml +++ b/tools/filters/maf/maf_to_bed.xml @@ -6,7 +6,7 @@ - + diff --git a/tools/filters/maf/maf_to_fasta.xml b/tools/filters/maf/maf_to_fasta.xml index c35e4d95fe4..c073f263880 100644 --- a/tools/filters/maf/maf_to_fasta.xml +++ b/tools/filters/maf/maf_to_fasta.xml @@ -15,7 +15,7 @@ - + @@ -23,7 +23,7 @@ - + @@ -184,5 +184,5 @@ will be converted to (**note** that the second MAF block, which does not have mm - +