diff --git a/tool_conf.xml.main b/tool_conf.xml.main index 62aad6782e5..abaf3d994f6 100644 --- a/tool_conf.xml.main +++ b/tool_conf.xml.main @@ -40,6 +40,7 @@ + diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index da52ee0cd9d..0349ce81980 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -75,6 +75,7 @@ + diff --git a/tools/maf/maf_reverse_complement.xml b/tools/maf/maf_reverse_complement.xml index 5602a0bc600..0eb16d0e81a 100644 --- a/tools/maf/maf_reverse_complement.xml +++ b/tools/maf/maf_reverse_complement.xml @@ -1,4 +1,4 @@ - + a MAF file maf_reverse_complement.py $input1 $out_file1 $species diff --git a/tools/maf/maf_to_interval.py b/tools/maf/maf_to_interval.py new file mode 100644 index 00000000000..52fb2862cbc --- /dev/null +++ b/tools/maf/maf_to_interval.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python + +""" +Read a maf and output intervals for specified list of species. +""" +import sys, os +from galaxy import eggs +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.align import maf +from galaxy.tools.util import maf_utilities + +assert sys.version_info[:2] >= ( 2, 4 ) + +def __main__(): + input_filename = sys.argv[1] + output_filename = sys.argv[2] + output_id = sys.argv[3] + #where to store files that become additional output + database_tmp_dir = sys.argv[4] + primary_spec = sys.argv[5] + species = sys.argv[6].split( ',' ) + all_species = sys.argv[7].split( ',' ) + partial = sys.argv[8] + keep_gaps = sys.argv[9] + out_files = {} + + if "None" in species: + species = [] + + if primary_spec not in species: + species.append( primary_spec ) + if primary_spec not in all_species: + all_species.append( primary_spec ) + + all_species.sort() + for spec in species: + if spec == primary_spec: + out_files[ spec ] = open( output_filename, 'wb+' ) + else: + out_files[ spec ] = open( os.path.join( database_tmp_dir, 'primary_%s_%s_visible_interval_%s' % ( output_id, spec, spec ) ), 'wb+' ) + out_files[ spec ].write( '#chrom\tstart\tend\tstrand\tscore\tname\t%s\n' % ( '\t'.join( all_species ) ) ) + num_species = len( all_species ) + + file_in = open( input_filename, 'r' ) + maf_reader = maf.Reader( file_in ) + + for i, m in enumerate( maf_reader ): + for j, block in enumerate( maf_utilities.iter_blocks_split_by_species( m ) ): + if len( block.components ) < num_species and partial == "partial_disallowed": continue + sequences = {} + for c in block.components: + spec, chrom = maf_utilities.src_split( c.src ) + if keep_gaps == 'remove_gaps': + sequences[ spec ] = c.text.replace( '-', '' ) + else: + sequences[ spec ] = c.text + sequences = '\t'.join( [ sequences.get( spec, '' ) for spec in all_species ] ) + for spec in species: + c = block.get_component_by_src_start( spec ) + if c is not None: + spec2, chrom = maf_utilities.src_split( c.src ) + assert spec2 == spec, Exception( 'Species name inconsistancy found in component: %s != %s' % ( spec, spec2 ) ) + out_files[ spec ].write( "%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % ( chrom, c.forward_strand_start, c.forward_strand_end, c.strand, m.score, "%s_%s_%s" % (spec, i, j), sequences ) ) + file_in.close() + for file_out in out_files.values(): + file_out.close() + +if __name__ == "__main__": __main__() diff --git a/tools/maf/maf_to_interval.xml b/tools/maf/maf_to_interval.xml new file mode 100644 index 00000000000..8704481c04e --- /dev/null +++ b/tools/maf/maf_to_interval.xml @@ -0,0 +1,127 @@ + + Converts a MAF formated file to the Interval format + maf_to_interval.py $input1 $out_file1 $out_file1.id $__new_file_path__ $input1.dbkey $species $input1.metadata.species $complete_blocks $remove_gaps + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +**What it does** + +This tool converts every MAF block to a set of genomic intervals describing the position of that alignment block within a corresponding genome. Sequences from aligning species are also included in the output. + +The interface for this tool contains several options: + + * **MAF file to convert**. Choose multiple alignments from history to be converted to BED format. + * **Choose species**. Choose additional species from the alignment to be included in the output + * **Exclude blocks which have a species missing**. if an alignment block does not contain any one of the species found in the alignment set and this option is set to **exclude blocks with missing species**, then coordinates of such a block **will not** be included in the output (see **Example 2** below). + * **Remove Gap characters from sequences**. Gaps can be removed from sequences before they are output. + + +----- + +**Example 1**: **Include only reference genome** (hg18 in this case) and **include blocks with missing species**: + +For the following alignment:: + + ##maf version=1 + a score=68686.000000 + s hg18.chr20 56827368 75 + 62435964 GACAGGGTGCATCTGGGAGGG---CCTGCCGGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- + s panTro2.chr20 56528685 75 + 62293572 GACAGGGTGCATCTGAGAGGG---CCTGCCAGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- + s rheMac2.chr10 89144112 69 - 94855758 GACAGGGTGCATCTGAGAGGG---CCTGCTGGGCCTTTG-TTCAAAACTAGATATGCCCCAACTCCAATTCTA------- + s mm8.chr2 173910832 61 + 181976762 AGAAGGATCCACCT------------TGCTGGGCCTCTGCTCCAGCAAGACCCACCTCCCAACTCAAATGCCC------- + s canFam2.chr24 46551822 67 + 50763139 CG------GCGTCTGTAAGGGGCCACCGCCCGGCCTGTG-CTCAAAGCTACAAATGACTCAACTCCCAACCGA------C + + a score=10289.000000 + s hg18.chr20 56827443 37 + 62435964 ATGTGCAGAAAATGTGATACAGAAACCTGCAGAGCAG + s panTro2.chr20 56528760 37 + 62293572 ATGTGCAGAAAATGTGATACAGAAACCTGCAGAGCAG + s rheMac2.chr10 89144181 37 - 94855758 ATGTGCGGAAAATGTGATACAGAAACCTGCAGAGCAG + +the tool will create **a single** history item containing the following (**note** the name field is numbered iteratively: hg18_0_0, hg18_1_0 etc. where the first number is the block number and the second number is the iteration through the block (if a species appears twice in a block, that interval will be repeated) and sequences for each species are included in the order specified in the header: the field is left empty when no sequence is available for that species):: + + #chrom start end strand score name canFam2 hg18 mm8 panTro2 rheMac2 + chr20 56827368 56827443 + 68686.0 hg18_0_0 CG------GCGTCTGTAAGGGGCCACCGCCCGGCCTGTG-CTCAAAGCTACAAATGACTCAACTCCCAACCGA------C GACAGGGTGCATCTGGGAGGG---CCTGCCGGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- AGAAGGATCCACCT------------TGCTGGGCCTCTGCTCCAGCAAGACCCACCTCCCAACTCAAATGCCC------- GACAGGGTGCATCTGAGAGGG---CCTGCCAGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- GACAGGGTGCATCTGAGAGGG---CCTGCTGGGCCTTTG-TTCAAAACTAGATATGCCCCAACTCCAATTCTA------- + chr20 56827443 56827480 + 10289.0 hg18_1_0 ATGTGCAGAAAATGTGATACAGAAACCTGCAGAGCAG ATGTGCAGAAAATGTGATACAGAAACCTGCAGAGCAG ATGTGCGGAAAATGTGATACAGAAACCTGCAGAGCAG + + +----- + +**Example 2**: **Include hg18 and mm8** and **exclude blocks with missing species**: + +For the following alignment:: + + ##maf version=1 + a score=68686.000000 + s hg18.chr20 56827368 75 + 62435964 GACAGGGTGCATCTGGGAGGG---CCTGCCGGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- + s panTro2.chr20 56528685 75 + 62293572 GACAGGGTGCATCTGAGAGGG---CCTGCCAGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- + s rheMac2.chr10 89144112 69 - 94855758 GACAGGGTGCATCTGAGAGGG---CCTGCTGGGCCTTTG-TTCAAAACTAGATATGCCCCAACTCCAATTCTA------- + s mm8.chr2 173910832 61 + 181976762 AGAAGGATCCACCT------------TGCTGGGCCTCTGCTCCAGCAAGACCCACCTCCCAACTCAAATGCCC------- + s canFam2.chr24 46551822 67 + 50763139 CG------GCGTCTGTAAGGGGCCACCGCCCGGCCTGTG-CTCAAAGCTACAAATGACTCAACTCCCAACCGA------C + + a score=10289.000000 + s hg18.chr20 56827443 37 + 62435964 ATGTGCAGAAAATGTGATACAGAAACCTGCAGAGCAG + s panTro2.chr20 56528760 37 + 62293572 ATGTGCAGAAAATGTGATACAGAAACCTGCAGAGCAG + s rheMac2.chr10 89144181 37 - 94855758 ATGTGCGGAAAATGTGATACAGAAACCTGCAGAGCAG + +the tool will create **two** history items (one for hg18 and one for mm8) containing the following (**note** that both history items contain only one line describing the first alignment block. The second MAF block is not included in the output because it does not contain mm8): + +History item **1** (for hg18):: + + #chrom start end strand score name canFam2 hg18 mm8 panTro2 rheMac2 + chr20 56827368 56827443 + 68686.0 hg18_0_0 CG------GCGTCTGTAAGGGGCCACCGCCCGGCCTGTG-CTCAAAGCTACAAATGACTCAACTCCCAACCGA------C GACAGGGTGCATCTGGGAGGG---CCTGCCGGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- AGAAGGATCCACCT------------TGCTGGGCCTCTGCTCCAGCAAGACCCACCTCCCAACTCAAATGCCC------- GACAGGGTGCATCTGAGAGGG---CCTGCCAGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- GACAGGGTGCATCTGAGAGGG---CCTGCTGGGCCTTTG-TTCAAAACTAGATATGCCCCAACTCCAATTCTA------- + + +History item **2** (for mm8):: + + #chrom start end strand score name canFam2 hg18 mm8 panTro2 rheMac2 + chr2 173910832 173910893 + 68686.0 mm8_0_0 CG------GCGTCTGTAAGGGGCCACCGCCCGGCCTGTG-CTCAAAGCTACAAATGACTCAACTCCCAACCGA------C GACAGGGTGCATCTGGGAGGG---CCTGCCGGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- AGAAGGATCCACCT------------TGCTGGGCCTCTGCTCCAGCAAGACCCACCTCCCAACTCAAATGCCC------- GACAGGGTGCATCTGAGAGGG---CCTGCCAGGCCTTTA-TTCAACACTAGATACGCCCCATCTCCAATTCTAATGGAC- GACAGGGTGCATCTGAGAGGG---CCTGCTGGGCCTTTG-TTCAAAACTAGATATGCCCCAACTCCAATTCTA------- + + +------- + +.. class:: infomark + +**About formats** + +**MAF format** multiple alignment format file. This format stores multiple alignments at the DNA level between entire genomes. + + - The .maf format is line-oriented. Each multiple alignment ends with a blank line. + - Each sequence in an alignment is on a single line. + - Lines starting with # are considered to be comments. + - Each multiple alignment is in a separate paragraph that begins with an "a" line and contains an "s" line for each sequence in the multiple alignment. + - Some MAF files may contain two optional line types: + + - An "i" line containing information about what is in the aligned species DNA before and after the immediately preceding "s" line; + - An "e" line containing information about the size of the gap between the alignments that span the current block. + + + + +