streampile

Forward-only pileups streamed from coordinate-sorted BAM and CRAM records, and a table of the alleles at every base.

Installation

pip install streampile

Quickstart

Building Pileups

A StreamingPileupBuilder reads records once and piles them up at the 0-based positions you ask for, moving forward only.

>>> from pysam import AlignmentFile
>>> from streampile import StreamingPileupBuilder
>>>
>>> with (
...     AlignmentFile("tests/data/reads.bam") as reads,
...     StreamingPileupBuilder(reads, min_base_quality=30) as builder,
... ):
...     first = builder.pileup("chr1", 10)
...     second = builder.pileup("chr1", 12)
>>>
>>> first.filtered_depth, first.bases
(4, ['A', 'T', 'G', 'A'])
>>> second.filtered_depth, second.bases
(3, ['G', 'G', 'G'])

Filter reads with read_filter, and count each template once with without_overlaps():

>>> with (
...     AlignmentFile("tests/data/reads.bam") as reads,
...     StreamingPileupBuilder(reads, read_filter=lambda read: read.is_paired) as builder,
... ):
...     pileup = builder.pileup("chr1", 50)
>>>
>>> pileup.bases, pileup.without_overlaps().bases
(['T', 'T'], ['T'])

Pass tap=writer.write to receive every record, in input order, once the builder has moved past it.

Sweeping a Territory

columns yields the pileup at every position of a span.

>>> with AlignmentFile("tests/data/reads.bam") as reads, StreamingPileupBuilder(reads) as builder:
...     sum(pileup.unfiltered_depth for pileup in builder.columns("chr1", 0, 20))
73

Tabulating Alleles

tabulate counts the reads of every allele at every base of a bedspec Territory. Alleles are normalized VCF alleles at 1-based positions, so they match a VCF by CHROM, POS, REF, and ALT.

>>> from bedspec import Bed3
>>> from bedspec import Territory
>>> from pysam import FastaFile
>>> from streampile import tabulate
>>>
>>> territory = Territory([Bed3("chr1", start=9, end=12)])
>>> with (
...     AlignmentFile("tests/data/reads.bam") as reads,
...     FastaFile("tests/data/reference.fa") as reference,
... ):
...     bases = list(tabulate(reads, reference, territory, min_base_quality=30))
>>>
>>> for base in bases:
...     print(base.pos, base.ref, base.depth, base.alts, base.alt_reads)
10 A 4 () ()
11 A 4 ('T', 'GG') (1, 1)
12 C 4 () ()

Each base also splits its reads by strand and counts no-calls apart from its depth.

Reading a Table

>>> from streampile import TabulationReader
>>>
>>> for base in TabulationReader.from_path("tests/data/counts.tsv"):
...     print(base.pos, base.depth, base.alts, base.alt_reads)
21 5 () ()
22 5 () ()
23 4 ('C',) (1,)
24 4 () ()

Command Line

streampile tabulate \
    --bam tests/data/reads.bam \
    --ref tests/data/reference.fa \
    --intervals tests/data/territory.bed \
    --min-base-quality 30 \
    --min-mapping-quality 20 \
    --out counts.tsv

##streampile-tabulation=1
##streampile-version=0.1.0
##bam=tests/data/reads.bam
##reference=tests/data/reference.fa
##territory=tests/data/territory.bed
##min_base_quality=30
##min_mapping_quality=20
##exclude_flags=0xf00
#contig	pos	ref	depth	no_calls	ref_reads	ref_fwd	ref_rev	alt_refs	alts	alt_reads	alt_fwd	alt_rev
chr1	21	T	5	0	5	3	2					
chr1	22	G	5	0	5	3	2					
chr1	23	C	4	0	3	2	1	CA	C	1	0	1
chr1	24	A	4	0	3	2	1					

Write to a .gz path with --index tbi to compress and index the table, then read a region back:

for base in TabulationReader.query("counts.tsv.gz", "chr1", 20, 24):
    print(base.pos, base.depth, base.alts)

Development and Testing

See the contributing guide for more information.

The streaming design follows the StreamingPileupBuilder of fgbio; see NOTICE.