Homology search with MMseqs2#

MMseqs2 is a software suite for fast searching and clustering of large sequence sets. Its interface is the MMseqsApp class from the biotite.application_v2.mmseqs subpackage. In contrast to the MSA programs from the previous chapter, MMseqs2 is not a single command, but a collection of subcommands that read and write databases. Accordingly, MMseqsApp has one method per supported subcommand and each of them returns a Future that resolves to a Database object.

Creating databases#

We start with some homologous sequences to search in: a few globins from different species, plus insulin, which is unrelated to the globins. Furthermore, we need a query sequence, in our case the human hemoglobin subunit delta. Both are fetched from NCBI Entrez as FASTA files.

from tempfile import NamedTemporaryFile, gettempdir
import biotite.database.entrez as entrez

globin_uids = [
    "P69905", "P01942", "P01966",  # hemoglobin subunit alpha
    "P68871", "P02088", "P02070",  # hemoglobin subunit beta
    "P02144", "P02185",            # myoglobin
    "P01308",                      # insulin
]
temp_file = NamedTemporaryFile(suffix=".fasta", delete=False)
target_path = entrez.fetch_single_file(
    globin_uids, temp_file.name, "protein", "fasta"
)
query_path = entrez.fetch("P02042", gettempdir(), "fa", "protein", "fasta")

MMseqs2 does not work on FASTA files directly, but on its own database format. Hence, the first step is to convert the files into a SequenceDatabase via MMseqsApp.create_db().

import biotite.application_v2.mmseqs as mmseqs

app = mmseqs.MMseqsApp()
target_db = app.create_db(target_path).result()
query_db = app.create_db(query_path).result()
print(type(target_db).__name__)

SequenceDatabase

A Database is a directory containing the files MMseqs2 wrote. By default, this directory is temporary and is deleted when the object is discarded. If you want to keep a database, for example to search it again in a later session, use Database.copy() or Database.move() to put it into a persistent directory.

Searching a database#

Now we search the target database for sequences homologous to the query with MMseqsApp.search(). By default, MMseqs2 only computes the score of each hit. The a option enables backtracing, i.e. the alignment path is stored as well, which we will need later on to get the actual alignments. Like in the previous chapter, the command line options are passed as keyword arguments.

alignment_db = app.search(query_db, target_db, a=True).result()
print(type(alignment_db).__name__)

AlignmentDatabase

The hits are again stored in a database, an AlignmentDatabase. To make them human-readable, MMseqs2 provides the subcommand to convert them into a BLAST-like table, available via MMseqsApp.convert_alignments(). If no output path is given, the table is written into a temporary file.

table_file = app.convert_alignments(
    alignment_db,
    format_mode=mmseqs.AlignmentFormatMode.BLAST_TABLE_WITH_HEADERS,
).result()
with open(table_file.name) as file:
    print(file.read())
table_file.close()

query	target	fident	alnlen	mismatch	gapopen	qstart	qend	tstart	tend	evalue	bits
P02042.2	P68871.2	0.931	147	10	0	1	147	1	147	4.966E-95	291
P02042.2	P02070.1	0.833	144	24	0	4	147	2	145	3.184E-82	254
P02042.2	P02088.2	0.782	147	32	0	1	147	1	147	1.770E-79	246
P02042.2	P69905.2	0.434	145	74	3	4	146	3	141	2.229E-34	116
P02042.2	P01942.2	0.400	145	79	2	4	146	3	141	4.727E-32	109
P02042.2	P01966.2	0.420	145	76	3	4	146	3	141	6.477E-32	109

Only the hemoglobin subunits are found: The sequence identity of the query to the beta chains is high, to the alpha chains still moderate, but the myoglobins are too dissimilar for a sequence search with default parameters. Also note that the identifiers are the accessions MMseqs2 parsed from the FASTA headers.

Working with Biotite objects#

So far we have used the subcommands as they are: with files as input and output. However, usually you have Sequence objects at hand and want to get Alignment objects back. The biotite.application_v2.mmseqs subpackage provides utility functions for this purpose. To begin with, create_database_from_sequences() creates a SequenceDatabase from a dictionary of sequences. This way we also get to choose the identifiers ourselves.

import biotite.sequence.io.fasta as fasta

fasta_file = fasta.FastaFile.read(target_path)
target_sequences = {
    # Use the UniProt entry name from the header as identifier,
    # e.g. 'HBA_HUMAN'
    header.split("|")[2].split()[0]: sequence
    for header, sequence in fasta.get_sequences(fasta_file).items()
}
query_sequence = fasta.get_sequence(fasta.FastaFile.read(query_path))

target_db = mmseqs.create_database_from_sequences(app, target_sequences)
query_db = mmseqs.create_database_from_sequences(
    app, {"HBD_HUMAN": query_sequence}
)
alignment_db = app.search(query_db, target_db, a=True).result()

Instead of parsing the BLAST-like table, get_alignment_table() extracts the columns of interest.

table = mmseqs.get_alignment_table(
    app, alignment_db, ["query", "target", "pident", "evalue"]
)
for target, identity, e_value in zip(
    table["target"], table["pident"], table["evalue"]
):
    print(f"{target:12}{identity:>8}{e_value:>12}")

HBB_HUMAN     93.100   4.966E-95
HBB_BOVIN     83.300   3.184E-82
HBB1_MOUSE    78.200   1.770E-79
HBA_HUMAN     43.400   2.229E-34
HBA_MOUSE     40.000   4.727E-32
HBA_BOVIN     42.000   6.477E-32

Finally, get_alignments() reconstructs an Alignment for each hit, which is the reason we enabled backtracing above. The alignments are ordered by increasing E-value, so the first one is the best hit.

alignments = mmseqs.get_alignments(app, alignment_db)
print(list(alignments.keys()))
print()
best_alignment = alignments["HBD_HUMAN", "HBB_HUMAN"]
print(best_alignment)

[('HBD_HUMAN', 'HBB_HUMAN'), ('HBD_HUMAN', 'HBB_BOVIN'), ('HBD_HUMAN', 'HBB1_MOUSE'), ('HBD_HUMAN', 'HBA_HUMAN'), ('HBD_HUMAN', 'HBA_MOUSE'), ('HBD_HUMAN', 'HBA_BOVIN')]

MVHLTPEEKTAVNALWGKVNVDAVGGEALGRLLVVYPWTQRFFESFGDLSSPDAVMGNPKVKAHGKKVLG
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLG

AFSDGLAHLDNLKGTFSQLSELHCDKLHVDPENFRLLGNVLVCVLARNFGKEFTPQMQAAYQKVVAGVAN
AFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVAN

ALAHKYH
ALAHKYH

Clustering#

Besides searching, the other main purpose of MMseqs2 is clustering sequences by their similarity, in order to reduce redundancy in a sequence set. This is done by MMseqsApp.cluster(), that returns a ClusterDatabase, and get_clusters(), that maps the representative sequence of each cluster to its members. Here we cluster the globins at a minimum sequence identity of 50 %.

cluster_db = app.cluster(target_db, min_seq_id=0.5).result()
clusters = mmseqs.get_clusters(app, cluster_db)
for representative, members in clusters.items():
    print(f"{representative}: {members}")

HBA_BOVIN: ['HBA_BOVIN', 'HBA_HUMAN', 'HBA_MOUSE']
HBB_HUMAN: ['HBB_HUMAN', 'HBB1_MOUSE', 'HBB_BOVIN']
MYG_PHYMC: ['MYG_PHYMC', 'MYG_HUMAN']
INS_HUMAN: ['INS_HUMAN']

As expected, the alpha chains, the beta chains and the myoglobins form separate clusters, while insulin is on its own.

Increasing the sensitivity with profiles#

Remember that the myoglobins were not found by the sequence search. A common solution for detecting such distant homologs is a profile search: The hits from the first search are summarized into a sequence profile via MMseqsApp.result_to_profile(), that captures the conservation of each position of the query. The resulting ProfileDatabase can be used as query for another search in place of the sequence database. Combined with a high sensitivity, requested via the s option, the myoglobins are found this time.

profile_db = app.result_to_profile(alignment_db).result()
profile_alignment_db = app.search(
    profile_db, target_db, a=True, s=7.5
).result()
table = mmseqs.get_alignment_table(
    app, profile_alignment_db, ["target", "pident", "evalue"]
)
for target, identity, e_value in zip(
    table["target"], table["pident"], table["evalue"]
):
    print(f"{target:12}{identity:>8}{e_value:>12}")

HBB_HUMAN     93.100   3.303E-86
HBB1_MOUSE    78.200   5.220E-84
HBB_BOVIN     83.300   2.321E-82
HBA_MOUSE     40.000   2.403E-72
HBA_HUMAN     42.700   3.295E-72
HBA_BOVIN     41.300   1.062E-70
MYG_HUMAN     26.200   3.773E-15
MYG_PHYMC     24.100   1.205E-11

If you are interested in the profile itself, get_matrix_from_profile() converts it into a position-specific SubstitutionMatrix.

Structure search with Foldseek#

Foldseek applies the same concepts to protein structures: It converts structures into sequences of structural states and then uses the machinery of MMseqs2 to search and cluster them. Consequently, its interface FoldseekApp provides the same methods and databases, only that create_database_from_structures() creates the database from AtomArray objects instead of sequences.

On this page
Edit on GitHub
Show Source