Memory-mapped access to an uncompressed, faidx-indexed genome FASTA.
The uncompressed file is memory-mapped, and (chromosome, position) is translated to a
byte offset with the .fai index, so nothing is loaded up front. The operating system
page cache serves the touched pages, and those pages are reclaimable rather than
anonymous memory. Alleles are compared as prefixes of the reference starting at the
variant position, ignoring case (soft-masked bases are lowercase).
Classes:
-
FaiEntry
–
-
IndexedFasta
–
An uncompressed FASTA and its .fai entries, keyed by gwaslab chromosome code.
Functions:
-
contig_to_gwaslab_code
–
gwaslab numeric chromosome code for a main contig name (chr1, 1, chrX, chrM, MT), else None.
-
gwaslab_code_to_contig_name
–
Bare contig name (1-22, X, Y, MT) for a gwaslab numeric chromosome code.
-
reference_is_acgt
–
Boolean array: the lengths[i] reference bases from 1-based positions[i] are in bounds
-
reference_matches
–
Boolean array: allele i equals the reference sequence starting at 1-based positions[i].
Attributes:
module-attribute
DEFAULT_MAX_GATHER_BYTES = 64 * 1024 * 1024
module-attribute
module-attribute
FASTA_FILENAME = 'genome.fa'
Attributes:
instance-attribute
instance-attribute
instance-attribute
instance-attribute
An uncompressed FASTA and its .fai entries, keyed by gwaslab chromosome code.
Methods:
Attributes:
instance-attribute
entries: Mapping[int, FaiEntry]
instance-attribute
classmethod
open(directory: Path) -> IndexedFasta
Source code in mecfs_bio/build_system/task/genome_reference_harmonization/fasta.py
| @classmethod
def open(cls, directory: Path) -> "IndexedFasta":
fasta_path = directory / FASTA_FILENAME
fai_path = directory / (FASTA_FILENAME + FAI_SUFFIX)
assert fasta_path.is_file(), f"missing {fasta_path}"
assert fai_path.is_file(), f"missing {fai_path}"
entries: dict[int, FaiEntry] = {}
for line in fai_path.read_text().splitlines():
name, length, offset, line_bases, line_bytes = line.split("\t")[:5]
code = contig_to_gwaslab_code(name)
if code is None:
continue
assert code not in entries, f"two FASTA contigs map to chromosome {code}"
entries[code] = FaiEntry(
length=int(length),
offset=int(offset),
line_bases=int(line_bases),
line_bytes=int(line_bytes),
)
return cls(fasta_path=fasta_path, entries=entries)
|
contig_to_gwaslab_code(name: str) -> int | None
gwaslab numeric chromosome code for a main contig name (chr1, 1, chrX, chrM, MT), else None.
Source code in mecfs_bio/build_system/task/genome_reference_harmonization/fasta.py
| def contig_to_gwaslab_code(name: str) -> int | None:
"""gwaslab numeric chromosome code for a main contig name (chr1, 1, chrX, chrM, MT), else None."""
match = _MAIN_CONTIG.match(name)
if match is None:
return None
token = match.group(1)
if token == "M":
token = "MT"
if token in GWASLAB_CHROM_CODE_FOR_NAME:
return GWASLAB_CHROM_CODE_FOR_NAME[token]
return int(token)
|
gwaslab_code_to_contig_name(code: int) -> str
Bare contig name (1-22, X, Y, MT) for a gwaslab numeric chromosome code.
Source code in mecfs_bio/build_system/task/genome_reference_harmonization/fasta.py
| def gwaslab_code_to_contig_name(code: int) -> str:
"""Bare contig name (1-22, X, Y, MT) for a gwaslab numeric chromosome code."""
return GWASLAB_CHROM_NAME_FOR_CODE.get(code, str(code))
|
reference_is_acgt(
fasta: IndexedFasta,
chrom: int,
positions: ndarray,
lengths: ndarray,
max_gather_bytes: int = DEFAULT_MAX_GATHER_BYTES,
) -> np.ndarray
Boolean array: the lengths[i] reference bases from 1-based positions[i] are in bounds
and each is A, C, G or T in either case. False where the span holds N or another IUPAC
ambiguity code.
Source code in mecfs_bio/build_system/task/genome_reference_harmonization/fasta.py
| def reference_is_acgt(
fasta: IndexedFasta,
chrom: int,
positions: np.ndarray,
lengths: np.ndarray,
max_gather_bytes: int = DEFAULT_MAX_GATHER_BYTES,
) -> np.ndarray:
"""Boolean array: the lengths[i] reference bases from 1-based positions[i] are in bounds
and each is A, C, G or T in either case. False where the span holds N or another IUPAC
ambiguity code."""
assert (
positions.ndim == 1 and lengths.ndim == 1 and len(positions) == len(lengths)
), "positions and lengths must be one-dimensional and the same length"
assert chrom in fasta.entries, f"chromosome {chrom} is not in {fasta.fasta_path}"
assert max_gather_bytes > 0
result = np.zeros(len(positions), dtype=bool)
if len(positions) == 0:
return result
assert (lengths > 0).all(), "spans must be non-empty"
entry = fasta.entries[chrom]
genome = np.memmap(fasta.fasta_path, dtype=np.uint8, mode="r")
for piece in _length_slices(lengths, max_gather_bytes):
reference = _gather_reference(
genome, entry, positions[piece.rows], piece.length
)
result[piece.rows] = reference.in_bounds & np.isin(
reference.bases, _ACGT_CODES
).all(axis=1)
return result
|
reference_matches(
fasta: IndexedFasta,
chrom: int,
positions: ndarray,
alleles: Series,
max_gather_bytes: int = DEFAULT_MAX_GATHER_BYTES,
) -> np.ndarray
Boolean array: allele i equals the reference sequence starting at 1-based positions[i].
Source code in mecfs_bio/build_system/task/genome_reference_harmonization/fasta.py
| def reference_matches(
fasta: IndexedFasta,
chrom: int,
positions: np.ndarray,
alleles: pl.Series,
max_gather_bytes: int = DEFAULT_MAX_GATHER_BYTES,
) -> np.ndarray:
"""Boolean array: allele i equals the reference sequence starting at 1-based positions[i]."""
assert positions.ndim == 1 and len(positions) == len(alleles), (
"positions and alleles must be one-dimensional and the same length"
)
assert chrom in fasta.entries, f"chromosome {chrom} is not in {fasta.fasta_path}"
assert max_gather_bytes > 0
result = np.zeros(len(positions), dtype=bool)
if len(positions) == 0:
return result
entry = fasta.entries[chrom]
genome = np.memmap(fasta.fasta_path, dtype=np.uint8, mode="r")
lengths = alleles.str.len_bytes().to_numpy()
assert (lengths > 0).all(), "alleles must be non-empty"
for piece in _length_slices(lengths, max_gather_bytes):
reference = _gather_reference(
genome, entry, positions[piece.rows], piece.length
)
observed = (
np.frombuffer(
"".join(alleles.gather(piece.rows).to_list()).encode("ascii"),
dtype=np.uint8,
).reshape(len(piece.rows), piece.length)
& _ASCII_UPPERCASE_MASK
)
result[piece.rows] = reference.in_bounds & (reference.bases == observed).all(
axis=1
)
return result
|