# Copyright 2016 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. """Bio.Align.AlignInfo related tests.""" import math import unittest from Bio import AlignIO from Bio import BiopythonDeprecationWarning from Bio.Align import MultipleSeqAlignment from Bio.Align.AlignInfo import SummaryInfo from Bio.Data import IUPACData from Bio.motifs import Motif from Bio.Seq import Seq from Bio.SeqRecord import SeqRecord class AlignInfoTests(unittest.TestCase): """Test basic usage.""" def assertAlmostEqualList(self, list1, list2, **kwargs): self.assertEqual(len(list1), len(list2)) for v1, v2 in zip(list1, list2): self.assertAlmostEqual(v1, v2, **kwargs) def test_nucleotides(self): filename = "GFF/multi.fna" fmt = "fasta" msa = AlignIO.read(filename, fmt) summary = SummaryInfo(msa) alignment = msa.alignment motif = Motif("ACGT", alignment) with self.assertWarns(BiopythonDeprecationWarning): c = summary.dumb_consensus(threshold=0.1, ambiguous="N") # dumb_consensus uses ambiguous if multiple letters have the same score self.assertEqual(c, "ANGNCCCC") c = motif.counts.calculate_consensus(identity=0.1) # Instead, EMBOSS uses the first letter it encounters self.assertEqual(c, "AaGcCCCC") with self.assertWarns(BiopythonDeprecationWarning): c = summary.dumb_consensus(ambiguous="N") self.assertEqual(c, "NNNNNNNN") c = motif.counts.calculate_consensus(identity=0.7) self.assertEqual(c, "NNNNNNNN") with self.assertWarns(BiopythonDeprecationWarning): c = summary.gap_consensus(ambiguous="N") self.assertEqual(c, "NNNNNNNN") expected = {"A": 0.25, "G": 0.25, "T": 0.25, "C": 0.25} with self.assertWarns(BiopythonDeprecationWarning): m = summary.pos_specific_score_matrix(chars_to_ignore=["-"], axis_seq=c) counts = motif.counts for i in range(alignment.length): for letter in "ACGT": self.assertAlmostEqual(counts[letter][i], m[i][letter]) self.assertEqual( str(m), """ A C G T N 2.0 0.0 1.0 0.0 N 1.0 1.0 1.0 0.0 N 1.0 0.0 2.0 0.0 N 0.0 1.0 1.0 1.0 N 1.0 2.0 0.0 0.0 N 0.0 2.0 1.0 0.0 N 1.0 2.0 0.0 0.0 N 0.0 2.0 1.0 0.0 """, ) # provide the frequencies and chars to ignore explicitly. with self.assertWarns(BiopythonDeprecationWarning): ic = summary.information_content( e_freq_table=expected, chars_to_ignore=["-"] ) self.assertAlmostEqual(ic, 7.32029999423075) ic = sum(motif.relative_entropy) self.assertAlmostEqual(ic, 7.32029999423075) def test_proteins(self): letters = IUPACData.protein_letters a = MultipleSeqAlignment( [ SeqRecord(Seq("MHQAIFIYQIGYP*LKSGYIQSIRSPEYDNW-"), id="ID001"), SeqRecord(Seq("MH--IFIYQIGYAYLKSGYIQSIRSPEY-NW*"), id="ID002"), SeqRecord(Seq("MHQAIFIYQIGYPYLKSGYIQSIRSPEYDNW*"), id="ID003"), ] ) self.assertEqual(32, a.get_alignment_length()) s = SummaryInfo(a) alignment = a.alignment motif = Motif(letters + "*", alignment) counts = motif.counts with self.assertWarns(BiopythonDeprecationWarning): dumb_consensus = s.dumb_consensus() self.assertEqual(dumb_consensus, "MHQAIFIYQIGYXXLKSGYIQSIRSPEYDNW*") consensus = counts.calculate_consensus(identity=0.7) self.assertEqual(consensus, dumb_consensus) with self.assertWarns(BiopythonDeprecationWarning): c = s.gap_consensus(ambiguous="X") self.assertEqual(c, "MHXXIFIYQIGYXXLKSGYIQSIRSPEYXNWX") with self.assertWarns(BiopythonDeprecationWarning): m = s.pos_specific_score_matrix(chars_to_ignore=["-", "*"], axis_seq=c) j = 0 all_letters = s._get_all_letters() for i in range(alignment.length): for letter in letters: count = counts[letter][i] if letter in all_letters: self.assertAlmostEqual(count, m[j][letter]) else: self.assertAlmostEqual(count, 0.0) j += 1 self.assertEqual( str(m), """ A D E F G H I K L M N P Q R S W Y M 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 H 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 X 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 X 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 I 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 F 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 I 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Y 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 Q 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 I 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 G 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Y 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 X 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0 X 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 L 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 K 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 S 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 G 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Y 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 I 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Q 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 S 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 I 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 R 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 S 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 P 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 E 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Y 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 X 0.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 N 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 W 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 X 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 """, ) base_freq = 1.0 / len(letters) e_freq_table = {letter: base_freq for letter in letters} with self.assertWarns(BiopythonDeprecationWarning): ic = s.information_content( e_freq_table=e_freq_table, chars_to_ignore=["-", "*"] ) self.assertAlmostEqual(ic, 133.061475107) motif = Motif(letters, alignment) ic = sum(motif.relative_entropy) self.assertAlmostEqual(ic, 133.061475107) def test_pseudo_count(self): # use example from # http://biologie.univ-mrs.fr/upload/p202/01.4.PSSM_theory.pdf msa = MultipleSeqAlignment( [ SeqRecord(Seq("AACCACGTTTAA"), id="ID001"), SeqRecord(Seq("CACCACGTGGGT"), id="ID002"), SeqRecord(Seq("CACCACGTTCGC"), id="ID003"), SeqRecord(Seq("GCGCACGTGGGG"), id="ID004"), SeqRecord(Seq("TCGCACGTTGTG"), id="ID005"), SeqRecord(Seq("TGGCACGTGTTT"), id="ID006"), SeqRecord(Seq("TGACACGTGGGA"), id="ID007"), SeqRecord(Seq("TTACACGTGCGC"), id="ID008"), ] ) summary = SummaryInfo(msa) expected = {"A": 0.325, "G": 0.175, "T": 0.325, "C": 0.175} with self.assertWarns(BiopythonDeprecationWarning): ic = summary.information_content( e_freq_table=expected, log_base=math.exp(1), pseudo_count=1 ) self.assertAlmostEqual(ic, 7.546369561463767) ic_vector = [ 0.11112361, 0.08677812, 0.35598044, 1.29445419, 0.80272907, 1.29445419, 1.29445419, 0.80272907, 0.60929642, 0.39157892, 0.46539767, 0.03739368, ] self.assertAlmostEqualList(summary.ic_vector, ic_vector) # One more time, now using a new-style Alignment object: alignment = msa.alignment motif = Motif("ACGT", alignment) motif.background = expected motif.pseudocounts = expected self.assertAlmostEqualList(motif.relative_entropy * math.log(2), ic_vector) ic = sum(ic_vector) self.assertAlmostEqual(ic, 7.546369561463767) if __name__ == "__main__": runner = unittest.TextTestRunner(verbosity=2) unittest.main(testRunner=runner)