Genome Toolkit Series

Genome Toolkit. Part 4.5: Adding a Minimal Streaming FASTA Loader

rebelCoder

Welcome back to the Genome Toolkit series!

In Part 4.4, we built Genome Toolkit’s first loading layer. We moved our original sequence out of application.py, created the neutral SequenceRecord, added text.get_sequence(), loaded a sequence from sample.txt, converted that parsed record into validated DNA, and kept both of our original k-mer algorithms unchanged.

Our current data flow is:

sample.txt
    ↓
text.get_sequence()
    ↓
SequenceRecord
    ↓
DNA
    ↓
dna.sequence
    ↓
existing k-mer algorithms

That means most of the loading architecture is already in place.

In Part 4.5, we do not need another record model or another biological model. We only need to extend the loading layer so Genome Toolkit can understand a much more useful biological sequence format: FASTA.

We will add one new loader module:

load/fasta.py

and extend our existing loading API.

The new input side becomes:

plain text ─┐
            ├→ SequenceRecord → DNA → dna.sequence → algorithms
FASTA ──────┘

Most importantly, we are still not changing DNA or either k-mer algorithm. The algorithms still accept a normal Python string, and our application will continue passing them:

seq = dna.sequence

Part 4.6 is where we will finally upgrade the algorithm boundary itself.

Starting With Our Working Project

We continue from the exact final state of Part 4.4:

genome_toolkit/
├── .git/
├── .gitignore
├── .venv/
├── README.md
├── application.py
├── pyproject.toml
├── uv.lock
├── samples/
│   └── sample.txt
└── src/
    └── genome_toolkit/
        ├── __init__.py
        ├── py.typed
        ├── algorithms/
        │   ├── __init__.py
        │   └── kmer.py
        ├── sequence/
        │   ├── __init__.py
        │   ├── base.py
        │   └── dna.py
        └── load/
            ├── __init__.py
            ├── records.py
            └── text.py

Our existing loading boundary is already established:

record = text.get_sequence(
    SAMPLES_DIR / "sample.txt"
)

dna = DNA.model_validate(
    record,
    from_attributes=True,
)

And the final bridge into our algorithms is still:

seq = dna.sequence

so the scientific part of application.py remains:

print(f"Sequence: {seq}")
print(f"k-mer: {kmer}")
print(f"Repeats found: {count_kmer(seq, kmer)}")
print(f"Most frequent k-mer: {find_most_frequent_kmers(seq, k_len)}")

Running:

uv run python application.py

still gives us:

Sequence: AATTTTAAAAC
k-mer: AA
Repeats found: 4
Most frequent k-mer: ['TTT', 'AAA']

That output is our checkpoint again. We are going to change the source format without changing the biological validation or scientific calculations.

What Is FASTA and Why Is It So Common?

FASTA is one of the most widely used formats for storing and exchanging biological sequences. It has been used in bioinformatics for decades and has effectively become a de facto standard because it is simple, human-readable, easy for software to parse, and supported by a huge range of biological databases and analysis tools.

FASTA can store different biological sequence types, including:

DNA
RNA
protein / amino-acid sequences

The FASTA format itself does not decide what those symbols mean. It gives us a simple record structure: a header followed by sequence data. Our biological models still decide whether the sequence is valid DNA, RNA, protein, or something else.

FASTA is also one of those formats we will keep running into again and again in bioinformatics. If we download a bacterial genome, a viral genome, a chromosome, an assembled contig, a transcript sequence, or a protein sequence from a biological database, FASTA is one of the most common formats we will see.

There is one useful distinction to remember. Reference and assembled sequences are very commonly distributed as FASTA, while raw sequencing reads are commonly distributed as FASTQ, because FASTQ stores a quality score for every sequenced symbol in addition to the sequence itself.

So as Genome Toolkit moves from small examples toward real biological datasets and experiments, being able to read FASTA is not an optional convenience. It gives us access to one of the standard ways biological sequence data is stored and exchanged.

If FASTA is completely new to you, I recommend reading the Wikipedia overview before continuing. It gives a useful introduction to the format, its history, and the basic record structure:

https://en.wikipedia.org/wiki/FASTA_format

If you followed our earlier DNA Toolkit series, FASTA should already look familiar. We used FASTA files there to load biological sequences and work with them further, including extracting protein sequences from the loaded DNA.

This time, however, we are building FASTA support directly into Genome Toolkit itself. We want the loader to do more than simply read one small example file: it should be able to scan FASTA records one line at a time, find the record we need, and avoid keeping unrelated sequences in memory. We will see exactly how that works later in this article when we build get_sequence().

At its simplest, one FASTA record looks like this:

>identifier optional description
SEQUENCE

A single FASTA file can contain many records:

>original_example Original Genome Toolkit sequence
AATTTT
AAAAC

>second_example Another small sequence
ACGT
ACGT

For the first record, we want to extract:

identifier  → original_example
description → Original Genome Toolkit sequence
sequence    → AATTTTAAAAC

Notice that the sequence is wrapped across two lines in the file, but biologically it is still one continuous sequence.

Compared with our plain-text loader, FASTA adds:

headers
identifiers
optional descriptions
wrapped sequence lines
multiple records per file

So even though Part 4.5 changes only a few files, fasta.py itself contains enough new logic to deserve a proper walkthrough.

Our FASTA Scope for Part 4.5

In this part, we are going to add support for the common FASTA structure used by NCBI and many other sequence resources:

>identifier optional description
SEQUENCE

The first whitespace-separated value after > will become the identifier. Everything after it will be preserved as the optional description.

Our new loader will understand this general FASTA structure, but it will not try to decode provider-specific metadata inside the header.

UniProt is one of the major resources for protein sequences and protein annotation. Its FASTA headers can include information such as protein accessions, entry names, organism names, gene names, and other identifiers connected to the protein record.

Ensembl is a major genome annotation resource. It provides reference genomes together with annotated genes, transcripts, proteins, and related genomic information for many species. Its FASTA headers can also include extra identifiers and metadata describing the sequence record.

So although both resources can provide FASTA files, the text inside their headers can carry more information than the simple:

>identifier optional description

structure we need to understand for Genome Toolkit right now.

We do not need to decode all of those provider-specific details yet. For Part 4.5, our goal is simply:

read a FASTA file
discover records
select one record
return SequenceRecord

This simple FASTA loader will be more than enough for Genome Toolkit and many of the real genome-data experiments we are planning next.

At the same time, it is useful to see where this loading layer could grow in the future.

There are two different kinds of extensions we may eventually need:

actual biological file formats
+
provider-specific interpretation of data inside those formats

For example, our loading package could one day grow into something like this:

src/genome_toolkit/
└── load/
    ├── records.py
    │
    ├── fasta.py       # generic FASTA: DNA / RNA / protein sequences
    ├── fastq.py       # sequence + per-base quality scores
    ├── genbank.py     # sequence + rich biological annotations/features
    ├── embl.py        # EMBL sequence records + annotations
    │
    ├── uniprot.py     # interpret UniProt-specific FASTA metadata
    └── ensembl.py     # interpret Ensembl-specific FASTA metadata

The important distinction is:

FASTA / FASTQ / GenBank / EMBL
→ actual biological data formats

UniProt / Ensembl / NCBI
→ data providers that may use those formats
  while adding their own metadata and conventions

So a generic FASTA record might still look like:

record = SequenceRecord(
    identifier="M57671.1",
    description="Octodon degus insulin mRNA, complete cds",
    sequence="AATTTTAAAAC",
)

while a future UniProt-aware loader could extract more protein-specific metadata from a UniProt FASTA header:

protein = Protein(
    identifier="P12345",
    entry_name="ABC_HUMAN",
    organism="Homo sapiens",
    taxonomy_id=9606,
    gene="ABC1",
    sequence="MKTLLVAG...",
)

A future FASTQ loader would solve a different problem because FASTQ stores sequence reads together with per-base quality information:

read = SequenceRead(
    identifier="read_001",
    sequence="ACGTACGT",
    quality="IIIIIIII",
)

And a GenBank loader could preserve rich biological annotations and features together with the sequence:

record = GenBankRecord(
    identifier="M57671.1",
    organism="Octodon degus",
    sequence="AATTTTAAAAC...",
    features=[...],
)

Other biological data, such as:

GFF / GTF annotations
VCF variants

would probably deserve their own focused modules rather than being forced into the same sequence-loader shape.

That is the long-term direction: Genome Toolkit can gradually become a more robust biological data-processing package by adding focused support for the formats and metadata we actually encounter in our work.

But we are not going to build all of that now.

For the experiments coming next, a simple generic FASTA loader is more than enough. We will add new loaders, richer metadata handling, and other format-specific functionality only when a real experiment gives us a reason to need them.

What Are We Adding?

We need one new module:

src/genome_toolkit/load/fasta.py

with three functions:

_parse_header()     shared FASTA header parsing
get_headers()       discover the records in a FASTA file
get_sequence()      stream through the file and return one record

That is the whole FASTA layer for Part 4.5.

The existing boundary from Part 4.4 stays the same:

FASTA
  ↓
SequenceRecord
  ↓
DNA

The FASTA loader reads the file and returns neutral sequence data. Our DNA model still decides whether those symbols are valid DNA.

Creating fasta.py

Inside our existing loading package, create:

src/genome_toolkit/load/fasta.py

Our loading package becomes:

load/
├── __init__.py
├── fasta.py                   # <-- NEW
├── records.py
└── text.py

Start with:

"""Streaming helpers for reading FASTA files."""

from pathlib import Path

from .records import SequenceRecord

Both imports are familiar from Part 4.4. Path gives us a consistent way to work with the FASTA filepath, and SequenceRecord is still the neutral object returned by our loaders.

Parsing One FASTA Header

Both public FASTA functions need to understand a header such as:

>M57671.1 Octodon degus insulin mRNA, complete cds

Add one small helper:

def _parse_header(line: str) -> tuple[str, str | None]:
    """Parse one FASTA header into its identifier and optional description.

    Args:
        line: FASTA header line including the leading `>` marker.

    Returns:
        A tuple containing the sequence identifier and its optional
        description. The description is `None` when the header contains
        only an identifier.

    Raises:
        ValueError: If the header does not contain a sequence identifier.
    """
    # Remove the leading ">" marker and surrounding whitespace.
    header = line[1:].strip()

    # Reject a header with no identifier.
    if not header:
        raise ValueError("A FASTA header must contain a sequence identifier.")

    # Split the identifier from the optional description.
    parts = header.split(maxsplit=1)

    # Extract the identifier and optional description.
    identifier = parts[0]
    description = parts[1] if len(parts) == 2 else None

    # Return both values.
    return identifier, description

The leading underscore in:

_parse_header()

marks it as an internal helper. Our public functions will be get_headers() and get_sequence().

The first line:

header = line[1:].strip()

removes the leading > and surrounding whitespace.

Then:

parts = header.split(maxsplit=1)

splits the header only once.

For:

M57671.1 Octodon degus insulin mRNA, complete cds

we get:

identifier  → M57671.1
description → Octodon degus insulin mRNA, complete cds

Using maxsplit=1 is important because the description can contain many words.

Finally:

return identifier, description

returns both values as:

tuple[str, str | None]

The description is optional, so a header such as:

>M57671.1

simply returns:

("M57671.1", None)

Discovering FASTA Records With get_headers()

Now let us add our first public FASTA function.

Its job is simple:

Scan the file and tell us which FASTA records are available.

Add:

def get_headers(
    filepath: str | Path,
) -> list[tuple[str, str | None]]:
    """Read FASTA headers and return the available records in file order.

    The function scans the FASTA source without storing sequence data.
    Each returned item contains the record identifier and its optional
    description.

    Args:
        filepath: Path to the FASTA file as a string or `Path` object.

    Returns:
        A list of `(identifier, description)` tuples in file order.
        `description` is `None` when a header contains only an identifier.

    Raises:
        ValueError: If the source contains no FASTA records, sequence data
            appears before the first header, or a FASTA header does not
            contain a sequence identifier.
    """
    path = Path(filepath)
    headers: list[tuple[str, str | None]] = []

As in Part 4.4, we normalize the filepath with:

path = Path(filepath)

Then we create an empty list for the parsed headers.

Now scan the file:

with path.open("r", encoding="utf-8") as file:
    for raw_line in file:
        line = raw_line.strip()

        if not line:
            continue

        if line.startswith(">"):
            headers.append(_parse_header(line))
            continue

        if not headers:
            raise ValueError(
                f"'{path}' is not a valid FASTA source: "
                "sequence data appears before the first header."
            )

The logic is compact.

Blank lines are ignored.

If the line starts with:

>

we parse it and append the resulting:

(identifier, description)

tuple to headers.

If we encounter sequence data before we have seen even one header:

if not headers:

the file does not have the FASTA structure we expect, so we stop with a clear error.

After the scan:

if not headers:
    raise ValueError(f"'{path}' does not contain any FASTA records.")

return headers

An empty file, or a file containing only blank lines, therefore cannot silently return an empty record list.

The complete function is:

def get_headers(
    filepath: str | Path,
) -> list[tuple[str, str | None]]:
    """Read FASTA headers and return the available records in file order.

    The function scans the FASTA source without storing sequence data.
    Each returned item contains the record identifier and its optional
    description.

    Args:
        filepath: Path to the FASTA file as a string or `Path` object.

    Returns:
        A list of `(identifier, description)` tuples in file order.
        `description` is `None` when a header contains only an identifier.

    Raises:
        ValueError: If the source contains no FASTA records, sequence data
            appears before the first header, or a FASTA header does not
            contain a sequence identifier.
    """
    # Convert the supplied filepath into a Path.
    path = Path(filepath)

    # Prepare an empty list for discovered headers.
    headers: list[tuple[str, str | None]] = []

    # Open the FASTA file.
    with path.open("r", encoding="utf-8") as file:
        # Read the file one line at a time.
        for raw_line in file:
            # Remove surrounding whitespace.
            line = raw_line.strip()

            # Skip empty lines.
            if not line:
                continue

            # If this is a header, parse and save it.
            if line.startswith(">"):
                headers.append(_parse_header(line))
                continue

            # Reject sequence data before the first header.
            if not headers:
                raise ValueError(
                    f"'{path}' is not a valid FASTA source: "
                    "sequence data appears before the first header."
                )

    # Reject a file containing no FASTA headers.
    if not headers:
        raise ValueError(f"'{path}' does not contain any FASTA records.")

    # Return all discovered headers.
    return headers

Notice that get_headers() does not store the sequence data. It scans the file and collects only the identifiers and descriptions we need for record discovery.

Selecting One FASTA Record

Now we can build the main loader:

get_sequence()

We want to select one record either by identifier:

fasta.get_sequence(
    "samples/sample.fasta",
    identifier="M57671.1",
)

or by its zero-based position:

fasta.get_sequence(
    "samples/sample.fasta",
    index=0,
)

Both calls will return one SequenceRecord.

Creating the get_sequence() Boundary

Start with:

def get_sequence(
    filepath: str | Path,
    *,
    identifier: str | None = None,
    index: int | None = None,
) -> SequenceRecord:
    """Return one FASTA record by identifier or zero-based index.

    Args:
        filepath: Path to the FASTA file.
        identifier: Identifier of the record to load.
        index: Zero-based position of the record to load.

    Returns:
        The selected FASTA record as a neutral `SequenceRecord`.

    Raises:
        ValueError: If selector arguments are invalid, sequence data appears
            before the first header, the requested record cannot be found, or
            the selected record contains no sequence data.
    """

The * makes identifier and index keyword-only arguments.

That gives us clear calls such as:

identifier="M57671.1"

or:

index=0

instead of relying on positional values whose meaning is harder to see.

Requiring One Selector

We want exactly one of these:

identifier
index

So add:

if (identifier is None) == (index is None):
    raise ValueError("Provide exactly one of 'identifier' or 'index'.")

The rule is:

identifier only → valid
index only      → valid
both            → invalid
neither         → invalid

Then clean and check the supplied value:

if identifier is not None:
    identifier = identifier.strip()

    if not identifier:
        raise ValueError("'identifier' must not be empty.")

if index is not None and index < 0:
    raise ValueError("'index' is zero-based and must be non-negative.")

Our FASTA indexes are zero-based:

first record  → 0
second record → 1
third record  → 2

Preparing the Record Search

Now add the small amount of state we need while moving through the file:

path = Path(filepath)
current_index = 0
seen_header = False
selected_identifier: str | None = None
selected_description: str | None = None
sequence_parts: list[str] = []

current_index tracks the current FASTA record.

We start directly at:

0

because FASTA record indexes are zero-based. That gives us a simple reading order:

check record 0
advance to 1
check record 1
advance to 2
...

seen_header lets us distinguish a real FASTA record from sequence data that appears before the first header, without using a special negative index value.

selected_identifier and selected_description stay empty until we find the requested record.

sequence_parts will store only the wrapped sequence lines belonging to that selected record.

Streaming Through the FASTA File

Now scan the file:

with path.open("r", encoding="utf-8") as file:
    for raw_line in file:
        line = raw_line.strip()

        if not line:
            continue

        if line.startswith(">"):
            if selected_identifier is not None:
                break

            current_identifier, current_description = _parse_header(line)

            matches = (
                current_identifier == identifier
                if identifier is not None
                else current_index == index
            )

            if matches:
                selected_identifier = current_identifier
                selected_description = current_description

            current_index += 1
            continue

Every header starts another FASTA record.

Once:

selected_identifier is not None

the target record has already been found. The next header therefore means that selected sequence is complete, so:

break

stops reading the file.

For each header, we first parse and check the current record using its existing zero-based index. Only after that header has been processed do we advance:

current_index += 1

So the order stays easy to follow:

check record 0 → advance to 1
check record 1 → advance to 2
check record 2 → advance to 3

Matching by Identifier or Index

This part decides whether the current record is the one we want:

matches = (
    current_identifier == identifier
    if identifier is not None
    else current_index == index
)

If the caller supplied an identifier, we compare identifiers.

Otherwise, we compare the zero-based index.

So:

identifier="original_example"

means:

current_identifier == "original_example"

while:

index=4

means:

current_index == 4

If the record matches:

if matches:
    selected_identifier = current_identifier
    selected_description = current_description

we remember its metadata.

Collecting Only the Selected Sequence

For a sequence line:

if not seen_header:
    raise ValueError(
        f"'{path}' is not a valid FASTA source: "
        "sequence data appears before the first header."
    )

if selected_identifier is not None:
    sequence_parts.append("".join(line.split()))

seen_header remains False until the first FASTA header appears, so sequence data before that point is invalid FASTA structure for our loader.

Once the requested record has been selected, its sequence lines are added to:

sequence_parts

For:

AATTTT
AAAAC

we collect:

["AATTTT", "AAAAC"]

Records before the target are scanned, but their sequence lines are not stored.

When the next header appears after the selected record, the loop stops.

Finishing the Record

After the scan, we handle three simple cases.

If we never saw a header:

if not seen_header:
    raise ValueError(f"'{path}' does not contain any FASTA records.")

If the requested identifier or index was never found:

if selected_identifier is None:
    target = identifier if identifier is not None else f"index {index}"

    raise ValueError(
        f"No sequence found matching '{target}' in '{path}'."
    )

Then join the wrapped sequence lines:

sequence = "".join(sequence_parts)

If the selected record contains no sequence:

if not sequence:
    raise ValueError(
        f"Sequence '{selected_identifier}' in '{path}' "
        "does not contain sequence data."
    )

Finally, return our existing neutral record:

return SequenceRecord(
    identifier=selected_identifier,
    description=selected_description,
    sequence=sequence,
)

The FASTA-specific work is finished.

The Complete get_sequence() Function

Here is the complete function we just built:

def get_sequence(
    filepath: str | Path,
    *,
    identifier: str | None = None,
    index: int | None = None,
) -> SequenceRecord:
    """Return one FASTA record by identifier or zero-based index.

    Args:
        filepath: Path to the FASTA file.
        identifier: Identifier of the record to load.
        index: Zero-based position of the record to load.

    Returns:
        The selected FASTA record as a neutral `SequenceRecord`.

    Raises:
        ValueError: If selector arguments are invalid, sequence data appears
            before the first header, the requested record cannot be found, or
            the selected record contains no sequence data.
    """
    # Require exactly one selector: identifier OR index.
    if (identifier is None) == (index is None):
        raise ValueError("Provide exactly one of 'identifier' or 'index'.")

    # Validate the identifier if one was provided.
    if identifier is not None:
        identifier = identifier.strip()

        if not identifier:
            raise ValueError("'identifier' must not be empty.")

    # Validate the zero-based index if one was provided.
    if index is not None and index < 0:
        raise ValueError("'index' is zero-based and must be non-negative.")

    # Convert the supplied filepath into a Path.
    path = Path(filepath)

    # Start with the first FASTA record at index 0.
    current_index = 0

    # Track whether we have seen at least one FASTA header.
    seen_header = False

    # Prepare storage for the selected record.
    selected_identifier: str | None = None
    selected_description: str | None = None
    sequence_parts: list[str] = []

    # Open the FASTA file.
    with path.open("r", encoding="utf-8") as file:
        # Read the file one line at a time.
        for raw_line in file:
            # Remove surrounding whitespace.
            line = raw_line.strip()

            # Skip empty lines.
            if not line:
                continue

            # If this is a FASTA header:
            if line.startswith(">"):
                # Stop if the selected record has already ended.
                if selected_identifier is not None:
                    break

                # Mark that a FASTA header has been seen.
                seen_header = True

                # Parse the current header.
                current_identifier, current_description = _parse_header(line)

                # Check whether this record matches the requested ID or index.
                matches = (
                    current_identifier == identifier
                    if identifier is not None
                    else current_index == index
                )

                # Remember the selected identifier and description.
                if matches:
                    selected_identifier = current_identifier
                    selected_description = current_description

                # Move to the next record number after processing this header.
                current_index += 1

                # Continue to the next line.
                continue

            # Reject sequence data before the first header.
            if not seen_header:
                raise ValueError(
                    f"'{path}' is not a valid FASTA source: "
                    "sequence data appears before the first header."
                )

            # If this is our selected record, save the sequence line.
            if selected_identifier is not None:
                sequence_parts.append("".join(line.split()))

    # Reject a file containing no FASTA records.
    if not seen_header:
        raise ValueError(f"'{path}' does not contain any FASTA records.")

    # Reject the request if the selected record was never found.
    if selected_identifier is None:
        target = identifier if identifier is not None else f"index {index}"

        raise ValueError(
            f"No sequence found matching '{target}' in '{path}'."
        )

    # Join the selected sequence lines.
    sequence = "".join(sequence_parts)

    # Reject a selected record containing no sequence data.
    if not sequence:
        raise ValueError(
            f"Sequence '{selected_identifier}' in '{path}' "
            "does not contain sequence data."
        )

    # Return the selected data as a SequenceRecord.
    return SequenceRecord(
        identifier=selected_identifier,
        description=selected_description,
        sequence=sequence,
    )

What Makes This Streaming?

The key line is:

for raw_line in file:

Python gives us one line at a time instead of loading the entire FASTA file into another large string first.

For get_sequence():

records before target
    scan, do not store their sequence

selected record
    store its sequence lines

next header
    stop reading

records after target
    never read

So memory mainly depends on the selected sequence rather than every sequence in the FASTA file.

We still scan from the beginning until we reach the target. This is not random-access indexing, but it is exactly the behavior we need for our current experiments.

Why Use _parse_header() as a Helper?

Both public functions need to interpret FASTA headers in exactly the same way.

Without a shared helper:

get_headers()     header parsing
get_sequence()    same header parsing again

With _parse_header():

get_headers()  ─┐
                ├→ _parse_header()
get_sequence() ─┘

we keep that format rule in one place.

get_headers() can focus on discovering the available records, while get_sequence() can focus on finding and returning one selected record.

Biological validation is still separate. The FASTA loader understands the file structure; DNA decides whether the selected sequence contains valid DNA symbols.

The Complete fasta.py

Our complete new module is now:

"""Streaming helpers for reading FASTA files."""

from pathlib import Path

from .records import SequenceRecord


def _parse_header(line: str) -> tuple[str, str | None]:
    """Parse one FASTA header into its identifier and optional description.

    Args:
        line: FASTA header line including the leading `>` marker.

    Returns:
        A tuple containing the sequence identifier and its optional
        description. The description is `None` when the header contains
        only an identifier.

    Raises:
        ValueError: If the header does not contain a sequence identifier.
    """
    # Remove the leading ">" marker and surrounding whitespace.
    header = line[1:].strip()

    # Reject a header with no identifier.
    if not header:
        raise ValueError("A FASTA header must contain a sequence identifier.")

    # Split the identifier from the optional description.
    parts = header.split(maxsplit=1)

    # Extract the identifier and optional description.
    identifier = parts[0]
    description = parts[1] if len(parts) == 2 else None

    # Return both values.
    return identifier, description


def get_headers(
    filepath: str | Path,
) -> list[tuple[str, str | None]]:
    """Read FASTA headers and return the available records in file order.

    The function scans the FASTA source without storing sequence data.
    Each returned item contains the record identifier and its optional
    description.

    Args:
        filepath: Path to the FASTA file as a string or `Path` object.

    Returns:
        A list of `(identifier, description)` tuples in file order.
        `description` is `None` when a header contains only an identifier.

    Raises:
        ValueError: If the source contains no FASTA records, sequence data
            appears before the first header, or a FASTA header does not
            contain a sequence identifier.
    """
    # Convert the supplied filepath into a Path.
    path = Path(filepath)

    # Prepare an empty list for discovered headers.
    headers: list[tuple[str, str | None]] = []

    # Open the FASTA file.
    with path.open("r", encoding="utf-8") as file:
        # Read the file one line at a time.
        for raw_line in file:
            # Remove surrounding whitespace.
            line = raw_line.strip()

            # Skip empty lines.
            if not line:
                continue

            # If this is a header, parse and save it.
            if line.startswith(">"):
                headers.append(_parse_header(line))
                continue

            # Reject sequence data before the first header.
            if not headers:
                raise ValueError(
                    f"'{path}' is not a valid FASTA source: "
                    "sequence data appears before the first header."
                )

    # Reject a file containing no FASTA headers.
    if not headers:
        raise ValueError(f"'{path}' does not contain any FASTA records.")

    # Return all discovered headers.
    return headers


def get_sequence(
    filepath: str | Path,
    *,
    identifier: str | None = None,
    index: int | None = None,
) -> SequenceRecord:
    """Return one FASTA record by identifier or zero-based index.

    Args:
        filepath: Path to the FASTA file.
        identifier: Identifier of the record to load.
        index: Zero-based position of the record to load.

    Returns:
        The selected FASTA record as a neutral `SequenceRecord`.

    Raises:
        ValueError: If selector arguments are invalid, sequence data appears
            before the first header, the requested record cannot be found, or
            the selected record contains no sequence data.
    """
    # Require exactly one selector: identifier OR index.
    if (identifier is None) == (index is None):
        raise ValueError("Provide exactly one of 'identifier' or 'index'.")

    # Validate the identifier if one was provided.
    if identifier is not None:
        identifier = identifier.strip()

        if not identifier:
            raise ValueError("'identifier' must not be empty.")

    # Validate the zero-based index if one was provided.
    if index is not None and index < 0:
        raise ValueError("'index' is zero-based and must be non-negative.")

    # Convert the supplied filepath into a Path.
    path = Path(filepath)

    # Start with the first FASTA record at index 0.
    current_index = 0

    # Track whether we have seen at least one FASTA header.
    seen_header = False

    # Prepare storage for the selected record.
    selected_identifier: str | None = None
    selected_description: str | None = None
    sequence_parts: list[str] = []

    # Open the FASTA file.
    with path.open("r", encoding="utf-8") as file:
        # Read the file one line at a time.
        for raw_line in file:
            # Remove surrounding whitespace.
            line = raw_line.strip()

            # Skip empty lines.
            if not line:
                continue

            # If this is a FASTA header:
            if line.startswith(">"):
                # Stop if the selected record has already ended.
                if selected_identifier is not None:
                    break

                # Mark that a FASTA header has been seen.
                seen_header = True

                # Parse the current header.
                current_identifier, current_description = _parse_header(line)

                # Check whether this record matches the requested ID or index.
                matches = (
                    current_identifier == identifier
                    if identifier is not None
                    else current_index == index
                )

                # Remember the selected identifier and description.
                if matches:
                    selected_identifier = current_identifier
                    selected_description = current_description

                # Move to the next record number after processing this header.
                current_index += 1

                # Continue to the next line.
                continue

            # Reject sequence data before the first header.
            if not seen_header:
                raise ValueError(
                    f"'{path}' is not a valid FASTA source: "
                    "sequence data appears before the first header."
                )

            # If this is our selected record, save the sequence line.
            if selected_identifier is not None:
                sequence_parts.append("".join(line.split()))

    # Reject a file containing no FASTA records.
    if not seen_header:
        raise ValueError(f"'{path}' does not contain any FASTA records.")

    # Reject the request if the selected record was never found.
    if selected_identifier is None:
        target = identifier if identifier is not None else f"index {index}"

        raise ValueError(
            f"No sequence found matching '{target}' in '{path}'."
        )

    # Join the selected sequence lines.
    sequence = "".join(sequence_parts)

    # Reject a selected record containing no sequence data.
    if not sequence:
        raise ValueError(
            f"Sequence '{selected_identifier}' in '{path}' "
            "does not contain sequence data."
        )

    # Return the selected data as a SequenceRecord.
    return SequenceRecord(
        identifier=selected_identifier,
        description=selected_description,
        sequence=sequence,
    )

Exporting FASTA Beside Plain Text

We have added the implementation, so now update:

src/genome_toolkit/load/__init__.py

Our current file is:

"""Sequence loaders and neutral records."""

from . import text
from .records import SequenceRecord

__all__ = ["SequenceRecord", "text"]

Add the new source:

"""Sequence loaders and neutral records."""

from . import fasta, text
from .records import SequenceRecord

__all__ = ["SequenceRecord", "fasta", "text"]

Now our public loading API grows naturally:

text.get_sequence(...)
fasta.get_headers(...)
fasta.get_sequence(...)

We do not need names such as:

get_text_sequence()
get_fasta_sequence()
get_fasta_headers()

in one flat namespace. The source module already tells us which format we are working with.

Adding a Repeatable FASTA Sample

For Part 4.5, we will use a multi-record FASTA sample rather than a file created only around our short test sequence.

Create:

samples/sample.fasta

with:

>M57671.1 Octodon degus insulin mRNA, complete cds
TGCGTTAGGCTAAACCTTGGGCCCGGCTTAGGCTAAACCTTGGGCCCGGCTTGGGCCCGGCTT
AGGCTAAACCTTGGGCCCGGCTTAGGCTAAACCTTGGGCCCGGCTTAGGCTAAACCTTGGGCC

>ALPHA_GENE_4582 one of the alpha proteins
TGCGTTAGGCTAAACCTTGGGCCCGGCTTAGGCTAAACCTTGGGCCCGGCTTGGGCCCGGCTT
AGGCTAAACCTTGGGCCCGGCTTAGGCTAAACCTTGGGCCCGGCTTAGGCTAAACCTTGGGCC

>BETA_REGULATOR_991 a regulatory sequence from the beta cluster
CCGGTTAAGCTTCCGGTTAAGCTTCCGGTTAAGCTTCCGGTTAAGCTTCCGGTTAAGCTTCC
GGTTAAGCTTCCGGTTAAGCTTCCGGTTAAGCTTCCGGTTAAGCTTCCGGTTAAGCTTCCGG

>GAMMA_OPERON_SEQ3 a short sequence from the gamma operon region
GATTACAGATTACA

>original_example Original Genome Toolkit test sequence
AATTTT
AAAAC

This gives us a much better FASTA sample for the loader because the file contains several records, different identifiers and descriptions, and sequences of different lengths.

At the end, we keep our original Genome Toolkit test sequence:

AATTTTAAAAC

as its own FASTA record:

>original_example Original Genome Toolkit test sequence
AATTTT
AAAAC

That lets us demonstrate a realistic multi-record FASTA file while still preserving the exact scientific checkpoint we have used throughout the series.

Our samples directory is now:

samples/
├── sample.txt
└── sample.fasta              # <-- NEW

Checking the Available FASTA Records

Before changing our main application, let us verify get_headers() by itself.

Use:

from genome_toolkit.load import fasta


headers = fasta.get_headers("samples/sample.fasta")
print(headers)

Running it should give us:

[
    ('M57671.1', 'Octodon degus insulin mRNA, complete cds'),
    ('ALPHA_GENE_4582', 'one of the alpha proteins'),
    ('BETA_REGULATOR_991', 'a regulatory sequence from the beta cluster'),
    ('GAMMA_OPERON_SEQ3', 'a short sequence from the gamma operon region'),
    ('original_example', 'Original Genome Toolkit test sequence')
]

In a normal terminal, Python will usually print that list on one line. The important part is that Genome Toolkit discovers all five records in file order and keeps each identifier together with its description.

The list position also gives us the zero-based indexes:

index 0 → M57671.1
index 1 → ALPHA_GENE_4582
index 2 → BETA_REGULATOR_991
index 3 → GAMMA_OPERON_SEQ3
index 4 → original_example

For our main application, we will select our existing checkpoint by identifier:

identifier="original_example"

That makes the intention clearer than depending on where the record happens to appear in the file.

Updating Our Existing application.py

Part 4.4 used the plain-text source:

from genome_toolkit.load import text

and:

record = text.get_sequence(
    SAMPLES_DIR / "sample.txt"
)

Part 4.5 changes only that source-specific part.

Replace the import with:

from genome_toolkit.load import fasta

and load our original Genome Toolkit test sequence by identifier:

record = fasta.get_sequence(
    SAMPLES_DIR / "sample.fasta",
    identifier="original_example",
)

Everything after the returned SequenceRecord stays exactly the same.

The complete application.py is:

from pathlib import Path

from genome_toolkit.algorithms import (
    count_kmer,
    find_most_frequent_kmers,
)
from genome_toolkit.load import fasta
from genome_toolkit.sequence import DNA


SAMPLES_DIR = Path(__file__).resolve().parent / "samples"


record = fasta.get_sequence(
    SAMPLES_DIR / "sample.fasta",
    identifier="original_example",
)

dna = DNA.model_validate(
    record,
    from_attributes=True,
)

seq = dna.sequence
kmer = "AA"
k_len = 3


print(f"Sequence: {seq}")
print(f"k-mer: {kmer}")
print(f"Repeats found: {count_kmer(seq, kmer)}")
print(f"Most frequent k-mer: {find_most_frequent_kmers(seq, k_len)}")

Our workflow is now:

sample.fasta
     ↓
fasta.get_sequence()
     ↓
SequenceRecord
     ↓
DNA.model_validate(...)
     ↓
validated DNA
     ↓
dna.sequence
     ↓
existing k-mer algorithms

Notice again what did not change:

src/genome_toolkit/sequence/
src/genome_toolkit/algorithms/

We have not changed the algorithm signatures, loops, or return values.

Running the Complete Part 4.5 Application

Run:

uv run python application.py

We should still see:

Sequence: AATTTTAAAAC
k-mer: AA
Repeats found: 4
Most frequent k-mer: ['TTT', 'AAA']

The biological input now came from:

samples/sample.fasta

instead of:

samples/sample.txt

but the same validated sequence reached the same algorithms and produced the same result.

That is exactly the kind of continuity we want during this refactor.

Our New Project Structure

At the end of Part 4.5, our project looks like this:

genome_toolkit/
├── .git/
├── .gitignore
├── .venv/
├── README.md
├── application.py              # <-- UPDATED
├── pyproject.toml
├── uv.lock
├── samples/
│   ├── sample.txt
│   └── sample.fasta            # <-- NEW
└── src/
    └── genome_toolkit/
        ├── __init__.py
        ├── py.typed
        ├── algorithms/
        │   ├── __init__.py
        │   └── kmer.py
        ├── sequence/
        │   ├── __init__.py
        │   ├── base.py
        │   └── dna.py
        └── load/
            ├── __init__.py      # <-- UPDATED
            ├── fasta.py         # <-- NEW
            ├── records.py
            └── text.py

No dependency was added, so:

pyproject.toml
uv.lock

remain unchanged.

The only package implementation files changed in Part 4.5 are:

load/fasta.py        # NEW
load/__init__.py     # UPDATED

plus:

sample.fasta         # NEW
application.py       # UPDATED

That is a small structural change, even though the FASTA parser itself contains several new pieces of logic.

Summary

In Part 4.5, we extended the loading layer we built in Part 4.4 and added our first FASTA loader.

The new fasta.py module gives us three focused functions:

_parse_header()
get_headers()
get_sequence()

Together they let Genome Toolkit:

parse FASTA headers
discover record identifiers and descriptions
select one record by identifier or index
join wrapped sequence lines
read the file incrementally
return SequenceRecord

Both loading paths now reach the same neutral object:

plain text ─┐
            ├→ SequenceRecord
FASTA ──────┘

and the rest of our workflow remains unchanged:

SequenceRecord
      ↓
DNA
      ↓
dna.sequence
      ↓
existing k-mer algorithms

This simple FASTA support is enough for many of the real genome-data experiments we are planning next. As our experiments grow, Genome Toolkit can grow with them by adding new loaders and richer biological-data support only when we actually need them.

New Concepts We Learned

  • FASTA — One of the most common formats for storing and exchanging DNA, RNA, and protein sequences.
  • FASTA header — A line beginning with > that contains a record identifier and may also contain a description.
  • Generic FASTA loader — Our loader understands the common FASTA record structure without trying to decode every provider-specific header convention.
  • Private helper function — _parse_header() keeps the shared header parsing rule in one place for both public FASTA functions.
  • split(maxsplit=1) — Splits the identifier from the rest of the header while keeping a multi-word description together.
  • Tuple return values and unpacking — _parse_header() returns (identifier, description), which we can assign directly to two variables.
  • Keyword-only arguments — The * in get_sequence() gives us clear calls such as identifier="..." or index=0.
  • Zero-based indexing — The first FASTA record is index 0, the second is index 1, and so on.
  • Conditional expression — matches chooses the identifier comparison when an identifier was supplied, otherwise it compares the record index.
  • continue and break — continue moves to the next loop iteration, while break stops the scan once the selected FASTA record is complete.
  • Streaming record selection — get_sequence() reads the file line by line, ignores sequence data from records we do not need, stores only the selected sequence, and stops at the next header.

The most important architectural idea remains the same:

The loader understands the file format. The biological model decides whether the selected symbols are valid DNA.

What is Next?

Our input side now supports both plain-text and FASTA sequence files:

plain text ─┐
            ├→ SequenceRecord → DNA
FASTA ──────┘

Our two k-mer functions still accept a plain string:

sequence: str

and return simple Python values:

count_kmer()               → int
find_most_frequent_kmers() → list[str]

In Part 4.6, we will improve how those algorithms receive data and how they return scientific results.

The actual k-mer calculations will stay the same.

The full source code for Genome Toolkit is available here:

https://github.com/rebelC0der/Genome_Toolkit

I hope adding our first FASTA loader and learning how to discover and stream biological sequence records through Genome Toolkit was useful for your bioinformatics and programming journey! If you found this article valuable and want to help us continue building rebelScience, please consider supporting our project. You can explore various ways to contribute here:

https://rebelscience.club/cryptocurrency-donations/

Until next time, rebelCoder, signing out.

References

Video version of this article:

Related Posts

Genome Toolkit Series

Genome Toolkit. Part 4.4: Loading Sequences From Plain Text

rebelCoder
Genome Toolkit Series

Genome Toolkit. Part 4.3: Building Validated Biological Models

rebelCoder

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.