From 1788c4a611d93be6fb97be021df024b269b39a76 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 5 Oct 2026 08:58:59 +0000 Subject: [PATCH 1/4] Fix parser/CLI edge cases, align Python with R, strengthen CI Python: - read_fasta/read_fastq: a bare ">" or "@" header raised an uncaught IndexError (traceback from the CLI); now a clear ValueError. FASTQ readers tolerate blank lines between records. - CLI `gc` printed NA-less output and aborted the whole file on an all-N record; it now prints NA for that record and continues. - CLI `stats` on an empty FASTA reported "max() arg is an empty sequence"; now "no sequences found". - gc_content(ignore_n=True) now drops every non-ACGT symbol from the denominator, matching the R package (was N only). - kmer_count is case-insensitive and skips k-mers with ambiguous bases by default (skip_ambiguous=False restores the old counting). - __version__ is read from package metadata (single source of truth); version bumped to 0.2.0 for the behaviour changes. Ship py.typed. - Doctests run under pytest; ruff lint config added. R: - reverse_complement() accepts IUPAC ambiguity codes, like Python. - count_motif() errors on an empty motif instead of returning nchar(seq). CI: - Run on pull_request as well as push. - Python: test 3.10-3.13 (requires-python is >=3.10) and lint with ruff. - R: full R CMD check instead of devtools::test(), so stale Rd files, NAMESPACE problems and broken examples are caught. Docs: fix the placeholder install URLs in python/README.md, add CI badges and usage examples to the top-level README. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WqrAUWjyEmFz36zphqfoBs --- .github/workflows/python-tests.yml | 19 ++++++++++++++----- .github/workflows/r-tests.yml | 24 +++++++++++++++--------- README.md | 29 ++++++++++++++++++++++++++--- python/README.md | 11 ++++++----- python/pyproject.toml | 17 ++++++++++++++--- python/src/rushtools/__init__.py | 7 ++++++- python/src/rushtools/cli.py | 10 ++++++++-- python/src/rushtools/parsers.py | 12 ++++++++++-- python/src/rushtools/py.typed | 0 python/src/rushtools/seq.py | 27 ++++++++++++++++++++------- python/tests/test_cli.py | 30 ++++++++++++++++++++++++++++++ python/tests/test_parsers.py | 28 ++++++++++++++++++++++++++++ python/tests/test_seq.py | 11 +++++++++++ r/NEWS.md | 7 +++++++ r/R/sequences.R | 12 ++++++++---- r/man/reverse_complement.Rd | 5 +++-- r/tests/testthat/test-motif.R | 5 +++++ r/tests/testthat/test-sequences.R | 5 +++++ 18 files changed, 216 insertions(+), 43 deletions(-) create mode 100644 python/src/rushtools/py.typed diff --git a/.github/workflows/python-tests.yml b/.github/workflows/python-tests.yml index 68aa6fe..f5a6bad 100644 --- a/.github/workflows/python-tests.yml +++ b/.github/workflows/python-tests.yml @@ -1,10 +1,16 @@ name: Python tests -on: push +on: + push: + pull_request: jobs: test: runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.10", "3.11", "3.12", "3.13"] defaults: run: working-directory: python @@ -13,12 +19,15 @@ jobs: - uses: actions/setup-python@v5 with: - python-version: "3.12" + python-version: ${{ matrix.python-version }} - - name: Install package and pytest + - name: Install package with dev extras run: | python -m pip install --upgrade pip - python -m pip install -e . pytest + python -m pip install -e ".[dev]" - - name: Run pytest + - name: Lint with ruff + run: ruff check . + + - name: Run pytest (unit tests + doctests) run: pytest diff --git a/.github/workflows/r-tests.yml b/.github/workflows/r-tests.yml index 3f22cd1..ba30a28 100644 --- a/.github/workflows/r-tests.yml +++ b/.github/workflows/r-tests.yml @@ -1,13 +1,14 @@ -name: R tests +name: R CMD check -on: push +on: + push: + pull_request: jobs: - test: + check: runs-on: ubuntu-latest - defaults: - run: - working-directory: r + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - uses: actions/checkout@v4 @@ -18,7 +19,12 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: working-directory: r - extra-packages: any::devtools + extra-packages: any::rcmdcheck + needs: check - - name: Run devtools::test() - run: Rscript -e 'devtools::test(stop_on_failure = TRUE)' + # Runs the tests and also catches stale man/ pages, NAMESPACE problems + # and broken examples, which devtools::test() alone would miss. + - uses: r-lib/actions/check-r-package@v2 + with: + working-directory: r + upload-snapshots: true diff --git a/README.md b/README.md index 9a6b9b9..593470d 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,8 @@ # rushtools +[![Python tests](https://github.com/halireena/rushtools/actions/workflows/python-tests.yml/badge.svg)](https://github.com/halireena/rushtools/actions/workflows/python-tests.yml) +[![R CMD check](https://github.com/halireena/rushtools/actions/workflows/r-tests.yml/badge.svg)](https://github.com/halireena/rushtools/actions/workflows/r-tests.yml) + My personal function libraries, in two languages: | Folder | What it is | @@ -24,7 +27,26 @@ pip install "git+https://github.com/halireena/rushtools#subdirectory=python" ``` ```python -import rushtools +from rushtools import gc_content, read_fasta + +for rid, seq in read_fasta("genes.fasta"): + print(rid, round(gc_content(seq), 3)) +``` + +It also installs a `rushtools` command: + +```bash +rushtools stats genes.fasta # sequence count, total bp, longest, N50 +rushtools gc genes.fasta # GC content per sequence (NA for all-N records) +rushtools orfs genes.fasta --min-len 90 +``` + +R usage: + +```r +gc_content(c(a = "ATGC", b = "GGCC")) +reverse_complement("ATGCRY") +tidy_de_table(res, lfc_col = "logFC", p_col = "P.Value") ``` ## Run the tests locally @@ -38,6 +60,7 @@ devtools::test() Python (from the `python/` folder): ```bash -pip install -e . pytest -pytest +pip install -e ".[dev]" +pytest # unit tests + docstring examples +ruff check . # lint ``` diff --git a/python/README.md b/python/README.md index 3c8f0c6..9ef9cf7 100644 --- a/python/README.md +++ b/python/README.md @@ -5,17 +5,18 @@ Small, tested Python helpers for everyday plant bioinformatics, built while lear ## Install ```bash -pip install git+https://github.com/USER/rushtools.git # latest version -pip install git+https://github.com/USER/rushtools.git@v0.1.0 # a specific tagged version +pip install "git+https://github.com/halireena/rushtools#subdirectory=python" # latest version +pip install "git+https://github.com/halireena/rushtools@v0.1.0#subdirectory=python" # a tagged version ``` For development (editable install + tests): ```bash -git clone https://github.com/USER/rushtools.git -cd rushtools +git clone https://github.com/halireena/rushtools.git +cd rushtools/python pip install -e ".[dev]" -pytest +pytest # unit tests + the examples in the docstrings +ruff check . # lint ``` ## Use in Python diff --git a/python/pyproject.toml b/python/pyproject.toml index 7399132..a37684c 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "rushtools" -version = "0.1.0" +version = "0.2.0" description = "Small, tested helpers for everyday plant bioinformatics: sequences, FASTA/FASTQ, assembly stats." readme = "README.md" requires-python = ">=3.10" @@ -14,7 +14,7 @@ keywords = ["bioinformatics", "DNA", "FASTA", "teaching"] dependencies = [] [project.optional-dependencies] -dev = ["pytest>=7"] +dev = ["pytest>=7", "ruff>=0.5"] [project.scripts] rushtools = "rushtools.cli:main" @@ -22,5 +22,16 @@ rushtools = "rushtools.cli:main" [tool.setuptools.packages.find] where = ["src"] +[tool.setuptools.package-data] +rushtools = ["py.typed"] + [tool.pytest.ini_options] -testpaths = ["tests"] +testpaths = ["tests", "src"] +addopts = "--doctest-modules" + +[tool.ruff] +line-length = 110 +target-version = "py310" + +[tool.ruff.lint] +select = ["E", "F", "W", "I", "B", "UP"] diff --git a/python/src/rushtools/__init__.py b/python/src/rushtools/__init__.py index 07d19aa..1c84e78 100644 --- a/python/src/rushtools/__init__.py +++ b/python/src/rushtools/__init__.py @@ -3,6 +3,8 @@ The most useful functions are imported here, so users can write ``from rushtools import gc_content`` instead of ``from rushtools.seq import gc_content``. """ +from importlib.metadata import PackageNotFoundError, version + from rushtools.assembly import n50 from rushtools.parsers import read_fasta, read_fastq, write_fasta from rushtools.seq import ( @@ -15,7 +17,10 @@ translate, ) -__version__ = "0.1.0" +try: + __version__ = version("rushtools") +except PackageNotFoundError: # running from a source tree without installing + __version__ = "0.0.0" __all__ = [ "CODON_TABLE", diff --git a/python/src/rushtools/cli.py b/python/src/rushtools/cli.py index 76a2654..5f742fa 100644 --- a/python/src/rushtools/cli.py +++ b/python/src/rushtools/cli.py @@ -13,11 +13,17 @@ def _cmd_gc(args: argparse.Namespace) -> None: print("id\tlength\tgc") for rid, seq in read_fasta(args.fasta): - print(f"{rid}\t{len(seq)}\t{gc_content(seq):.{args.digits}f}") + try: + gc = f"{gc_content(seq):.{args.digits}f}" + except ValueError: # no A/C/G/T bases, e.g. an all-N record + gc = "NA" + print(f"{rid}\t{len(seq)}\t{gc}") def _cmd_stats(args: argparse.Namespace) -> None: lengths = [len(seq) for _, seq in read_fasta(args.fasta)] + if not lengths: + raise ValueError(f"{args.fasta}: no sequences found") print(f"sequences\t{len(lengths)}") print(f"total_bp\t{sum(lengths)}") print(f"longest\t{max(lengths)}") @@ -65,7 +71,7 @@ def main(argv: list[str] | None = None) -> int: args = build_parser().parse_args(argv) try: args.func(args) - except (FileNotFoundError, ValueError) as err: + except (OSError, ValueError) as err: print(f"rushtools: error: {err}", file=sys.stderr) return 1 return 0 diff --git a/python/src/rushtools/parsers.py b/python/src/rushtools/parsers.py index 9a1983f..00f364c 100644 --- a/python/src/rushtools/parsers.py +++ b/python/src/rushtools/parsers.py @@ -24,7 +24,10 @@ def read_fasta(path) -> Iterator[tuple[str, str]]: if line.startswith(">"): if rid is not None: yield rid, "".join(chunks) - rid, chunks = line[1:].split()[0], [] + words = line[1:].split() + if not words: + raise ValueError(f"{path}: empty FASTA header") + rid, chunks = words[0], [] else: if rid is None: raise ValueError(f"{path}: sequence before the first '>' header") @@ -40,12 +43,17 @@ def read_fastq(path) -> Iterator[tuple[str, str, list[int]]]: header = fh.readline() if not header: return + if not header.strip(): + continue # tolerate blank lines between records seq = fh.readline().rstrip("\n") plus = fh.readline() qual = fh.readline().rstrip("\n") if not header.startswith("@") or not plus.startswith("+") or len(seq) != len(qual): raise ValueError(f"{path}: malformed FASTQ record near {header.strip()!r}") - yield header[1:].split()[0], seq, [ord(c) - 33 for c in qual] + words = header[1:].split() + if not words: + raise ValueError(f"{path}: empty FASTQ header") + yield words[0], seq, [ord(c) - 33 for c in qual] def write_fasta(records: Iterable[tuple[str, str]], path, width: int = 60) -> int: diff --git a/python/src/rushtools/py.typed b/python/src/rushtools/py.typed new file mode 100644 index 0000000..e69de29 diff --git a/python/src/rushtools/seq.py b/python/src/rushtools/seq.py index 22dfd7a..683cdd9 100644 --- a/python/src/rushtools/seq.py +++ b/python/src/rushtools/seq.py @@ -20,7 +20,8 @@ def gc_content(seq: str, ignore_n: bool = True) -> float: Args: seq: DNA sequence, any case. - ignore_n: leave N bases out of the denominator. + ignore_n: leave N and any other non-ACGT symbols out of the denominator + (matches the R package). Raises: ValueError: if the sequence has no usable bases. @@ -30,7 +31,7 @@ def gc_content(seq: str, ignore_n: bool = True) -> float: """ s = seq.upper() gc = s.count("G") + s.count("C") - denom = len(s) - s.count("N") if ignore_n else len(s) + denom = sum(s.count(b) for b in "ACGT") if ignore_n else len(s) if denom == 0: raise ValueError("sequence has no A/C/G/T bases") return gc / denom @@ -99,13 +100,25 @@ def find_orfs(seq: str, min_len: int = 30) -> list[dict]: return sorted(orfs, key=lambda o: (o["start"], o["strand"])) -def kmer_count(seq: str, k: int) -> dict[str, int]: - """Count overlapping k-mers. Returns {kmer: count}.""" +def kmer_count(seq: str, k: int, *, skip_ambiguous: bool = True) -> dict[str, int]: + """Count overlapping k-mers, case-insensitively. Returns {kmer: count}. + + Args: + seq: DNA sequence, any case. + k: k-mer length. + skip_ambiguous: leave out k-mers containing anything other than A/C/G/T. + + >>> kmer_count("ACGNacg", 2) + {'AC': 2, 'CG': 2} + """ if k < 1: raise ValueError("k must be >= 1") + s = seq.upper() counts: dict[str, int] = {} - for i in range(len(seq) - k + 1): - kmer = seq[i:i + k] + for i in range(len(s) - k + 1): + kmer = s[i:i + k] + if skip_ambiguous and kmer.strip("ACGT"): + continue counts[kmer] = counts.get(kmer, 0) + 1 return counts @@ -114,4 +127,4 @@ def hamming(a: str, b: str) -> int: """Number of positions at which two equal-length sequences differ.""" if len(a) != len(b): raise ValueError(f"sequences differ in length: {len(a)} vs {len(b)}") - return sum(x != y for x, y in zip(a, b)) + return sum(x != y for x, y in zip(a, b, strict=True)) diff --git a/python/tests/test_cli.py b/python/tests/test_cli.py index 2929fa0..c7a1229 100644 --- a/python/tests/test_cli.py +++ b/python/tests/test_cli.py @@ -18,3 +18,33 @@ def test_cli_stats(tmp_path, capsys): def test_cli_missing_file(capsys): assert main(["gc", "does_not_exist.fasta"]) == 1 assert "error" in capsys.readouterr().err + + +def test_cli_gc_reports_na_for_all_n_record(tmp_path, capsys): + fasta = tmp_path / "n.fasta" + write_fasta([("allN", "NNNN"), ("b", "ACGT")], fasta) + assert main(["gc", str(fasta)]) == 0 + lines = capsys.readouterr().out.splitlines() + assert lines[1:] == ["allN\t4\tNA", "b\t4\t0.500"] + + +def test_cli_stats_empty_file(tmp_path, capsys): + fasta = tmp_path / "empty.fasta" + fasta.write_text("") + assert main(["stats", str(fasta)]) == 1 + assert "no sequences found" in capsys.readouterr().err + + +def test_cli_empty_header_is_clean_error(tmp_path, capsys): + fasta = tmp_path / "bad.fasta" + fasta.write_text(">\nACGT\n") + assert main(["gc", str(fasta)]) == 1 + assert "empty FASTA header" in capsys.readouterr().err + + +def test_cli_orfs(tmp_path, capsys): + fasta = tmp_path / "orf.fasta" + write_fasta([("g", "CC" + "ATG" + "GCT" * 10 + "TAA" + "CC")], fasta) + assert main(["orfs", str(fasta), "--min-len", "30"]) == 0 + lines = capsys.readouterr().out.splitlines() + assert lines[1] == "g\t+\t3\t38\t36\t11" diff --git a/python/tests/test_parsers.py b/python/tests/test_parsers.py index 3989365..cb01b44 100644 --- a/python/tests/test_parsers.py +++ b/python/tests/test_parsers.py @@ -39,3 +39,31 @@ def test_n50(): assert n50([2, 2, 2, 3, 3, 4, 8, 8]) == 8 with pytest.raises(ValueError): n50([]) + + +def test_fasta_empty_header_raises(tmp_path): + path = tmp_path / "bad.fasta" + path.write_text(">\nACGT\n") + with pytest.raises(ValueError, match="empty FASTA header"): + list(read_fasta(path)) + + +def test_fastq_empty_header_raises(tmp_path): + path = tmp_path / "bad.fastq" + path.write_text("@\nACGT\n+\nIIII\n") + with pytest.raises(ValueError, match="empty FASTQ header"): + list(read_fastq(path)) + + +def test_fastq_tolerates_blank_lines_and_gzip(tmp_path): + path = tmp_path / "r.fastq.gz" + with gzip.open(path, "wt") as fh: + fh.write("@r1\nAC\n+\nII\n\n@r2\nGT\n+\n!!\n\n") + assert [r[0] for r in read_fastq(path)] == ["r1", "r2"] + + +def test_fastq_malformed_raises(tmp_path): + path = tmp_path / "bad.fastq" + path.write_text("@r1\nACGT\n+\nII\n") + with pytest.raises(ValueError, match="malformed"): + list(read_fastq(path)) diff --git a/python/tests/test_seq.py b/python/tests/test_seq.py index 8dc4a73..ed506db 100644 --- a/python/tests/test_seq.py +++ b/python/tests/test_seq.py @@ -69,3 +69,14 @@ def test_hamming(): assert hamming("GATTACA", "GACTATA") == 2 with pytest.raises(ValueError): hamming("ATG", "AT") + + +def test_gc_content_ignores_all_ambiguity_codes(): + # matches the R package: anything other than A/C/G/T leaves the denominator + assert gc_content("GCRYSS") == 1.0 + assert gc_content("GCRYSS", ignore_n=False) == pytest.approx(2 / 6) + + +def test_kmer_count_case_insensitive_and_skips_ambiguous(): + assert kmer_count("acgNN", 2) == {"AC": 1, "CG": 1} + assert kmer_count("acgNN", 2, skip_ambiguous=False) == {"AC": 1, "CG": 1, "GN": 1, "NN": 1} diff --git a/r/NEWS.md b/r/NEWS.md index 441887d..d4ce97c 100644 --- a/r/NEWS.md +++ b/r/NEWS.md @@ -1,3 +1,10 @@ +# rushtools (development version) + +* `reverse_complement()` now accepts IUPAC ambiguity codes (R, Y, S, W, K, + M, B, D, H, V), matching the Python package. +* `count_motif()` now errors on an empty motif instead of returning a + meaningless count. + # rushtools 0.2.0 * New `count_motif()` counts motif occurrences in DNA sequences. diff --git a/r/R/sequences.R b/r/R/sequences.R index fd9b1eb..b55cac1 100644 --- a/r/R/sequences.R +++ b/r/R/sequences.R @@ -24,22 +24,23 @@ gc_content <- function(seqs, ignore_n = TRUE) { #' Reverse complement of DNA sequences #' -#' @param seqs Character vector of DNA sequences containing only A, C, G, T -#' or N (any case). +#' @param seqs Character vector of DNA sequences (any case) using A, C, G, T, +#' N or the other IUPAC ambiguity codes (R, Y, S, W, K, M, B, D, H, V). #' @return Character vector of upper-case reverse complements, named like `seqs`. #' @examples #' reverse_complement("ATGCCC") #' reverse_complement(c(fwd = "GGATCCNA")) +#' reverse_complement("ACGTRY") # IUPAC codes: R <-> Y #' @export reverse_complement <- function(seqs) { stopifnot("`seqs` must be character" = is.character(seqs)) seqs <- toupper(seqs) - bad <- grepl("[^ACGTN]", seqs) + bad <- grepl("[^ACGTRYSWKMBDHVN]", seqs) if (any(bad)) { stop("Non-DNA characters found in sequence(s) ", paste(which(bad), collapse = ", "), call. = FALSE) } - comp <- chartr("ACGTN", "TGCAN", seqs) + comp <- chartr("ACGTRYSWKMBDHVN", "TGCAYRSWMKVHDBN", seqs) out <- vapply(strsplit(comp, ""), function(b) paste(rev(b), collapse = ""), FUN.VALUE = character(1)) names(out) <- names(seqs) @@ -58,6 +59,9 @@ reverse_complement <- function(seqs) { #' @export count_motif <- function(seqs, motif) { stopifnot(is.character(seqs), is.character(motif), length(motif) == 1) + if (is.na(motif) || !nzchar(motif)) { + stop("`motif` must be a non-empty string", call. = FALSE) + } hits <- gregexpr(toupper(motif), toupper(seqs), fixed = TRUE) out <- vapply(hits, function(h) sum(h > 0), FUN.VALUE = integer(1)) names(out) <- names(seqs) diff --git a/r/man/reverse_complement.Rd b/r/man/reverse_complement.Rd index 9fef630..7cde374 100644 --- a/r/man/reverse_complement.Rd +++ b/r/man/reverse_complement.Rd @@ -7,8 +7,8 @@ reverse_complement(seqs) } \arguments{ -\item{seqs}{Character vector of DNA sequences containing only A, C, G, T -or N (any case).} +\item{seqs}{Character vector of DNA sequences (any case) using A, C, G, T, +N or the other IUPAC ambiguity codes (R, Y, S, W, K, M, B, D, H, V).} } \value{ Character vector of upper-case reverse complements, named like \code{seqs}. @@ -19,4 +19,5 @@ Reverse complement of DNA sequences \examples{ reverse_complement("ATGCCC") reverse_complement(c(fwd = "GGATCCNA")) +reverse_complement("ACGTRY") # IUPAC codes: R <-> Y } diff --git a/r/tests/testthat/test-motif.R b/r/tests/testthat/test-motif.R index 62e3352..18ecf42 100644 --- a/r/tests/testthat/test-motif.R +++ b/r/tests/testthat/test-motif.R @@ -5,3 +5,8 @@ test_that("count_motif counts non-overlapping matches", { expect_equal(count_motif("AAAA", "AA"), 2L) expect_named(count_motif(c(x = "AA"), "A"), "x") }) + +test_that("count_motif rejects an empty motif", { + expect_error(count_motif("ACGT", ""), "non-empty") + expect_error(count_motif("ACGT", NA_character_), "non-empty") +}) diff --git a/r/tests/testthat/test-sequences.R b/r/tests/testthat/test-sequences.R index 770a56a..5a3a195 100644 --- a/r/tests/testthat/test-sequences.R +++ b/r/tests/testthat/test-sequences.R @@ -22,6 +22,11 @@ test_that("reverse_complement is correct and reversible", { expect_equal(reverse_complement(reverse_complement(s)), s) }) +test_that("reverse_complement handles IUPAC ambiguity codes", { + expect_equal(reverse_complement("ACGTRYKMBV"), "BVKMRYACGT") + expect_equal(reverse_complement("SWN"), "NWS") +}) + test_that("reverse_complement rejects RNA and other letters", { expect_error(reverse_complement("AUGC"), "Non-DNA") }) From 902fd0b2ea05e729a4d14c6161ab7d4aadd0289c Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 5 Oct 2026 09:11:01 +0000 Subject: [PATCH 2/4] Fix NA in R reverse_complement() and two silent CLI failures R: - reverse_complement(NA) returned the string "NA" instead of NA. Python CLI: - `rushtools gc --digits -1` printed NA for every record (the format error was mistaken for "no A/C/G/T bases"); negative values are now rejected by argparse with a clear message. - `rushtools revcomp HELLO` printed "OLLED" without complaint; it now errors on characters that are not DNA/IUPAC codes. The library function reverse_complement() is unchanged. Tests for each, plus a regression test that FASTQ files with Windows line endings read cleanly. pytest now also works before installing (pythonpath = ["src"]). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WqrAUWjyEmFz36zphqfoBs --- python/pyproject.toml | 1 + python/src/rushtools/cli.py | 21 +++++++++++++++++---- python/tests/test_cli.py | 30 ++++++++++++++++++++++++++++++ python/tests/test_parsers.py | 7 +++++++ r/NEWS.md | 2 ++ r/R/sequences.R | 2 ++ r/man/reverse_complement.Rd | 1 + r/tests/testthat/test-sequences.R | 5 +++++ 8 files changed, 65 insertions(+), 4 deletions(-) diff --git a/python/pyproject.toml b/python/pyproject.toml index a37684c..ebd438a 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -27,6 +27,7 @@ rushtools = ["py.typed"] [tool.pytest.ini_options] testpaths = ["tests", "src"] +pythonpath = ["src"] # so `pytest` also works before `pip install -e .` addopts = "--doctest-modules" [tool.ruff] diff --git a/python/src/rushtools/cli.py b/python/src/rushtools/cli.py index 5f742fa..7d60604 100644 --- a/python/src/rushtools/cli.py +++ b/python/src/rushtools/cli.py @@ -9,15 +9,25 @@ from rushtools.parsers import read_fasta from rushtools.seq import find_orfs, gc_content, reverse_complement +_IUPAC_DNA = set("ACGTRYSWKMBDHVNacgtryswkmbdhvn") + + +def _non_negative_int(text: str) -> int: + value = int(text) + if value < 0: + raise argparse.ArgumentTypeError(f"must be 0 or more, got {value}") + return value + def _cmd_gc(args: argparse.Namespace) -> None: print("id\tlength\tgc") for rid, seq in read_fasta(args.fasta): try: - gc = f"{gc_content(seq):.{args.digits}f}" + gc = gc_content(seq) except ValueError: # no A/C/G/T bases, e.g. an all-N record - gc = "NA" - print(f"{rid}\t{len(seq)}\t{gc}") + print(f"{rid}\t{len(seq)}\tNA") + else: + print(f"{rid}\t{len(seq)}\t{gc:.{args.digits}f}") def _cmd_stats(args: argparse.Namespace) -> None: @@ -31,6 +41,9 @@ def _cmd_stats(args: argparse.Namespace) -> None: def _cmd_revcomp(args: argparse.Namespace) -> None: + bad = sorted(set(args.sequence) - _IUPAC_DNA) + if bad: + raise ValueError(f"not a DNA sequence (unexpected characters: {''.join(bad)})") print(reverse_complement(args.sequence)) @@ -48,7 +61,7 @@ def build_parser() -> argparse.ArgumentParser: p = sub.add_parser("gc", help="GC content of every sequence in a FASTA file") p.add_argument("fasta", help="input FASTA (plain or .gz)") - p.add_argument("-d", "--digits", type=int, default=3, help="decimal places (default: 3)") + p.add_argument("-d", "--digits", type=_non_negative_int, default=3, help="decimal places (default: 3)") p.set_defaults(func=_cmd_gc) p = sub.add_parser("stats", help="number of sequences, total length, N50") diff --git a/python/tests/test_cli.py b/python/tests/test_cli.py index c7a1229..8ab2265 100644 --- a/python/tests/test_cli.py +++ b/python/tests/test_cli.py @@ -1,3 +1,5 @@ +import pytest + from rushtools.cli import main from rushtools.parsers import write_fasta @@ -48,3 +50,31 @@ def test_cli_orfs(tmp_path, capsys): assert main(["orfs", str(fasta), "--min-len", "30"]) == 0 lines = capsys.readouterr().out.splitlines() assert lines[1] == "g\t+\t3\t38\t36\t11" + + +def test_cli_gc_rejects_negative_digits(tmp_path, capsys): + # used to print "NA" for every record instead of complaining + fasta = tmp_path / "x.fasta" + write_fasta([("a", "ACGT")], fasta) + with pytest.raises(SystemExit) as exc: + main(["gc", str(fasta), "--digits", "-1"]) + assert exc.value.code == 2 + assert "must be 0 or more" in capsys.readouterr().err + + +def test_cli_gc_digits(tmp_path, capsys): + fasta = tmp_path / "x.fasta" + write_fasta([("a", "GGC")], fasta) + assert main(["gc", str(fasta), "-d", "0"]) == 0 + assert capsys.readouterr().out.splitlines()[1] == "a\t3\t1" + + +def test_cli_revcomp_rejects_non_dna(capsys): + # used to print "OLLEH" -> "OLLED" without any warning + assert main(["revcomp", "HELLO"]) == 1 + assert "not a DNA sequence" in capsys.readouterr().err + + +def test_cli_revcomp_accepts_iupac_and_lowercase(capsys): + assert main(["revcomp", "acgtRY"]) == 0 + assert capsys.readouterr().out.strip() == "RYacgt" diff --git a/python/tests/test_parsers.py b/python/tests/test_parsers.py index cb01b44..bc8a89c 100644 --- a/python/tests/test_parsers.py +++ b/python/tests/test_parsers.py @@ -67,3 +67,10 @@ def test_fastq_malformed_raises(tmp_path): path.write_text("@r1\nACGT\n+\nII\n") with pytest.raises(ValueError, match="malformed"): list(read_fastq(path)) + + +def test_read_fastq_windows_line_endings(tmp_path): + # text mode turns \r\n into \n, so files saved on Windows read cleanly + path = tmp_path / "crlf.fastq" + path.write_bytes(b"@r1\r\nACGT\r\n+\r\nII5!\r\n") + assert list(read_fastq(path)) == [("r1", "ACGT", [40, 40, 20, 0])] diff --git a/r/NEWS.md b/r/NEWS.md index d4ce97c..7b2ab4d 100644 --- a/r/NEWS.md +++ b/r/NEWS.md @@ -4,6 +4,8 @@ M, B, D, H, V), matching the Python package. * `count_motif()` now errors on an empty motif instead of returning a meaningless count. +* `reverse_complement()` keeps `NA` as `NA`; it used to return the text + `"NA"`. # rushtools 0.2.0 diff --git a/r/R/sequences.R b/r/R/sequences.R index b55cac1..9bec1a3 100644 --- a/r/R/sequences.R +++ b/r/R/sequences.R @@ -27,6 +27,7 @@ gc_content <- function(seqs, ignore_n = TRUE) { #' @param seqs Character vector of DNA sequences (any case) using A, C, G, T, #' N or the other IUPAC ambiguity codes (R, Y, S, W, K, M, B, D, H, V). #' @return Character vector of upper-case reverse complements, named like `seqs`. +#' `NA` inputs stay `NA`. #' @examples #' reverse_complement("ATGCCC") #' reverse_complement(c(fwd = "GGATCCNA")) @@ -43,6 +44,7 @@ reverse_complement <- function(seqs) { comp <- chartr("ACGTRYSWKMBDHVN", "TGCAYRSWMKVHDBN", seqs) out <- vapply(strsplit(comp, ""), function(b) paste(rev(b), collapse = ""), FUN.VALUE = character(1)) + out[is.na(seqs)] <- NA_character_ # otherwise NA comes back as the text "NA" names(out) <- names(seqs) out } diff --git a/r/man/reverse_complement.Rd b/r/man/reverse_complement.Rd index 7cde374..ee6e7b6 100644 --- a/r/man/reverse_complement.Rd +++ b/r/man/reverse_complement.Rd @@ -12,6 +12,7 @@ N or the other IUPAC ambiguity codes (R, Y, S, W, K, M, B, D, H, V).} } \value{ Character vector of upper-case reverse complements, named like \code{seqs}. +\code{NA} inputs stay \code{NA}. } \description{ Reverse complement of DNA sequences diff --git a/r/tests/testthat/test-sequences.R b/r/tests/testthat/test-sequences.R index 5a3a195..4de844d 100644 --- a/r/tests/testthat/test-sequences.R +++ b/r/tests/testthat/test-sequences.R @@ -30,3 +30,8 @@ test_that("reverse_complement handles IUPAC ambiguity codes", { test_that("reverse_complement rejects RNA and other letters", { expect_error(reverse_complement("AUGC"), "Non-DNA") }) + +test_that("reverse_complement keeps NA as NA (not the text \"NA\")", { + expect_identical(reverse_complement(c(a = "ATG", b = NA)), + c(a = "CAT", b = NA_character_)) +}) From fb50ade545b24e1fe5e4cce5ab7b9596a1330c38 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 5 Oct 2026 09:14:09 +0000 Subject: [PATCH 3/4] Make the READMEs runnable for newcomers: example data, outputs, glossary - examples/: tiny made-up genes.fasta and reads.fastq, so every example that used to point at a non-existent "genes.fasta" now runs as written. - Top-level README: prerequisites, "Getting started in 5 minutes" with a virtual environment, the real output under every command, R install from the clone (no extra packages needed), a glossary (FASTA, FASTQ, Phred, GC, IUPAC, ORF, k-mer, N50, fold change, BH/padj, ...) and a troubleshooting table. The Python loop example now handles the all-N record instead of crashing on it. - python/README: function table, outputs, a working "branch or tag" install line (there is no v0.1.0 tag yet), venv steps. - r/README: the install line pointed at "yourGitHubName/rushtools" without subdir = "r"; fixed. Examples for every function with output; how to keep DESeq2 gene IDs (row names) through tidy_de_table(); how to run the tests with testthat only. - CITATION.cff, CONTRIBUTING.md, python/CHANGELOG.md. - python/LICENSE copyright holder now matches the R licence; R version 0.2.0.9000 to match the "development version" heading in NEWS.md. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WqrAUWjyEmFz36zphqfoBs --- CITATION.cff | 18 ++++ CONTRIBUTING.md | 22 +++++ README.md | 228 ++++++++++++++++++++++++++++++++++++++----- examples/README.md | 15 +++ examples/genes.fasta | 10 ++ examples/reads.fastq | 12 +++ python/CHANGELOG.md | 31 ++++++ python/LICENSE | 2 +- python/README.md | 103 ++++++++++++++++--- r/DESCRIPTION | 2 +- r/README.md | 68 ++++++++++++- 11 files changed, 468 insertions(+), 43 deletions(-) create mode 100644 CITATION.cff create mode 100644 CONTRIBUTING.md create mode 100644 examples/README.md create mode 100644 examples/genes.fasta create mode 100644 examples/reads.fastq create mode 100644 python/CHANGELOG.md diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..8f37c7f --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,18 @@ +cff-version: 1.2.0 +message: "If you use rushtools, please cite it as below." +type: software +title: "rushtools: small, tested helpers for everyday plant bioinformatics in R and Python" +authors: + - given-names: Halireena Rushdiha + family-names: Mohomed +version: 0.2.0 +license: MIT +repository-code: "https://github.com/halireena/rushtools" +keywords: + - bioinformatics + - DNA + - FASTA + - FASTQ + - differential expression + - R + - Python diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..485f333 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,22 @@ +# Contributing + +This is a small personal project, but bug reports, questions and +suggestions are very welcome: please open an issue on GitHub with the +command or code you ran, what you expected and what happened (a few lines +of an example file help a lot). + +If you'd like to send a fix: + +1. Fork and clone the repository, then make a branch. +2. Python: `cd python && pip install -e ".[dev]"`, then make sure + `pytest` and `ruff check .` both pass. +3. R: from `r/`, run `Rscript -e 'testthat::test_local(".")'`. If you + change a roxygen comment (`#'`), regenerate the help pages with + `Rscript -e 'roxygen2::roxygenise()'`. +4. Add a test for any bug you fix, and a line to `python/CHANGELOG.md` or + `r/NEWS.md`. +5. Keep the two packages in step where a function exists in both + (e.g. `gc_content()` and `reverse_complement()`). + +The example files in `examples/` are used by the READMEs; if you change +them, update the outputs shown there. diff --git a/README.md b/README.md index 593470d..637cc54 100644 --- a/README.md +++ b/README.md @@ -3,64 +3,242 @@ [![Python tests](https://github.com/halireena/rushtools/actions/workflows/python-tests.yml/badge.svg)](https://github.com/halireena/rushtools/actions/workflows/python-tests.yml) [![R CMD check](https://github.com/halireena/rushtools/actions/workflows/r-tests.yml/badge.svg)](https://github.com/halireena/rushtools/actions/workflows/r-tests.yml) -My personal function libraries, in two languages: +My personal function libraries for everyday plant bioinformatics, in two languages: | Folder | What it is | |-------------|----------------------------------------------| -| `r/` | R package `rushtools` (`R/`, `tests/`, `DESCRIPTION`) | -| `python/` | Python package `rushtools` (`src/` layout, `pyproject.toml`, `tests/`) | +| `python/` | Python package `rushtools` (`src/` layout, `pyproject.toml`, `tests/`): DNA sequence helpers, FASTA/FASTQ readers, assembly stats and a `rushtools` command | +| `r/` | R package `rushtools` (`R/`, `tests/`, `DESCRIPTION`): GC content, reverse complements, fold changes, z-scores, p-value adjustment, tidying differential expression tables | +| `examples/` | Tiny made-up FASTA and FASTQ files so every example below runs as written | Tests for both run automatically on every push (see the Actions tab). +New to some of the words? See the [glossary](#glossary) at the bottom. -## Install the R package +## What you need first -```r -# install.packages("remotes") -remotes::install_github("halireena/rushtools", subdir = "r") -library(rushtools) +- **git**, to download the code (`git --version` to check) +- **Python 3.10 or newer** for the Python package (`python --version`; on macOS/Linux it may be called `python3`) +- **R 4.0 or newer** for the R package (`R --version`) + +You only need the language you plan to use. + +## Getting started in 5 minutes (Python) + +Run these in a terminal. They download the code, make a private Python +environment (a "virtual environment", so nothing clashes with other projects), +install the package and try it on the example files. + +```bash +git clone https://github.com/halireena/rushtools.git +cd rushtools +python -m venv .venv +source .venv/bin/activate # Windows: .venv\Scripts\activate +pip install -e "./python[dev]" +rushtools --version ``` -## Install the Python package +```text +rushtools 0.2.0 +``` + +Now ask it about the example FASTA file. Each command prints a +tab-separated table you can paste into a spreadsheet: ```bash -pip install "git+https://github.com/halireena/rushtools#subdirectory=python" +rushtools stats examples/genes.fasta # sequence count, total bp, longest, N50 +``` + +```text +sequences 4 +total_bp 204 +longest 125 +n50 125 +``` + +```bash +rushtools gc examples/genes.fasta # GC content per sequence (NA for all-N records) ``` +```text +id length gc +geneA 125 0.496 +geneB 32 0.875 +geneC 27 0.095 +scaffold_gap 20 NA +``` + +`scaffold_gap` is only `N`s (unknown bases), so it has no GC content to report. + +```bash +rushtools orfs examples/genes.fasta --min-len 90 +``` + +```text +id strand start end length_nt protein_aa +geneA + 6 113 108 35 +``` + +One open reading frame in `geneA`, on the forward (`+`) strand, from base 6 to +base 113: 108 nucleotides including the stop codon, which code for 35 amino acids. + +```bash +rushtools revcomp ATGCCC +``` + +```text +GGGCAT +``` + +And the same from Python (save as a `.py` file and run it from the `rushtools` folder, +or paste it into `python`): + ```python from rushtools import gc_content, read_fasta -for rid, seq in read_fasta("genes.fasta"): - print(rid, round(gc_content(seq), 3)) +for rid, seq in read_fasta("examples/genes.fasta"): + try: + print(rid, round(gc_content(seq), 3)) + except ValueError as err: # the all-N record has no real bases to count + print(rid, "skipped:", err) ``` -It also installs a `rushtools` command: +```text +geneA 0.496 +geneB 0.875 +geneC 0.095 +scaffold_gap skipped: sequence has no A/C/G/T bases +``` + +Reading sequencing reads with their quality scores: + +```python +from rushtools import read_fastq + +for rid, seq, quals in read_fastq("examples/reads.fastq"): + print(rid, seq, "mean quality:", round(sum(quals) / len(quals), 1)) +``` + +```text +read1 ATGCCGGGGTACCGTAGTTT mean quality: 40.0 +read2 GTCTGCGACCTCACAAAATT mean quality: 30.5 +read3 TGGGCATGCNTTTATCCCGA mean quality: 12.5 +``` + +More Python functions and options: [`python/README.md`](python/README.md). +To use your own data, swap `examples/genes.fasta` for the path to your file +(plain or gzipped, e.g. `my_genes.fa.gz`). + +When you are done, `deactivate` leaves the virtual environment; next time, +`cd rushtools` and `source .venv/bin/activate` again. + +### Install without cloning + +If you just want the package and not the example files or tests: ```bash -rushtools stats genes.fasta # sequence count, total bp, longest, N50 -rushtools gc genes.fasta # GC content per sequence (NA for all-N records) -rushtools orfs genes.fasta --min-len 90 +pip install "git+https://github.com/halireena/rushtools#subdirectory=python" ``` -R usage: +## Getting started (R) + +Start R **from the top `rushtools` folder** you cloned above and install the +package straight from the `r/` folder. This needs no extra packages: + +```r +install.packages("r", repos = NULL, type = "source") +library(rushtools) +``` + +Or install from GitHub without cloning: + +```r +install.packages("remotes") # once +remotes::install_github("halireena/rushtools", subdir = "r") +library(rushtools) +``` + +Then: ```r gc_content(c(a = "ATGC", b = "GGCC")) -reverse_complement("ATGCRY") +#> a b +#> 0.5 1.0 + +reverse_complement("ATGCRY") # R and Y are IUPAC codes for "A or G" and "C or T" +#> [1] "RYGCAT" + +# a differential expression result, e.g. from limma +res <- data.frame(gene = c("g1", "g2", "g3"), + logFC = c(2.1, -1.5, 0.2), + P.Value = c(0.0001, 0.003, 0.7)) tidy_de_table(res, lfc_col = "logFC", p_col = "P.Value") +#> gene log2fc pvalue padj direction +#> 1 g1 2.1 1e-04 0.0003 up +#> 2 g2 -1.5 3e-03 0.0045 down +#> 3 g3 0.2 7e-01 0.7000 ns ``` +`tidy_de_table()` renames the columns to `log2fc`, `pvalue` and `padj` +(BH-adjusted p-value), and calls a gene `up` or `down` when `padj < 0.05` and +the fold change is at least 2x (`|log2fc| >= 1`). Every function has a help +page, e.g. `?tidy_de_table`. Full list: [`r/README.md`](r/README.md). + ## Run the tests locally -R (from the `r/` folder): +Python (from the `python/` folder, with the virtual environment active): -```r -devtools::test() +```bash +pytest # unit tests + the examples in the docstrings +ruff check . # lint (style and common-mistake checks) ``` -Python (from the `python/` folder): +R (from the `r/` folder; needs `install.packages("testthat")` once): ```bash -pip install -e ".[dev]" -pytest # unit tests + docstring examples -ruff check . # lint +Rscript -e 'testthat::test_local(".")' ``` + +or, if you use devtools, `devtools::test()` inside R. + +## Glossary + +| Term | Meaning | +|---|---| +| **FASTA** | Plain-text file of sequences. Each record is a `>` header line (the first word is the ID) followed by one or more lines of sequence. See [`examples/genes.fasta`](examples/genes.fasta). | +| **FASTQ** | Plain-text file of sequencing reads, 4 lines per read: `@ID`, sequence, `+`, and one quality character per base. See [`examples/reads.fastq`](examples/reads.fastq). | +| **Phred quality** | How sure the sequencer is about a base: Q = -10 log10(error probability). Q20 = 1 error in 100, Q30 = 1 in 1000, Q40 = 1 in 10,000. In FASTQ files each score is stored as the character with code Q + 33 (`I` = 40, `5` = 20, `!` = 0). | +| **GC content** | Fraction of bases that are G or C (0 to 1). Differs between genomes and genes; very high or low values can affect PCR and sequencing. | +| **N / IUPAC codes** | Letters for uncertain bases: `N` = any base, `R` = A or G, `Y` = C or T, `S` = G or C, `W` = A or T, `K` = G or T, `M` = A or C, `B`/`D`/`H`/`V` = "not A/C/G/T" respectively. | +| **Reverse complement** | The sequence of the other DNA strand, read in its own 5'->3' direction: swap A<->T and C<->G, then reverse. | +| **Strand / frame** | DNA is read in triplets (codons), so each strand has 3 possible reading frames: 6 in total. `+` is the sequence as written, `-` its reverse complement. | +| **ORF** | Open reading frame: a stretch that starts with the start codon `ATG` and runs to a stop codon (`TAA`, `TAG`, `TGA`). A candidate protein-coding region. | +| **k-mer** | Every substring of length k. `ATATG` has the 2-mers `AT`, `TA`, `AT`, `TG`. | +| **Hamming distance** | Number of positions at which two equal-length sequences differ. | +| **Contig / N50** | Contigs are the pieces a genome assembly is made of. N50 is the contig length L such that contigs of length L or longer hold at least half of all bases; bigger usually means a more complete assembly. | +| **log2 fold change** | How much a gene's expression changes: log2(treated / control). 1 = doubled, -1 = halved, 0 = no change. | +| **Pseudocount** | A small number (default 1) added before taking logs so that zero counts don't give `-Inf`. | +| **z-score** | (value - mean) / standard deviation: how many SDs a value is from the average. | +| **Multiple testing / BH / padj** | Testing thousands of genes at p < 0.05 gives many false positives by chance. The Benjamini-Hochberg (BH) adjustment turns p-values into `padj` values that control the false discovery rate (the expected share of false hits among the genes you call significant). | +| **Differential expression (DE)** | Finding genes whose expression differs between conditions, usually with DESeq2, limma or edgeR in R. | + +## Troubleshooting + +| Problem | Fix | +|---|---| +| `python: command not found` | Use `python3` instead (macOS/Linux), or install Python 3.10+ from python.org. | +| `ERROR: Package 'rushtools' requires a different Python` | Your Python is older than 3.10. Check with `python --version`. | +| `rushtools: command not found` | The virtual environment is not active: run `source .venv/bin/activate` (Windows: `.venv\Scripts\activate`) from the `rushtools` folder, then try again. | +| `No such file or directory: 'examples/genes.fasta'` | Run the command from the top `rushtools` folder (the one with `examples/` in it), or give the full path to the file. | +| `ValueError: sequence has no A/C/G/T bases` | `gc_content()` was given a sequence of only `N`s (or an empty one). Skip it with `try`/`except`, as in the example above; the `rushtools gc` command prints `NA` instead. | +| `error: sequence before the first '>' header` | The file is not FASTA. If its first line starts with `@`, it is FASTQ: use `read_fastq()`. | +| `ModuleNotFoundError: No module named 'rushtools'` when running `pytest` | Make sure the virtual environment is active and you ran `pip install -e "./python[dev]"`. | +| R: `there is no package called 'remotes'` | Run `install.packages("remotes")` first, or install from the clone with `install.packages("r", repos = NULL, type = "source")`. | +| R: `invalid package 'r'` / `no packages specified` | Start R from the top `rushtools` folder (or `setwd("path/to/rushtools")`), not from inside `r/`. | +| R: `Non-DNA characters found in sequence(s) 1` | The sequence contains something other than A/C/G/T/N/IUPAC codes, e.g. `U` (RNA) or a space. | + +## Citing, licence and contributing + +- If rushtools helps your work, please cite it: GitHub shows a "Cite this repository" button generated from [`CITATION.cff`](CITATION.cff). +- MIT licence ([`python/LICENSE`](python/LICENSE), [`r/LICENSE.md`](r/LICENSE.md)). +- What changed between versions: [`python/CHANGELOG.md`](python/CHANGELOG.md) and [`r/NEWS.md`](r/NEWS.md). +- Bug reports and suggestions are welcome; see [`CONTRIBUTING.md`](CONTRIBUTING.md). diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..44bb241 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,15 @@ +# Example data + +Tiny, **made-up** files so every example in the READMEs runs as written. +They are not real plant sequences. + +| File | What it is | +|---|---| +| `genes.fasta` | 4 DNA sequences in FASTA format: one with a 108-nt open reading frame (`geneA`), a GC-rich one, an AT-rich one with some `N`s, and one that is only `N`. | +| `reads.fastq` | 3 sequencing reads in FASTQ format, with good, falling and poor quality scores. | + +Run the examples from the top folder of the repository, e.g. + +```bash +rushtools stats examples/genes.fasta +``` diff --git a/examples/genes.fasta b/examples/genes.fasta new file mode 100644 index 0000000..2a09fcb --- /dev/null +++ b/examples/genes.fasta @@ -0,0 +1,10 @@ +>geneA made-up gene with one open reading frame +CTTCAATGCCGGGGTACCGTAGTTTGTCTGCGACCTCACAAAATTTGGGCATGCTTTTAT +CCCGACGCTCTCTATCCACACGATTGGCGACGTCGGGGCTCAAGAAGAATTAAGTTAGCT +AACGT +>geneB made-up GC-rich fragment +GCGGCCGCTAGCGGCCGGGCGCATGCCGCGGC +>geneC made-up AT-rich fragment with a run of unknown bases (N) +ATTATANNNNNNTATAAATTGCATTTA +>scaffold_gap made-up record that is only N (no real bases) +NNNNNNNNNNNNNNNNNNNN diff --git a/examples/reads.fastq b/examples/reads.fastq new file mode 100644 index 0000000..0ce3b87 --- /dev/null +++ b/examples/reads.fastq @@ -0,0 +1,12 @@ +@read1 made-up good-quality read +ATGCCGGGGTACCGTAGTTT ++ +IIIIIIIIIIIIIIIIIIII +@read2 made-up read whose quality drops towards the end +GTCTGCGACCTCACAAAATT ++ +IIIIIIIIIIII55555+++ +@read3 made-up poor read with an uncalled base (N) +TGGGCATGCNTTTATCCCGA ++ +555555+++!++++++++++ diff --git a/python/CHANGELOG.md b/python/CHANGELOG.md new file mode 100644 index 0000000..8df652c --- /dev/null +++ b/python/CHANGELOG.md @@ -0,0 +1,31 @@ +# Changelog (Python package) + +## Unreleased + +- `rushtools gc --digits` with a negative number now gives a clear error; + it used to print `NA` for every record. +- `rushtools revcomp` now refuses input that is not DNA (e.g. `HELLO`) + instead of printing nonsense. The `reverse_complement()` function itself + is unchanged. +- `pytest` works from `python/` even before `pip install -e .`. + +## 0.2.0 + +- `read_fasta()` / `read_fastq()`: a bare `>` or `@` header is a clear + `ValueError` instead of an `IndexError`; FASTQ blank lines between + records are tolerated. +- `rushtools gc` prints `NA` for records with no A/C/G/T bases and carries + on with the rest of the file. +- `rushtools stats` on an empty file says "no sequences found". +- `gc_content(ignore_n=True)` leaves every non-ACGT symbol out of the + denominator, matching the R package (it used to drop only `N`). +- `kmer_count()` is case-insensitive and skips k-mers with ambiguous bases + (`skip_ambiguous=False` for the old behaviour). +- `__version__` comes from the package metadata; ships `py.typed`. + +## 0.1.0 + +- First version: `gc_content()`, `reverse_complement()`, `translate()`, + `find_orfs()`, `kmer_count()`, `hamming()`, `n50()`, FASTA/FASTQ readers, + `write_fasta()` and the `rushtools` command (`gc`, `stats`, `revcomp`, + `orfs`). diff --git a/python/LICENSE b/python/LICENSE index a81570f..47ad064 100644 --- a/python/LICENSE +++ b/python/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2026 Rush +Copyright (c) 2026 Halireena Rushdiha Mohomed Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/python/README.md b/python/README.md index 9ef9cf7..2268791 100644 --- a/python/README.md +++ b/python/README.md @@ -1,46 +1,127 @@ # rushtools Small, tested Python helpers for everyday plant bioinformatics, built while learning Python. +No dependencies beyond the Python standard library. Needs Python 3.10 or newer. ## Install ```bash -pip install "git+https://github.com/halireena/rushtools#subdirectory=python" # latest version -pip install "git+https://github.com/halireena/rushtools@v0.1.0#subdirectory=python" # a tagged version +pip install "git+https://github.com/halireena/rushtools#subdirectory=python" # latest version +pip install "git+https://github.com/halireena/rushtools@main#subdirectory=python" # a branch or tag after @ ``` -For development (editable install + tests): +(`pip install git+...` needs git installed.) + +For development (editable install + tests), in a virtual environment: ```bash git clone https://github.com/halireena/rushtools.git cd rushtools/python +python -m venv .venv +source .venv/bin/activate # Windows: .venv\Scripts\activate pip install -e ".[dev]" pytest # unit tests + the examples in the docstrings ruff check . # lint ``` +The examples below use the tiny made-up files in +[`../examples/`](../examples/) and are run from the top `rushtools` folder +(`cd ..` if you are in `python/`). New words (ORF, N50, Phred...) are explained in the +[glossary](../README.md#glossary). + ## Use in Python ```python from rushtools import gc_content, reverse_complement, translate, read_fasta, n50 -gc_content("ATGCGC") # 0.6667 -reverse_complement("ATGC") # 'GCAT' -translate("ATGTGGTAA") # 'MW*' -for rid, seq in read_fasta("genes.fasta"): +gc_content("ATGCGC") # 0.6666666666666666 (4 of 6 bases are G or C) +reverse_complement("ATGC") # 'GCAT' +translate("ATGTGGTAA") # 'MW*' (* = stop codon) +translate("ATGTGGTAA", to_stop=True) # 'MW' +n50([100, 200, 300, 400, 500]) # 400 + +for rid, seq in read_fasta("examples/genes.fasta"): print(rid, len(seq)) ``` +```text +geneA 125 +geneB 32 +geneC 27 +scaffold_gap 20 +``` + +| Function | What it does | +|---|---| +| `gc_content(seq, ignore_n=True)` | GC fraction (0-1). N and other non-ACGT letters are left out unless `ignore_n=False`. Raises `ValueError` if there are no A/C/G/T bases. | +| `reverse_complement(seq)` | Reverse complement; keeps upper/lower case and handles IUPAC codes. | +| `translate(seq, to_stop=False, frame=0)` | DNA or RNA to protein with the standard genetic code; `*` marks a stop, `X` a codon with N. | +| `find_orfs(seq, min_len=30)` | ATG-to-stop open reading frames in all six frames, as dicts with `strand`, `frame`, `start`, `end` (1-based, on the forward strand), `length` (nt, stop included) and `protein`. | +| `kmer_count(seq, k)` | Counts of every overlapping k-mer, e.g. `kmer_count("ATATG", 2)` gives `{'AT': 2, 'TA': 1, 'TG': 1}`. | +| `hamming(a, b)` | Number of mismatches between two equal-length sequences. | +| `read_fasta(path)` | Yields `(id, sequence)` per record; the id is the first word of the header. Plain or `.gz`. | +| `read_fastq(path)` | Yields `(id, sequence, qualities)` per read, with qualities as Phred scores (`I` = 40). Plain or `.gz`; 4-line records. | +| `write_fasta(records, path, width=60)` | Writes `(id, sequence)` pairs; returns how many were written. | +| `n50(lengths)` | N50 of a list of contig lengths. | +| `CODON_TABLE` | The standard genetic code as `{"ATG": "M", ...}`. | + +ORFs in the example gene: + +```python +from rushtools import find_orfs, read_fasta + +seqs = dict(read_fasta("examples/genes.fasta")) +for orf in find_orfs(seqs["geneA"], min_len=90): + print(orf["strand"], orf["start"], orf["end"], orf["length"], orf["protein"]) +``` + +```text ++ 6 113 108 MPGYRSLSATSQNLGMLLSRRSLSTRLATSGLKKN +``` + ## Use on the command line +Every command prints a tab-separated table (or a single sequence) and +`rushtools COMMAND --help` explains its options. + ```bash rushtools --help -rushtools gc genes.fasta -rushtools stats genes.fasta +rushtools gc examples/genes.fasta +rushtools stats examples/genes.fasta rushtools revcomp ATGCCC -rushtools orfs genes.fasta --min-len 90 +rushtools orfs examples/genes.fasta --min-len 90 +``` + +Output of the last four: + +```text +id length gc +geneA 125 0.496 +geneB 32 0.875 +geneC 27 0.095 +scaffold_gap 20 NA +``` + +```text +sequences 4 +total_bp 204 +longest 125 +n50 125 +``` + +```text +GGGCAT ``` +```text +id strand start end length_nt protein_aa +geneA + 6 113 108 35 +``` + +`gc` prints `NA` for a record with no A/C/G/T bases; `-d 2` changes the +number of decimals. In `orfs`, `length_nt` includes the stop codon and +`protein_aa` does not. To save a table: `rushtools gc examples/genes.fasta > gc.tsv`. + ## Licence -MIT +MIT. Changes between versions: [CHANGELOG.md](CHANGELOG.md). diff --git a/r/DESCRIPTION b/r/DESCRIPTION index 26e69da..84166e6 100644 --- a/r/DESCRIPTION +++ b/r/DESCRIPTION @@ -1,6 +1,6 @@ Package: rushtools Title: Bioinformatics Helper Functions for Sequences and Expression Data -Version: 0.2.0 +Version: 0.2.0.9000 Authors@R: person("Halireena Rushdiha", "Mohomed", , "217553546+halireena@users.noreply.github.com", role = c("aut", "cre")) Description: Small, tested helper functions for everyday plant diff --git a/r/README.md b/r/README.md index c9a91c3..0f01756 100644 --- a/r/README.md +++ b/r/README.md @@ -1,30 +1,88 @@ # rushtools Small, tested helper functions for everyday plant bioinformatics in R. +Needs R 4.0 or newer; no dependencies beyond base R. ## Installation +From GitHub: + +```r +install.packages("remotes") # once +remotes::install_github("halireena/rushtools", subdir = "r") +``` + +Or from a clone of the repository, with R started in the top `rushtools` folder: + ```r -# install.packages("remotes") -remotes::install_github("yourGitHubName/rushtools") +install.packages("r", repos = NULL, type = "source") ``` ## Functions | Function | What it does | |---|---| -| `gc_content()` | GC fraction of DNA sequences | -| `reverse_complement()` | reverse complement of DNA sequences | +| `gc_content()` | GC fraction of DNA sequences (`NA` if a sequence has no A/C/G/T) | +| `reverse_complement()` | reverse complement of DNA sequences, IUPAC codes allowed | +| `count_motif()` | count motif occurrences in DNA sequences (e.g. a restriction site) | | `log2_fold_change()` | log2 fold change between two groups (raw or log2 input) | | `z_score()` | z-score standardisation | | `adjust_p()` | multiple-testing correction (BH by default) with input checks | | `tidy_de_table()` | standardise DESeq2 / limma / edgeR result tables | -| `count_motif()` | count motif occurrences in DNA sequences | + +Each has a help page with examples, e.g. `?count_motif`. Unfamiliar terms +(BH, IUPAC, fold change...) are explained in the +[glossary](../README.md#glossary). ## Example ```r library(rushtools) + gc_content(c("ATGC", "GGCC")) +#> [1] 0.5 1.0 + reverse_complement("ATGCCC") +#> [1] "GGGCAT" + +count_motif(c(s1 = "GAATTCGAATTC", s2 = "ATGC"), "GAATTC") # EcoRI site +#> s1 s2 +#> 2 0 + +log2_fold_change(c(400, 420, 380), c(50, 55, 45)) # about 8x higher +#> [1] 2.975033 + +z_score(c(2, 4, 6, 8)) +#> [1] -1.1618950 -0.3872983 0.3872983 1.1618950 + +adjust_p(c(a = 0.01, b = 0.02, c = 0.03, d = 0.5)) +#> a b c d +#> 0.04 0.04 0.04 0.50 + +res <- data.frame(gene = c("g1", "g2", "g3"), logFC = c(2.1, -1.5, 0.2), + P.Value = c(0.0001, 0.003, 0.7)) +tidy_de_table(res, lfc_col = "logFC", p_col = "P.Value") +#> gene log2fc pvalue padj direction +#> 1 g1 2.1 1e-04 0.0003 up +#> 2 g2 -1.5 3e-03 0.0045 down +#> 3 g3 0.2 7e-01 0.7000 ns +``` + +For a DESeq2 result the default column names already match. DESeq2 keeps +gene IDs in the row names, which `tidy_de_table()` resets, so move them into +a column first: + +```r +de <- as.data.frame(DESeq2::results(dds)) +tidy_de_table(cbind(gene = rownames(de), de), padj_col = "padj") ``` + +## Tests + +From this `r/` folder (needs `install.packages("testthat")` once): + +```bash +Rscript -e 'testthat::test_local(".")' +``` + +Changes between versions: [NEWS.md](NEWS.md). From 6b424f8914986ffb805740e0e39555cf732c98e9 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 5 Oct 2026 09:16:32 +0000 Subject: [PATCH 4/4] Make the test and example commands copy-pasteable from the top folder Spell out the `cd python` / `cd r` / `cd ..` steps instead of "from the X folder", so each block runs as written straight after the previous one; note where to get the example file if you installed without cloning. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WqrAUWjyEmFz36zphqfoBs --- README.md | 8 ++++++-- python/README.md | 11 +++++++++-- r/README.md | 3 ++- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 637cc54..9b679ca 100644 --- a/README.md +++ b/README.md @@ -185,17 +185,21 @@ page, e.g. `?tidy_de_table`. Full list: [`r/README.md`](r/README.md). ## Run the tests locally -Python (from the `python/` folder, with the virtual environment active): +Python (from the top `rushtools` folder, with the virtual environment active): ```bash +cd python pytest # unit tests + the examples in the docstrings ruff check . # lint (style and common-mistake checks) +cd .. ``` -R (from the `r/` folder; needs `install.packages("testthat")` once): +R (needs `install.packages("testthat")` once): ```bash +cd r Rscript -e 'testthat::test_local(".")' +cd .. ``` or, if you use devtools, `devtools::test()` inside R. diff --git a/python/README.md b/python/README.md index 2268791..42a400d 100644 --- a/python/README.md +++ b/python/README.md @@ -25,8 +25,15 @@ ruff check . # lint ``` The examples below use the tiny made-up files in -[`../examples/`](../examples/) and are run from the top `rushtools` folder -(`cd ..` if you are in `python/`). New words (ORF, N50, Phred...) are explained in the +[`../examples/`](../examples/) and are run from the top `rushtools` folder: + +```bash +cd .. # from rushtools/python back to rushtools, where examples/ lives +``` + +(If you installed without cloning, download +[`examples/genes.fasta`](../examples/genes.fasta) or use your own FASTA file.) +New words (ORF, N50, Phred...) are explained in the [glossary](../README.md#glossary). ## Use in Python diff --git a/r/README.md b/r/README.md index 0f01756..b0cf211 100644 --- a/r/README.md +++ b/r/README.md @@ -79,9 +79,10 @@ tidy_de_table(cbind(gene = rownames(de), de), padj_col = "padj") ## Tests -From this `r/` folder (needs `install.packages("testthat")` once): +From the top `rushtools` folder (needs `install.packages("testthat")` once): ```bash +cd r Rscript -e 'testthat::test_local(".")' ```