Reads¶
For a depth check without retaining a BAM, use
data.count_reads(file, regions, max_records=1_000_000). The result has
records and an acquisition receipt. It counts records in the indexed region
union, including original duplicates, and does not fetch mates. Exceeding the cap
raises RecordLimitError; complete counts and overflow results are cached for
offline reuse. The standalone equivalent is osteosarc.count_reads(...).
Copy just the reads you need out of the remote BAMs into small local BAMs: test data for a unit test in seconds, without downloading any BAM whole. You need SAMtools 1.21 or newer on your PATH; reopening reads you already fetched doesn't.
Make test data¶
From a terminal, name a sample and a variant:
osteosarc reads T0_tumor --assay rna-seq --variant DYNC1H1-chr14-101980529 --padding 100 --to tests/data
A sample stands for each of its indexed BAMs, here its bulk RNA-seq ones, read
four at a time. Osteosarc streams only the reads within 100 bases of the variant, saves each result
in tests/data as a small indexed BAM named for its source and the variant, and
prints the paths. It skips a BAM whose genome build doesn't match the variant's,
saying so. Give a file's key instead of a sample to read one BAM, and repeat
--variant for several variants.
In Python:
from osteosarc import Dataset
data = Dataset.open(offline=False)
source = data.file("rna-seq/reprocessed/BG003082/BG003082.Aligned.sortedByCoord.out.md.bam")
targets = data.variants(ids=["DYNC1H1-chr14-101980529"], status="ready")
subset = data.extract_reads(source, variants=targets, padding=100, to="tests/data")
print(subset.path) # tests/data/BG003082.Aligned.sortedByCoord.out.md.DYNC1H1-chr14-101980529.bam
Without a folder to put it in, the BAM stays in the cache. Asking again reuses it,
even offline. To find a BAM, see samples and files;
data.samples["T0_tumor"] lists a sample's BAMs.
Regions instead of variants¶
from osteosarc import Region
regions = [Region("chr14", 101980528, 101980530, "GRCh38")]
subset = data.extract_reads(source, regions)
with subset.open() as bam:
print(bam.count("chr14", 101980528, 101980530))
In Python, regions are zero-based and half-open. On the command line they're
one-based and inclusive, like SAMtools, and in the BAM's own genome build unless you
say otherwise with --assembly (a sample's BAMs may differ, so with a sample you must):
osteosarc reads rna-seq/reprocessed/BG003082/BG003082.Aligned.sortedByCoord.out.md.bam chr14:101980529-101980530
Region.from_samtools("chr14:101980529-101980530", assembly="GRCh38") converts one;
a Region's assembly can be None, for the build of whatever BAM it's read from.
Overlapping regions are merged. Reads come back exactly as stored, duplicates and
tags included; a region with no reads gives an empty BAM.
Filter reads¶
from osteosarc import ReadFilter
filtered = data.extract_reads(source, regions, filters=ReadFilter(min_mapq=20, exclude_flags=0x500))
Flag 0x100 marks secondary alignments and 0x400 duplicates. A filter can also keep
only some cell barcodes, as in ReadFilter(barcodes=("AAACCTG-1",)), or some read
names. On the command line, use --min-mapq and --exclude-flags.
Mates and split reads¶
--fetch-pairs also fetches each read's mate, even outside your regions, except a
mate that didn't align and has no position of its own, as STAR stores them: finding
those means reading every such read at the end of the BAM, and they have a sequence
but no alignment. That's much quicker whenever the regions hold reads whose mates
have no position, as RNA-seq regions of any size usually do: six regions of one
RNA-seq BAM took 40 s instead of 68 s. Otherwise it reads the BAM a second time for
the mates, which takes about a second longer. --unplaced-mates fetches those mates
too. In Python, fetch_pairs=True fetches every mate, as samtools view
--fetch-pairs does, and adding unplaced_mates=False leaves those out. Bundles you
make leave them out (test data).
To follow split reads as well, pass a recovery policy (or --recover-linked):
from osteosarc import RecoveryPolicy
linked = data.extract_reads(source, regions, recovery=RecoveryPolicy(max_rounds=2))
print(linked.receipt["status"], linked.receipt["records"])
Recovery fetches each mate and split-read piece from the same BAM and keeps a
partner only if its name, position and strand match the link. It never invents
one: a partner it can't find stays missing, and the receipt lists it. It stops at
limits on rounds, intervals, bases and records. Running out of rounds marks the
result truncated, and too many records is an error. With
RecoveryPolicy(on_timeout="incomplete"), a query that times out gives what was
found so far, marked incomplete, and the next call retries it.
For large seed regions, partner_batch_size=256 divides partner locations into
separate indexed queries. With on_timeout="incomplete", a slow batch is marked
missing while later batches continue; successfully retrieved partners survive.
max_partner_queries=64 caps the number of attempted partner queries, including
timeouts, and partner_timeout=30 sets their timeout in seconds independently of
seed extraction. These optional limits leave the default request strategy and
cache identities unchanged. The receipt records failed batches, visited regions
and exhausted limits; a failed query never establishes absence.
Ordinary seed context is streamed rather than retained as alignment objects in memory. All seed records still survive, with their exact stored values and multiplicities. Only names with an outgoing mate/SA pointer enter partner queries; partners already present in the seed can satisfy those pointers even without pointers of their own. An unchanged recovered BAM shares the cached seed file. Counts and per-record provenance still consume memory proportional to the number of distinct records. Recovery requests use schema 3, so earlier recovery derivatives remain intact and are not reused for the changed queries.
Genome builds¶
info = data.inspect_alignment(source)
print(info.assembly)
Osteosarc reads the genome build from the chromosome lengths in the BAM header,
and refuses a request whose regions or variants come from another build; it never
converts coordinates. GRCh37 mitochondrial regions also need the chromosome's
length, as in Region("MT", 12990, 13000, "GRCh37", reference_length=16569),
because the two GRCh37 references differ there. CRAM files need a local reference
FASTA, given with reference=.
What's cached¶
Each result is cached by its request, with a receipt of its regions, filters,
checksums, header, tool versions and read count. Osteosarc checks a remote BAM's
HTTP headers before and after reading, to catch a file that changed mid-read.
osteosarc downloads lists every extract with its local path.
When the source is a File from a snapshot, header inspection and regional
extraction also compare its inventory size and numeric modification timestamp
with the server's headers. A date mismatch raises IntegrityError before
reading the alignment or fetching its index. Cached headers and subsets are
checked against their recorded source dates too, including older receipts,
without network access or rewriting the saved provenance.
New receipts record inventory_modification: matched means both dates were
available and agree within the existing one-second tolerance. HTTP dates have
second precision (RFC 9110).
inventory_timestamp_unavailable or last_modified_unavailable means that date
check could not be performed; other source identity and byte checks still run.
Malformed server dates fail when a numeric inventory timestamp is present.
An arbitrary URL without a snapshot timestamp has no inventory date to verify.
Local BAMs¶
The same extraction works on any local indexed BAM:
from osteosarc import extract_reads
regional = extract_reads(str(subset.path), regions)
print(regional.path, regional.receipt["records"])
For a bundle of test data with a balanced set of reads at each variant, and the shared test data the OpenVax libraries use, see test data.