def test_to_phylip_no_positions(self):
        d1 = DNASequence('', id="d1")
        d2 = DNASequence('', id="d2")
        a = Alignment([d1, d2])

        with self.assertRaises(SequenceCollectionError):
            a.to_phylip()

    def test_update_ids_sequence_attributes_propagated(self):
        # 1 seq
        exp_sc = Alignment([
            DNA('ACGT', id="abc", description='desc', quality=range(4))
        ])
        exp_id_map = {'abc': 'seq1'}

        obj = Alignment([
            DNA('ACGT', id="seq1", description='desc', quality=range(4))
        ])

        obs_sc, obs_id_map = obj.update_ids(ids=('abc',))
        self._assert_sequence_collections_equal(obs_sc, exp_sc)
        self.assertEqual(obs_id_map, exp_id_map)

        # 2 seqs
        exp_sc = Alignment([
            DNA('ACGT', id="abc", description='desc1', quality=range(4)),
            DNA('TGCA', id="def", description='desc2', quality=range(4)[::-1])
        ])
        exp_id_map = {'abc': 'seq1', 'def': 'seq2'}

        obj = Alignment([
            DNA('ACGT', id="seq1", description='desc1', quality=(0, 1, 2, 3)),
            DNA('TGCA', id="seq2", description='desc2', quality=(3, 2, 1, 0))
        ])

        obs_sc, obs_id_map = obj.update_ids(ids=('abc', 'def'))
        self._assert_sequence_collections_equal(obs_sc, exp_sc)
        self.assertEqual(obs_id_map, exp_id_map)

    def test_omit_gap_positions(self):
        expected = self.a2
        self.assertEqual(self.a2.omit_gap_positions(1.0), expected)
        self.assertEqual(self.a2.omit_gap_positions(0.51), expected)

        r1 = RNA('UUAU', metadata={'id': "r1"})
        r2 = RNA('ACGU', metadata={'id': "r2"})
        expected = Alignment([r1, r2])
        self.assertEqual(self.a2.omit_gap_positions(0.49), expected)

        r1 = RNA('UUAU', metadata={'id': "r1"})
        r2 = RNA('ACGU', metadata={'id': "r2"})
        expected = Alignment([r1, r2])
        self.assertEqual(self.a2.omit_gap_positions(0.0), expected)

        self.assertEqual(self.empty.omit_gap_positions(0.0), self.empty)
        self.assertEqual(self.empty.omit_gap_positions(0.49), self.empty)
        self.assertEqual(self.empty.omit_gap_positions(1.0), self.empty)

        # Test to ensure floating point precision bug isn't present. See the
        # tests for Alignment.position_frequencies for more details.
        seqs = []
        for i in range(33):
            seqs.append(DNA('-.', metadata={'id': str(i)}))
        aln = Alignment(seqs)
        self.assertEqual(aln.omit_gap_positions(1 - np.finfo(float).eps),
                         Alignment([DNA('', metadata={'id': str(i)})
                                    for i in range(33)]))

def compute_distance_matrix(msa_file, csvfile="distance_mat.csv"):
    """
    load up some aligned sequences, and compute a distance matrix 
    compute distances between the sequences using the hamming function

    see also: 
    scipy.spatial.distance.hamming

    @args msa_file: multiple sequence alignment in fasta format 
    @type msa_file: str 
    @args csvfile: output distance matrix file in csv format 
    @type csvfile: str 
    """

    records = []
    for rec in SeqIO.parse(msa_file, "fasta"):
        records.append(RNA(rec.seq, rec.id))

    aln = Alignment(records)
    master_dm = aln.distances()

    ## writing the result to a csv file
    csv_header_row = [header for header in master_dm.ids]

    ## result as a list of list
    with open(csvfile, "w") as output:
        writer = csv.writer(output, lineterminator="\n")

        writer.writerows([csv_header_row])
        writer.writerows(master_dm)

    output.close()

    def test_to_phylip_unequal_sequence_lengths(self):
        d1 = DNASequence('A-CT', id="d1")
        d2 = DNASequence('TTA', id="d2")
        d3 = DNASequence('.-AC', id="d3")
        a = Alignment([d1, d2, d3])

        with self.assertRaises(SequenceCollectionError):
            a.to_phylip()

    def test_to_phylip_map_labels(self):
        """to_phylip functions as expected with label mapping
        """
        d1 = DNASequence("..ACC-GTTGG..", id="d1")
        d2 = DNASequence("TTACCGGT-GGCC", id="d2")
        d3 = DNASequence(".-ACC-GTTGC--", id="d3")
        a = Alignment([d1, d2, d3])

        phylip_str, id_map = a.to_phylip(map_labels=True, label_prefix="s")
        self.assertEqual(id_map, {"s1": "d1", "s3": "d3", "s2": "d2"})
        expected = "\n".join(["3 13", "s1 ..ACC-GTTGG..", "s2 TTACCGGT-GGCC", "s3 .-ACC-GTTGC--"])
        self.assertEqual(phylip_str, expected)

    def test_to_phylip(self):
        """to_phylip functions as expected
        """
        d1 = DNASequence("..ACC-GTTGG..", id="d1")
        d2 = DNASequence("TTACCGGT-GGCC", id="d2")
        d3 = DNASequence(".-ACC-GTTGC--", id="d3")
        a = Alignment([d1, d2, d3])

        phylip_str, id_map = a.to_phylip(map_labels=False)
        self.assertEqual(id_map, {"d1": "d1", "d3": "d3", "d2": "d2"})
        expected = "\n".join(["3 13", "d1 ..ACC-GTTGG..", "d2 TTACCGGT-GGCC", "d3 .-ACC-GTTGC--"])
        self.assertEqual(phylip_str, expected)

    def test_majority_consensus(self):
        d1 = DNASequence('TTT', id="d1")
        d2 = DNASequence('TT-', id="d2")
        d3 = DNASequence('TC-', id="d3")
        a1 = Alignment([d1, d2, d3])
        self.assertTrue(a1.majority_consensus().equals(DNASequence('TT-')))

        d1 = DNASequence('T', id="d1")
        d2 = DNASequence('A', id="d2")
        a1 = Alignment([d1, d2])
        self.assertTrue(a1.majority_consensus() in
                        [DNASequence('T'), DNASequence('A')])

        self.assertEqual(self.empty.majority_consensus(), '')

    def test_majority_consensus(self):
        """majority_consensus functions as expected
        """
        d1 = DNASequence("TTT", id="d1")
        d2 = DNASequence("TT-", id="d2")
        d3 = DNASequence("TC-", id="d3")
        a1 = Alignment([d1, d2, d3])
        self.assertEqual(a1.majority_consensus(), DNASequence("TT-"))

        d1 = DNASequence("T", id="d1")
        d2 = DNASequence("A", id="d2")
        a1 = Alignment([d1, d2])
        self.assertTrue(a1.majority_consensus() in [DNASequence("T"), DNASequence("A")])

        self.assertEqual(self.empty.majority_consensus(), "")

    def test_generate_lane_mask(self):

        sample_alignment = """>1
        AAAAT
        >2
        AAAGG
        >3
        AACCC
        >4
        A----""".split('\n')
        aln = Alignment.from_fasta_records(parse_fasta(sample_alignment), DNA)

        actual_lanemask = generate_lane_mask(aln, 0.00)
        self.assertEqual(actual_lanemask, "11111")
        actual_lanemask = generate_lane_mask(aln, 0.10)
        self.assertEqual(actual_lanemask, "11100")
        actual_lanemask = generate_lane_mask(aln, 0.20)
        self.assertEqual(actual_lanemask, "11100")
        actual_lanemask = generate_lane_mask(aln, 0.40)
        self.assertEqual(actual_lanemask, "11000")
        actual_lanemask = generate_lane_mask(aln, 0.60)
        self.assertEqual(actual_lanemask, "11000")
        actual_lanemask = generate_lane_mask(aln, 0.80)
        self.assertEqual(actual_lanemask, "10000")
        actual_lanemask = generate_lane_mask(aln, 1.00)
        self.assertEqual(actual_lanemask, "00000")

def align_two_alignments(aln1_fp, aln2_fp, moltype, params=None):
    """Returns an Alignment object from two existing Alignments.

    Parameters
    ----------
    aln1_fp : string
        file path of 1st alignment
    aln2_fp : string
        file path of 2nd alignment
    params : dict of parameters to pass in to the Mafft app controller.

    Returns
    -------
        The aligned sequences.
    """

    # Create Mafft app.
    app = Mafft(InputHandler='_input_as_paths',
                params=params,
                SuppressStderr=False)
    app._command = 'mafft-profile'

    # Get results using int_map as input to app
    res = app([aln1_fp, aln2_fp])

    return Alignment.read(res['StdOut'], constructor=moltype)

    def test_omit_gap_sequences(self):
        expected = self.a2
        self.assertEqual(self.a2.omit_gap_sequences(1.0), expected)
        self.assertEqual(self.a2.omit_gap_sequences(0.20), expected)

        expected = Alignment([self.r2])
        self.assertEqual(self.a2.omit_gap_sequences(0.19), expected)

        self.assertEqual(self.empty.omit_gap_sequences(0.0), self.empty)
        self.assertEqual(self.empty.omit_gap_sequences(0.2), self.empty)
        self.assertEqual(self.empty.omit_gap_sequences(1.0), self.empty)

        # Test to ensure floating point precision bug isn't present. See the
        # tests for Alignment.position_frequencies for more details.
        aln = Alignment([DNA('.' * 33, id='abc'), DNA('-' * 33, id='def')])
        self.assertEqual(aln.omit_gap_sequences(1 - np.finfo(float).eps),
                         Alignment([]))

    def setUp(self):
        self.d1 = DNASequence("..ACC-GTTGG..", id="d1")
        self.d2 = DNASequence("TTACCGGT-GGCC", id="d2")
        self.d3 = DNASequence(".-ACC-GTTGC--", id="d3")

        self.r1 = RNASequence("UUAU-", id="r1")
        self.r2 = RNASequence("ACGUU", id="r2")

        self.seqs1 = [self.d1, self.d2, self.d3]
        self.seqs2 = [self.r1, self.r2]

        self.seqs1_t = [("d1", "..ACC-GTTGG.."), ("d2", "TTACCGGT-GGCC"), ("d3", ".-ACC-GTTGC--")]
        self.seqs2_t = [("r1", "UUAU-"), ("r2", "ACGUU")]

        self.a1 = Alignment(self.seqs1)
        self.a2 = Alignment(self.seqs2)
        self.a3 = Alignment(self.seqs2, score=42.0, start_end_positions=[(0, 3), (5, 9)])
        self.a4 = Alignment(self.seqs2, score=-42.0, start_end_positions=[(1, 4), (6, 10)])
        self.empty = Alignment([])

 def test_position_frequencies_floating_point_precision(self):
     # Test that a position with no variation yields a frequency of exactly
     # 1.0. Note that it is important to use self.assertEqual here instead
     # of self.assertAlmostEqual because we want to test for exactly 1.0. A
     # previous implementation of Alignment.position_frequencies added
     # (1 / sequence_count) for each occurrence of a character in a position
     # to compute the frequencies (see
     # https://github.com/biocore/scikit-bio/issues/801). In certain cases,
     # this yielded a frequency slightly less than 1.0 due to roundoff
     # error. The test case here uses an alignment of 10 sequences with no
     # variation at a position. This test case exposes the roundoff error
     # present in the previous implementation because 1/10 added 10 times
     # yields a number slightly less than 1.0. This occurs because 1/10
     # cannot be represented exactly as a floating point number.
     seqs = []
     for i in range(10):
         seqs.append(DNA('A', metadata={'id': str(i)}))
     aln = Alignment(seqs)
     self.assertEqual(aln.position_frequencies(),
                      [defaultdict(float, {'A': 1.0})])

    def test_majority_consensus(self):
        # empty cases
        self.assertEqual(
            self.empty.majority_consensus(), Sequence(''))
        self.assertEqual(
            self.no_positions.majority_consensus(), RNA(''))

        # alignment where all sequences are the same
        aln = Alignment([DNA('AG', metadata={'id': 'a'}),
                         DNA('AG', metadata={'id': 'b'})])
        self.assertEqual(aln.majority_consensus(), DNA('AG'))

        # no ties
        d1 = DNA('TTT', metadata={'id': "d1"})
        d2 = DNA('TT-', metadata={'id': "d2"})
        d3 = DNA('TC-', metadata={'id': "d3"})
        a1 = Alignment([d1, d2, d3])
        self.assertEqual(a1.majority_consensus(), DNA('TT-'))

        # ties
        d1 = DNA('T', metadata={'id': "d1"})
        d2 = DNA('A', metadata={'id': "d2"})
        a1 = Alignment([d1, d2])
        self.assertTrue(a1.majority_consensus() in
                        [DNA('T'), DNA('A')])

    def test_majority_consensus(self):
        # empty cases
        self.assertTrue(
            self.empty.majority_consensus().equals(Sequence('')))
        self.assertTrue(
            self.no_positions.majority_consensus().equals(RNA('')))

        # alignment where all sequences are the same
        aln = Alignment([DNA('AG', id='a'),
                         DNA('AG', id='b')])
        self.assertTrue(aln.majority_consensus().equals(DNA('AG')))

        # no ties
        d1 = DNA('TTT', id="d1")
        d2 = DNA('TT-', id="d2")
        d3 = DNA('TC-', id="d3")
        a1 = Alignment([d1, d2, d3])
        self.assertTrue(a1.majority_consensus().equals(DNA('TT-')))

        # ties
        d1 = DNA('T', id="d1")
        d2 = DNA('A', id="d2")
        a1 = Alignment([d1, d2])
        self.assertTrue(a1.majority_consensus() in
                        [DNA('T'), DNA('A')])

    def setUp(self):
        self.d1 = DNASequence('..ACC-GTTGG..', id="d1")
        self.d2 = DNASequence('TTACCGGT-GGCC', id="d2")
        self.d3 = DNASequence('.-ACC-GTTGC--', id="d3")

        self.r1 = RNASequence('UUAU-', id="r1")
        self.r2 = RNASequence('ACGUU', id="r2")

        self.seqs1 = [self.d1, self.d2, self.d3]
        self.seqs2 = [self.r1, self.r2]

        self.seqs1_t = [('d1', '..ACC-GTTGG..'), ('d2', 'TTACCGGT-GGCC'),
                        ('d3', '.-ACC-GTTGC--')]
        self.seqs2_t = [('r1', 'UUAU-'), ('r2', 'ACGUU')]

        self.a1 = Alignment(self.seqs1)
        self.a2 = Alignment(self.seqs2)
        self.a3 = Alignment(self.seqs2, score=42.0,
                            start_end_positions=[(0, 3), (5, 9)])
        self.a4 = Alignment(self.seqs2, score=-42.0,
                            start_end_positions=[(1, 4), (6, 10)])

        # no sequences
        self.empty = Alignment([])

        # sequences, but no positions
        self.no_positions = Alignment([RNA('', id='a'), RNA('', id='b')])

def align_unaligned_seqs(seqs_fp, moltype=DNA, params=None, accurate=False):
    """Aligns unaligned sequences

    Parameters
    ----------
    seqs_fp : string
        file path of the input fasta file
    moltype : {skbio.DNA, skbio.RNA, skbio.Protein}
    params : dict-like type
        It pass the additional parameter settings to the application.
        Default is None.
    accurate : boolean
        Perform accurate alignment or not. It will sacrifice performance
        if set to True. Default is False.

    Returns
    -------
    Alignment object
        The aligned sequences.

    See Also
    --------
    skbio.Alignment
    skbio.DNA
    skbio.RNA
    skbio.Protein
    """
    # Create Mafft app.
    app = Mafft(InputHandler='_input_as_path', params=params)

    # Turn on correct sequence type
    app.Parameters[MOLTYPE_MAP[moltype]].on()

    # Do not report progress
    app.Parameters['--quiet'].on()

    # More accurate alignment, sacrificing performance.
    if accurate:
        app.Parameters['--globalpair'].on()
        app.Parameters['--maxiterate'].Value = 1000

    # Get results using int_map as input to app
    res = app(seqs_fp)

    # Get alignment as dict out of results
    alignment = Alignment.read(res['StdOut'], constructor=moltype)

    # Clean up
    res.cleanUp()

    return alignment

def parse_deblur_output(seqs_fp, derep_clusters):
    """ Parse deblur output file into an OTU map.

    Parameters
    ----------
    seqs_fp: string
        file path to deblurred sequences
    derep_clusters: dictionary
        dictionary of dereplicated sequences map

    Returns
    -------
    clusters: dictionary
        dictionary of clusters including dereplicated sequence labels

    Notes
    -----
    For each deblurred sequence in seqs_fp, use the sequence label to
    obtain all dereplicated sequence labels belonging to it
    (from derep_clusters) to create entries in a new dictionary where the keys
    are actual sequences (not the labels). Note not all sequences
    in derep_clusters will be in seqs_fp since they could have been removed in
    the artifact filtering step.
    """
    clusters = {}
    # Replace representative sequence name with actual sequence in cluster
    msa_fa = Alignment.read(seqs_fp, format='fasta')
    for label, seq in Alignment.iteritems(msa_fa):
        cluster_id = label.split(';')[0]
        seq2 = str(seq.degap())
        if seq2 not in clusters:
            clusters[seq2] = []
        if cluster_id not in derep_clusters:
            raise ValueError(
                'Seed ID %s does not exist in .uc file' % cluster_id)
        else:
            clusters[seq2].extend(derep_clusters[cluster_id])
    return clusters

    def test_update_ids_sequence_attributes_propagated(self):
        # 1 seq
        exp_sc = Alignment([
            DNA('ACGT', metadata={'id': "abc", 'description': 'desc'},
                positional_metadata={'quality': range(4)})
        ])
        exp_id_map = {'abc': 'seq1'}

        obj = Alignment([
            DNA('ACGT', metadata={'id': "seq1", 'description': 'desc'},
                positional_metadata={'quality': range(4)})
        ])

        obs_sc, obs_id_map = obj.update_ids(ids=('abc',))
        self.assertEqual(obs_sc, exp_sc)
        self.assertEqual(obs_id_map, exp_id_map)

        # 2 seqs
        exp_sc = Alignment([
            DNA('ACGT', metadata={'id': "abc", 'description': 'desc1'},
                positional_metadata={'quality': range(4)}),
            DNA('TGCA', metadata={'id': "def", 'description': 'desc2'},
                positional_metadata={'quality': range(4)[::-1]})
        ])
        exp_id_map = {'abc': 'seq1', 'def': 'seq2'}

        obj = Alignment([
            DNA('ACGT', metadata={'id': "seq1", 'description': 'desc1'},
                positional_metadata={'quality': (0, 1, 2, 3)}),
            DNA('TGCA', metadata={'id': "seq2", 'description': 'desc2'},
                positional_metadata={'quality': (3, 2, 1, 0)})
        ])

        obs_sc, obs_id_map = obj.update_ids(ids=('abc', 'def'))
        self.assertEqual(obs_sc, exp_sc)
        self.assertEqual(obs_id_map, exp_id_map)