# HG changeset patch
# User peterjc
# Date 1394114293 18000
# Node ID 0c24e4e2177d89c86643d9789c959bd8b763f568

Uploaded v0.0.3, first stable release.

diff -r 000000000000 -r 0c24e4e2177d test-data/demo_nuc_align.fasta
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/demo_nuc_align.fasta	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,6 @@
+>Alpha
+GATGAGGAACGA
+>Beta
+GATGAG---CGU
+>Gamma
+GAT------CGG
diff -r 000000000000 -r 0c24e4e2177d test-data/demo_nucs.fasta
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/demo_nucs.fasta	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,6 @@
+>Alpha
+GATGAGGAACGA
+>Beta
+GATGAGCGU
+>Gamma
+GATCGG
diff -r 000000000000 -r 0c24e4e2177d test-data/demo_nucs_trailing_stop.fasta
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/demo_nucs_trailing_stop.fasta	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,6 @@
+>Alpha
+GATGAGGAACGA
+>Beta
+GATGAGCGU
+>Gamma
+GATCGGTAA
diff -r 000000000000 -r 0c24e4e2177d test-data/demo_prot_align.fasta
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/demo_prot_align.fasta	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,6 @@
+>Alpha
+DEER
+>Beta
+DE-R
+>Gamma
+D--R
diff -r 000000000000 -r 0c24e4e2177d tools/align_back_trans/README.rst
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/align_back_trans/README.rst	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,108 @@
+Galaxy tool to back-translate a protein alignment to nucleotides
+================================================================
+
+This tool is copyright 2012-2014 by Peter Cock, The James Hutton Institute
+(formerly SCRI, Scottish Crop Research Institute), UK. All rights reserved.
+See the licence text below (MIT licence).
+
+This tool is a short Python script (using Biopython library functions) to
+load a protein alignment, and matching nucleotide FASTA file of unaligned
+sequences, which are threaded onto the protein alignment in order to produce
+a codon aware nucleotide alignment - which can be viewed as a back translation.
+
+This tool is available from the Galaxy Tool Shed at:
+
+* http://toolshed.g2.bx.psu.edu/view/peterjc/align_back_trans
+
+
+Automated Installation
+======================
+
+This should be straightforward using the Galaxy Tool Shed, which should be
+able to automatically install the dependency on Biopython, and then install
+this tool and run its unit tests.
+
+
+Manual Installation
+===================
+
+There are just two files to install to use this tool from within Galaxy:
+
+* ``align_back_trans.py`` (the Python script)
+* ``align_back_trans.xml`` (the Galaxy tool definition)
+
+The suggested location is in a dedicated ``tools/align_back_trans`` folder.
+
+You will also need to modify the ``tools_conf.xml`` file to tell Galaxy to offer
+the tool. One suggested location is in the multiple alignments section. Simply
+add the line::
+
+    <tool file="align_back_trans/align_back_trans.xml" />
+
+You will also need to install Biopython 1.62 or later. If you want to run
+the unit tests, include this line in ``tools_conf.xml.sample`` and the sample
+FASTA files under the ``test-data`` directory. Then::
+
+    ./run_functional_tests.sh -id align_back_trans
+
+That's it.
+
+
+History
+=======
+
+======= ======================================================================
+Version Changes
+------- ----------------------------------------------------------------------
+v0.0.1  - Initial version, based on a previously written Python script
+v0.0.2  - Optionally check the translation is consistent
+======= ======================================================================
+
+
+Developers
+==========
+
+This script was initially developed on this repository:
+https://github.com/peterjc/picobio/blob/master/align/align_back_trans.py
+
+With the addition of a Galaxy wrapper, developement moved here:
+https://github.com/peterjc/pico_galaxy/tree/master/tools/align_back_trans
+
+For making the "Galaxy Tool Shed" http://toolshed.g2.bx.psu.edu/ tarball use
+the following command from the Galaxy root folder::
+
+    $ tar -czf align_back_trans.tar.gz tools/align_back_trans/README.rst tools/align_back_trans/align_back_trans.py tools/align_back_trans/align_back_trans.xml tools/align_back_trans/tool_dependencies.xml test-data/demo_nucs.fasta test-data/demo_nucs_trailing_stop.fasta test-data/demo_prot_align.fasta test-data/demo_nuc_align.fasta
+
+Check this worked::
+
+    $ tar -tzf align_back_trans.tar.gz
+    tools/align_back_trans/README.rst
+    tools/align_back_trans/align_back_trans.py
+    tools/align_back_trans/align_back_trans.xml
+    tools/align_back_trans/tool_dependencies.xml
+    test-data/demo_nucs.fasta
+    test-data/demo_nucs_trailing_stop.fasta
+    test-data/demo_prot_align.fasta
+    test-data/demo_nuc_align.fasta
+
+
+Licence (MIT)
+=============
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff -r 000000000000 -r 0c24e4e2177d tools/align_back_trans/align_back_trans.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/align_back_trans/align_back_trans.py	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,193 @@
+#!/usr/bin/env python
+"""Back-translate a protein alignment to nucleotides
+
+This tool is a short Python script (using Biopython library functions) to
+load a protein alignment, and matching nucleotide FASTA file of unaligned
+sequences, in order to produce a codon aware nucleotide alignment - which
+can be viewed as a back translation.
+
+The development repository for this tool is here:
+
+* https://github.com/peterjc/pico_galaxy/tree/master/tools/align_back_trans  
+
+This tool is available with a Galaxy wrapper from the Galaxy Tool Shed at:
+
+* http://toolshed.g2.bx.psu.edu/view/peterjc/align_back_trans
+
+See accompanying text file for licence details (MIT licence).
+
+This is version 0.0.3 of the script.
+"""
+
+import sys
+from Bio.Seq import Seq
+from Bio.Alphabet import generic_dna, generic_protein
+from Bio.Align import MultipleSeqAlignment
+from Bio import SeqIO
+from Bio import AlignIO
+from Bio.Data.CodonTable import ambiguous_generic_by_id
+
+if "-v" in sys.argv or "--version" in sys.argv:
+    print "v0.0.3"
+    sys.exit(0)
+
+def stop_err(msg, error_level=1):
+    """Print error message to stdout and quit with given error level."""
+    sys.stderr.write("%s\n" % msg)
+    sys.exit(error_level)
+
+def check_trans(identifier, nuc, prot, table):
+    """Returns nucleotide sequence if works (can remove trailing stop)"""
+    if len(nuc) % 3:
+        stop_err("Nucleotide sequence for %s is length %i (not a multiple of three)"
+                 % (identifier, nuc))
+
+    p = str(prot).upper().replace("*", "X")
+    t = str(nuc.translate(table)).upper().replace("*", "X")
+    if len(t) == len(p) + 1:
+        if str(nuc)[-3:].upper() in ambiguous_generic_by_id[table].stop_codons:
+            #Allow this...
+            t = t[:-1]
+            nuc  = nuc[:-3] #edit return value
+    if len(t) != len(p) and p in t:
+        stop_err("%s translation matched but only as subset of nucleotides, "
+                 "wrong start codon?" % identifier)
+    if len(t) != len(p) and p[1:] in t:
+        stop_err("%s translation matched (ignoring first base) but only "
+                 "as subset of nucleotides, wrong start codon?" % identifier)
+    if len(t) != len(p):
+        stop_err("Inconsistent lengths for %s, ungapped protein %i, "
+                 "tripled %i vs ungapped nucleotide %i" %
+                 (identifier,
+                  len(p),
+                  len(p) * 3,
+                  len(nuc)))
+
+
+    if t == p:
+        return nuc
+    elif p.startswith("M") and "M" + t[1:] == p:
+        #Close, was there a start codon?
+        if str(nuc[0:3]).upper() in ambiguous_generic_by_id[table].start_codons:
+            return nuc
+        else:
+            stop_err("Translation check failed for %s\n"
+                     "Would match if %s was a start codon (check correct table used)\n"
+                     % (identifier, nuc[0:3].upper()))
+    else:
+        #Allow * vs X here? e.g. internal stop codons
+        m = "".join("." if x==y else "!" for (x,y) in zip(p,t))
+        if len(prot) < 70:
+            sys.stderr.write("Protein:     %s\n" % p)
+            sys.stderr.write("             %s\n" % m)
+            sys.stderr.write("Translation: %s\n" % t)
+        else:
+            for offset in range(0, len(p), 60):
+                sys.stderr.write("Protein:     %s\n" % p[offset:offset+60])
+                sys.stderr.write("             %s\n" % m[offset:offset+60])
+                sys.stderr.write("Translation: %s\n\n" % t[offset:offset+60])
+        stop_err("Translation check failed for %s\n" % identifier)
+
+def sequence_back_translate(aligned_protein_record, unaligned_nucleotide_record, gap, table=0):
+    #TODO - Separate arguments for protein gap and nucleotide gap?
+    if not gap or len(gap) != 1:
+        raise ValueError("Please supply a single gap character")
+
+    alpha = unaligned_nucleotide_record.seq.alphabet
+    if hasattr(alpha, "gap_char"):
+        gap_codon = alpha.gap_char * 3
+        assert len(gap_codon) == 3
+    else:
+        from Bio.Alphabet import Gapped
+        alpha = Gapped(alpha, gap)
+        gap_codon = gap*3
+
+    ungapped_protein = aligned_protein_record.seq.ungap(gap)
+    ungapped_nucleotide = unaligned_nucleotide_record.seq
+    if table:
+        ungapped_nucleotide = check_trans(aligned_protein_record.id, ungapped_nucleotide, ungapped_protein, table)
+    elif len(ungapped_protein) * 3 != len(ungapped_nucleotide):
+        stop_err("Inconsistent lengths for %s, ungapped protein %i, "
+                 "tripled %i vs ungapped nucleotide %i" %
+                 (aligned_protein_record.id,
+                  len(ungapped_protein),
+                  len(ungapped_protein) * 3,
+                  len(ungapped_nucleotide)))
+
+    seq = []
+    nuc = str(ungapped_nucleotide)
+    for amino_acid in aligned_protein_record.seq:
+        if amino_acid == gap:
+            seq.append(gap_codon)
+        else:
+            seq.append(nuc[:3])
+            nuc = nuc[3:]
+    assert not nuc, "Nucleotide sequence for %r longer than protein %r" \
+        % (unaligned_nucleotide_record.id, aligned_protein_record.id)
+
+    aligned_nuc = unaligned_nucleotide_record[:] #copy for most annotation
+    aligned_nuc.letter_annotation = {} #clear this
+    aligned_nuc.seq = Seq("".join(seq), alpha) #replace this
+    assert len(aligned_protein_record.seq) * 3 == len(aligned_nuc)
+    return aligned_nuc
+
+def alignment_back_translate(protein_alignment, nucleotide_records, key_function=None, gap=None, table=0):
+    """Thread nucleotide sequences onto a protein alignment."""
+    #TODO - Separate arguments for protein and nucleotide gap characters?
+    if key_function is None:
+        key_function = lambda x: x
+    if gap is None:
+        gap = "-"
+
+    aligned = []
+    for protein in protein_alignment:
+        try:
+            nucleotide = nucleotide_records[key_function(protein.id)]
+        except KeyError:
+            raise ValueError("Could not find nucleotide sequence for protein %r" \
+                             % protein.id)
+        aligned.append(sequence_back_translate(protein, nucleotide, gap, table))
+    return MultipleSeqAlignment(aligned)
+
+
+if len(sys.argv) == 4:
+    align_format, prot_align_file, nuc_fasta_file = sys.argv[1:]
+    nuc_align_file = sys.stdout
+    table = 0
+elif len(sys.argv) == 5:
+    align_format, prot_align_file, nuc_fasta_file, nuc_align_file = sys.argv[1:]
+    table = 0
+elif len(sys.argv) == 6:
+    align_format, prot_align_file, nuc_fasta_file, nuc_align_file, table = sys.argv[1:]
+else:
+    stop_err("""This is a Python script for 'back-translating' a protein alignment,
+
+It requires three, four or five arguments:
+- alignment format (e.g. fasta, clustal),
+- aligned protein file (in specified format),
+- unaligned nucleotide file (in fasta format).
+- aligned nucleotiode output file (in same format), optional.
+- NCBI translation table (0 for none), optional
+
+The nucleotide alignment is printed to stdout if no output filename is given.
+
+Example usage:
+
+$ python align_back_trans.py fasta demo_prot_align.fasta demo_nucs.fasta demo_nuc_align.fasta
+
+Warning: If the output file already exists, it will be overwritten.
+
+This script is available with sample data and a Galaxy wrapper here:
+https://github.com/peterjc/pico_galaxy/tree/master/tools/align_back_trans
+http://toolshed.g2.bx.psu.edu/view/peterjc/align_back_trans
+""")
+
+try:
+    table = int(table)
+except:
+    stop_err("Bad table argument %r" % table)
+
+prot_align = AlignIO.read(prot_align_file, align_format, alphabet=generic_protein)
+nuc_dict = SeqIO.index(nuc_fasta_file, "fasta")
+nuc_align = alignment_back_translate(prot_align, nuc_dict, gap="-", table=table)
+AlignIO.write(nuc_align, nuc_align_file, align_format)
diff -r 000000000000 -r 0c24e4e2177d tools/align_back_trans/align_back_trans.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/align_back_trans/align_back_trans.xml	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,131 @@
+<tool id="align_back_trans" name="Thread nucleotides onto a protein alignment (back-translation)" version="0.0.3">
+    <description>Gives a codon aware alignment</description>
+    <requirements>
+        <requirement type="package" version="1.63">biopython</requirement>
+        <requirement type="python-module">Bio</requirement>
+    </requirements>
+    <version_command interpreter="python">align_back_trans.py --version</version_command>
+    <command interpreter="python">
+align_back_trans.py $prot_align.ext "$prot_align" "$nuc_file" "$out_nuc_align" "$table"
+    </command>
+    <stdio>
+        <!-- Anything other than zero is an error -->
+        <exit_code range="1:" />
+        <exit_code range=":-1" />
+    </stdio>
+    <inputs>
+        <param name="prot_align" type="data" format="fasta,muscle,clustal" label="Aligned protein file" help="Mutliple sequence file in FASTA, ClustalW or PHYLIP format." />
+        <param name="table" type="select" label="Genetic code" help="Tables from the NCBI, these determine the start and stop codons">
+            <option value="1">1. Standard</option>
+            <option value="2">2. Vertebrate Mitochondrial</option>
+            <option value="3">3. Yeast Mitochondrial</option>
+            <option value="4">4. Mold, Protozoan, Coelenterate Mitochondrial and Mycoplasma/Spiroplasma</option>
+            <option value="5">5. Invertebrate Mitochondrial</option>
+            <option value="6">6. Ciliate Macronuclear and Dasycladacean</option>
+            <option value="9">9. Echinoderm Mitochondrial</option>
+            <option value="10">10. Euplotid Nuclear</option>
+            <option value="11">11. Bacterial</option>
+            <option value="12">12. Alternative Yeast Nuclear</option>
+            <option value="13">13. Ascidian Mitochondrial</option>
+            <option value="14">14. Flatworm Mitochondrial</option>
+            <option value="15">15. Blepharisma Macronuclear</option>
+            <option value="16">16. Chlorophycean Mitochondrial</option>
+            <option value="21">21. Trematode Mitochondrial</option>
+            <option value="22">22. Scenedesmus obliquus</option>
+            <option value="23">23. Thraustochytrium Mitochondrial</option>
+            <option value="0">Don't check the translation</option>
+        </param>
+        <param name="nuc_file" type="data" format="fasta" label="Unaligned nucleotide sequences" help="FASTA format, using same identifiers as your protein alignment" />
+    </inputs>
+    <outputs>
+        <data name="out_nuc_align" format="fasta" label="${prot_align.name} (back-translated)">
+            <!-- TODO - Replace this with format="input:prot_align" if/when that works -->
+            <change_format>
+                <when input_dataset="prot_align" attribute="extension" value="clustal" format="clustal" />
+                <when input_dataset="prot_align" attribute="extension" value="phylip" format="phylip" />
+            </change_format>
+        </data>
+    </outputs>
+    <tests>
+        <test>
+            <param name="prot_align" value="demo_prot_align.fasta" />
+            <param name="nuc_file" value="demo_nucs.fasta" />
+            <param name="table" value="0" />
+            <output name="out_nuc_align" file="demo_nuc_align.fasta" />
+        </test>
+        <test>
+            <param name="prot_align" value="demo_prot_align.fasta" />
+            <param name="nuc_file" value="demo_nucs_trailing_stop.fasta" />
+            <param name="table" value="11" />
+            <output name="out_nuc_align" file="demo_nuc_align.fasta" />
+        </test>
+    </tests>
+    <help>
+**What it does**
+
+Takes an input file of aligned protein sequences (typically FASTA or Clustal
+format), and a matching file of unaligned nucleotide sequences (FASTA format,
+using the same identifiers), and threads the nucleotide sequences onto the
+protein alignment to produce a codon aware nucleotide alignment - which can
+be viewed as a back translation.
+
+If you specify one of the standard NCBI genetic codes (recommended), then the
+translation is verified. This will allow fuzzy matching if stop codons in the
+protein sequence have been reprented as X, and will allow for a trailing stop
+codon present in the nucleotide sequences but not the protein.
+
+Note - the protein and nucleotide sequences must use the same identifers.
+
+Note - If no translation table is specified, the provided nucleotide sequences
+should be exactly three times the length of the protein sequences (exluding the gaps).
+
+Note - the nucleotide FASTA file may contain extra sequences not in the
+protein alignment, they will be ignored. This can be useful if for example
+you have a nucleotide FASTA file containing all the genes in an organism,
+while the protein alignment is for a specific gene family.
+
+**Example**
+
+Given this protein alignment in FASTA format::
+
+    >Alpha
+    DEER
+    >Beta
+    DE-R
+    >Gamma
+    D--R
+
+and this matching unaligned nucleotide FASTA file::
+
+    >Alpha
+    GATGAGGAACGA
+    >Beta
+    GATGAGCGU
+    >Gamma
+    GATCGG
+
+the tool would return this nucleotide alignment::
+
+    >Alpha
+    GATGAGGAACGA
+    >Beta
+    GATGAG---CGU
+    >Gamma
+    GAT------CGG
+
+Notice that all the gaps are multiples of three in length.
+
+
+**Citation**
+
+This tool uses Biopython, so if you use this Galaxy tool in work leading to a
+scientific publication please cite the following paper:
+
+Cock et al (2009). Biopython: freely available Python tools for computational
+molecular biology and bioinformatics. Bioinformatics 25(11) 1422-3.
+http://dx.doi.org/10.1093/bioinformatics/btp163 pmid:19304878.
+
+This tool is available to install into other Galaxy Instances via the Galaxy
+Tool Shed at http://toolshed.g2.bx.psu.edu/view/peterjc/align_back_trans
+    </help>
+</tool>
diff -r 000000000000 -r 0c24e4e2177d tools/align_back_trans/tool_dependencies.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/align_back_trans/tool_dependencies.xml	Thu Mar 06 08:58:13 2014 -0500
@@ -0,0 +1,6 @@
+<?xml version="1.0"?>
+<tool_dependency>
+    <package name="biopython" version="1.63">
+        <repository changeset_revision="a5c49b83e983" name="package_biopython_1_63" owner="biopython" toolshed="http://toolshed.g2.bx.psu.edu" />
+    </package>
+</tool_dependency>