Multiple sequence alignments

The biotite.application_v2 subpackage provides interfaces to various multiple sequence alignments (MSAs) programs. For our example we choose the software MUSCLE: The subpackage biotite.application_v2.muscle contains the class Muscle3App that does the job. First we get some homologous input sequences from NCBI Entrez.

from tempfile import NamedTemporaryFile
import biotite.sequence.io.fasta as fasta
import biotite.database.entrez as entrez

# Use cyclotide sequences again, but this time more than two
query = (
    entrez.SimpleQuery("Cyclotide") &
    entrez.SimpleQuery("cter") &
    entrez.SimpleQuery("srcdb_swiss-prot", field="Properties") ^
    entrez.SimpleQuery("Precursor")
)
uids = entrez.search(query, "protein")
temp_file = NamedTemporaryFile(suffix=".fasta", delete=False)
fasta_file = fasta.FastaFile.read(
    entrez.fetch_single_file(uids, temp_file.name, "protein", "fasta")
)
sequences = list(fasta.get_sequences(fasta_file).values())

We simply input the sequences to Muscle3App.run(), wait for the application to finish and get the resulting Alignment object from the Muscle3Result.

import biotite.application_v2.muscle as muscle

app = muscle.Muscle3App()
result = app.run(sequences).result()
alignment = result.alignment
print(alignment)

G-IPCGESCVFIPCTVTALLGCSCKDKVCYKN
G-IPCGESCVFIPC-ISTVIGCSCKNKVCYRN
G-IPCGESCVFIPC-ITAAIGCSCKSKVCYRN
G-IPCGESCVFIPC-ITGIAGCSCKSKVCYRN
GSAFCGETCVLGTC---YTPDCSCTALVCLKN
H-EPCGESCVFIPC-ITTVVGCSCKNKVCYD-
H-EPCGESCVFIPC-ITTVVGCSCKNKVCYN-
GTVPCGESCVFIPC-ITGIAGCSCKNKVCYID
G-LPCGESCVFIPC-ITTVVGCSCKNKVCYND
G-LPCGESCVFIPC-ITTVVGCSCKNKVCYNN
G-IPCGESCVFIPC-ISSVVGCSCKSKVCYLD
G-IPCAESCVWIPCTVTALLGCSCKDKVCYLD
G-IPCAESCVWIPCTVTALLGCSCKDKVCYLN
G-VPCAESCVWIPCTVTALLGCSCKDKVCYLD

In most MSA software even more information than the mere alignment can be extracted. For instance the guide tree, that was used for the alignment, is also part of the result.

import matplotlib.pyplot as plt
import biotite.sequence.graphics as graphics

tree = result.guide_tree_identity

fig, ax = plt.subplots(figsize=(6.0, 3.0), constrained_layout=True)
graphics.plot_dendrogram(ax, tree, orientation="left", show_distance=True)
_ = ax.set_xlabel("Distance")

If you have MUSCLE version 5 installed, use Muscle5App instead, which is used in the same way.

Variety of MSA software

The alternatives to MUSCLE are Clustal Omega and MAFFT. To use them, simply replace Muscle3App with ClustalOmegaApp or MafftApp.

import biotite.application_v2.clustalo as clustalo

alignment = clustalo.ClustalOmegaApp().run(sequences).result().alignment
print(alignment)

-GIPCGESCVFIPCTVTALLGCSCKDKVCYKN
-GIPCGESCVFIPCI-STVIGCSCKNKVCYRN
-GIPCGESCVFIPCI-TAAIGCSCKSKVCYRN
-GIPCGESCVFIPCI-TGIAGCSCKSKVCYRN
GSAFCGETCVLGTCY---TPDCSCTALVCLKN
-HEPCGESCVFIPCI-TTVVGCSCKNKVCYD-
-HEPCGESCVFIPCI-TTVVGCSCKNKVCYN-
GTVPCGESCVFIPCI-TGIAGCSCKNKVCYID
-GLPCGESCVFIPCI-TTVVGCSCKNKVCYND
-GLPCGESCVFIPCI-TTVVGCSCKNKVCYNN
-GIPCGESCVFIPCI-SSVVGCSCKSKVCYLD
-GIPCAESCVWIPCTVTALLGCSCKDKVCYLD
-GIPCAESCVWIPCTVTALLGCSCKDKVCYLN
-GVPCAESCVWIPCTVTALLGCSCKDKVCYLD

As shown in the output, the alignment with Clustal Omega slightly differs from the one performed with MUSCLE.

Custom substitution matrices and sequence types

If the MSA software supports protein sequence alignment and custom substitution matrices, e.g. MUSCLE and MAFFT, almost any type of sequence can be aligned. Let’s show this on the example of a nonsense alphabet.

import numpy as np
import biotite.application_v2.mafft as mafft
import biotite.sequence as seq
import biotite.sequence.align as align

alphabet = seq.Alphabet(("foo", "bar", 42))
sequences = [
    seq.GeneralSequence(alphabet, symbols) for symbols in [
        ["foo", "bar", 42, "foo", "foo", 42, 42],
        ["foo", 42, "foo", "bar", "foo", 42, 42],
    ]
]
matrix = align.SubstitutionMatrix(
    alphabet, alphabet, np.array([
        [ 100, -100, -100],
        [-100,  100, -100],
        [-100, -100,  100]
    ])
)
result = mafft.MafftApp().run(sequences, matrix=matrix).result()
alignment = result.alignment
# As the alphabet does not have characters as symbols
# the alignment cannot be directly printed
# However, we can print the trace
print(alignment.trace)

[[ 0  0]
 [ 1 -1]
 [ 2  1]
 [ 3  2]
 [-1  3]
 [ 4  4]
 [ 5  5]
 [ 6  6]]

This works, because internally the sequences and the matrix are converted into protein sequences/matrix. Then the masquerading sequences are aligned via the software and finally the sequences are mapped back into the original sequence type.