# Copyright 2009 by Peter Cock.  All rights reserved.
# This code is part of the Biopython distribution and governed by its
# license.  Please see the LICENSE file that should have been included
# as part of this package.
"""Runs a few EMBOSS tools to check our wrappers and parsers."""

import os
import sys
import unittest
import subprocess
from Bio._py3k import StringIO

from Bio.Emboss.Applications import WaterCommandline, NeedleCommandline
from Bio.Emboss.Applications import SeqretCommandline, SeqmatchallCommandline
from Bio import SeqIO
from Bio import AlignIO
from Bio import MissingExternalDependencyError
from Bio.Application import _escape_filename
from Bio.Alphabet import generic_protein, generic_dna, generic_nucleotide
from Bio.Seq import Seq, translate
from Bio.SeqRecord import SeqRecord
#from Bio.Data.IUPACData import ambiguous_dna_letters

#################################################################

#Try to avoid problems when the OS is in another language
os.environ['LANG'] = 'C'

exes_wanted = ["water", "needle", "seqret", "transeq", "seqmatchall",
               "embossversion"]
exes = dict()  # Dictionary mapping from names to exe locations

if "EMBOSS_ROOT" in os.environ:
    #Windows default installation path is C:\mEMBOSS which contains the exes.
    #EMBOSS also sets an environment variable which we will check for.
    path = os.environ["EMBOSS_ROOT"]
    if os.path.isdir(path):
        for name in exes_wanted:
            if os.path.isfile(os.path.join(path, name+".exe")):
                exes[name] = os.path.join(path, name+".exe")
        del name
    else:
        raise MissingExternalDependencyError(
                  "$EMBOSS_ROOT=%r which does not exist!" % path)
    del path
if sys.platform!="win32":
    from Bio._py3k import getoutput
    for name in exes_wanted:
        #This will "just work" if installed on the path as normal on Unix
        output = getoutput("%s -help" % name)
        if "not found" not in output and "not recognized" not in output:
            exes[name] = name
        del output
    del name

if len(exes) < len(exes_wanted):
    raise MissingExternalDependencyError(
        "Install EMBOSS if you want to use Bio.Emboss.")


def get_emboss_version():
    """Returns a tuple of three ints, e.g. (6,1,0)"""
    #Windows and Unix versions of EMBOSS seem to differ in
    #which lines go to stdout and stderr - so merge them.
    child = subprocess.Popen(_escape_filename(exes["embossversion"]),
                             stdout=subprocess.PIPE,
                             stderr=subprocess.STDOUT,
                             universal_newlines=True,
                             shell=(sys.platform!="win32"))
    stdout, stderr = child.communicate()
    child.stdout.close()  # This is both stdout and stderr
    del child
    assert stderr is None  # Send to stdout instead
    for line in stdout.split("\n"):
        if line.strip()=="Report the current EMBOSS version number":
            # e.g.
            # $ embossversion
            # Report the current EMBOSS version number
            # 6.5.7.0
            pass
        elif line.strip()=="Reports the current EMBOSS version number":
            # e.g.
            # $ embossversion
            # Reports the current EMBOSS version number
            # 6.3.1
            pass
        elif line.startswith("Writes the current EMBOSS version number"):
            pass
        elif line.count(".")==2:
            return tuple(int(v) for v in line.strip().split("."))
        elif line.count(".")==3:
            #e.g. I installed mEMBOSS-6.2.0.1-setup.exe
            #which reports 6.2.0.1 - for this return (6,2,0)
            return tuple(int(v) for v in line.strip().split("."))[:3]
        else:
            #Either we can't understand the output, or this is really
            #an error message not caught earlier (e.g. not in English)
            raise MissingExternalDependencyError(
                "Install EMBOSS if you want to use Bio.Emboss (%s)."
                % line)
    #In case there was no output at all...
    raise MissingExternalDependencyError("Could not get EMBOSS version")


#To avoid confusing known errors from old versions of EMBOSS ...
emboss_version = get_emboss_version()
if emboss_version < (6, 1, 0):
    raise MissingExternalDependencyError(
        "Test requires EMBOSS 6.1.0 patch 3 or later.")

#################################################################


#Top level function as this makes it easier to use for debugging:
def emboss_convert(filename, old_format, new_format):
    """Run seqret, returns handle."""
    #Setup, this assumes for all the format names used
    #Biopython and EMBOSS names are consistent!
    cline = SeqretCommandline(exes["seqret"],
                              sequence = filename,
                              sformat = old_format,
                              osformat = new_format,
                              auto = True,  # no prompting
                              stdout = True)
    #Run the tool,
    child = subprocess.Popen(str(cline),
                             stdin=subprocess.PIPE,
                             stdout=subprocess.PIPE,
                             stderr=subprocess.PIPE,
                             universal_newlines=True,
                             shell=(sys.platform!="win32"))
    child.stdin.close()
    child.stderr.close()
    return child.stdout


#Top level function as this makes it easier to use for debugging:
def emboss_piped_SeqIO_convert(records, old_format, new_format):
    """Run seqret, returns records (as a generator)."""
    #Setup, this assumes for all the format names used
    #Biopython and EMBOSS names are consistent!
    cline = SeqretCommandline(exes["seqret"],
                              sformat = old_format,
                              osformat = new_format,
                              auto = True,  # no prompting
                              filter = True)
    #Run the tool,
    child = subprocess.Popen(str(cline),
                             stdin=subprocess.PIPE,
                             stdout=subprocess.PIPE,
                             stderr=subprocess.PIPE,
                             universal_newlines=True,
                             shell=(sys.platform!="win32"))
    SeqIO.write(records, child.stdin, old_format)
    child.stdin.close()
    child.stderr.close()
    #TODO - Is there a nice way to return an iterator AND
    #automatically close the handle?
    records = list(SeqIO.parse(child.stdout, new_format))
    child.stdout.close()
    return records


#Top level function as this makes it easier to use for debugging:
def emboss_piped_AlignIO_convert(alignments, old_format, new_format):
    """Run seqret, returns alignments (as a generator)."""
    #Setup, this assumes for all the format names used
    #Biopython and EMBOSS names are consistent!
    cline = SeqretCommandline(exes["seqret"],
                              sformat = old_format,
                              osformat = new_format,
                              auto = True,  # no prompting
                              filter = True)
    #Run the tool,
    child = subprocess.Popen(str(cline),
                             stdin=subprocess.PIPE,
                             stdout=subprocess.PIPE,
                             stderr=subprocess.PIPE,
                             universal_newlines=True,
                             shell=(sys.platform!="win32"))
    try:
        AlignIO.write(alignments, child.stdin, old_format)
    except Exception as err:
        child.stdin.close()
        child.stderr.close()
        child.stdout.close()
        raise
    child.stdin.close()
    child.stderr.close()
    #TODO - Is there a nice way to return an iterator AND
    #automatically close the handle?
    try:
        aligns = list(AlignIO.parse(child.stdout, new_format))
    except Exception as err:
        child.stdout.close()
        raise
    child.stdout.close()
    return aligns


#Top level function as this makes it easier to use for debugging:
def compare_records(old_list, new_list):
    """Check two lists of SeqRecords agree, raises a ValueError if mismatch."""
    if len(old_list) != len(new_list):
        raise ValueError("%i vs %i records" % (len(old_list), len(new_list)))
    for old, new in zip(old_list, new_list):
        #Note the name matching is a bit fuzzy, e.g. truncation and
        #no spaces in PHYLIP files.
        if old.id != new.id and old.name != new.name \
        and (old.id not in new.id) and (new.id not in old.id) \
        and (old.id.replace(" ", "_") != new.id.replace(" ", "_")):
            raise ValueError("'%s' or '%s' vs '%s' or '%s' records"
                             % (old.id, old.name, new.id, new.name))
        if len(old.seq) != len(new.seq):
            raise ValueError("%i vs %i" % (len(old.seq), len(new.seq)))
        if str(old.seq).upper() != str(new.seq).upper():
            if str(old.seq).replace("X", "N")==str(new.seq) :
                raise ValueError("X -> N (protein forced into nucleotide?)")
            if len(old.seq) < 200:
                raise ValueError("'%s' vs '%s'" % (old.seq, new.seq))
            else:
                raise ValueError("'%s...%s' vs '%s...%s'"
                                 % (old.seq[:60], old.seq[-10:],
                                    new.seq[:60], new.seq[-10:]))
        if old.features and new.features \
        and len(old.features) != len(new.features):
            raise ValueError("%i vs %i features"
                             % (len(old.features, len(new.features))))
        #TODO - check annotation
    return True


#Top level function as this makes it easier to use for debugging:
def compare_alignments(old_list, new_list):
    """Check two lists of Alignments agree, raises a ValueError if mismatch."""
    if len(old_list) != len(new_list):
        raise ValueError("%i vs %i alignments" % (len(old_list), len(new_list)))
    for old, new in zip(old_list, new_list):
        if len(old) != len(new):
            raise ValueError("Alignment with %i vs %i records"
                             % (len(old), len(new)))
        compare_records(old, new)
    return True


class SeqRetSeqIOTests(unittest.TestCase):
    """Check EMBOSS seqret against Bio.SeqIO for converting files."""

    def tearDown(self):
        clean_up()

    def check_SeqIO_to_EMBOSS(self, in_filename, in_format, skip_formats=[],
                              alphabet=None):
        """Can Bio.SeqIO write files seqret can read back?"""
        if alphabet:
            records = list(SeqIO.parse(in_filename, in_format, alphabet))
        else:
            records = list(SeqIO.parse(in_filename, in_format))
        for temp_format in ["genbank", "embl", "fasta"]:
            if temp_format in skip_formats:
                continue
            new_records = list(emboss_piped_SeqIO_convert(records, temp_format, "fasta"))
            try:
                self.assertTrue(compare_records(records, new_records))
            except ValueError as err:
                raise ValueError("Disagree on file %s %s in %s format: %s"
                                 % (in_format, in_filename, temp_format, err))

    def check_EMBOSS_to_SeqIO(self, filename, old_format,
                              skip_formats=[]):
        """Can Bio.SeqIO read seqret's conversion of the file?"""
        #TODO: Why can't we read EMBOSS's swiss output?
        self.assertTrue(os.path.isfile(filename))
        old_records = list(SeqIO.parse(filename, old_format))
        for new_format in ["genbank", "fasta", "pir", "embl", "ig"]:
            if new_format in skip_formats:
                continue
            handle = emboss_convert(filename, old_format, new_format)
            new_records = list(SeqIO.parse(handle, new_format))
            handle.close()
            try:
                self.assertTrue(compare_records(old_records, new_records))
            except ValueError as err:
                raise ValueError("Disagree on %s file %s in %s format: %s"
                                 % (old_format, filename, new_format, err))

    def check_SeqIO_with_EMBOSS(self, filename, old_format, skip_formats=[],
                                alphabet=None):
        #Check EMBOSS can read Bio.SeqIO output...
        self.check_SeqIO_to_EMBOSS(filename, old_format, skip_formats,
                                   alphabet)
        #Check Bio.SeqIO can read EMBOSS seqret output...
        self.check_EMBOSS_to_SeqIO(filename, old_format, skip_formats)

    def test_abi(self):
        """SeqIO agrees with EMBOSS' Abi to FASTQ conversion."""
        #This lets use check the id, sequence, and quality scores
        for filename in ["Abi/3730.ab1", "Abi/empty.ab1"]:
            old = SeqIO.read(filename, "abi")
            handle = emboss_convert(filename, "abi", "fastq-sanger")
            new = SeqIO.read(handle, "fastq-sanger")
            handle.close()
            if emboss_version == (6, 4, 0) and new.id == "EMBOSS_001":
                #Avoid bug in EMBOSS 6.4.0 (patch forthcoming)
                pass
            else:
                self.assertEqual(old.id, new.id)
            self.assertEqual(str(old.seq), str(new.seq))
            if emboss_version < (6, 3, 0) and new.letter_annotations["phred_quality"] == [1]*len(old):
                #Apparent bug in EMBOSS 6.2.0.1 on Windows
                pass
            else:
                self.assertEqual(old.letter_annotations, new.letter_annotations)

    def test_genbank(self):
        """SeqIO & EMBOSS reading each other's conversions of a GenBank file."""
        self.check_SeqIO_with_EMBOSS("GenBank/cor6_6.gb", "genbank")

    def test_genbank2(self):
        """SeqIO & EMBOSS reading each other's conversions of another GenBank file."""
        self.check_SeqIO_with_EMBOSS("GenBank/NC_000932.gb", "genbank")

    def test_embl(self):
        """SeqIO & EMBOSS reading each other's conversions of an EMBL file."""
        self.check_SeqIO_with_EMBOSS("EMBL/U87107.embl", "embl")

    def test_ig(self):
        """SeqIO & EMBOSS reading each other's conversions of an ig file."""
        #NOTE - EMBOSS considers "genbank" to be for nucleotides only,
        #and will turn "X" into "N" for GenBank output.
        self.check_SeqIO_to_EMBOSS("IntelliGenetics/VIF_mase-pro.txt", "ig",
                                   alphabet=generic_protein,
                                   skip_formats=["genbank", "embl"])
        #TODO - What does a % in an ig sequence mean?
        #e.g. "IntelliGenetics/vpu_nucaligned.txt"
        #and  "IntelliGenetics/TAT_mase_nuc.txt"
        #EMBOSS seems to ignore them.

    def test_pir(self):
        """SeqIO & EMBOSS reading each other's conversions of a PIR file."""
        #Skip genbank here, EMBOSS mangles the LOCUS line:
        self.check_SeqIO_with_EMBOSS("NBRF/clustalw.pir", "pir",
                               skip_formats=["genbank"])
        #Skip EMBL here, EMBOSS mangles the ID line
        #Skip GenBank, EMBOSS 6.0.1 on Windows won't output proteins as GenBank
        self.check_SeqIO_with_EMBOSS("NBRF/DMB_prot.pir", "pir",
                               skip_formats=["embl", "genbank"])

    def test_clustalw(self):
        """SeqIO & EMBOSS reading each other's conversions of a Clustalw file."""
        self.check_SeqIO_with_EMBOSS("Clustalw/hedgehog.aln", "clustal",
                                   skip_formats=["embl", "genbank"])
        self.check_SeqIO_with_EMBOSS("Clustalw/opuntia.aln", "clustal",
                                   skip_formats=["embl", "genbank"])


class SeqRetAlignIOTests(unittest.TestCase):
    """Check EMBOSS seqret against Bio.SeqIO for converting files."""

    def tearDown(self):
        clean_up()

    def check_EMBOSS_to_AlignIO(self, filename, old_format,
                              skip_formats=[]):
        """Can AlignIO read seqret's conversion of the file?"""
        self.assertTrue(os.path.isfile(filename), filename)
        old_aligns = list(AlignIO.parse(filename, old_format))
        formats = ["clustal", "phylip", "ig"]
        if len(old_aligns) == 1:
            formats.extend(["fasta", "nexus"])
        for new_format in formats:
            if new_format in skip_formats:
                continue
            handle = emboss_convert(filename, old_format, new_format)
            try:
                new_aligns = list(AlignIO.parse(handle, new_format))
            except:
                handle.close()
                raise ValueError("Can't parse %s file %s in %s format."
                                 % (old_format, filename, new_format))
            handle.close()
            try:
                self.assertTrue(compare_alignments(old_aligns, new_aligns))
            except ValueError as err:
                raise ValueError("Disagree on %s file %s in %s format: %s"
                                 % (old_format, filename, new_format, err))

    def check_AlignIO_to_EMBOSS(self, in_filename, in_format, skip_formats=[],
                                alphabet=None):
        """Can Bio.AlignIO write files seqret can read back?"""
        if alphabet:
            old_aligns = list(AlignIO.parse(in_filename, in_format, alphabet))
        else:
            old_aligns = list(AlignIO.parse(in_filename, in_format))

        formats = ["clustal", "phylip"]
        if len(old_aligns) == 1:
            formats.extend(["fasta", "nexus"])
        for temp_format in formats:
            if temp_format in skip_formats:
                continue
            #PHYLIP is a simple format which explicitly supports
            #multiple alignments (unlike FASTA).
            try:
                new_aligns = list(emboss_piped_AlignIO_convert(old_aligns,
                                                               temp_format,
                                                               "phylip"))
            except ValueError as e:
                #e.g. ValueError: Need a DNA, RNA or Protein alphabet
                #from writing Nexus files...
                continue
            try:
                self.assertTrue(compare_alignments(old_aligns, new_aligns))
            except ValueError as err:
                raise ValueError("Disagree on file %s %s in %s format: %s"
                                 % (in_format, in_filename, temp_format, err))

    def check_AlignIO_with_EMBOSS(self, filename, old_format, skip_formats=[],
                                  alphabet=None):
        #Check EMBOSS can read Bio.AlignIO output...
        self.check_AlignIO_to_EMBOSS(filename, old_format, skip_formats,
                                   alphabet)
        #Check Bio.AlignIO can read EMBOSS seqret output...
        self.check_EMBOSS_to_AlignIO(filename, old_format, skip_formats)

    def test_align_clustalw(self):
        """AlignIO & EMBOSS reading each other's conversions of a ClustalW file."""
        self.check_AlignIO_with_EMBOSS("Clustalw/hedgehog.aln", "clustal")
        self.check_AlignIO_with_EMBOSS("Clustalw/opuntia.aln", "clustal")
        self.check_AlignIO_with_EMBOSS("Clustalw/odd_consensus.aln", "clustal",
                               skip_formats=["nexus"])  # TODO - why not nexus?
        self.check_AlignIO_with_EMBOSS("Clustalw/protein.aln", "clustal")
        self.check_AlignIO_with_EMBOSS("Clustalw/promals3d.aln", "clustal")

    def test_clustalw(self):
        """AlignIO & EMBOSS reading each other's conversions of a PHYLIP file."""
        self.check_AlignIO_with_EMBOSS("Phylip/horses.phy", "phylip")
        self.check_AlignIO_with_EMBOSS("Phylip/hennigian.phy", "phylip")
        self.check_AlignIO_with_EMBOSS("Phylip/reference_dna.phy", "phylip")
        self.check_AlignIO_with_EMBOSS("Phylip/reference_dna2.phy", "phylip")
        self.check_AlignIO_with_EMBOSS("Phylip/interlaced.phy", "phylip")
        self.check_AlignIO_with_EMBOSS("Phylip/interlaced2.phy", "phylip")
        self.check_AlignIO_with_EMBOSS("Phylip/random.phy", "phylip")


class PairwiseAlignmentTests(unittest.TestCase):
    """Run pairwise alignments with water and needle, and parse them."""

    def tearDown(self):
        clean_up()

    def pairwise_alignment_check(self, query_seq,
                                 targets, alignments,
                                 local=True):
        """Check pairwise alignment data is sane."""
        #The datasets should be small, so making iterators into lists is OK
        targets = list(targets)
        alignments = list(alignments)
        self.assertEqual(len(targets), len(alignments))
        for target, alignment in zip(targets, alignments):
            self.assertEqual(len(alignment), 2)
            #self.assertEqual(target.id, alignment[1].id) #too strict
            if alignment[1].id not in target.id \
            and alignment[1].id not in target.name:
                raise AssertionError("%s vs %s or %s"
                                     % (alignment[1].id, target.id, target.name))
            if local:
                #Local alignment
                self.assertTrue(str(alignment[0].seq).replace("-", "")
                             in query_seq)
                self.assertTrue(str(alignment[1].seq).replace("-", "").upper()
                             in str(target.seq).upper())
            else:
                #Global alignment
                self.assertEqual(str(query_seq), str(alignment[0].seq).replace("-", ""))
                self.assertEqual(str(target.seq).upper(),
                                 str(alignment[1].seq).replace("-", "").upper())
        return True

    def run_water(self, cline):
        #Run the tool,
        stdout, stderr = cline()
        self.assertTrue(stderr.strip().startswith("Smith-Waterman local alignment"),
                        stderr)
        if cline.outfile:
            self.assertEqual(stdout.strip(), "")
            self.assertTrue(os.path.isfile(cline.outfile),
                            "Missing output file %r from:\n%s" % (cline.outfile, cline))
        else :
            #Don't use this yet... could return stdout handle instead?
            return stdout

    def test_water_file(self):
        """water with the asis trick, output to a file."""
        #Setup, try a mixture of keyword arguments and later additions:
        cline = WaterCommandline(cmd=exes["water"],
                                 gapopen="10", gapextend="0.5")
        #Try using both human readable names, and the literal ones:
        cline.set_parameter("asequence", "asis:ACCCGGGCGCGGT")
        cline.set_parameter("-bsequence", "asis:ACCCGAGCGCGGT")
        #Try using a property set here:
        cline.outfile = "Emboss/temp with space.water"
        self.assertEqual(str(eval(repr(cline))), str(cline))
        #Run the tool,
        self.run_water(cline)
        #Check we can parse the output...
        align = AlignIO.read(cline.outfile, "emboss")
        self.assertEqual(len(align), 2)
        self.assertEqual(str(align[0].seq), "ACCCGGGCGCGGT")
        self.assertEqual(str(align[1].seq), "ACCCGAGCGCGGT")
        #Clean up,
        os.remove(cline.outfile)

    def test_water_piped(self):
        """water with asis trick, output piped to stdout."""
        cline = WaterCommandline(cmd=exes["water"],
                                 asequence="asis:ACCCGGGCGCGGT",
                                 bsequence="asis:ACCCGAGCGCGGT",
                                 gapopen=10,
                                 gapextend=0.5,
                                 auto=True, filter=True)
        self.assertEqual(str(cline),
                         exes["water"] + " -auto -filter"
                         + " -asequence=asis:ACCCGGGCGCGGT"
                         + " -bsequence=asis:ACCCGAGCGCGGT"
                         + " -gapopen=10 -gapextend=0.5")
        #Run the tool,
        child = subprocess.Popen(str(cline),
                                 stdin=subprocess.PIPE,
                                 stdout=subprocess.PIPE,
                                 stderr=subprocess.PIPE,
                                 universal_newlines=True,
                                 shell=(sys.platform!="win32"))
        child.stdin.close()
        #Check we could read it's output
        align = AlignIO.read(child.stdout, "emboss")
        self.assertEqual(len(align), 2)
        self.assertEqual(str(align[0].seq), "ACCCGGGCGCGGT")
        self.assertEqual(str(align[1].seq), "ACCCGAGCGCGGT")
        #Check no error output:
        self.assertEqual(child.stderr.read(), "")
        self.assertEqual(0, child.wait())
        child.stdout.close()
        child.stderr.close()

    def test_needle_file(self):
        """needle with the asis trick, output to a file."""
        #Setup,
        cline = NeedleCommandline(cmd=exes["needle"])
        cline.set_parameter("-asequence", "asis:ACCCGGGCGCGGT")
        cline.set_parameter("-bsequence", "asis:ACCCGAGCGCGGT")
        cline.set_parameter("-gapopen", "10")
        cline.set_parameter("-gapextend", "0.5")
        #EMBOSS would guess this, but let's be explicit:
        cline.set_parameter("-snucleotide", "True")
        cline.set_parameter("-outfile", "Emboss/temp with space.needle")
        self.assertEqual(str(eval(repr(cline))), str(cline))
        #Run the tool,
        stdout, stderr = cline()
        #Check it worked,
        self.assertTrue(stderr.strip().startswith("Needleman-Wunsch global alignment"), stderr)
        self.assertEqual(stdout.strip(), "")
        filename = cline.outfile
        self.assertTrue(os.path.isfile(filename),
                        "Missing output file %r from:\n%s" % (filename, cline))
        #Check we can parse the output...
        align = AlignIO.read(filename, "emboss")
        self.assertEqual(len(align), 2)
        self.assertEqual(str(align[0].seq), "ACCCGGGCGCGGT")
        self.assertEqual(str(align[1].seq), "ACCCGAGCGCGGT")
        #Clean up,
        os.remove(filename)

    def test_needle_piped(self):
        """needle with asis trick, output piped to stdout."""
        cline = NeedleCommandline(cmd=exes["needle"],
                                 asequence="asis:ACCCGGGCGCGGT",
                                 bsequence="asis:ACCCGAGCGCGGT",
                                 gapopen=10,
                                 gapextend=0.5,
                                 auto=True, filter=True)
        self.assertEqual(str(cline),
                         exes["needle"] + " -auto -filter"
                         + " -asequence=asis:ACCCGGGCGCGGT"
                         + " -bsequence=asis:ACCCGAGCGCGGT"
                         + " -gapopen=10 -gapextend=0.5")
        #Run the tool,
        child = subprocess.Popen(str(cline),
                                 stdin=subprocess.PIPE,
                                 stdout=subprocess.PIPE,
                                 stderr=subprocess.PIPE,
                                 universal_newlines=True,
                                 shell=(sys.platform!="win32"))
        child.stdin.close()
        #Check we could read it's output
        align = AlignIO.read(child.stdout, "emboss")
        self.assertEqual(len(align), 2)
        self.assertEqual(str(align[0].seq), "ACCCGGGCGCGGT")
        self.assertEqual(str(align[1].seq), "ACCCGAGCGCGGT")
        #Check no error output:
        self.assertEqual(child.stderr.read(), "")
        self.assertEqual(0, child.wait())
        child.stdout.close()
        child.stderr.close()

    def test_water_file2(self):
        """water with the asis trick and nucleotide FASTA file, output to a file."""
        #Setup,
        query = "ACACACTCACACACACTTGGTCAGAGATGCTGTGCTTCTTGGAAGCAAGGNCTCAAAGGCAAGGTGCACGCAGAGGGACGTTTGAGTCTGGGATGAAGCATGTNCGTATTATTTATATGATGGAATTTCACGTTTTTATG"
        out_file = "Emboss/temp_test2.water"
        in_file = "Fasta/f002"
        self.assertTrue(os.path.isfile(in_file))
        if os.path.isfile(out_file):
            os.remove(out_file)
        cline = WaterCommandline(cmd=exes["water"])
        cline.set_parameter("-asequence", "asis:%s" % query)
        cline.set_parameter("-bsequence", in_file)
        cline.set_parameter("-gapopen", "10")
        cline.set_parameter("-gapextend", "0.5")
        cline.set_parameter("-outfile", out_file)
        self.assertEqual(str(eval(repr(cline))), str(cline))
        #Run the tool,
        self.run_water(cline)
        #Check we can parse the output and it is sensible...
        self.pairwise_alignment_check(query,
                                      SeqIO.parse(in_file, "fasta"),
                                      AlignIO.parse(out_file, "emboss"),
                                      local=True)
        #Clean up,
        os.remove(out_file)

    def test_water_file3(self):
        """water with the asis trick and GenBank file, output to a file."""
        #Setup,
        query = "TGTTGTAATGTTTTAATGTTTCTTCTCCCTTTAGATGTACTACGTTTGGA"
        out_file = "Emboss/temp_test3.water"
        in_file = "GenBank/cor6_6.gb"
        self.assertTrue(os.path.isfile(in_file))
        if os.path.isfile(out_file):
            os.remove(out_file)
        cline = WaterCommandline(cmd=exes["water"])
        cline.set_parameter("asequence", "asis:%s" % query)
        cline.set_parameter("bsequence", in_file)
        #TODO - Tell water this is a GenBank file!
        cline.set_parameter("gapopen", "1")
        cline.set_parameter("gapextend", "0.5")
        cline.set_parameter("outfile", out_file)
        self.assertEqual(str(eval(repr(cline))), str(cline))
        #Run the tool,
        self.run_water(cline)
        #Check we can parse the output and it is sensible...
        self.pairwise_alignment_check(query,
                                      SeqIO.parse(in_file, "genbank"),
                                      AlignIO.parse(out_file, "emboss"),
                                      local=True)
        #Clean up,
        os.remove(out_file)

    def test_water_file4(self):
        """water with the asis trick and SwissProt file, output to a file."""
        #Setup,
        query = "DVCTGKALCDPVTQNIKTYPVKIENLRVMI"
        out_file = "Emboss/temp_test4.water"
        in_file = "SwissProt/sp004"
        self.assertTrue(os.path.isfile(in_file))
        if os.path.isfile(out_file):
            os.remove(out_file)
        cline = WaterCommandline(cmd=exes["water"])
        cline.set_parameter("-asequence", "asis:%s" % query)
        cline.set_parameter("-bsequence", in_file)
        #EMBOSS should work this out, but let's be explicit:
        cline.set_parameter("-sprotein", True)
        #TODO - Tell water this is a SwissProt file!
        cline.set_parameter("-gapopen", "20")
        cline.set_parameter("-gapextend", "5")
        cline.set_parameter("-outfile", out_file)
        self.assertEqual(str(eval(repr(cline))), str(cline))
        #Run the tool,
        self.run_water(cline)
        #Check we can parse the output and it is sensible...
        self.pairwise_alignment_check(query,
                                      SeqIO.parse(in_file, "swiss"),
                                      AlignIO.parse(out_file, "emboss"),
                                      local=True)
        #Clean up,
        os.remove(out_file)

    def test_needle_piped2(self):
        """needle with asis trick, and nucleotide FASTA file, output piped to stdout."""
        #TODO - Support needle in Bio.Emboss.Applications
        #(ideally with the -auto and -filter arguments)
        #Setup,
        query = "ACACACTCACACACACTTGGTCAGAGATGCTGTGCTTCTTGGAA"
        cline = exes["needle"]
        cline += " -asequence asis:" + query
        cline += " -bsequence Fasta/f002"
        cline += " -auto"  # no prompting
        cline += " -filter"  # use stdout
        #Run the tool,
        child = subprocess.Popen(str(cline),
                                 stdin=subprocess.PIPE,
                                 stdout=subprocess.PIPE,
                                 stderr=subprocess.PIPE,
                                 universal_newlines=True,
                                 shell=(sys.platform!="win32"))
        child.stdin.close()
        #Check we can parse the output and it is sensible...
        self.pairwise_alignment_check(query,
                                      SeqIO.parse("Fasta/f002", "fasta"),
                                      AlignIO.parse(child.stdout, "emboss"),
                                      local=False)
        #Check no error output:
        self.assertEqual(child.stderr.read(), "")
        self.assertEqual(0, child.wait())
        child.stdout.close()
        child.stderr.close()

    def test_water_needs_output(self):
        """water without output file or stdout/filter should give error."""
        cline = WaterCommandline(cmd=exes["water"],
                                 asequence="asis:ACCCGGGCGCGGT",
                                 bsequence="asis:ACCCGAGCGCGGT",
                                 gapopen=10,
                                 gapextend=0.5,
                                 auto=True)
        self.assertTrue(cline.auto)
        self.assertTrue(not cline.stdout)
        self.assertTrue(not cline.filter)
        self.assertEqual(cline.outfile, None)
        self.assertRaises(ValueError, str, cline)

    def test_needle_needs_output(self):
        """needle without output file or stdout/filter should give error."""
        cline = NeedleCommandline(cmd=exes["needle"],
                                 asequence="asis:ACCCGGGCGCGGT",
                                 bsequence="asis:ACCCGAGCGCGGT",
                                 gapopen=10,
                                 gapextend=0.5,
                                 auto=True)
        self.assertTrue(cline.auto)
        self.assertTrue(not cline.stdout)
        self.assertTrue(not cline.filter)
        self.assertEqual(cline.outfile, None)
        self.assertRaises(ValueError, str, cline)

    def test_seqtmatchall_piped(self):
        """seqmatchall with pair output piped to stdout."""
        cline = SeqmatchallCommandline(cmd=exes["seqmatchall"],
                                       sequence="Fasta/f002",
                                       aformat="pair", wordsize=9,
                                       auto=True, stdout=True)
        self.assertEqual(str(cline),
                         exes["seqmatchall"] + " -auto -stdout"
                         + " -sequence=Fasta/f002"
                         + " -wordsize=9 -aformat=pair")
        #Run the tool,
        child = subprocess.Popen(str(cline),
                                 stdin=subprocess.PIPE,
                                 stdout=subprocess.PIPE,
                                 stderr=subprocess.PIPE,
                                 universal_newlines=True,
                                 shell=(sys.platform!="win32"))
        child.stdin.close()
        #Check we could read it's output
        for align in AlignIO.parse(child.stdout, "emboss") :
            self.assertEqual(len(align), 2)
            self.assertEqual(align.get_alignment_length(), 9)
        #Check no error output:
        self.assertEqual(child.stderr.read(), "")
        self.assertEqual(0, child.wait())
        child.stdout.close()
        child.stderr.close()


#Top level function as this makes it easier to use for debugging:
def emboss_translate(sequence, table=None, frame=None):
    """Call transeq, returns protein sequence as string."""
    #TODO - Support transeq in Bio.Emboss.Applications?
    #(doesn't seem worthwhile as Biopython can do translations)

    if not sequence:
        raise ValueError(sequence)

    #Setup,
    cline = exes["transeq"]

    if len(sequence) < 100:
        filename = None
        cline += " -sequence asis:%s" % sequence
    else:
        #There are limits on command line string lengths...
        #use a temp file instead.
        filename = "Emboss/temp_transeq.txt"
        SeqIO.write(SeqRecord(sequence, id="Test"), filename, "fasta")
        cline += " -sequence %s" % filename

    cline += " -auto"  # no prompting
    cline += " -filter"  # use stdout
    if table is not None:
        cline += " -table %s" % str(table)
    if frame is not None:
        cline += " -frame %s" % str(frame)
    #Run the tool,
    child = subprocess.Popen(str(cline),
                             stdin=subprocess.PIPE,
                             stdout=subprocess.PIPE,
                             stderr=subprocess.PIPE,
                             universal_newlines=True,
                             shell=(sys.platform!="win32"))
    out, err = child.communicate()
    #Check no error output:
    if err != "":
        raise ValueError(str(cline) + "\n" + err)

    #Check we could read it's output
    record = SeqIO.read(StringIO(out), "fasta")

    if 0 != child.wait():
        raise ValueError(str(cline))

    if filename:
        os.remove(filename)
        if not record.id.startswith("Test"):
            raise ValueError(str(cline))
    else:
        if not record.id.startswith("asis"):
            raise ValueError(str(cline))
    return str(record.seq)


#Top level function as this makes it easier to use for debugging:
def check_translation(sequence, translation, table=None):
    if table is None:
        t = 1
    else:
        t = table
    if translation != str(sequence.translate(t)) \
    or translation != str(translate(sequence, t)) \
    or translation != translate(str(sequence), t):
        #More details...
        for i, amino in enumerate(translation):
            codon = sequence[i*3:i*3+3]
            if amino != str(codon.translate(t)):
                raise ValueError("%s -> %s not %s (table %s)"
                         % (codon, amino, codon.translate(t), t))
        #Shouldn't reach this line:
        raise ValueError("%s -> %s (table %s)"
                         % (sequence, translation, t))
    return True


class TranslationTests(unittest.TestCase):
    """Run pairwise alignments with water and needle, and parse them."""

    def tearDown(self):
        clean_up()

    def test_simple(self):
        """transeq vs Bio.Seq for simple translations (including alt tables)."""

        examples = [Seq("ACGTGACTGACGTAGCATGCCACTAGG"),
                    #Unamibguous TA? codons:
                    Seq("TAATACTATTAG", generic_dna),
                    #Most of the ambiguous TA? codons:
                    Seq("TANTARTAYTAMTAKTAHTABTADTAV", generic_dna),
                    #Problem cases,
                    #
                    #Seq("TAW", generic_dna),
                    #W = A or T, but EMBOSS does TAW -> X
                    #TAA -> Y, TAT ->Y, so in Biopython TAW -> Y
                    #
                    #Seq("TAS", generic_dna),
                    #S = C or G, but EMBOSS does TAS -> Y
                    #TAG -> *, TAC ->Y, so in Biopython TAS -> X (Y or *)
                    #
                    #Seq("AAS", generic_dna),
                    #On table 9, EMBOSS gives N, we give X.
                    #S = C or G, so according to my reading of
                    #table 9 on the NCBI page, AAC=N, AAG=K
                    #suggesting this is a bug in EMBOSS.
                    #
                    Seq("ACGGGGGGGGTAAGTGGTGTGTGTGTAGT", generic_dna),
                    ]

        for sequence in examples:
            #EMBOSS treats spare residues differently... avoid this issue
            if len(sequence) % 3 != 0:
                sequence = sequence[:-(len(sequence)%3)]
            self.assertEqual(len(sequence) % 3, 0)
            self.assertTrue(len(sequence) > 0)
            self.check(sequence)

    def check(self, sequence):
        """Compare our translation to EMBOSS's using all tables.

        Takes a Seq object (and a filename containing it)."""
        translation = emboss_translate(sequence)
        self.assertTrue(check_translation(sequence, translation))

        for table in [1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16, 21, 22, 23]:
            translation = emboss_translate(sequence, table)
            self.assertTrue(check_translation(sequence, translation, table))
        return True

    def translate_all_codons(self, letters):
        sequence = Seq("".join(c1+c3+c3
                               for c1 in letters
                               for c2 in letters
                               for c3 in letters),
                       generic_nucleotide)
        self.check(sequence)

    #def test_all_ambig_dna_codons(self):
    #    """transeq vs Bio.Seq on ambiguous DNA codons (inc. alt tables)."""
    #    self.translate_all_codons(ambiguous_dna_letters)

    def test_all_unambig_dna_codons(self):
        """transeq vs Bio.Seq on unambiguous DNA codons (inc. alt tables)."""
        self.translate_all_codons("ATCGatcg")

    def test_all_unambig_rna_codons(self):
        """transeq vs Bio.Seq on unambiguous RNA codons (inc. alt tables)."""
        self.translate_all_codons("AUCGaucg")

    def test_mixed_unambig_rna_codons(self):
        """transeq vs Bio.Seq on unambiguous DNA/RNA codons (inc. alt tables)."""
        self.translate_all_codons("ATUCGatucg")


def clean_up():
    """Fallback clean up method to remove temp files."""
    for filename in os.listdir("Emboss"):
        if filename.startswith("temp_"):
            try:
                os.remove(filename)
            except:
                pass

if __name__ == "__main__":
    runner = unittest.TextTestRunner(verbosity = 2)
    unittest.main(testRunner=runner)
    clean_up()