From 449d2c99c8b9edf83533863eba0e903486bdee81 Mon Sep 17 00:00:00 2001 From: Cade Mirchandani Date: Tue, 21 Jul 2026 12:52:11 -0700 Subject: [PATCH 1/2] feat!: rebuild py-ft as a pure-python toolkit over pysam + pyMA Drops the Rust/pyo3 bindings (Fiberbam, Fiberdata, Fiberwriter, center) and the maturin build. py-ft is now plain Python: Fiberdata wraps a pysam record's molecular_annotation and exposes the fiber-seq accessors (m6a, cpg, nuc, msp, fire) over iter_type. fetch() iterates mapped and unmapped reads (so unaligned fiber-seq BAMs work), and Fiberdata.is_mapped lets callers filter. Plotting moves to an optional [viz] extra. region_to_* and the fiberplot CLI, which drove the old Rust Fiberbam, are stubbed with NotImplementedError until they're ported to fetch(); plot.py is unchanged. Adds a root release workflow for py-ft (py-ft-v* tags) using PyPI trusted publishing; the pyMA release workflow lives with the pyMA changes. --- .github/workflows/release-python-pyft.yml | 40 ++ py-ft/.github/workflows/CI.yml | 120 ------ py-ft/Cargo.toml | 21 - py-ft/pyproject.toml | 67 ++- py-ft/python/pyft/__init__.py | 16 +- py-ft/python/pyft/fiberdata.py | 144 +++++++ py-ft/python/pyft/fiberplot.py | 11 +- py-ft/python/pyft/tests/test_fiberdata.py | 77 ++++ py-ft/python/pyft/utils.py | 58 +-- py-ft/src/fiberdata.rs | 498 ---------------------- py-ft/src/lib.rs | 25 -- 11 files changed, 346 insertions(+), 731 deletions(-) create mode 100644 .github/workflows/release-python-pyft.yml delete mode 100644 py-ft/.github/workflows/CI.yml delete mode 100644 py-ft/Cargo.toml create mode 100644 py-ft/python/pyft/fiberdata.py create mode 100644 py-ft/python/pyft/tests/test_fiberdata.py delete mode 100644 py-ft/src/fiberdata.rs delete mode 100644 py-ft/src/lib.rs diff --git a/.github/workflows/release-python-pyft.yml b/.github/workflows/release-python-pyft.yml new file mode 100644 index 00000000..c6bd9c3e --- /dev/null +++ b/.github/workflows/release-python-pyft.yml @@ -0,0 +1,40 @@ +# Build and publish py-ft to PyPI. +# +# Triggered by a tag like `py-ft-v0.4.0`. Since py-ft is now a pure-Python +# package (over pysam + molecular-annotation), this is a single universal +# wheel + sdist — no per-platform matrix. Bump the version in +# py-ft/pyproject.toml, then push the tag. Publish pyMA first: py-ft depends on +# `molecular-annotation` from PyPI. +# +# Uses PyPI trusted publishing (OIDC) — no API token. Requires a matching +# trusted publisher on PyPI for project `pyft` (owner/repo +# fiberseq/fibertools-rs, workflow release-python-pyft.yml, no environment). +name: release py-ft +on: + push: + tags: ["py-ft-v*"] + workflow_dispatch: + +permissions: + contents: read + +jobs: + build-publish: + runs-on: ubuntu-latest + permissions: + id-token: write + contents: read + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: Build sdist + wheel + run: | + python -m pip install --upgrade build + python -m build py-ft --outdir dist + - name: Publish to PyPI (trusted publishing) + if: startsWith(github.ref, 'refs/tags/') + uses: pypa/gh-action-pypi-publish@release/v1 + with: + skip-existing: true diff --git a/py-ft/.github/workflows/CI.yml b/py-ft/.github/workflows/CI.yml deleted file mode 100644 index 8268eaf6..00000000 --- a/py-ft/.github/workflows/CI.yml +++ /dev/null @@ -1,120 +0,0 @@ -# This file is autogenerated by maturin v1.1.0 -# To update, run -# -# maturin generate-ci github -# -name: CI - -on: - push: - branches: - - main - - master - tags: - - '*' - pull_request: - workflow_dispatch: - -permissions: - contents: read - -jobs: - linux: - runs-on: ubuntu-latest - strategy: - matrix: - target: [x86_64, x86, aarch64, armv7, s390x, ppc64le] - steps: - - uses: actions/checkout@v3 - - uses: actions/setup-python@v4 - with: - python-version: '3.10' - - name: Build wheels - uses: PyO3/maturin-action@v1 - with: - target: ${{ matrix.target }} - args: --release --out dist --find-interpreter - sccache: 'true' - manylinux: auto - - name: Upload wheels - uses: actions/upload-artifact@v3 - with: - name: wheels - path: dist - - windows: - runs-on: windows-latest - strategy: - matrix: - target: [x64, x86] - steps: - - uses: actions/checkout@v3 - - uses: actions/setup-python@v4 - with: - python-version: '3.10' - architecture: ${{ matrix.target }} - - name: Build wheels - uses: PyO3/maturin-action@v1 - with: - target: ${{ matrix.target }} - args: --release --out dist --find-interpreter - sccache: 'true' - - name: Upload wheels - uses: actions/upload-artifact@v3 - with: - name: wheels - path: dist - - macos: - runs-on: macos-latest - strategy: - matrix: - target: [x86_64, aarch64] - steps: - - uses: actions/checkout@v3 - - uses: actions/setup-python@v4 - with: - python-version: '3.10' - - name: Build wheels - uses: PyO3/maturin-action@v1 - with: - target: ${{ matrix.target }} - args: --release --out dist --find-interpreter - sccache: 'true' - - name: Upload wheels - uses: actions/upload-artifact@v3 - with: - name: wheels - path: dist - - sdist: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v3 - - name: Build sdist - uses: PyO3/maturin-action@v1 - with: - command: sdist - args: --out dist - - name: Upload sdist - uses: actions/upload-artifact@v3 - with: - name: wheels - path: dist - - release: - name: Release - runs-on: ubuntu-latest - if: "startsWith(github.ref, 'refs/tags/')" - needs: [linux, windows, macos, sdist] - steps: - - uses: actions/download-artifact@v3 - with: - name: wheels - - name: Publish to PyPI - uses: PyO3/maturin-action@v1 - env: - MATURIN_PYPI_TOKEN: ${{ secrets.PYPI_API_TOKEN }} - with: - command: upload - args: --skip-existing * diff --git a/py-ft/Cargo.toml b/py-ft/Cargo.toml deleted file mode 100644 index 63807984..00000000 --- a/py-ft/Cargo.toml +++ /dev/null @@ -1,21 +0,0 @@ -[package] -edition = "2021" -name = "pyft" -version = "0.3.0" -readme = "./README.md" - -[lib] -crate-type = ["cdylib"] -name = "pyft" - -[dependencies] -anyhow = "1.0.58" -colored = "2.0.0" -env_logger = "0.9.0" -fibertools-rs = { path = "../" } -itertools = "0.10.5" -lazy_static = "1.4.0" -log = "0.4" -pyo3 = "0.19.1" -regex = "1.9.1" -rust-htslib = "0.46" diff --git a/py-ft/pyproject.toml b/py-ft/pyproject.toml index 37fba7bc..9d80e6a4 100644 --- a/py-ft/pyproject.toml +++ b/py-ft/pyproject.toml @@ -1,31 +1,72 @@ [build-system] -build-backend = "maturin" -requires = ["maturin>=1.1,<2.0"] +requires = ["hatchling"] +build-backend = "hatchling.build" [project] -dynamic = ["version"] +name = "pyft" +version = "0.4.0" +description = "Inspection and visualization toolkit for fiber-seq BAMs" +readme = "README.md" +requires-python = ">=3.9" +license = { text = "MIT" } +keywords = ["bioinformatics", "genomics", "fiber-seq", "BAM", "pysam"] classifiers = [ - "Programming Language :: Rust", - "Programming Language :: Python :: Implementation :: CPython", - "Programming Language :: Python :: Implementation :: PyPy", + "Programming Language :: Python :: 3", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Bio-Informatics", ] -name = "pyft" -requires-python = ">=3.7" dependencies = [ + "pysam>=0.20", + "molecular-annotation>=0.1", +] + +[project.optional-dependencies] +# Plotting layer (pyft.plot / pyft.utils, the fiberplot CLI). Kept optional so +# the inspection API stays a light install. +viz = [ "pandas>=2.0", "altair==5.5", "vegafusion[embed]", - "tqdm", - "typing-extensions", "polars>=1.0,<2.0", + "tqdm", ] +dev = ["pytest>=7"] [project.scripts] fiberplot = "pyft.fiberplot:main" -[tool.maturin] -features = ["pyo3/extension-module"] -python-source = "python" +[tool.hatch.build.targets.wheel] +packages = ["python/pyft"] + +[tool.hatch.build.targets.wheel.force-include] +# nothing extra; tests are excluded from the wheel by default + +[tool.pytest.ini_options] +testpaths = ["python/pyft/tests"] [tool.ruff] extend-include = ["*.ipynb"] + +# --------------------------------------------------------------------------- +# Pixi dev environment (https://pixi.sh). Run from `py-ft/`: +# pixi run test +# molecular-annotation is built editable from the sibling crate so py-ft always +# tests against the local pyMA (mirrors the Rust path dependency). +# --------------------------------------------------------------------------- +[tool.pixi.workspace] +name = "pyft" +channels = ["conda-forge", "bioconda"] +platforms = ["osx-arm64", "osx-64", "linux-64"] + +[tool.pixi.dependencies] +python = ">=3.10,<3.13" +pip = "*" +pytest = ">=7" +pysam = ">=0.20" + +[tool.pixi.pypi-dependencies] +pyft = { path = ".", editable = true } +molecular-annotation = { path = "../molecular-annotation/python", editable = true } + +[tool.pixi.tasks] +test = "pytest python/pyft/tests -v" diff --git a/py-ft/python/pyft/__init__.py b/py-ft/python/pyft/__init__.py index 7e062abb..c2fd6c50 100644 --- a/py-ft/python/pyft/__init__.py +++ b/py-ft/python/pyft/__init__.py @@ -1,8 +1,10 @@ -# this imports all the rust functions -from .pyft import * -from . import utils -from . import plot +"""pyft: an inspection and visualization toolkit for fiber-seq BAMs. -__doc__ = pyft.__doc__ -if hasattr(pyft, "__all__"): - __all__ = pyft.__all__ +Reads fiber-seq annotations (nuc/msp/fire and m6A/5mC base mods) via pysam and +the molecular_annotation library. The plotting helpers live in `pyft.plot` and +`pyft.utils` and are imported on demand (they pull in pandas/altair). +""" + +from .fiberdata import Feature, Fiberdata, fetch + +__all__ = ["Feature", "Fiberdata", "fetch"] diff --git a/py-ft/python/pyft/fiberdata.py b/py-ft/python/pyft/fiberdata.py new file mode 100644 index 00000000..1bbaad4b --- /dev/null +++ b/py-ft/python/pyft/fiberdata.py @@ -0,0 +1,144 @@ +"""A thin fiber-seq read model over pysam + the molecular_annotation library. + +`Fiberdata` wraps a single BAM record's molecular annotations and exposes the +fiber-seq vocabulary (m6a, 5mC, nuc, msp, fire) on top of the generic +`molecular_annotation` types. + +Coordinates are 0-based half-open. `Feature.start/end` are in molecular +(forward) orientation; `Feature.ref_start/ref_end` are lifted to the reference +(None where the annotation falls outside an aligned block). +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Optional + +import pysam +from molecular_annotation import from_record as _ma_from_record + +if TYPE_CHECKING: + from collections.abc import Iterator + +# MSP precision at or above which an MSP is also a FIRE element. Only used when +# the BAM does not already carry a distinct "fire" annotation type. +FIRE_MIN_QUAL = 200 + + +@dataclass(frozen=True) +class Feature: + """A single annotation. Coordinates are 0-based half-open.""" + + start: int + end: int + ref_start: Optional[int] + ref_end: Optional[int] + qual: Optional[int] + + +class Fiberdata: + """A single fiber-seq read: record metadata plus typed annotation accessors.""" + + def __init__(self, record: "pysam.AlignedSegment", annotations): + self._rec = record + self._annot = annotations + + @classmethod + def from_record(cls, record: "pysam.AlignedSegment") -> "Fiberdata": + return cls(record, _ma_from_record(record)) + + # --- read metadata, straight off the pysam record --- + + @property + def name(self) -> str: + return self._rec.query_name + + @property + def strand(self) -> str: + return "-" if self._rec.is_reverse else "+" + + @property + def is_mapped(self) -> bool: + return not self._rec.is_unmapped + + @property + def chrom(self) -> Optional[str]: + return None if self._rec.is_unmapped else self._rec.reference_name + + @property + def start(self) -> Optional[int]: + return None if self._rec.is_unmapped else self._rec.reference_start + + @property + def end(self) -> Optional[int]: + return None if self._rec.is_unmapped else self._rec.reference_end + + @property + def length(self) -> int: + return self._annot.read_length + + def tag(self, name: str, default=None): + try: + return self._rec.get_tag(name) + except KeyError: + return default + + # --- typed accessors: fiber-seq names -> molecular_annotation type strings --- + + def _features(self, type_name: str) -> list[Feature]: + items = self._annot.iter_type(type_name) + if not items: + return [] + # iter_type tuple: (q_start, q_end, f_start, f_end, ref_start, ref_end, quals, name) + return [ + Feature(f_start, f_end, ref_start, ref_end, quals[0] if quals else None) + for _qs, _qe, f_start, f_end, ref_start, ref_end, quals, _name in items + ] + + @property + def m6a(self) -> list[Feature]: + return self._features("a") + + @property + def cpg(self) -> list[Feature]: # 5mC + return self._features("m") + + @property + def nuc(self) -> list[Feature]: + return self._features("nuc") + + @property + def msp(self) -> list[Feature]: + return self._features("msp") + + @property + def fire(self) -> list[Feature]: + # FIRE is either its own annotation type, or the high-precision MSP subset. + if "fire" in self._annot.annotation_type_names(): + return self._features("fire") + return [f for f in self.msp if f.qual is not None and f.qual >= FIRE_MIN_QUAL] + + def __repr__(self) -> str: + return ( + f"Fiberdata({self.name} {self.chrom}:{self.start}-{self.end} {self.strand} " + f"m6a={len(self.m6a)} cpg={len(self.cpg)} nuc={len(self.nuc)} msp={len(self.msp)})" + ) + + +def fetch(bam_path, region=None) -> "Iterator[Fiberdata]": + """Yield `Fiberdata` over a BAM file. + + Without `region`, iterates every read in file order -- mapped **and** + unmapped -- so an unaligned fiber-seq BAM (raw PacBio/ONT output, before + alignment) works too; filter with `Fiberdata.is_mapped` if you only want + one. With `region`, uses pysam's indexed fetch, which returns the mapped + reads overlapping that region. + + Args: + bam_path: Path to a BAM/CRAM file (indexed if `region` is given). + region: Optional (chrom, start, end) tuple passed to `pysam.fetch`. + """ + with pysam.AlignmentFile(str(bam_path), "rb", check_sq=False) as bam: + source = bam.fetch(*region) if region is not None else bam.fetch(until_eof=True) + for record in source: + yield Fiberdata.from_record(record) diff --git a/py-ft/python/pyft/fiberplot.py b/py-ft/python/pyft/fiberplot.py index 83b8cfe5..550f2d3d 100644 --- a/py-ft/python/pyft/fiberplot.py +++ b/py-ft/python/pyft/fiberplot.py @@ -88,17 +88,14 @@ def create_fiberplot( if output_file is None: output_file = "fiberplot.html" - # Load the BAM file - print(f"Loading BAM file: {bam_file}", file=sys.stderr) - fiberbam = pyft.Fiberbam(str(bam_file)) - - # Fetch and convert to dataframe + # Fetch and convert to dataframe. NOTE: region_to_* are being ported to the + # pysam + molecular_annotation backend and currently raise NotImplementedError. print(f"Fetching region: {region[0]}:{region[1]}-{region[2]}", file=sys.stderr) if centered: print(f"Creating centered plot (strand={strand})", file=sys.stderr) df = utils.region_to_centered_df( - fiberbam, + bam_file, region, strand=strand, max_flank=max_flank, @@ -107,7 +104,7 @@ def create_fiberplot( chart = plot.centered_chart(df, width=width, height=height) else: print("Creating region plot", file=sys.stderr) - df = utils.region_to_df(fiberbam, region, min_basemod_qual=min_basemod_qual) + df = utils.region_to_df(bam_file, region, min_basemod_qual=min_basemod_qual) # Add centering columns required by centered_chart df["centering_position"] = region[1] df["centering_strand"] = strand diff --git a/py-ft/python/pyft/tests/test_fiberdata.py b/py-ft/python/pyft/tests/test_fiberdata.py new file mode 100644 index 00000000..43f80fa3 --- /dev/null +++ b/py-ft/python/pyft/tests/test_fiberdata.py @@ -0,0 +1,77 @@ +"""Tests for pyft.fiberdata.""" + +from pathlib import Path + +import pytest + +from pyft.fiberdata import Feature, Fiberdata, fetch + +# Reuse the fibertools-rs repo fixture (all.bam); absent outside the repo. +_BAM = Path(__file__).parents[4] / "tests" / "data" / "all.bam" + + +class TestRealFixture: + """Against a real fiber-seq BAM (all 22 reads are mapped).""" + + pytestmark = pytest.mark.skipif( + not _BAM.exists(), + reason="fiber-seq BAM fixture only present inside the fibertools-rs repo", + ) + + def _first(self) -> Fiberdata: + for fd in fetch(_BAM): + return fd + raise AssertionError("no records in fixture") + + def test_fetch_yields_fibers(self): + assert sum(1 for _ in fetch(_BAM)) == 22 + + def test_first_fiber_metadata_and_counts(self): + fd = self._first() + # Known values for m54329U_210323_190418/5048829/ccs (reverse-aligned). + assert fd.name == "m54329U_210323_190418/5048829/ccs" + assert fd.is_mapped + assert fd.chrom == "ptg000001l" + assert fd.strand == "-" + assert fd.length == 15524 + assert len(fd.m6a) == 1541 + assert len(fd.cpg) == 135 + assert len(fd.nuc) == 76 + assert len(fd.msp) == 77 + + def test_features_are_typed_and_lift_to_reference(self): + fd = self._first() + m6a = fd.m6a + assert all(isinstance(f, Feature) for f in m6a) + assert all(0 <= f.start < fd.length for f in m6a) + assert len([f for f in m6a if f.ref_start is not None]) > 0 + assert all(0 <= f.qual <= 255 for f in m6a) + + +class TestUnmappedRead: + """Fiberdata works on unmapped reads (e.g. unaligned fiber-seq BAMs).""" + + def _unmapped_record(self, pysam): + import array + + header = pysam.AlignmentHeader.from_dict({"HD": {"VN": "1.0"}}) + rec = pysam.AlignedSegment(header) + rec.query_name = "unaligned_read" + rec.query_sequence = "ACAGAA" # A at 0,2,4,5 + rec.flag = 4 # unmapped + rec.set_tag("MA", "6") # read length only, no structural annotations + rec.set_tag("MM", "A+a,1,0,0;") # m6A at 2,4,5 + rec.set_tag("ML", array.array("B", [200, 150, 100])) + return rec + + def test_unmapped_read_parses_without_reference(self): + pysam = pytest.importorskip("pysam") + fd = Fiberdata.from_record(self._unmapped_record(pysam)) + + assert fd.is_mapped is False + assert fd.chrom is None and fd.start is None and fd.end is None + assert fd.strand == "+" + # base mods still decode from the molecular sequence... + assert [f.start for f in fd.m6a] == [2, 4, 5] + # ...but there is nothing to lift to, so reference coords are None. + assert all(f.ref_start is None for f in fd.m6a) diff --git a/py-ft/python/pyft/utils.py b/py-ft/python/pyft/utils.py index 009b82d3..769efbef 100644 --- a/py-ft/python/pyft/utils.py +++ b/py-ft/python/pyft/utils.py @@ -329,45 +329,23 @@ def base_mod_end(list_of_starts): data_dict["qual"].append(fiber.cpg.ml) -def region_to_centered_df( - fiberbam, region, strand="+", max_flank=None, min_basemod_qual=125 -): - """ - Takes a fiberbam and a region and returns a pandas dataframe with reference centered positions - """ - data_dict = empty_data_dict() - for fiber in fiberbam.center( - region[0], start=region[1], end=region[2], strand=strand - ): - _add_fiber_to_data_dict(fiber, data_dict) - df = pd.DataFrame.from_dict(data_dict).explode(WIDE_COLUMNS) - df["centering_position"] = region[1] - if strand == "-": - df["centering_position"] = region[2] - 1 - df["centering_strand"] = strand - - # trim the dataframe to only include fibers that overlap - if max_flank is not None: - df = df[(df.start < +max_flank) & (df.end > -max_flank)] - - is_basemod = df.type.isin(["m6a", "5mC"]) - df = df[~(is_basemod & (df.qual < min_basemod_qual))] - return df +# NOTE: region_to_centered_df / region_to_df built their dataframes from the +# now-removed Rust `Fiberbam`. They are stubbed while py-ft is ported to the +# pysam + molecular_annotation backend; the replacements will loop over +# `pyft.fetch(bam_path, region)` and `Fiberdata` accessors (see project notes). -def region_to_df(fiberbam, region, min_basemod_qual=125): - """ - Takes a fiberbam and a region and returns a pandas dataframe with fibers - that overlap the region - """ - data_dict = empty_data_dict() - for fiber in fiberbam.fetch( - region[0], - start=region[1], - end=region[2], - ): - _add_fiber_to_data_dict(fiber, data_dict) - df = pd.DataFrame.from_dict(data_dict).explode(WIDE_COLUMNS) - is_basemod = df.type.isin(["m6a", "5mC"]) - df = df[~(is_basemod & (df.qual < min_basemod_qual))] - return df +def region_to_centered_df(bam_path, region, strand="+", max_flank=None, min_basemod_qual=125): + """Reference-centered dataframe for a region (not yet ported).""" + raise NotImplementedError( + "region_to_centered_df is being ported to the pysam + molecular_annotation " + "backend; use pyft.fetch() for now." + ) + + +def region_to_df(bam_path, region, min_basemod_qual=125): + """Per-fiber dataframe for a region (not yet ported).""" + raise NotImplementedError( + "region_to_df is being ported to the pysam + molecular_annotation " + "backend; use pyft.fetch() for now." + ) diff --git a/py-ft/src/fiberdata.rs b/py-ft/src/fiberdata.rs deleted file mode 100644 index faeb56da..00000000 --- a/py-ft/src/fiberdata.rs +++ /dev/null @@ -1,498 +0,0 @@ -use fibertools_rs::utils::bamlift; -//use fibertools_rs::bio_io; -use fibertools_rs::fiber::FiberseqData; -use fibertools_rs::utils::input_bam::FiberFilters; -use pyo3::iter::IterNextOutput; -use pyo3::prelude::*; -use rust_htslib::{bam, bam::ext::BamRecordExtensions, bam::record::Aux, bam::Read}; -use std::collections::HashMap; -use std::time; -use std::vec::IntoIter; -#[pyclass] -/// Class for fiberseq data. This class corresponds to a single record in the bam file. -pub struct Fiberdata { - /// Number of c_s - #[pyo3(get, set)] - pub ec: i64, - /// Name of the read - #[pyo3(get, set)] - pub qname: String, - /// SAM flag - #[pyo3(get, set)] - pub sam_flag: u16, - /// Read quality - #[pyo3(get, set)] - pub rq: String, - /// Haplotype tag - #[pyo3(get, set)] - pub hp: String, - /// Read sequence - #[pyo3(get, set)] - pub seq: String, - /// Chromosome - #[pyo3(get, set)] - pub chrom: String, - /// Chromosome start - #[pyo3(get, set)] - pub start: i64, - /// Chromosome end - #[pyo3(get, set)] - pub end: i64, - /// Strand - #[pyo3(get, set)] - pub strand: char, - #[pyo3(get, set)] - /// Read group - pub rg: String, - /// m6a :class:`~pyft.Basemods` object - #[pyo3(get, set)] - pub m6a: Basemods, - /// cpg :class:`~pyft.Basemods` object - #[pyo3(get, set)] - pub cpg: Basemods, - /// :class:`~pyft.Ranges` object for msp features - #[pyo3(get, set)] - pub msp: Ranges, - /// :class:`~pyft.Ranges` object for nuc features - #[pyo3(get, set)] - pub nuc: Ranges, - /// Aligned block pairs from the bam record - #[pyo3(get)] - pub aligned_block_pairs: Vec<([i64; 2], [i64; 2])>, - /// offset for centering the fiber - #[pyo3(get, set)] - offset: Option, - /// record - fiber: FiberseqData, -} -#[pymethods] -impl Fiberdata { - pub fn __str__(&self) -> String { - format!( - "fiber: {}\t\ - chrom: {}\tstart: {}\tend {}\t\ - num m6a: {}\t num cpg: {}\t\ - num nuc: {}\t num msp: {}", - self.qname, - self.chrom, - self.start, - self.end, - self.m6a.starts.len(), - self.cpg.starts.len(), - self.nuc.starts.len(), - self.msp.starts.len() - ) - } - - /// return the length of the sequence - pub fn get_seq_length(&self) -> usize { - self.seq.len() - } - - /// liftover query (fiber) positions to the reference - pub fn lift_query_positions(&self, positions: Vec) -> Vec> { - bamlift::lift_query_positions(&self.aligned_block_pairs, &positions) - .expect("Positions must be sorted before calling liftover") - } - - /// liftover reference positions to the query (fiber) - pub fn lift_reference_positions(&self, positions: Vec) -> Vec> { - bamlift::lift_reference_positions(&self.aligned_block_pairs, &positions) - .expect("Positions must be sorted before calling liftover") - } -} - -impl Fiberdata { - fn new(fiber: FiberseqData, offset: Option) -> Self { - // PB features - let ec = fiber.ec.round() as i64; - - // record features - let qname = std::str::from_utf8(fiber.record.qname()).unwrap(); - let sam_flag = fiber.record.flags(); - let rq = match fiber.get_rq() { - Some(x) => format!("{}", x), - None => ".".to_string(), - }; - let hp = fiber.get_hp(); - let seq = String::from_utf8_lossy(&fiber.record.seq().as_bytes()).to_string(); - let (chrom, start, end, strand): (&str, i64, i64, char) = if fiber.record.is_unmapped() { - (".", 0, 0, '.') - } else { - let strand = if fiber.record.is_reverse() { '-' } else { '+' }; - ( - &fiber.target_name, - fiber.record.reference_start(), - fiber.record.reference_end(), - strand, - ) - }; - let rg = if let std::result::Result::Ok(Aux::String(f)) = fiber.record.aux(b"RG") { - log::trace!("{f}"); - f - } else { - "." - }; - - // fiberseq features - let m6a = Basemods::new( - fiber.m6a.starts().into_iter().map(Some).collect(), - fiber.m6a.reference_starts(), - fiber.m6a.qual(), - ); - let cpg = Basemods::new( - fiber.cpg.starts().into_iter().map(Some).collect(), - fiber.cpg.reference_starts(), - fiber.cpg.qual(), - ); - - let nuc = Ranges { - starts: fiber.nuc.starts().into_iter().map(Some).collect(), - ends: fiber.nuc.ends().into_iter().map(Some).collect(), - lengths: fiber.nuc.lengths().into_iter().map(Some).collect(), - qual: fiber.nuc.qual(), - reference_starts: fiber.nuc.reference_starts(), - reference_ends: fiber.nuc.reference_ends(), - reference_lengths: fiber.nuc.reference_lengths(), - }; - let msp = Ranges { - starts: fiber.msp.starts().into_iter().map(Some).collect(), - ends: fiber.msp.ends().into_iter().map(Some).collect(), - lengths: fiber.msp.lengths().into_iter().map(Some).collect(), - qual: fiber.msp.qual(), - reference_starts: fiber.msp.reference_starts(), - reference_ends: fiber.msp.reference_ends(), - reference_lengths: fiber.msp.reference_lengths(), - }; - - let aligned_block_pairs = fiber.record.aligned_block_pairs().collect(); - Self { - ec, - qname: qname.to_string(), - sam_flag, - rq, - hp, - seq, - chrom: chrom.to_string(), - start, - end, - strand, - rg: rg.to_string(), - m6a, - cpg, - msp, - nuc, - aligned_block_pairs, - offset, - fiber, - } - } -} - -#[pyclass] -pub struct Fiberwriter { - writer: bam::Writer, -} - -#[pymethods] -impl Fiberwriter { - #[new] - pub fn new(output_bam_path: &str, template_bam_path: &str) -> Self { - let header = bam::Header::from_template( - bam::Reader::from_path(template_bam_path) - .expect("unable to open template bam file") - .header(), - ); - let writer = bam::Writer::from_path(output_bam_path, &header, bam::Format::Bam) - .expect("unable to open output bam file"); - Self { writer } - } - /// Write a Fiberdata object to a bam file. - #[pyo3(signature = (fiberdata))] - pub fn write(&mut self, fiberdata: &Fiberdata) -> PyResult<()> { - self.writer - .write(&fiberdata.fiber.record) - .expect("unable to write record"); - pyo3::prelude::PyResult::Ok(()) - } -} - -#[pyclass] -/// Open a fiberseq bam file. Must have an index. -pub struct Fiberbam { - bam: bam::IndexedReader, - reader: bam::Reader, - header: bam::Header, - target_dict: HashMap, - start: time::Instant, -} - -/// pure rust functions for Fiberbam -impl Fiberbam { - fn _fetch_helper( - &mut self, - chrom: &str, - start: Option, - end: Option, - ) -> Vec { - let fetch_args = match (start, end) { - (Some(start), Some(end)) => bam::FetchDefinition::from((chrom, start, end)), - _ => bam::FetchDefinition::from(chrom), - }; - - let head_view = bam::HeaderView::from_header(&self.header); - - self.bam.fetch(fetch_args).expect("unable to fetch region"); - let records: Vec = self.bam.records().map(|r| r.unwrap()).collect(); - - log::info!( - "{} records fetched in {:.2}s", - records.len(), - self._time_from_last() - ); - let fiberdata = FiberseqData::from_records(records, &head_view, &FiberFilters::default()); - log::info!( - "Fiberdata made for {} records in {:.2}s", - fiberdata.len(), - self._time_from_last() - ); - fiberdata - } -} - -#[pymethods] -impl Fiberbam { - #[new] - #[pyo3(signature = (bam_file, threads = 8))] - fn new(bam_file: &str, threads: usize) -> Self { - let mut bam = - bam::IndexedReader::from_path(bam_file).expect("unable to open indexed bam file"); - bam.set_threads(threads).unwrap(); - - let mut reader = bam::Reader::from_path(bam_file).expect("unable to open bam file"); - reader.set_threads(threads).unwrap(); - - let header = bam::Header::from_template(bam.header()); - let head_view = bam::HeaderView::from_header(&header); - let target_dict = FiberseqData::dict_from_head_view(&head_view); - let start = time::Instant::now(); - Self { - bam, - reader, - header, - target_dict, - start, - } - } - - /// Returns an iterator over :class:`~pyft.Fiberdata` objects for the selected region. - /// Arguments for the region to fetch can be provided in multiple ways, e.g.: - /// - /// fiberdata.fetch("chr1", 100, 200) - /// - /// fiberdata.fetch("chr1", start=100, end=200) - /// - /// fiberdata.fetch("chr1:100-200") - /// - /// fiberdata.fetch("chr1") - #[pyo3(signature = (chrom, start = None, end=None))] - pub fn fetch(&mut self, chrom: &str, start: Option, end: Option) -> Fiberiter { - Fiberiter::build_fiberdata_iter(self._fetch_helper(chrom, start, end)) - } - - /// Returns an iterator over :class:`~pyft.Fiberdata` objects; however, the data is centered around the region that has been fetched (`should` work the same as **ft center**). - #[pyo3(signature = (chrom, start, end, strand = '+'))] - pub fn center(&mut self, chrom: &str, start: i64, end: i64, strand: char) -> Fiberiter { - let position = if strand == '-' { end - 1 } else { start }; - let center_position = fibertools_rs::subcommands::center::CenterPosition { - chrom: chrom.to_string(), - position, - strand, - }; - let fiberdata = self._fetch_helper(chrom, Some(position), Some(position + 1)); - let fiberdata: Vec = fiberdata - .into_iter() - .filter_map(|fiber| fiber.center(¢er_position)) - .collect(); - log::info!( - "Fiberdata centered for {} records in {:.2}s", - fiberdata.len(), - self._time_from_last() - ); - Fiberiter::build_fiberdata_iter(fiberdata) - } - - fn _time_from_last(&mut self) -> f64 { - let elapsed = self.start.elapsed().as_secs_f64(); - self.start = time::Instant::now(); - elapsed - } - - fn __next__(&mut self) -> IterNextOutput { - let data = self.reader.records().next(); - match data { - Some(record) => { - let record = record.unwrap(); - let target_name = self.target_dict.get(&record.tid()); - let fiber = FiberseqData::new(record, target_name, &FiberFilters::default()); - IterNextOutput::Yield(Fiberdata::new(fiber, None)) - } - None => IterNextOutput::Return("Ended"), - } - } - - fn __iter__(slf: PyRef<'_, Self>) -> PyRef<'_, Self> { - slf - } -} - -#[pyclass] -/// An iterator over :class:`~pyft.Fiberdata` objects. -pub struct Fiberiter { - fiberdata: IntoIter, - length: usize, -} - -/// rust only functions for Fiberiter -impl Fiberiter { - fn build_fiberdata_iter(fiberdata: Vec) -> Self { - let length = fiberdata.len(); - Self { - fiberdata: fiberdata.into_iter(), - length, - } - } -} - -#[pymethods] -impl Fiberiter { - fn __next__(&mut self) -> IterNextOutput { - let data = self.fiberdata.next(); - match data { - Some(fiber) => IterNextOutput::Yield(Fiberdata::new(fiber, None)), - None => IterNextOutput::Return("Ended"), - } - } - - fn __iter__(slf: PyRef<'_, Self>) -> PyRef<'_, Self> { - slf - } - - fn __len__(&self) -> usize { - self.length - } - - fn len(&self) -> usize { - self.length - } -} - -/// Class for describing base modifications within fiberseq data. -/// For example, m6a and 5mC (cpg) features. -#[pyclass] -#[derive(Clone)] -pub struct Basemods { - /// Basemod starts - #[pyo3(get)] - starts: Vec>, - /// Basemod reference starts - #[pyo3(get)] - reference_starts: Vec>, - /// Basemod ML - #[pyo3(get)] - ml: Vec, -} - -#[pymethods] -impl Basemods { - #[new] - pub fn new(starts: Vec>, reference_starts: Vec>, ml: Vec) -> Self { - Self { - starts, - reference_starts, - ml, - } - } - - /// return reference ends, start + 1 - pub fn get_reference_ends(&self) -> Vec> { - self.reference_starts - .iter() - .map(|x| x.as_ref().map(|x| x + 1)) - .collect() - } - - /// return ends, start + 1 - pub fn get_ends(&self) -> Vec> { - self.starts - .iter() - .map(|x| x.as_ref().map(|x| x + 1)) - .collect() - } -} - -/// Class for describing ranges within fiberseq data. -/// For example, nucleosomes and msp features. -#[pyclass] -#[derive(Clone)] -pub struct Ranges { - /// Range starts - #[pyo3(get, set)] - pub starts: Vec>, - /// Range ends - #[pyo3(get, set)] - pub ends: Vec>, - /// Range lengths - #[pyo3(get, set)] - pub lengths: Vec>, - /// quals - #[pyo3(get, set)] - pub qual: Vec, - /// Reference starts - #[pyo3(get, set)] - pub reference_starts: Vec>, - /// Reference ends - #[pyo3(get, set)] - pub reference_ends: Vec>, - /// Reference lengths - #[pyo3(get, set)] - pub reference_lengths: Vec>, -} - -#[pymethods] -impl Ranges { - #[new] - pub fn new( - starts: Vec>, - lengths: Vec>, - qual: Vec, - reference_starts: Vec>, - reference_lengths: Vec>, - ) -> Self { - let ends = starts - .iter() - .zip(lengths.iter()) - .map(|(s, l)| match (s, l) { - (Some(s), Some(l)) => Some(s + l), - _ => None, - }) - .collect(); - let reference_ends = reference_starts - .iter() - .zip(reference_lengths.iter()) - .map(|(s, l)| match (s, l) { - (Some(s), Some(l)) => Some(s + l), - _ => None, - }) - .collect(); - Self { - starts, - ends, - lengths, - qual, - reference_starts, - reference_ends, - reference_lengths, - } - } -} diff --git a/py-ft/src/lib.rs b/py-ft/src/lib.rs deleted file mode 100644 index c401fe09..00000000 --- a/py-ft/src/lib.rs +++ /dev/null @@ -1,25 +0,0 @@ -/// A module for using fibertools-rs -mod fiberdata; - -use env_logger::{Builder, Target}; -use log::LevelFilter; -use pyo3::prelude::*; - -/// A python module for using rust to access data from fiberseq BAM files. -#[pymodule] -fn pyft(_py: Python, m: &PyModule) -> PyResult<()> { - Builder::new() - .format_timestamp_secs() - .target(Target::Stderr) - .filter(None, LevelFilter::Info) - .init(); - - m.add("__version__", env!("CARGO_PKG_VERSION"))?; - m.add_class::()?; - //m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - Ok(()) -} From b2c1526922a52725cd80a9556fe3c9f85e0a32e4 Mon Sep 17 00:00:00 2001 From: Cade Mirchandani Date: Tue, 21 Jul 2026 15:55:25 -0700 Subject: [PATCH 2/2] docs: build py-ft docs on the pysam + pyMA API The docs still built the old Rust/maturin package (rust toolchain + compiling the crate), which is why the Read the Docs build failed. py-ft is pure Python now, so: - .readthedocs.yaml drops rust/cmake; installing the package just pulls pysam + molecular-annotation from PyPI so autodoc can import pyft. - conf.py builds without nbsphinx; api.rst documents pyft.fiberdata via autodoc (Napoleon for the Google-style docstrings). - add __version__ to the package (conf.py reads it). - replace the old Fiberbam notebooks with a quickstart vignette on fetch() and the Fiberdata accessors; refresh index.rst. --- .readthedocs.yaml | 8 +- py-ft/docs/api.rst | 13 +- py-ft/docs/conf.py | 14 +- py-ft/docs/index.rst | 43 +- py-ft/docs/requirements.txt | 7 - .../docs/vignettes/basic-usage-examples.ipynb | 484 ------------------ py-ft/docs/vignettes/centered-plot.ipynb | 442 ---------------- py-ft/docs/vignettes/index.rst | 5 +- py-ft/docs/vignettes/quickstart.rst | 94 ++++ py-ft/python/pyft/__init__.py | 9 +- 10 files changed, 149 insertions(+), 970 deletions(-) delete mode 100644 py-ft/docs/vignettes/basic-usage-examples.ipynb delete mode 100644 py-ft/docs/vignettes/centered-plot.ipynb create mode 100644 py-ft/docs/vignettes/quickstart.rst diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 0643bb76..81bc60fd 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -8,16 +8,14 @@ sphinx: build: os: "ubuntu-22.04" - apt_packages: - - cmake tools: python: "3.10" - rust: "latest" +# py-ft is now a pure-Python package (pysam + molecular-annotation), so the docs +# build no longer needs Rust/cmake. Installing the package pulls its runtime +# deps from PyPI so autodoc can import pyft. python: install: - requirements: py-ft/docs/requirements.txt - method: pip path: ./py-ft -#conda: -# environment: py-ft/docs/environment.yml diff --git a/py-ft/docs/api.rst b/py-ft/docs/api.rst index 1a063d81..f28963f7 100644 --- a/py-ft/docs/api.rst +++ b/py-ft/docs/api.rst @@ -1,16 +1,11 @@ - API Reference -================= +============= -.. automodule:: pyft - :members: - :undoc-members: - :show-inheritance: - :member-order: bysource +``pyft.fiberdata`` +------------------ -.. automodule:: utils +.. automodule:: pyft.fiberdata :members: :undoc-members: :show-inheritance: :member-order: bysource - diff --git a/py-ft/docs/conf.py b/py-ft/docs/conf.py index cd079375..b9f25059 100644 --- a/py-ft/docs/conf.py +++ b/py-ft/docs/conf.py @@ -6,10 +6,9 @@ import re import sys -# import sphinx -# import sphinx.ext.autosummary as autosummary -# sys.path.insert(0, os.path.abspath("../")) -sys.path.insert(0, os.path.abspath("../python/pyft")) +# Put the package source on the path so autodoc can import `pyft` even when the +# package isn't pip-installed (Read the Docs installs it; local builds may not). +sys.path.insert(0, os.path.abspath("../python")) import pyft # -- Project information ----------------------------------------------------- @@ -27,16 +26,13 @@ extensions = [ "sphinx.ext.autodoc", - # "sphinx_autodoc_typehints", + "sphinx.ext.napoleon", "sphinx.ext.viewcode", - # "sphinx.ext.napoleon", "sphinx_rtd_theme", "sphinx.ext.intersphinx", - # "edit_on_github", - "nbsphinx", ] -source_suffix = [".rst", ".py"] +source_suffix = [".rst"] templates_path = ["_templates"] exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] diff --git a/py-ft/docs/index.rst b/py-ft/docs/index.rst index f4af18b0..ec289775 100644 --- a/py-ft/docs/index.rst +++ b/py-ft/docs/index.rst @@ -1,8 +1,8 @@ Getting started ---------------- +=============== -pyft: python bindings for fibertools-rs ---------------------------------------- +pyft: a fiber-seq inspection toolkit for Python +----------------------------------------------- .. image:: https://readthedocs.org/projects/py-ft/badge/?version=latest :target: https://py-ft.readthedocs.io/en/latest/?badge=latest @@ -12,7 +12,14 @@ pyft: python bindings for fibertools-rs :target: https://badge.fury.io/py/pyft -**pyft** provides a python API for the rust library `fibertools-rs `_. The inspiration for this API is to make analysis in python easier and faster; therefore, only extraction of data from a fiberseq bam is supported and not writing. +**pyft** is a lightweight toolkit for inspecting and visualizing `fiber-seq +`_ BAMs from Python. It reads a +record with `pysam `_ and decodes its fiber-seq +annotations through the `molecular-annotation +`_ library, then exposes them +under the fiber-seq vocabulary: ``m6a``, ``cpg`` (5mC), ``nuc``, ``msp`` and +``fire``. It is read/inspection-focused — for producing or editing fiber-seq +annotations, use `fibertools-rs `_. Install @@ -20,14 +27,32 @@ Install .. code-block:: bash pip install pyft - + +The plotting helpers (``pyft.plot`` / ``pyft.utils`` and the ``fiberplot`` CLI) +require extra dependencies: + +.. code-block:: bash + + pip install "pyft[viz]" + + +Quick example +============= +.. code-block:: python + + from pyft import fetch + + for fiber in fetch("fiberseq.bam"): + print(fiber.name, fiber.strand, "m6A:", len(fiber.m6a), "MSPs:", len(fiber.msp)) + Vignettes ========= -The `vignettes `_ are a good place to start to understand the capabilities of pyft. +The :doc:`vignettes ` are a good place to start. + Indices and tables -================== +=================== * :ref:`genindex` * :ref:`modindex` @@ -38,9 +63,7 @@ Indices and tables :hidden: :maxdepth: 2 :caption: pyft - + self api.rst vignettes/index.rst - - diff --git a/py-ft/docs/requirements.txt b/py-ft/docs/requirements.txt index f1dc4b42..d64916c5 100644 --- a/py-ft/docs/requirements.txt +++ b/py-ft/docs/requirements.txt @@ -1,11 +1,4 @@ -sphinx_bootstrap_theme sphinx<7.0.0 -sphinx-autodoc-typehints -#==1.23.0 sphinx-rtd-theme==1.2.2 -lxml-html-clean docutils==0.18.1 setuptools -nbsphinx -pypandoc_binary -#sphinxawesome-theme diff --git a/py-ft/docs/vignettes/basic-usage-examples.ipynb b/py-ft/docs/vignettes/basic-usage-examples.ipynb deleted file mode 100644 index e2f6c526..00000000 --- a/py-ft/docs/vignettes/basic-usage-examples.ipynb +++ /dev/null @@ -1,484 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Basic use cases for pyft" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] 6 records fetched in 0.00s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] Fiberdata made for 6 records in 0.02s\n", - "100%|██████████| 6/6 [00:00<00:00, 509.30it/s]" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " chrom motif_start motif_end strand n_spanning_fibers n_spanning_msps \\\n", - "0 chr11 5204946 5204981 + 181 92 \n", - "0 chr11 5204946 5204981 + 181 92 \n", - "0 chr11 5204946 5204981 + 181 92 \n", - "0 chr11 5204946 5204981 + 181 92 \n", - "0 chr11 5204946 5204981 + 181 92 \n", - ".. ... ... ... ... ... ... \n", - "16 chr19 45817350 45817385 + 136 124 \n", - "16 chr19 45817350 45817385 + 136 124 \n", - "16 chr19 45817350 45817385 + 136 124 \n", - "16 chr19 45817350 45817385 + 136 124 \n", - "16 chr19 45817350 45817385 + 136 124 \n", - "\n", - " n_overlapping_nucs module:0-8 module:8-16 module:16-23 module:23-29 \\\n", - "0 85 False False False False \n", - "0 85 False False False False \n", - "0 85 False False True True \n", - "0 85 False False False False \n", - "0 85 False False False False \n", - ".. ... ... ... ... ... \n", - "16 8 False True True True \n", - "16 8 False True True True \n", - "16 8 False False False True \n", - "16 8 True False True True \n", - "16 8 False False False False \n", - "\n", - " module:29-35 fire_qual fiber_name n_modules \\\n", - "0 False 247 m64076_211222_124721/148505307/ccs 5 \n", - "0 False -1 m64076_211222_124721/51053256/ccs 5 \n", - "0 False 246 m64076_211222_124721/62391018/ccs 5 \n", - "0 False -1 m64076_211222_124721/97191992/ccs 5 \n", - "0 False -1 m64076_211222_124721/99419016/ccs 5 \n", - ".. ... ... ... ... \n", - "16 True 0 m64076_211222_124721/157222001/ccs 5 \n", - "16 True 246 m64076_211222_124721/65339699/ccs 5 \n", - "16 False 0 m64076_211222_124721/6882497/ccs 5 \n", - "16 False 243 m64076_211222_124721/31394454/ccs 5 \n", - "16 False -1 m64076_211222_124721/100926481/ccs 5 \n", - "\n", - " has_spanning_msp \n", - "0 True \n", - "0 False \n", - "0 True \n", - "0 False \n", - "0 False \n", - ".. ... \n", - "16 True \n", - "16 True \n", - "16 True \n", - "16 True \n", - "16 False \n", - "\n", - "[2065 rows x 16 columns]\n", - " chrom motif_start motif_end strand fire_qual \\\n", - "0 chr11 5204946 5204981 + 247 \n", - "1 chr11 5204946 5204981 + -1 \n", - "2 chr11 5204946 5204981 + 246 \n", - "3 chr11 5204946 5204981 + -1 \n", - "4 chr11 5204946 5204981 + -1 \n", - "... ... ... ... ... ... \n", - "10320 chr19 45817350 45817385 + 0 \n", - "10321 chr19 45817350 45817385 + 246 \n", - "10322 chr19 45817350 45817385 + 0 \n", - "10323 chr19 45817350 45817385 + 243 \n", - "10324 chr19 45817350 45817385 + -1 \n", - "\n", - " fiber_name has_spanning_msp footprinted \\\n", - "0 m64076_211222_124721/148505307/ccs True False \n", - "1 m64076_211222_124721/51053256/ccs False False \n", - "2 m64076_211222_124721/62391018/ccs True False \n", - "3 m64076_211222_124721/97191992/ccs False False \n", - "4 m64076_211222_124721/99419016/ccs False False \n", - "... ... ... ... \n", - "10320 m64076_211222_124721/157222001/ccs True True \n", - "10321 m64076_211222_124721/65339699/ccs True True \n", - "10322 m64076_211222_124721/6882497/ccs True False \n", - "10323 m64076_211222_124721/31394454/ccs True False \n", - "10324 m64076_211222_124721/100926481/ccs False False \n", - "\n", - " start end centering_position centering_strand type \n", - "0 0 8 5204946 + not-footprinted \n", - "1 0 8 5204946 + not-footprinted \n", - "2 0 8 5204946 + not-footprinted \n", - "3 0 8 5204946 + not-footprinted \n", - "4 0 8 5204946 + not-footprinted \n", - "... ... ... ... ... ... \n", - "10320 29 35 45817350 + footprinted \n", - "10321 29 35 45817350 + footprinted \n", - "10322 29 35 45817350 + not-footprinted \n", - "10323 29 35 45817350 + not-footprinted \n", - "10324 29 35 45817350 + not-footprinted \n", - "\n", - "[10325 rows x 13 columns]\n", - " chrom fiber_start fiber_end fiber_name strand \\\n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - ".. ... ... ... ... ... \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "\n", - " type start end qual \n", - "0 msp 26333672 26333727 0 \n", - "0 msp 26333848 26333890 0 \n", - "0 msp 26334056 26334094 0 \n", - "0 msp 26334254 26334319 0 \n", - "0 msp 26334561 26334565 0 \n", - ".. ... ... ... ... \n", - "23 5mC 26365739 26365740 213 \n", - "23 5mC 26366886 26366887 255 \n", - "23 5mC 26367221 26367222 172 \n", - "23 5mC 26367226 26367227 246 \n", - "23 5mC 26367254 26367255 252 \n", - "\n", - "[9111 rows x 9 columns]\n", - " chrom fiber_start fiber_end fiber_name strand \\\n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - "0 chr22 26333471 26371209 m64076_210328_012155/35587949/ccs + \n", - ".. ... ... ... ... ... \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "23 chr22 26354168 26367283 m54329U_210326_192251/160237619/ccs - \n", - "\n", - " type start end qual centering_position centering_strand \n", - "0 msp -16802 -16734 0 26354169 - \n", - "0 msp -16489 -16476 0 26354169 - \n", - "0 msp -16230 -16185 0 26354169 - \n", - "0 msp -16044 -16004 0 26354169 - \n", - "0 msp -15865 -15810 0 26354169 - \n", - ".. ... ... ... ... ... ... \n", - "23 5mC -708 -707 173 26354169 - \n", - "23 5mC -667 -666 224 26354169 - \n", - "23 5mC -591 -590 135 26354169 - \n", - "23 5mC -95 -94 228 26354169 - \n", - "23 5mC -61 -60 178 26354169 - \n", - "\n", - "[9111 rows x 11 columns]\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] 6 records fetched in 0.02s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] Fiberdata made for 6 records in 0.02s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] Fiberdata centered for 6 records in 0.00s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] 6 records fetched in 0.04s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] Fiberdata made for 6 records in 0.02s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] 6 records fetched in 0.01s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] Fiberdata made for 6 records in 0.02s\n", - "[2024-05-28T23:34:23Z INFO pyft::fiberdata] Fiberdata centered for 6 records in 0.00s\n" - ] - } - ], - "source": [ - "# import pyft\n", - "# from pyft import pyft\n", - "import pyft\n", - "import tqdm\n", - "\n", - "bam_f = \"../../../tests/data/center.bam\"\n", - "fiberbam = pyft.Fiberbam(bam_f)\n", - "out_fiberbam = pyft.Fiberwriter(\"test.bam\", bam_f)\n", - "rgn = [\"chr22\", 26_354_169, 26_354_170]\n", - "for fiber in tqdm.tqdm(fiberbam.fetch(*rgn)):\n", - " # the number of ccs passes\n", - " fiber.ec\n", - " # the mps start positions\n", - " fiber.msp.starts\n", - " # the fire quality scores of the MSPs\n", - " fiber.msp.qual\n", - " # print the nuc reference starts\n", - " fiber.nuc.reference_starts\n", - " # lift query (fiber) positions to reference positions\n", - " fiber.lift_query_positions([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])\n", - " # lift reference positions to query (fiber) positions\n", - " fiber.lift_reference_positions([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])\n", - "\n", - " out_fiberbam.write(fiber)\n", - "\n", - "\n", - "for fiber in fiberbam.center(rgn[0], start=rgn[1], end=rgn[2], strand=\"-\"):\n", - " # returns the same fiber object as above; however, all the positions have been modified to be relative to the region fetched\n", - " # print(fiber.msp.reference_starts)\n", - " continue\n", - "\n", - "\n", - "# example of reading in a footprinting table\n", - "df = pyft.utils.read_footprint_table(\n", - " \"../../../tests/data/ctcf-footprints.bed.gz\", long=True\n", - ")\n", - "print(df)\n", - "\n", - "# read in a footprinting table and center the positions\n", - "df = pyft.utils.read_and_center_footprint_table(\n", - " \"../../../tests/data/ctcf-footprints.bed.gz\"\n", - ")\n", - "print(df)\n", - "\n", - "# read a region of a fiberbam into a pandas dataframe\n", - "df = pyft.utils.region_to_df(fiberbam, rgn)\n", - "print(df)\n", - "\n", - "# read a region of a fiberbam into a pandas dataframe and center the positions\n", - "df = pyft.utils.region_to_centered_df(fiberbam, rgn, strand=\"-\")\n", - "print(df)" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "fiber: m54329U_210814_130637/103874956/ccs\tchrom: .\tstart: 5506049\tend 5532904\tnum m6a: 1908\t num cpg: 379\tnum nuc: 141\t num msp: 142\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n", - "hi\n" - ] - } - ], - "source": [ - "bam_f = \"../../../tests/data/center.bam\"\n", - "fiberbam = pyft.Fiberbam(bam_f)\n", - "\n", - "# iterate over a fiberbam one fiber at a time\n", - "for idx, fiber in enumerate(fiberbam):\n", - " if idx > 10:\n", - " break\n", - " print(fiber)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": ".venv", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.1" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/py-ft/docs/vignettes/centered-plot.ipynb b/py-ft/docs/vignettes/centered-plot.ipynb deleted file mode 100644 index ff3ce641..00000000 --- a/py-ft/docs/vignettes/centered-plot.ipynb +++ /dev/null @@ -1,442 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Plotting with `pyft`" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "DataTransformerRegistry.enable('vegafusion')" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "import pyft\n", - "import pandas as pd\n", - "import altair as alt\n", - "\n", - "alt.data_transformers.enable(\"vegafusion\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Read in the results of a ft-footprint calculation and plot the results using `pyft`." - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
chrommotif_startmotif_endstrandfootprint_codesfire_qualfiber_namehas_spanning_mspfootprintedstartendcentering_positioncentering_strandtype
0chr1152049465204981+3247m64076_211222_124721/148505307/ccsTrueTrue015204946+footprinted
1chr1152049465204981+2-1m64076_211222_124721/51053256/ccsFalseTrue015204946+not-footprinted
\n", - "
" - ], - "text/plain": [ - " chrom motif_start motif_end strand footprint_codes fire_qual \\\n", - "0 chr11 5204946 5204981 + 3 247 \n", - "1 chr11 5204946 5204981 + 2 -1 \n", - "\n", - " fiber_name has_spanning_msp footprinted start \\\n", - "0 m64076_211222_124721/148505307/ccs True True 0 \n", - "1 m64076_211222_124721/51053256/ccs False True 0 \n", - "\n", - " end centering_position centering_strand type \n", - "0 1 5204946 + footprinted \n", - "1 1 5204946 + not-footprinted " - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dfm = pyft.utils.read_and_center_footprint_table(\n", - " \"../../../tests/data/ctcf-footprints.bed.gz\"\n", - ")\n", - "dfm.head(2)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Read in fiber data centered on the footprint locations. " - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "[2024-11-12T22:29:48Z INFO pyft::fiberdata] 181 records fetched in 0.01s\n", - "[2024-11-12T22:29:48Z INFO pyft::fiberdata] Fiberdata made for 181 records in 0.11s\n", - "[2024-11-12T22:29:48Z INFO pyft::fiberdata] Fiberdata centered for 181 records in 0.02s\n", - "[2024-11-12T22:29:48Z INFO pyft::fiberdata] 172 records fetched in 0.11s\n", - "[2024-11-12T22:29:48Z INFO pyft::fiberdata] Fiberdata made for 172 records in 0.10s\n", - "[2024-11-12T22:29:48Z INFO pyft::fiberdata] Fiberdata centered for 172 records in 0.07s\n" - ] - } - ], - "source": [ - "rgns = pd.read_csv(\"../../../tests/data/ctcf.bed.gz\", sep=\"\\t\", header=None, nrows=2)\n", - "rgns.columns = [\"chrom\", \"start\", \"end\", \"name\", \"score\", \"strand\", \"name2\"]\n", - "fiberbam = pyft.Fiberbam(\"../../../tests/data/ctcf.bam\")\n", - "centers = []\n", - "z = None\n", - "for idx, rgn in rgns.iterrows():\n", - " region = (rgn[\"chrom\"], rgn[\"start\"], rgn[\"end\"])\n", - " z = pyft.utils.region_to_centered_df(\n", - " fiberbam, region, strand=rgn[\"strand\"], max_flank=250\n", - " )\n", - " centers.append(z)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Combine the footprinting results with the fiber data centered around the footprints. " - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
chromfiber_startfiber_endfiber_namestrandtypestartendqualcentering_positioncentering_strandmotif_startmotif_endfootprint_codesfire_qualhas_spanning_mspfootprinted
0chr115184260.05205600.0m64076_211222_124721/148505307/ccs+msp-225-16005204946+NaNNaNNaNNaNNaNNaN
1chr115184260.05205600.0m64076_211222_124721/148505307/ccs+msp-571352475204946+NaNNaNNaNNaNNaNNaN
\n", - "
" - ], - "text/plain": [ - " chrom fiber_start fiber_end fiber_name strand \\\n", - "0 chr11 5184260.0 5205600.0 m64076_211222_124721/148505307/ccs + \n", - "1 chr11 5184260.0 5205600.0 m64076_211222_124721/148505307/ccs + \n", - "\n", - " type start end qual centering_position centering_strand motif_start \\\n", - "0 msp -225 -160 0 5204946 + NaN \n", - "1 msp -57 135 247 5204946 + NaN \n", - "\n", - " motif_end footprint_codes fire_qual has_spanning_msp footprinted \n", - "0 NaN NaN NaN NaN NaN \n", - "1 NaN NaN NaN NaN NaN " - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "both_dfs = pd.concat(centers + [dfm], axis=0).reset_index(drop=True)\n", - "both_dfs.head(2)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Show the chart within the notebook." - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "\n", - "\n", - "
\n", - "" - ], - "text/plain": [ - "alt.LayerChart(...)" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "chart = pyft.plot.centered_chart(both_dfs, width=400, height=200)\n", - "chart" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Save the chart as a html file and open it in your browser.\n" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "chart.save(\"/Users/mrvollger/Desktop/chart.html\")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": ".env", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.13.0" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/py-ft/docs/vignettes/index.rst b/py-ft/docs/vignettes/index.rst index 087b191c..6b809960 100644 --- a/py-ft/docs/vignettes/index.rst +++ b/py-ft/docs/vignettes/index.rst @@ -1,11 +1,10 @@ Vignettes ========= -Here are some examples of `pyft` in action: +Examples of ``pyft`` in action: .. toctree:: :maxdepth: 1 :titlesonly: - basic-usage-examples.py - centered-plot.ipynb + quickstart diff --git a/py-ft/docs/vignettes/quickstart.rst b/py-ft/docs/vignettes/quickstart.rst new file mode 100644 index 00000000..e7c543d7 --- /dev/null +++ b/py-ft/docs/vignettes/quickstart.rst @@ -0,0 +1,94 @@ +Quickstart +========== + +``pyft`` turns each record of a fiber-seq BAM into a :class:`~pyft.Fiberdata` +object that exposes the fiber-seq features under familiar names. + +Iterating a BAM +--------------- + +:func:`~pyft.fetch` yields one :class:`~pyft.Fiberdata` per read. Without a +region it walks the whole file, including **unmapped** reads, so it works on raw +(unaligned) fiber-seq BAMs as well as aligned ones: + +.. code-block:: python + + from pyft import fetch + + for fiber in fetch("fiberseq.bam"): + print( + fiber.name, + f"{fiber.chrom}:{fiber.start}-{fiber.end}", + fiber.strand, + "m6A:", len(fiber.m6a), + "5mC:", len(fiber.cpg), + "nuc:", len(fiber.nuc), + "MSP:", len(fiber.msp), + ) + +To restrict to a region of an indexed BAM, pass a ``(chrom, start, end)`` tuple +(this returns the mapped reads overlapping the region): + +.. code-block:: python + + for fiber in fetch("fiberseq.bam", region=("chr1", 100_000, 101_000)): + ... + +Annotation accessors +-------------------- + +Each accessor returns a list of :class:`~pyft.Feature`. Coordinates are 0-based +half-open; ``start``/``end`` are in molecular (read) orientation, while +``ref_start``/``ref_end`` are lifted to the reference (``None`` where the +annotation falls outside an aligned block). ``qual`` is the per-feature quality +(the ML probability for base mods, MSP precision for ``msp``/``fire``): + +.. code-block:: python + + fiber = next(fetch("fiberseq.bam")) + + # m6A calls that lifted to the reference + m6a_ref = [f.ref_start for f in fiber.m6a if f.ref_start is not None] + + # nucleosomes as reference intervals + nucs = [(f.ref_start, f.ref_end) for f in fiber.nuc] + + # high-confidence m6A only + strong = [f for f in fiber.m6a if f.qual >= 200] + +The available accessors are ``m6a``, ``cpg`` (5mC), ``nuc``, ``msp`` and +``fire`` (a distinct type if present, otherwise the high-precision MSP subset). + +Mapped vs. unmapped reads +------------------------- + +Unmapped reads still carry their molecular annotations; only the reference +coordinates are unavailable. Filter with :attr:`~pyft.Fiberdata.is_mapped`: + +.. code-block:: python + + mapped = [f for f in fetch("fiberseq.bam") if f.is_mapped] + +Plotting +-------- + +The accessors feed directly into any plotting library. This draws a browser-style +row per fiber — nucleosomes as blocks, m6A as ticks — using reference +coordinates (requires ``matplotlib``): + +.. code-block:: python + + import matplotlib.pyplot as plt + from matplotlib.patches import Rectangle + from pyft import fetch + + fig, ax = plt.subplots(figsize=(12, 6)) + for y, fiber in enumerate(fetch("fiberseq.bam", region=("chr1", 100_000, 101_000))): + for f in fiber.nuc: + if f.ref_start is not None: + ax.add_patch(Rectangle((f.ref_start, y - 0.2), + f.ref_end - f.ref_start, 0.4, color="0.7")) + xs = [f.ref_start for f in fiber.m6a if f.ref_start is not None] + ax.plot(xs, [y] * len(xs), "|", color="#800080", markersize=5) + ax.set(xlabel="chr1 position", ylabel="fiber") + fig.savefig("fibers.png", dpi=150) diff --git a/py-ft/python/pyft/__init__.py b/py-ft/python/pyft/__init__.py index c2fd6c50..3924a79c 100644 --- a/py-ft/python/pyft/__init__.py +++ b/py-ft/python/pyft/__init__.py @@ -5,6 +5,13 @@ `pyft.utils` and are imported on demand (they pull in pandas/altair). """ +from importlib.metadata import PackageNotFoundError, version + from .fiberdata import Feature, Fiberdata, fetch -__all__ = ["Feature", "Fiberdata", "fetch"] +try: + __version__ = version("pyft") +except PackageNotFoundError: # not installed (e.g. running from a source tree) + __version__ = "0.0.0" + +__all__ = ["Feature", "Fiberdata", "fetch", "__version__"]