From 1de9260af10bb57c4b59d0865720cf7a2cd287ca Mon Sep 17 00:00:00 2001 From: Elarwei Date: Wed, 24 Jun 2026 22:43:35 +0800 Subject: [PATCH] feat(alphafold): accept a user-provided custom MSA input (#52) Add an `msa` parameter (Python) / -msa, --msa flag (CLI) to gget alphafold that lets users supply a custom multiple sequence alignment (a3m or aligned FASTA) instead of running the internal jackhmmer search. The first sequence in the MSA must be the query. When a custom MSA is provided, gget skips the jackhmmer search and the genetic-database download entirely and builds the MSA features directly from the file. New helpers detect_msa_format()/parse_custom_msa()/ read_custom_msa() handle format detection and a3m/FASTA parsing (lowercase a3m insertions are folded into the deletion matrix, matching AlphaFold's own parser). Currently supported for single-sequence (monomer) predictions; clear errors are raised otherwise. Default behavior is unchanged (backward compatible). Resolves #52. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/en/alphafold.md | 2 + docs/src/en/updates.md | 1 + gget/gget_alphafold.py | 334 ++++++++++++++++++------ gget/main.py | 13 + tests/fixtures/test_alphafold_msa.a3m | 8 + tests/fixtures/test_alphafold_msa.fasta | 6 + tests/test_alphafold.py | 91 ++++++- 7 files changed, 367 insertions(+), 88 deletions(-) create mode 100644 tests/fixtures/test_alphafold_msa.a3m create mode 100644 tests/fixtures/test_alphafold_msa.fasta diff --git a/docs/src/en/alphafold.md b/docs/src/en/alphafold.md index c161c361a..a17b6c355 100644 --- a/docs/src/en/alphafold.md +++ b/docs/src/en/alphafold.md @@ -35,6 +35,8 @@ Path to folder to save prediction results in (str). Default: "./[date_time]_gget `-jhd` `--jackhmmer_savedir` Path to a parent directory in which to store the temporary jackhmmer files (str). By default, `gget alphafold` creates a "tmp" folder in your home directory (`~/tmp/jackhmmer/`), which can take up to ~2 GB of disk space. Use this argument to place these temporary files elsewhere, e.g. on a disk with more free space. Default: None. +`-msa` `--msa` +Path to a custom multiple sequence alignment (MSA) file to use instead of running the internal jackhmmer search (str). Accepts a3m (`.a3m`) or aligned FASTA (`.fasta`, `.fa`, `.afa`) files. The first sequence in the MSA must be the query (i.e. match the input `sequence`, ignoring gaps). When provided, `gget alphafold` skips the jackhmmer search entirely, so no genetic databases are downloaded. This is useful for submitting a manually curated MSA to improve folding predictions. Currently supported for single-sequence (monomer) predictions only. Default: None. **Flags** `-mfm` `--multimer_for_monomer` diff --git a/docs/src/en/updates.md b/docs/src/en/updates.md index 89be22281..d0dd98da1 100644 --- a/docs/src/en/updates.md +++ b/docs/src/en/updates.md @@ -5,6 +5,7 @@ #### *gget* officially became part of [*scverse*](https://scverse.org/) on June 9, 2026. 🥳🥳🥳 **Version ≥ 0.30.9** (XXX XX, 2026): +- [`gget alphafold`](alphafold.md): Added a new `msa` argument (`-msa`/`--msa` on the command line) that lets you provide a custom multiple sequence alignment (a3m or aligned FASTA) instead of running the internal jackhmmer search. The first sequence in the MSA must be the query. This enables folding from a manually curated MSA and skips the genetic database download entirely (currently supported for single-sequence/monomer predictions). Resolves [issue 52](https://github.com/scverse/gget/issues/52). **Version ≥ 0.30.8** (Jun 28, 2026): - [`gget g2p`](g2p.md): Either `gene` or `--uniprot_id` is now sufficient — whichever is missing is resolved via UniProt and cached. Gene→UniProt picks the canonical reviewed human Swiss-Prot entry; the resolution and its limitations are logged. The canonical pair is **always** prepended to the result as `gene_name` / `uniprot_id` columns (and stored on `df.attrs`), so the output schema is invariant regardless of input mode. Existing call sites continue to work. diff --git a/gget/gget_alphafold.py b/gget/gget_alphafold.py index 9e578ed3b..50f790d6c 100644 --- a/gget/gget_alphafold.py +++ b/gget/gget_alphafold.py @@ -23,6 +23,7 @@ import platform # noqa: E402 import random # noqa: E402 import shutil # noqa: E402 +import string # noqa: E402 import subprocess # noqa: E402 import sys # noqa: E402 from concurrent import futures # noqa: E402 @@ -159,6 +160,112 @@ def get_jackhmmer_dir(jackhmmer_savedir: str | None = None) -> str: return os.path.expanduser(os.path.join("~", "tmp", "jackhmmer", UUID)) +# Recognized custom MSA file extensions (a3m and aligned FASTA share the same parser; +# in plain aligned FASTA there are simply no lowercase insertion characters). +CUSTOM_MSA_EXTENSIONS = (".a3m", ".fasta", ".fa", ".afa") + + +def detect_msa_format(msa_path: str) -> str: + """Validate and return the format of a user-provided custom MSA file based on its extension. + + Args: + - msa_path Path to the custom MSA file (str). + + Returns "a3m" for '.a3m' files and "fasta" for aligned FASTA files ('.fasta', '.fa', '.afa'). + + Raises ValueError if the file extension is not recognized. + """ + lower = msa_path.lower() + if lower.endswith(".a3m"): + return "a3m" + if lower.endswith((".fasta", ".fa", ".afa")): + return "fasta" + raise ValueError( + f"Custom MSA file format not recognized for '{msa_path}'. " + f"gget alphafold accepts a3m ('.a3m') or aligned FASTA ('.fasta', '.fa', '.afa') files." + ) + + +def _parse_fasta_string(fasta_string: str) -> tuple[list[str], list[str]]: + """Parse a FASTA/a3m string into (sequences, descriptions).""" + sequences: list[str] = [] + descriptions: list[str] = [] + index = -1 + for line in fasta_string.splitlines(): + line = line.strip() + if line.startswith(">"): + index += 1 + descriptions.append(line[1:]) + sequences.append("") + continue + if not line: + continue + if index == -1: + raise ValueError("Custom MSA file must be in FASTA/a3m format (sequences preceded by '>' headers).") + sequences[index] += line + return sequences, descriptions + + +def parse_custom_msa(msa_string: str) -> tuple[list[str], list[list[int]], list[str]]: + """Parse a custom a3m/aligned-FASTA MSA string into AlphaFold MSA components. + + Lowercase characters in a3m sequences are treated as insertions relative to the query + and counted into the deletion matrix (matching AlphaFold's own a3m parsing). Aligned FASTA + files contain no lowercase insertions, so their deletion matrix is all zeros. + + Args: + - msa_string Contents of the custom MSA file (str). + + Returns a tuple of (aligned_sequences, deletion_matrix, descriptions). + + Raises ValueError if no sequences are found. + """ + sequences, descriptions = _parse_fasta_string(msa_string) + if not sequences or not any(sequences): + raise ValueError("No sequences found in the provided custom MSA.") + + deletion_matrix: list[list[int]] = [] + for msa_sequence in sequences: + deletion_vec: list[int] = [] + deletion_count = 0 + for residue in msa_sequence: + if residue.islower(): + deletion_count += 1 + else: + deletion_vec.append(deletion_count) + deletion_count = 0 + deletion_matrix.append(deletion_vec) + + # Build the aligned MSA matrix by removing the lowercase (insertion) characters. + deletion_table = str.maketrans("", "", string.ascii_lowercase) + aligned_sequences = [seq.translate(deletion_table) for seq in sequences] + return aligned_sequences, deletion_matrix, descriptions + + +def read_custom_msa(msa_path: str) -> Any: + """Read a user-provided custom MSA file and return an AlphaFold ``parsers.Msa`` object. + + Args: + - msa_path Path to the custom MSA file (a3m or aligned FASTA) (str). + + Returns an ``alphafold.data.parsers.Msa`` object. + """ + from alphafold.data import parsers + + # Validate the file extension up front (raises ValueError for unsupported formats). + detect_msa_format(msa_path) + + with open(msa_path) as f: + msa_string = f.read() + + aligned_sequences, deletion_matrix, descriptions = parse_custom_msa(msa_string) + return parsers.Msa( + sequences=aligned_sequences, + deletion_matrix=deletion_matrix, + descriptions=descriptions, + ) + + def clean_up(jackhmmer_dir: str | None = None) -> None: """Function to clean up temporary files after running gget alphafold. @@ -212,6 +319,7 @@ def alphafold( show_sidechains: bool = True, verbose: bool = True, jackhmmer_savedir: str | None = None, + msa: str | list[str] | None = None, ) -> None: """Predicts the structure of a protein using a slightly simplified version of AlphaFold v2.3.0 (https://doi.org/10.1038/s41586-021-03819-2). @@ -232,6 +340,12 @@ def alphafold( files (str). By default, gget creates a "tmp" folder in the home directory ("~/tmp/jackhmmer/"), which can take up to ~2 GB of disk space. Use this argument to place these temporary files elsewhere. Default: None. + - msa Path to a custom multiple sequence alignment (MSA) file to use instead of + running the internal jackhmmer search (str). Accepts a3m ('.a3m') or aligned + FASTA ('.fasta', '.fa', '.afa') files. The first sequence in the MSA must be the + query (i.e. match the input 'sequence', ignoring gaps). When provided, gget + alphafold skips the jackhmmer search entirely (no network/database download). + Currently supported for single-sequence (monomer) predictions only. Default: None. Saves the predicted aligned error (json) and the prediction (PDB) in the defined 'out' folder. @@ -481,55 +595,93 @@ class ModelType(enum.Enum): f"The accuracy of this algorithm has not been fully validated above 3000 residues, and you may experience long running times or run out of memory. Total sequence length is {total_sequence_length} residues." ) - ## Find the closest source - if verbose: - logger.info("Finding closest source for reference database.") - - ex = futures.ThreadPoolExecutor(3) - fs = [ex.submit(fetch, source) for source in ["", "-europe", "-asia"]] - source = None - for f in futures.as_completed(fs): - source = f.result() - ex.shutdown() - break - - DB_ROOT_PATH = f"https://storage.googleapis.com/alphafold-colab{source}/latest/" - MSA_DATABASES = [ - { - "db_name": "uniref90", - "db_path": f"{DB_ROOT_PATH}uniref90_2022_01.fasta", - "num_streamed_chunks": 62, - "z_value": 144_113_457, - }, - { - "db_name": "smallbfd", - "db_path": f"{DB_ROOT_PATH}bfd-first_non_consensus_sequences.fasta", - "num_streamed_chunks": 17, - "z_value": 65_984_053, - }, - { - "db_name": "mgnify", - "db_path": f"{DB_ROOT_PATH}mgy_clusters_2022_05.fasta", - "num_streamed_chunks": 120, - "z_value": 623_796_864, - }, - ] + ## Validate user-provided custom MSA input (https://github.com/scverse/gget/issues/52) + custom_msa_paths = None + if msa is not None: + # Normalize to a list of paths + custom_msa_paths = [msa] if isinstance(msa, str) else list(msa) - # Search UniProt and construct the all_seq features (only for heteromers, not homomers). - if model_type_to_use == ModelType.MULTIMER and len(set(sequences)) > 1: - MSA_DATABASES.extend( - [ - # Swiss-Prot and TrEMBL are concatenated together as UniProt - { - "db_name": "uniprot", - "db_path": f"{DB_ROOT_PATH}uniprot_2021_04.fasta", - "num_streamed_chunks": 101, - "z_value": 225_013_025 + 565_928, - }, - ] - ) + if len(seqs) > 1 or len(custom_msa_paths) > 1: + raise ValueError( + "Custom MSA input ('msa') is currently only supported for single-sequence (monomer) " + "predictions. Please provide a single sequence together with a single MSA file." + ) + + msa_path = custom_msa_paths[0] + if not os.path.isfile(msa_path): + raise FileNotFoundError(f"Custom MSA file not found: '{msa_path}'.") + + # Validate the file format (raises ValueError for unsupported extensions). + detect_msa_format(msa_path) + + # Parse the custom MSA and check that the first sequence matches the query. + with open(msa_path) as msa_file: + aligned_sequences, _, _ = parse_custom_msa(msa_file.read()) + query_in_msa = aligned_sequences[0].replace("-", "").upper() + if query_in_msa != sequences[0].upper(): + raise ValueError( + "The first sequence in the custom MSA must be the query sequence " + "(matching the input 'sequence', ignoring gaps).\n" + f"Query sequence: {sequences[0]}\n" + f"First MSA sequence: {query_in_msa}" + ) + + if verbose: + logger.info( + f"Using user-provided custom MSA from '{msa_path}' " + f"({len(aligned_sequences)} sequences); skipping the internal jackhmmer search." + ) + + ## Find the closest source (only needed for the internal jackhmmer search) + if msa is None: + if verbose: + logger.info("Finding closest source for reference database.") + + ex = futures.ThreadPoolExecutor(3) + fs = [ex.submit(fetch, source) for source in ["", "-europe", "-asia"]] + source = None + for f in futures.as_completed(fs): + source = f.result() + ex.shutdown() + break + + DB_ROOT_PATH = f"https://storage.googleapis.com/alphafold-colab{source}/latest/" + MSA_DATABASES = [ + { + "db_name": "uniref90", + "db_path": f"{DB_ROOT_PATH}uniref90_2022_01.fasta", + "num_streamed_chunks": 62, + "z_value": 144_113_457, + }, + { + "db_name": "smallbfd", + "db_path": f"{DB_ROOT_PATH}bfd-first_non_consensus_sequences.fasta", + "num_streamed_chunks": 17, + "z_value": 65_984_053, + }, + { + "db_name": "mgnify", + "db_path": f"{DB_ROOT_PATH}mgy_clusters_2022_05.fasta", + "num_streamed_chunks": 120, + "z_value": 623_796_864, + }, + ] + + # Search UniProt and construct the all_seq features (only for heteromers, not homomers). + if model_type_to_use == ModelType.MULTIMER and len(set(sequences)) > 1: + MSA_DATABASES.extend( + [ + # Swiss-Prot and TrEMBL are concatenated together as UniProt + { + "db_name": "uniprot", + "db_path": f"{DB_ROOT_PATH}uniprot_2021_04.fasta", + "num_streamed_chunks": 101, + "z_value": 225_013_025 + 565_928, + }, + ] + ) - TOTAL_JACKHMMER_CHUNKS = sum([cfg["num_streamed_chunks"] for cfg in MSA_DATABASES]) + TOTAL_JACKHMMER_CHUNKS = sum([cfg["num_streamed_chunks"] for cfg in MSA_DATABASES]) ### Search against existing databases # Resolve the temporary jackhmmer folder (optionally user-defined via jackhmmer_savedir) @@ -551,47 +703,57 @@ class ModelType(enum.Enum): for sequence_index, sequence in enumerate(sequences, start=1): # logger.info(f"Getting MSA for sequence {sequence_index}.") - ## Manage permissions to jackhmmer binary - command = f"chmod 755 {JACKHMMER_BINARY_PATH}" - with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process: - stderr = process.stderr.read().decode("utf-8") - # Exit system if the subprocess returned with an error - if process.wait() != 0: - if stderr: - # Log the standard error if it is not empty - sys.stderr.write(stderr) - logger.error("Giving chmod 755 permissions to jackhmmer binary failed.") - return - - # Save the target sequence in a fasta file - fasta_path = os.path.join(abs_out_path, f"target_{sequence_index}.fasta") - with open(fasta_path, "w") as f: - f.write(f">query\n{sequence}") - - # Don't do redundant work for multiple copies of the same chain in the multimer - if sequence not in raw_msa_results_for_sequence: - raw_msa_results = get_msa( - fasta_path=fasta_path, - msa_databases=MSA_DATABASES, - total_jackhmmer_chunks=TOTAL_JACKHMMER_CHUNKS, - ) - raw_msa_results_for_sequence[sequence] = raw_msa_results - else: - raw_msa_results = copy.deepcopy(raw_msa_results_for_sequence[sequence]) - - ## Extract the MSAs from the Stockholm files. - # NB: deduplication happens later in pipeline.make_msa_features. + # Build the single-chain MSA(s), either from a user-provided custom MSA or via jackhmmer. single_chain_msas = [] uniprot_msa = None - for db_name, db_results in raw_msa_results.items(): - merged_msa = notebook_utils.merge_chunked_msa(results=db_results, max_hits=MAX_HITS.get(db_name)) - if merged_msa.sequences and db_name != "uniprot": - single_chain_msas.append(merged_msa) - msa_size = len(set(merged_msa.sequences)) - if verbose: - logger.info(f"{msa_size} unique sequences found in {db_name} for sequence {sequence_index}.") - elif merged_msa.sequences and db_name == "uniprot": - uniprot_msa = merged_msa + + if msa is not None: + ## Use the user-provided custom MSA instead of running jackhmmer + custom_msa = read_custom_msa(custom_msa_paths[sequence_index - 1]) + single_chain_msas.append(custom_msa) + if verbose: + msa_size = len(set(custom_msa.sequences)) + logger.info(f"{msa_size} unique sequences found in the custom MSA for sequence {sequence_index}.") + else: + ## Manage permissions to jackhmmer binary + command = f"chmod 755 {JACKHMMER_BINARY_PATH}" + with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process: + stderr = process.stderr.read().decode("utf-8") + # Exit system if the subprocess returned with an error + if process.wait() != 0: + if stderr: + # Log the standard error if it is not empty + sys.stderr.write(stderr) + logger.error("Giving chmod 755 permissions to jackhmmer binary failed.") + return + + # Save the target sequence in a fasta file + fasta_path = os.path.join(abs_out_path, f"target_{sequence_index}.fasta") + with open(fasta_path, "w") as f: + f.write(f">query\n{sequence}") + + # Don't do redundant work for multiple copies of the same chain in the multimer + if sequence not in raw_msa_results_for_sequence: + raw_msa_results = get_msa( + fasta_path=fasta_path, + msa_databases=MSA_DATABASES, + total_jackhmmer_chunks=TOTAL_JACKHMMER_CHUNKS, + ) + raw_msa_results_for_sequence[sequence] = raw_msa_results + else: + raw_msa_results = copy.deepcopy(raw_msa_results_for_sequence[sequence]) + + ## Extract the MSAs from the Stockholm files. + # NB: deduplication happens later in pipeline.make_msa_features. + for db_name, db_results in raw_msa_results.items(): + merged_msa = notebook_utils.merge_chunked_msa(results=db_results, max_hits=MAX_HITS.get(db_name)) + if merged_msa.sequences and db_name != "uniprot": + single_chain_msas.append(merged_msa) + msa_size = len(set(merged_msa.sequences)) + if verbose: + logger.info(f"{msa_size} unique sequences found in {db_name} for sequence {sequence_index}.") + elif merged_msa.sequences and db_name == "uniprot": + uniprot_msa = merged_msa notebook_utils.show_msa_info(single_chain_msas=single_chain_msas, sequence_index=sequence_index) diff --git a/gget/main.py b/gget/main.py index 7a2944b09..827b54ed5 100644 --- a/gget/main.py +++ b/gget/main.py @@ -1341,6 +1341,18 @@ def main() -> None: "which can take up to ~2 GB of disk space. Use this argument to place them elsewhere." ), ) + parser_alphafold.add_argument( + "-msa", + "--msa", + type=str, + default=None, + required=False, + help=( + "Path to a custom multiple sequence alignment (MSA) file (a3m or aligned FASTA) to use\n" + "instead of running the internal jackhmmer search. The first sequence in the MSA must be\n" + "the query (matching the input sequence, ignoring gaps). Single-sequence predictions only." + ), + ) parser_alphafold.add_argument( "-q", "--quiet", @@ -3716,6 +3728,7 @@ def main() -> None: show_sidechains=False, verbose=args.quiet, jackhmmer_savedir=args.jackhmmer_savedir, + msa=args.msa, ) ## pdb return diff --git a/tests/fixtures/test_alphafold_msa.a3m b/tests/fixtures/test_alphafold_msa.a3m new file mode 100644 index 000000000..bbdcdf099 --- /dev/null +++ b/tests/fixtures/test_alphafold_msa.a3m @@ -0,0 +1,8 @@ +>query +MKVLAAGSTKDEFGHIKLMN +>sp|homolog1 +MKVLAAGSTKDEFGHIKLMN +>sp|homolog2 +MKVLAAGSTKDEF-HIKLMN +>sp|homolog3_with_insertion +MKVefLAAGSTKDEFGHIKLMN diff --git a/tests/fixtures/test_alphafold_msa.fasta b/tests/fixtures/test_alphafold_msa.fasta new file mode 100644 index 000000000..6b66e9885 --- /dev/null +++ b/tests/fixtures/test_alphafold_msa.fasta @@ -0,0 +1,6 @@ +>query +MKVLAAGSTKDEFGHIKLMN +>homolog1 +MKVLAAGSTKDEF-HIKLMN +>homolog2 +MKVLAAGSTKDEFGHIKLMN diff --git a/tests/test_alphafold.py b/tests/test_alphafold.py index 2c05786ed..fafe9d9fe 100644 --- a/tests/test_alphafold.py +++ b/tests/test_alphafold.py @@ -5,13 +5,27 @@ from unittest.mock import patch import gget.gget_alphafold as gget_alphafold -from gget.gget_alphafold import clean_up, get_jackhmmer_dir +from gget.gget_alphafold import ( + clean_up, + detect_msa_format, + get_jackhmmer_dir, + parse_custom_msa, +) from gget.gget_setup import UUID # AlphaFold requires heavy third-party dependencies (alphafold, openmm, jackhmmer, model # parameters) that are not available in the CI/test environment, so a full prediction run # cannot be exercised here. These tests validate the user-facing jackhmmer save-directory -# option (https://github.com/scverse/gget/issues/49) at the argument/path-handling level. +# option (https://github.com/scverse/gget/issues/49) at the argument/path-handling level and +# the user-provided custom MSA input feature (https://github.com/scverse/gget/issues/52) at +# the parsing/validation level (the logic that turns a custom a3m/FASTA file into AlphaFold +# MSA features). + +FIXTURE_DIR = os.path.join(os.path.dirname(__file__), "fixtures") +A3M_FIXTURE = os.path.join(FIXTURE_DIR, "test_alphafold_msa.a3m") +FASTA_FIXTURE = os.path.join(FIXTURE_DIR, "test_alphafold_msa.fasta") + +QUERY = "MKVLAAGSTKDEFGHIKLMN" class TestAlphafoldJackhmmerSavedir(unittest.TestCase): @@ -75,5 +89,78 @@ def test_cli_exposes_jackhmmer_savedir_flag(self): self.assertIn("--jackhmmer_savedir", result.stdout) +class TestAlphafoldCustomMSA(unittest.TestCase): + def test_detect_msa_format(self): + self.assertEqual(detect_msa_format("alignment.a3m"), "a3m") + self.assertEqual(detect_msa_format("/path/to/ALIGNMENT.A3M"), "a3m") + self.assertEqual(detect_msa_format("alignment.fasta"), "fasta") + self.assertEqual(detect_msa_format("alignment.fa"), "fasta") + self.assertEqual(detect_msa_format("alignment.afa"), "fasta") + + def test_detect_msa_format_unsupported(self): + with self.assertRaises(ValueError): + detect_msa_format("alignment.txt") + with self.assertRaises(ValueError): + detect_msa_format("alignment.sto") + + def test_parse_a3m_fixture(self): + with open(A3M_FIXTURE) as f: + aligned_sequences, deletion_matrix, descriptions = parse_custom_msa(f.read()) + + # Four sequences, the first of which is the query. + self.assertEqual(len(aligned_sequences), 4) + self.assertEqual(len(descriptions), 4) + self.assertEqual(aligned_sequences[0], QUERY) + + # After removing a3m insertions, every aligned row has the query's length. + for seq in aligned_sequences: + self.assertEqual(len(seq), len(QUERY)) + + # Deletion matrix rows align 1:1 with the aligned columns. + for seq, deletions in zip(aligned_sequences, deletion_matrix, strict=True): + self.assertEqual(len(deletions), len(seq)) + + # The query and the gapped homolog carry no insertions. + self.assertEqual(sum(deletion_matrix[0]), 0) + self.assertEqual(sum(deletion_matrix[2]), 0) + + # homolog3 has two lowercase insertion characters ("ef") after the third residue. + self.assertEqual(sum(deletion_matrix[3]), 2) + self.assertEqual(deletion_matrix[3][3], 2) + + def test_parse_fasta_fixture_has_no_insertions(self): + with open(FASTA_FIXTURE) as f: + aligned_sequences, deletion_matrix, _ = parse_custom_msa(f.read()) + + self.assertEqual(aligned_sequences[0], QUERY) + # Aligned FASTA contains no lowercase insertions -> deletion matrix is all zeros. + for deletions in deletion_matrix: + self.assertEqual(sum(deletions), 0) + + def test_query_matches_first_msa_sequence(self): + # Mirrors the validation gget performs: first MSA sequence (gaps removed) == query. + with open(A3M_FIXTURE) as f: + aligned_sequences, _, _ = parse_custom_msa(f.read()) + query_in_msa = aligned_sequences[0].replace("-", "").upper() + self.assertEqual(query_in_msa, QUERY) + + def test_parse_empty_msa_raises(self): + with self.assertRaises(ValueError): + parse_custom_msa("") + + def test_parse_non_fasta_raises(self): + with self.assertRaises(ValueError): + parse_custom_msa("MKVLAAG\nNOHEADERLINE\n") + + def test_cli_exposes_msa_flag(self): + """The command-line interface exposes the --msa option for gget alphafold.""" + result = subprocess.run( + [sys.executable, "-m", "gget", "alphafold", "--help"], + capture_output=True, + text=True, + ) + self.assertIn("--msa", result.stdout) + + if __name__ == "__main__": unittest.main()