diff --git a/.gitignore b/.gitignore index 86feca0..df2223b 100644 --- a/.gitignore +++ b/.gitignore @@ -25,3 +25,6 @@ share/python-wheels/ .installed.cfg *.egg MANIFEST + +# mkdocs build output (regenerate with: docs/build_docs.sh) +site/ diff --git a/README.md b/README.md index b072dfe..eb39d36 100755 --- a/README.md +++ b/README.md @@ -52,11 +52,16 @@ Easiest is to check [datasets](datasets) examples to see how the above files loo ## Documentation -* [PDF reference manual](https://github.com/bedapub/splicekit/raw/main/docs/splicekit_docs.pdf) -* [Google docs](https://docs.google.com/document/d/15ZRCeK8xyg3klLktZSHZ9k__Xw_BZRn_Q-J4W35JNnQ/edit?usp=sharing) of the above PDF (comment if you like) +Full documentation, including installation, a quick start guide, and reference pages for configuration, sample annotation, features, edgeR, motif/scanRBP analysis, juDGE plots, additional analyses, JBrowse2 and the command line, is available at: + +* [bedapub.github.io/splicekit](https://bedapub.github.io/splicekit/) ## Changelog +**Docs**: released in September 2026 + +* migrated documentation to a mkdocs-material site ([bedapub.github.io/splicekit](https://bedapub.github.io/splicekit/)), retiring the PDF/Google Docs manual + **v0.8.1**: released in July 2026 * added `bam_file` column support in `samples.tab` for per-sample BAM paths (subfolder layouts) diff --git a/docs/build_docs.sh b/docs/build_docs.sh new file mode 100755 index 0000000..dae0f87 --- /dev/null +++ b/docs/build_docs.sh @@ -0,0 +1,22 @@ +#!/bin/bash +# Builds the mkdocs documentation site (source: docs/, config: ../mkdocs.yml) +# into docs/site/. Uses mkdocs from the "pybio" micromamba environment +# directly by path, so this works whether or not that environment is +# currently activated. Run from anywhere -- it cd's to the repo root itself. +# +# Usage: +# docs/build_docs.sh # build the static site into docs/site/ +# docs/build_docs.sh serve # serve locally with live-reload for editing + +set -e + +MKDOCS=/home/gregor/micromamba/envs/pybio/bin/mkdocs +# this script lives in docs/, but mkdocs.yml lives one level up at the repo root +cd "$(dirname "$0")/.." + +if [ "$1" == "serve" ]; then + "$MKDOCS" serve +else + "$MKDOCS" build --strict + echo "Built docs/site/ -- open docs/site/index.html or run 'docs/build_docs.sh serve' to preview locally." +fi diff --git a/docs/pages/about.md b/docs/pages/about.md new file mode 100644 index 0000000..6e0ac0d --- /dev/null +++ b/docs/pages/about.md @@ -0,0 +1,15 @@ +# Authors + +[splicekit](https://github.com/bedapub/splicekit) is developed by [Gregor Rot](https://grexor.github.io/) and contributors at BEDA/Roche, together with [pybio](https://grexor.github.io/pybio/) and [scanRBP](https://grexor.github.io/scanRBP/). + +# Citing + +If you find **splicekit** useful in your work and research, please cite: + +Rot, G., Wehling, A., Schmucki, R., Berntenis, N., Zhang, J. D., & Ebeling, M. (2024)
+[splicekit: an integrative toolkit for splicing analysis from short-read RNA-seq](https://academic.oup.com/bioinformaticsadvances/article/4/1/vbae121/7735317)
+Bioinformatics Advances, 4(1). [doi:10.1093/bioadv/vbae121](https://doi.org/10.1093/bioadv/vbae121) + +# Issues and Suggestions + +Use the [GitHub repository issues page](https://github.com/bedapub/splicekit/issues) to report issues and leave suggestions, or write to [Gregor Rot](mailto:gregor.rot@gmail.com). diff --git a/docs/pages/additional-analyses.md b/docs/pages/additional-analyses.md new file mode 100644 index 0000000..2b2935f --- /dev/null +++ b/docs/pages/additional-analyses.md @@ -0,0 +1,60 @@ +# Additional analyses + +Beyond the core edgeR / motif / juDGE pipeline, splicekit ships several further analyses, each runnable as its own `splicekit` command and included in `splicekit process`. + +## Promiscuity: `splicekit promisc` + +```bash +splicekit promisc +``` + +Promiscuity analysis looks at genes and how many of their junctions change across conditions — genes with many independently-regulated junctions ("promiscuous" splicing) versus genes regulated through a single dominant junction. + +## Cluster of pairwise logFC: `splicekit clusterlogfc` + +```bash +splicekit clusterlogfc +``` + +Clusters comparisons by the logFC of their significant (FDR < 0.05) changes, computed separately at the junction, exon and gene level. This groups comparisons/treatments that produce similar splicing/expression signatures. + +## JUNE: junction-event analysis + +```bash +splicekit june +``` + +JUNE (junction-events) classifies pairs/sets of overlapping junctions into splicing event types — e.g. skipped exons and mutually exclusive exons — by comparing their shared and differing donor/acceptor coordinates, rather than looking at single junctions or exons in isolation. + +## rMATS + +```bash +splicekit rmats +``` + +Runs [rMATS](http://rnaseq-mats.sourceforge.net/) on the project's BAM files as an independent, complementary splicing-event caller (skipped exons, alternative 5'/3' splice sites, mutually exclusive exons, retained introns) alongside splicekit's own junction/exon-based analysis. + +## DEXSeq + +```bash +splicekit dexseq # all feature types +splicekit dexseq junctions # a single feature type +``` + +An alternative to `splicekit edgeR` for differential feature usage, using [DEXSeq](https://bioconductor.org/packages/release/bioc/html/DEXSeq.html) instead of edgeR. Requires `dexseq_scripts` to be set in `splicekit.config` to the path of the DEXSeq R scripts. + +## DonJuAn: `splicekit juan` + +```bash +splicekit juan +``` + +Merges donor/acceptor anchor edgeR results back into the junction results (`results/results_edgeR_junctions.tab`), so each junction's row also reports its anchors' differential usage statistics. Runs automatically as part of `splicekit process`, right after `splicekit edgeR`. + +## HTML report: `splicekit report` + +```bash +splicekit report +``` + +Assembles the project's results (edgeR, JUNE, and more) into a browsable HTML report under `report/`, served together with [JBrowse2](jbrowse2.md) when you run `splicekit web`. diff --git a/docs/pages/assets/judge.png b/docs/pages/assets/judge.png new file mode 100755 index 0000000..f0f46d7 Binary files /dev/null and b/docs/pages/assets/judge.png differ diff --git a/docs/pages/assets/splicekit_logo.png b/docs/pages/assets/splicekit_logo.png new file mode 100755 index 0000000..d5cae4c Binary files /dev/null and b/docs/pages/assets/splicekit_logo.png differ diff --git a/docs/pages/cli.md b/docs/pages/cli.md new file mode 100644 index 0000000..a75ff64 --- /dev/null +++ b/docs/pages/cli.md @@ -0,0 +1,62 @@ +# Command-line reference + +Every `splicekit` invocation prints its version first, then runs the requested command. Run `splicekit -help` for the built-in summary; this page is the fuller reference. + +## Processing + +| Command | Description | +|---|---| +| `splicekit process` | Run all available analyses in order: setup, annotation, features, edgeR, juan, judge, motifs, promisc, clusterlogfc, june, report, jbrowse2. | + +## Initialization + +| Command | Description | +|---|---| +| `splicekit setup` | Initialize the project folder structure. | +| `splicekit annotation` | Load `samples.tab` and build comparisons. See [Sample annotation](samples.md). | +| `splicekit features` | Create junction/exon/anchor/gene count tables from BAM files. See [Features & count tables](features.md). | + +## Splicing analyses + +| Command | Description | +|---|---| +| `splicekit edgeR [junctions\|exons\|anchors\|genes]` | Run edgeR differential usage analysis. All feature types if no sub-command is given. See [Differential splicing (edgeR)](edgeR.md). | +| `splicekit dexseq [junctions\|exons\|anchors\|genes]` | Run DEXSeq as an alternative to edgeR. See [Additional analyses](additional-analyses.md#dexseq). | +| `splicekit juan` | Merge donor/acceptor anchor edgeR results into the junction results. | +| `splicekit judge` | Generate juDGE plots (junction logFC vs. gene logFC). See [juDGE plots](judge.md). | +| `splicekit promisc` | Promiscuity analysis: genes and their junction changes across conditions. | +| `splicekit clusterlogfc` | Cluster comparisons by logFC of significant (FDR < 0.05) changes, at junction/exon/gene level. | +| `splicekit june` | JUNE junction-event analysis (e.g. skipped/mutually exclusive exons). | +| `splicekit rmats` | Process the project's BAM files with rMATS. | +| `splicekit cassettes` | Cassette exon analysis. | +| `splicekit patterns` | Sequence pattern analysis. | + +## Motif & scanRBP + +| Command | Description | +|---|---| +| `splicekit motifs [dreme\|scanrbp]` | Run motif logos, DREME and scanRBP analysis. All three if no sub-command is given. See [Motif & RNA-binding analysis](motifs.md). | + +## JBrowse2 & reporting + +| Command | Description | +|---|---| +| `splicekit jbrowse2 [process\|start]` | Process (build) JBrowse2 files and/or start the JBrowse2 web server. Both steps if no sub-command is given. | +| `splicekit jbrowse` | Alias for `splicekit jbrowse2`. | +| `splicekit web` | Start the local web server serving both the HTML report and JBrowse2. See [Exploring results](jbrowse2.md). | +| `splicekit report` | Generate the HTML report under `report/`. | + +## Other + +| Command | Description | +|---|---| +| `splicekit version` | Print the installed splicekit version. | + +## Global options + +| Option | Description | +|---|---| +| `-help` | Print the built-in usage summary (or a sub-command's usage, e.g. `splicekit edgeR -help`). | +| `-version` | Print the installed version and exit. | +| `-verbose` | Print more detailed progress output. | +| `-force` | Force recomputation of steps that would otherwise be skipped if their output already exists (used by `process` and `jbrowse2`). | diff --git a/docs/pages/configuration.md b/docs/pages/configuration.md new file mode 100644 index 0000000..b1f9a5b --- /dev/null +++ b/docs/pages/configuration.md @@ -0,0 +1,136 @@ +# Configuration + +splicekit reads parameters from two files in your project folder: + +- **`splicekit.config`** — a Python file (each line is `exec`'d) describing the study, sample annotation, genome and analysis parameters. Start from the [template](https://github.com/bedapub/splicekit/blob/main/splicekit/splicekit.config.template) or one of the [datasets](https://github.com/bedapub/splicekit/tree/main/datasets) examples. +- **`config.yaml`** — a Snakemake config file describing per-rule compute resources (cores/memory/time). Start from the [template config.yaml](https://github.com/bedapub/splicekit/blob/main/config.yaml). + +## splicekit.config + +### Core parameters + +**`study_name`**, default = `"descriptive short study name"` +Descriptive short study name, arbitrary string describing the study. + +**`library_type`**, default = `"paired-end"` +Either `"paired-end"` or `"single-end"`. + +**`library_strand`**, default = `"NONE"` +Possible values: + +- `"SECOND_READ_TRANSCRIPTION_STRAND"` +- `"FIRST_READ_TRANSCRIPTION_STRAND"` +- `"SINGLE_STRAND"` +- `"SINGLE_REVERSE"` +- `"NONE"` + +For unstranded data, use `"NONE"`. For paired-end stranded data, the most common value is `"SECOND_READ_TRANSCRIPTION_STRAND"`, meaning the second read of the pair maps in the transcript direction and the first read maps in the reverse direction. For stranded single-end sequencing, `"SINGLE_STRAND"` means the reads map in the transcript direction, and `"SINGLE_REVERSE"` means the reads map on the opposite strand of the transcripts. + +**`edgeR_FDR_thr`**, default = `0.05` +FDR threshold used to filter `results/results_edgeR_{feature_type}.tab`. See [Differential splicing (edgeR)](edgeR.md). + +**`dexseq_scripts`**, default = `""` +Path to DEXSeq R scripts, if you want to run `splicekit dexseq` as an alternative to `splicekit edgeR`. + +### Sample annotation parameters + +**`sample_column`**, default = `"sample_id"` +Which column in `samples.tab` holds the sample IDs. BAM files are then expected to be named `sample_id.bam` (see `bam_path`/`bam_column` below). + +**`treatment_column`**, default = `"treatment_id"` +Which column in `samples.tab` defines treatment (test and control labels). + +**`control_name`**, default = `"DMSO"` +The text identifying control samples in the treatment column. Non-control samples are compared against these. + +**`separate_column`**, default = `""` +Group samples by this column and only compare within a group (empty = don't separate). Useful when, say, samples come from two tissues and you only want within-tissue comparisons: add a tissue column to `samples.tab` and set `separate_column` to its name. + +**`group_column`**, default = `""` +Only compare a test sample against controls from the same "domain" (e.g. the same sequencing plate). + +### Genome parameters + +**`species`**, default = `None` +Genome species, resolved through `pybio`/Ensembl. Any species pybio knows about, or a custom species you registered with pybio via a FASTA+GTF pair. + +**`genome_version`**, default = `None` +Ensembl genome version, or your custom genome version. `None` takes the latest Ensembl release. + +**`bam_path`**, default = `""` +Folder where BAM files are stored, expected to contain `{sample_id}.bam` for every sample. Absolute, or relative to the project folder. + +**`bam_column`**, default = `"bam_file"` +Alternative to `bam_path`: the name of a column in `samples.tab` holding the full BAM path for each sample individually (useful when BAMs live in per-sample subfolders rather than one flat folder). At least one of `bam_path` or `bam_column` must resolve to something — splicekit exits with an error otherwise. + +### scanRBP parameters + +**`scanRBP`**, default = `True` +Should splicekit run the scanRBP RNA-protein binding analysis as part of `splicekit motifs`? See [Motif & RNA-binding analysis](motifs.md). + +**`protein`**, default = `"K562.TARDBP.0"` +ID of the protein PWM to scan regulated sequences against. List available IDs with `scanRBP search `, e.g.: + +``` +$ scanRBP search hnRNPA +scan_id protein tissue description +HNRNPA1.HepG2.00 HNRNPA1 HepG2 heterogeneous nuclear ribonucleoprotein A1 +HNRNPA1.HepG2.01 HNRNPA1 HepG2 heterogeneous nuclear ribonucleoprotein A1 +... +HNRNPA1.K562.00 HNRNPA1 K562 heterogeneous nuclear ribonucleoprotein A1 +``` + +**`protein_label`**, default = `"tdp43"` +Short descriptive label for the protein, used in file names and plot titles. + +### Processing parameters + +**`platform`**, default = `"desktop"` +`"desktop"`, `"cluster"` (LSF/`bsub`) or `"SLURM"`. When running under Snakemake, job submission is instead handled by `run_snakemake_local.sh`/`run_snakemake_slurm.sh` and `config.yaml`. + +**`container`**, default = `""` +Empty string assumes all software dependencies are installed on the machine/cluster. Set to `"singularity run docker://ghcr.io/bedapub/splicekit:main"` to run non-Python steps from the provided container image instead (pybio and scanRBP are installed via pip regardless). + +**`edgeR_memory`**, default = `"16GB"` +Memory reserved for edgeR cluster jobs (only relevant on `"cluster"`/`"SLURM"` platforms; under Snakemake, use `config.yaml` instead). + +### Visualization and labeling parameters + +**`short_names`**, default = `[]` +By default, no names are shortened or replaced in results or plots. Provide a list of triples to replace/shorten long strings: + +```python +# example splicekit.config short_names parameter +# replace cell_line_A with A, cell_line_B with B (only on an exact/"complete" match) +short_names = [("cell_line_A", "A", "complete"), ("cell_line_B", "B", "complete")] +``` + +Use `"partial"` instead of `"complete"` to also replace the string when it occurs inside a larger one (e.g. `...cell_line_A...` → `...A...`). + +## config.yaml + +`config.yaml` sizes the Snakemake rules — how many cores, how much memory and how much walltime each step gets, whether running locally or submitting to SLURM. + +```yaml +mapping: + perform_mapping: True # should splicekit map FASTQ -> BAM with pybio/STAR? (True/False) + alignIntronMax: 0 # STAR max intron size; 0 lets STAR derive it from its default window parameters + +defaults: + cores: 1 + mem: 4 # GB + time: "01:00:00" + +map_fastq_single: { cores: 8, mem: 8, time: "02:00:00" } +map_fastq_paired: { cores: 8, mem: 8, time: "02:00:00" } +bam_index: { cores: 8, mem: 2, time: "02:00:00" } +bam_bw: { cores: 8, mem: 2, time: "02:00:00" } +feature_counts: { cores: 8, mem: 2, time: "01:00:00" } +edgeR: { cores: 1, mem: 4, time: "04:00:00" } +edgeR_assemble: { cores: 1, mem: 16, time: "04:00:00" } +juan: { cores: 1, mem: 16, time: "01:00:00" } +``` + +Any rule without its own section falls back to `defaults`. Set `mapping.perform_mapping: False` if you're supplying your own BAM files and don't want splicekit to align FASTQs itself. + +After setting up `splicekit.config` and `config.yaml`, run the pipeline as described in [Quick Start](quickstart.md). diff --git a/docs/pages/coordinates.md b/docs/pages/coordinates.md new file mode 100644 index 0000000..d6092bb --- /dev/null +++ b/docs/pages/coordinates.md @@ -0,0 +1,13 @@ +# Genomic coordinates + +All genomic coordinates splicekit operates with are **0-based, left+right inclusive**. E.g. the range 100-103 includes coordinates 100, 101, 102 and 103; the first coordinate is 0. + +More specifically: + +- **Feature-specific**: all feature coordinates (junctions, anchors, exons) are given in numeric sort order regardless of strand — `feature_start` is always `<` `feature_stop`. Example: `chr1+_100_102` represents a junction spanning coordinates `[100, 101, 102]`; `chr1-_100_102` represents a junction spanning coordinates `[102, 101, 100]`. +- **Junction-specific**: junction coordinates cover/overlap 1 nucleotide of each adjoining exon. + - `chr1+_100_200` — a junction on chromosome 1 (`+` strand) from `[100..200]`. 100 is the last nucleotide of the upstream exon, 200 is the first nucleotide of the downstream exon. + - `chr1-_100_200` — a junction on chromosome 1 (`-` strand) from `[100..200]`. 200 is the last nucleotide of the upstream exon, 100 is the first nucleotide of the downstream exon. + +!!! important + RefSeq and Ensembl GTF files are 1-indexed. When splicekit reads files from RefSeq/Ensembl, it performs `coordinate -= 1` on every coordinate to keep them consistent with splicekit's internal 0-indexed structures. diff --git a/docs/pages/dependencies.md b/docs/pages/dependencies.md new file mode 100644 index 0000000..109008f --- /dev/null +++ b/docs/pages/dependencies.md @@ -0,0 +1,21 @@ +# Dependencies + +## Conda/micromamba environment (`splicekit.yaml`) + +Installed by `micromamba create -y -f splicekit.yaml`: + +[pigz](https://zlib.net/pigz/), [deeptools](https://deeptools.readthedocs.io/), [samtools](http://www.htslib.org/), [Snakemake](https://snakemake.readthedocs.io/), R (`r-base`, `r-locfit`), [STAR](https://github.com/alexdobin/STAR), [rMATS](http://rnaseq-mats.sourceforge.net/), [Node.js](https://nodejs.org/), Ghostscript, [subread](http://subread.sourceforge.net/) (`featureCounts`, > 2.0.6), [MEME](https://meme-suite.org/meme/) (> 5.5.1, for DREME), Perl's `cpanminus`, and the `snakemake-executor-plugin-cluster-generic` pip package (used by `run_snakemake_slurm.sh` for SLURM submission). + +## Installed by `install.sh` + +- R packages: `BiocManager`, `data.table`, `statmod`, `R.utils`, and Bioconductor's [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html). +- Perl modules for the JBrowse2/SOAP tooling (`XML::Compile::*`, `Log::Log4perl`, `Math::CDF`, `JSON`, and others). +- [`@jbrowse/cli`](https://www.npmjs.com/package/@jbrowse/cli) via `npm install -g`, used to set up the local [JBrowse2](jbrowse2.md) instance. + +## Python dependencies (installed by `pip install .`) + +[Levenshtein](https://pypi.org/project/Levenshtein/), [logomaker](https://logomaker.readthedocs.io/), [plotly](https://plotly.com/python/), `python-dateutil`, [pybio](https://grexor.github.io/pybio/), [scanRBP](https://grexor.github.io/scanRBP/), [pysam](https://pysam.readthedocs.io/), [numpy](https://numpy.org/), `psutil`, `beautifulsoup4`, `requests` and `rangehttpserver`. + +## Optional + +[Singularity](https://sylabs.io/singularity/) — only needed if you set `container = "singularity run docker://ghcr.io/bedapub/splicekit:main"` in `splicekit.config` instead of installing the conda environment directly. See [Installation: Container](installation.md#container). diff --git a/docs/pages/edgeR.md b/docs/pages/edgeR.md new file mode 100644 index 0000000..4432ef1 --- /dev/null +++ b/docs/pages/edgeR.md @@ -0,0 +1,60 @@ +# Differential splicing (edgeR) + +Running edgeR analysis on features (junctions, anchors, exons, genes) is a single command: + +```bash +splicekit edgeR # all feature types +splicekit edgeR junctions # a single feature type +splicekit edgeR exons +splicekit edgeR anchors +splicekit edgeR genes +``` + +Donor/acceptor anchor results are then merged back into the corresponding junction results by `splicekit juan` (part of `splicekit process`), so a junction's row also carries its anchors' edgeR statistics. + +## Results files + +Results are stored in `results/results_edgeR_{feature_type}.tab`, where `feature_type` is one of `genes`, `exons`, `junctions`, `donor_anchors`, `acceptor_anchors`. Only results with `FDR < splicekit.config.edgeR_FDR_thr` are reported (sorted by FDR), each linked to [JBrowse2](jbrowse2.md) via a URL. + +To explore all results without the FDR filter, use `results/results_edgeR_{feature_type}_all.tab`. + +### General columns + +| Column | Example | Description | +|---|---|---| +| `result_id` | `r1` | Integer result identifier, starting at 1. | +| `comparison` | `test_control` | Comparison name, from `annotation/comparisons.tab`. | +| `compound` | `treatment1` | Name of the treatment/compound tested. | +| `feature_id` | `chr1+_17741_17839` | ID of the reported feature: a gene/exon/junction/`[donor,acceptor]_anchor` ID. | +| `chr` | `1` | Chromosome of the feature. | +| `strand` | `+` | Strand of the feature (`+` or `-`). | +| `feature_start` | `17741` | Start of the feature (numerically, start < stop). See [Genomic coordinates](coordinates.md). | +| `feature_stop` | `17839` | Stop of the feature (numerically, stop > start). | +| `feature_length` | `250` | `feature_stop - feature_start + 1`. | +| `gene_id` | `ENSG00000120948` | Ensembl or RefSeq gene ID. | +| `gene_name` | `TARDBP` | Corresponding to `gene_id`. | +| `sum_feature_test` | `1000` | Sum of counts for this feature across all test samples. | +| `sum_feature_control` | `1000` | Sum of counts for this feature across all control samples. | +| `jbrowse_loc` | `3:342321..351243` | Genomic region shown in the JBrowse view. | +| `jbrowse_url` | | Link to the JBrowse view. | +| `logFC` | | Log fold change, from edgeR. | +| `exon.F` | | `exon.F` statistic, from edgeR. | +| `p_value` | | p-value, from edgeR. | +| `fdr` | | False discovery rate, from edgeR. | + +### Junction-specific (additional) columns + +| Column | Example | Description | +|---|---|---| +| `annotated` | `AA` | Two-letter code `AA`/`AN`/`NA`/`NN`: first letter for the donor site (5' of junction), second for the acceptor site (3' of junction); `A` = touches an annotated exon, `N` = does not. | +| `donor_anchor_id` | | ID of the donor anchor linked to this junction. | +| `acceptor_anchor_id` | | ID of the acceptor anchor linked to this junction. | +| `UTR` | | `first_exon_{start_pos}` if the junction touches any transcript's first exon of the gene. | + +### Exon-specific (additional) columns + +| Column | Example | Description | +|---|---|---| +| `delta_PSI` | | `test_PSI - control_PSI` (percentage spliced-in). | + +Next: [Motif & RNA-binding analysis](motifs.md) runs on the sequences around the regulated features found here. diff --git a/docs/pages/features.md b/docs/pages/features.md new file mode 100644 index 0000000..dceb33d --- /dev/null +++ b/docs/pages/features.md @@ -0,0 +1,24 @@ +# Features & count tables + +Running `splicekit features` creates count tables for junctions, anchors, exons and genes from the project's BAM files. + +## What are features? + +splicekit operates on 4 types of features: **junctions, anchors, exons and genes**. All feature IDs share the same format: `chrstrand_start_stop`, e.g. `chrX-_154371360_154374505`. See [Genomic coordinates](coordinates.md) for how coordinates are reported across splicekit. + +Junctions are detected directly from BAM files (independent of any pre-existing gene model), and reported in `reference/junctions.tab` together with a *donor anchor* and *acceptor anchor* — by default the 15nt regions flanking the junction's start and stop. These anchor regions are turned into `reference/donor_anchors.gtf` / `reference/acceptor_anchors.gtf` and quantified with **featureCounts**, alongside exon- and gene-level counts. See [File formats](fileformats.md) for the exact column layouts. + +## Feature data files + +Each individual sample gets one file (table) per feature type, under `data/sample_{feature_type}_data/`, listing every feature and its count in that sample. + +Example: `data/sample_exons_data/sample_99.tab` + +``` +GeneID Start End Length Symbol 1_test 2_test 3_control 4_control +1 58347029 58347353 325 A1BG 42 31 109 75 +1 58347640 58350370 2731 A1BG 0 0 3 1 +1 58350651 58351391 741 A1BG 0 0 10 1 +``` + +Next step: [Differential splicing (edgeR)](edgeR.md), which turns these per-sample count tables into per-comparison differential usage results. diff --git a/docs/pages/fileformats.md b/docs/pages/fileformats.md new file mode 100644 index 0000000..e2a0091 --- /dev/null +++ b/docs/pages/fileformats.md @@ -0,0 +1,33 @@ +# File formats + +## `reference/junctions.tab` + +Contains all junctions detected across every sample in the project. Only junctions that could be annotated to a gene are reported — including "novel" junctions that don't touch a RefSeq/Ensembl-annotated exon, as long as the junction's start and stop fall inside an annotated gene (see the `annotated` column). + +| Column | Example | Description | +|---|---|---| +| `junction_id` | `chr1+_17741_17839` | Unique ID: `chrstrand_start_stop`. | +| `donor_anchor_id` | `chr1+_17725_17740` | Matching donor anchor ID — by default the 15nt region upstream of the junction start. | +| `acceptor_anchor_id` | `chr1+_17840_17855` | Matching acceptor anchor ID — by default the 15nt region downstream of the junction stop. | +| `gene_id` | `ENSG00000120948` | Ensembl or RefSeq gene ID. A junction can be non-annotated (`annotated != "AA"`) but still assigned to a gene, meaning its start/stop fall inside the gene. | +| `gene_name` | `TARDBP` | Corresponding to `gene_id`. | +| `chr` | `1` | Chromosome. | +| `strand` | `+` | `+` or `-`. | +| `annotated` | `AA` | Two-letter code `AA`/`AN`/`NA`/`NN` — see [Genomic coordinates](coordinates.md) and [edgeR results](edgeR.md#junction-specific-additional-columns). | +| `count` | `553` | Raw read count across all samples in the project supporting this junction. | + +## `reference/donor_anchors.gtf` and `reference/acceptor_anchors.gtf` + +GTF files generated from all donor/acceptor anchors in `reference/junctions.tab`. Used by **featureCounts** to build anchor count tables across the project's samples. + +## `results/results_edgeR_{feature_type}.tab` + +See [Differential splicing (edgeR): Results files](edgeR.md#results-files) for the full column reference (general columns, plus junction- and exon-specific additions). + +## `data/sample_{feature_type}_data/*.tab` + +See [Features & count tables: Feature data files](features.md#feature-data-files). + +## `annotation/comparisons.tab` + +See [Sample annotation: Comparisons](samples.md#comparisons). diff --git a/docs/pages/index.md b/docs/pages/index.md new file mode 100644 index 0000000..b46cb8d --- /dev/null +++ b/docs/pages/index.md @@ -0,0 +1,30 @@ +# splicekit: splicing analysis from short-read RNA-seq + +![splicekit](assets/splicekit_logo.png) + +**splicekit** is a modular, integrative platform for splicing analysis of short-read RNA-seq data. Starting from aligned reads (BAM files) and a sample annotation table, it defines test-vs-control comparisons, builds per-feature count tables (junctions, anchors, exons, genes), and runs a battery of splicing analyses — differential feature usage with edgeR, motif and RNA-protein binding enrichment, junction-vs-gene expression comparisons and more — all self-contained in a single project folder. It integrates [pybio](https://github.com/grexor/pybio) for genome operations, [scanRBP](https://github.com/grexor/scanRBP) for RNA-protein binding, and ships its own [JBrowse2](https://jbrowse.org/jb2/) instance for browsing results. + +```bash +# create and activate the conda/micromamba environment, install splicekit +micromamba create -y -f splicekit.yaml +micromamba activate splicekit +./install.sh +pip install . + +# run the full pipeline with Snakemake +./run_snakemake_local.sh --configfile config.yaml +``` + +## What's included + +- **Comparisons from a sample sheet** — define test/control comparisons straight from `samples.tab`, optionally grouped or separated by extra columns. See [Sample annotation](samples.md). +- **Feature count tables** — junctions, anchors, exons and genes, built from BAM files. See [Features & count tables](features.md). +- **Differential splicing with edgeR** — per-feature differential usage results, linked directly to JBrowse2. See [Differential splicing (edgeR)](edgeR.md). +- **Motif & RNA-protein binding analysis** — donor/acceptor motif logos, DREME enrichment and scanRBP binding analysis. See [Motif & RNA-binding analysis](motifs.md). +- **juDGE plots** — junction logFC vs. gene logFC, distinguishing splicing modifiers from expression modifiers. See [juDGE plots](judge.md). +- **Promiscuity, clustering, JUNE and rMATS analyses** — see [Additional analyses](additional-analyses.md). +- **Integrated JBrowse2 + HTML report** — one local web server for both. See [Exploring results](jbrowse2.md). + +## Where to start + +New to splicekit? Read [Installation](installation.md) and then [Quick Start](quickstart.md) — together they take you from a fresh checkout to a running pipeline on example data in a few commands. Everything else in these docs is reference material for the individual analysis steps, configuration parameters and file formats. diff --git a/docs/pages/installation.md b/docs/pages/installation.md new file mode 100644 index 0000000..3fc3407 --- /dev/null +++ b/docs/pages/installation.md @@ -0,0 +1,45 @@ +# Installation + +Since v0.7, **splicekit** is a [Snakemake](https://snakemake.readthedocs.io/) pipeline with a Conda/micromamba environment. + +```bash +# clone the repository +git clone git@github.com:bedapub/splicekit.git +cd splicekit + +# create and activate the conda environment +micromamba create -y -f splicekit.yaml +micromamba activate splicekit + +# install remaining (non-conda) dependencies: R/edgeR, Perl modules, jbrowse-cli +./install.sh + +# install splicekit itself +pip install . +``` + +!!! note + `install.sh` installs R packages (`edgeR` via BiocManager), Perl SOAP/XML modules, and `@jbrowse/cli` via npm — these aren't packaged as conda dependencies, so run it once per environment. + +The `splicekit.yaml` environment brings in Snakemake, STAR, samtools, subread (`featureCounts`), MEME (for DREME), rMATS and the `snakemake-executor-plugin-cluster-generic` plugin used for SLURM submission. See [Dependencies](dependencies.md) for the full list. + +## Installing just the Python package + +If you already have the environment's tools on your `PATH` (or are only using splicekit's Python API / running individual `splicekit` CLI steps by hand), you can install the package on its own: + +```bash +pip install splicekit +``` + +or, from this repository directly: + +```bash +pip install git+https://github.com/bedapub/splicekit.git@main +``` + +!!! note + On some systems, **pip** installs the executable scripts under `~/.local/bin`. If this folder is not in your `PATH`, running `splicekit` will fail with `command not found`. Fix this with `export PATH="$PATH:~/.local/bin"` (add it to your `~/.profile` to persist across logins). Another option is to install inside a virtual environment (using [virtualenv](https://virtualenv.pypa.io/en/latest/)). + +## Container + +If you'd rather not install dependencies directly on the machine or cluster, set `container = "singularity run docker://ghcr.io/bedapub/splicekit:main"` in `splicekit.config`. splicekit will then run its non-Python steps through that imported Docker image (pybio and scanRBP are already installed as regular pip dependencies regardless of the container setting). See [Configuration](configuration.md#processing-parameters). diff --git a/docs/pages/jbrowse2.md b/docs/pages/jbrowse2.md new file mode 100644 index 0000000..01f804c --- /dev/null +++ b/docs/pages/jbrowse2.md @@ -0,0 +1,27 @@ +# Exploring results: web report & JBrowse2 + +To graphically explore results, splicekit provides an integrated [JBrowse2](https://jbrowse.org/jb2/) instance alongside its HTML report, both served from the same local web server. + +## Preparing and starting + +```bash +splicekit jbrowse2 process # process JBrowse2 files (genome, BAM tracks, etc.) +splicekit jbrowse2 start # start the local web server + +# equivalent to running both steps above: +splicekit jbrowse2 + +# also starts the same web server, alongside the HTML report: +splicekit web +``` + +`splicekit process` (the full pipeline) also runs the `jbrowse2` steps automatically. + +`splicekit jbrowse2 process` downloads and unpacks a local JBrowse2 web build the first time it runs, indexes the reference genome FASTA, and generates per-sample BigWig/CRAM tracks plus junction BED tracks from the project's BAM files. + +`splicekit web` / `splicekit jbrowse2 start` then serve everything from a single local HTTP server: + +- `http://:8007/report` — the HTML report generated by `splicekit report`. +- `http://:8007/jbrowse2/?config=splicekit_data/config.json` — the JBrowse2 genome browser, pre-configured with the project's genome and BAM/BigWig tracks. + +The [edgeR results tables](edgeR.md#results-files) link directly into this JBrowse2 instance via their `jbrowse_url` column, so you can jump straight from a significant junction/exon/gene to its genomic context. diff --git a/docs/pages/judge.md b/docs/pages/judge.md new file mode 100644 index 0000000..481a566 --- /dev/null +++ b/docs/pages/judge.md @@ -0,0 +1,13 @@ +# juDGE plots + +To characterize the effect of a treatment — is it mostly driving splicing changes, or general gene expression changes? — splicekit generates **juDGE plots**: junction logFC vs. gene logFC, one point per gene. + +```bash +splicekit judge +``` + +Plots are written to `results/judge/*` as PNG images, plus interactive Plotly HTML reports (hover to see the data point and gene name). + +![juDGE plot example](assets/judge.png) + +In a plot like the one above, a comparison where junctions (y axis) are perturbed much more than gene expression in general (x axis) is characterized as a **"splicing modifier"**. A comparison with more activity on the x axis — general gene expression — is labeled an **"expression modifier"**. diff --git a/docs/pages/motifs.md b/docs/pages/motifs.md new file mode 100644 index 0000000..690077e --- /dev/null +++ b/docs/pages/motifs.md @@ -0,0 +1,26 @@ +# Motif & RNA-binding analysis + +`splicekit motifs` analyzes the sequences around regulated splicing events found by [edgeR](edgeR.md). + +```bash +splicekit motifs # run all: motif logos, DREME, scanRBP +splicekit motifs dreme # only DREME +``` + +Motif analysis on donor site patterns (9nt sequences) runs on the top 100 hits of each comparison, producing motif logos and HTML reports under `results/motifs`. In addition to the logos, splicekit runs [DREME](https://meme-suite.org/meme/doc/dreme.html) on regulated sequences vs. control sequences to find enriched short motifs. + +## scanRBP: RNA-protein binding enrichment + +As part of the same motif analysis, splicekit identifies potential enrichment of RNA-protein binding at regulated sites (donor sites, acceptor sites and other regions), using [scanRBP](https://github.com/grexor/scanRBP). + +Once sets of control and regulated sequences are identified, scanRBP computes the log-odds of the binding signal for a chosen protein from its PWM. Bootstrapping the sequence labels estimates the probability that binding at regulated sequences differs from binding at controls (a log-FC of the binding signal). + +Configure which protein to scan in `splicekit.config`: + +```python +scanRBP = True # run the scanRBP step? (True/False) +protein = "K562.TARDBP.0" # PWM id, see: scanRBP search +protein_label = "tdp43" # short label used in file names and titles +``` + +See [Configuration](configuration.md#scanrbp-parameters) for the full parameter list, and the [scanRBP documentation](https://grexor.github.io/scanRBP/) for the standalone tool (`pip install scanRBP`) and its motif database. diff --git a/docs/pages/quickstart.md b/docs/pages/quickstart.md new file mode 100644 index 0000000..918bbe9 --- /dev/null +++ b/docs/pages/quickstart.md @@ -0,0 +1,48 @@ +# Quick Start + +To run splicekit you need: + +1. **A reference genome**, downloaded and processed with `pybio` (installed automatically as a splicekit dependency): + ```bash + pybio genome homo_sapiens # or: pybio search species, for other species + ``` +2. **Aligned reads in BAM format**, one file per sample. You can align FASTQ files yourself with STAR, or reuse a dataset's mapping script, e.g. [datasets/GSE221868/2_map.sh](https://github.com/bedapub/splicekit/blob/main/datasets/GSE221868/2_map.sh), which downloads the reference genome with pybio and aligns with STAR. +3. **`samples.tab`** — one line per sample, TAB delimited, connecting each `sample_id` to its `treatment_id`. See [Sample annotation](samples.md) and the [example samples.tab](https://github.com/bedapub/splicekit/blob/main/datasets/GSE126543/samples.tab). +4. **`splicekit.config`** — reference genome, BAM folder and the other core parameters. See [Configuration](configuration.md) and the [example splicekit.config](https://github.com/bedapub/splicekit/blob/main/datasets/GSE126543/splicekit.config). +5. **`config.yaml`** — per-rule Snakemake resources (cores/memory/time). Copy the [template config.yaml](https://github.com/bedapub/splicekit/blob/main/config.yaml) into your project folder and adjust it to your cluster/machine. + +The [datasets](https://github.com/bedapub/splicekit/tree/main/datasets) folder has four complete examples, each with its own scripts to download and process a public RNA-seq dataset from scratch. + +## Running the pipeline + +With `samples.tab`, `splicekit.config` and `config.yaml` in your project folder, run the whole pipeline with Snakemake: + +```bash +cd datasets/GSE126543 # example project folder +./1_download.sh # download sample FASTQs +pybio homo_sapiens # reference genome + +./run_snakemake_local.sh --configfile config.yaml # run locally +# or: +./run_snakemake_slurm.sh --configfile config.yaml # submit jobs to SLURM +``` + +`run_snakemake_slurm.sh` submits each Snakemake rule as its own SLURM job (via `snakemake-executor-plugin-cluster-generic`), sized per-rule from `config.yaml`. + +Once it finishes, explore the results: + +```bash +splicekit web +``` + +This starts a single local web server serving both the HTML report (`http://:8007/report`) and the [JBrowse2](jbrowse2.md) genome browser. + +!!! note + If you already have BAM files and want to skip Snakemake, you can run splicekit directly with `splicekit process` inside a folder containing `samples.tab` and `splicekit.config` — this runs the same analysis steps sequentially on a single machine. See [Command-line reference](cli.md). + +## Next steps + +- [Configuration](configuration.md) — every `splicekit.config` and `config.yaml` parameter. +- [Sample annotation](samples.md) — how `samples.tab` becomes `annotation/comparisons.tab`. +- [Features & count tables](features.md) — the four feature types and their count files. +- [Differential splicing (edgeR)](edgeR.md) — running and reading edgeR results. diff --git a/docs/pages/samples.md b/docs/pages/samples.md new file mode 100644 index 0000000..ec2bf9d --- /dev/null +++ b/docs/pages/samples.md @@ -0,0 +1,45 @@ +# Sample annotation + +The first step of the analysis, `splicekit annotation`, loads samples from `samples.tab` and builds the test-vs-control comparisons that every later step (features, edgeR, motifs, ...) works from. + +!!! note + `samples.tab` is **TAB delimited**. Lines starting with `#` are treated as comments. + +``` +sample_id treatment_id +sample1 control +sample2 control +sample3 test1 +sample4 test1 +sample5 test2 +sample6 test2 +``` + +splicekit expects a BAM file per sample, resolved either from `bam_path` (as `{bam_path}/{sample_id}.bam`) or from a per-sample path in a `bam_column` column of `samples.tab` — see [Configuration](configuration.md#genome-parameters). The `sample_column`, `treatment_column`, `control_name`, `separate_column` and `group_column` parameters in `splicekit.config` control how `samples.tab` is read and how comparisons are formed. + +## Comparisons + +Each non-control treatment (which may have several replicate samples) is compared against the samples matching `control_name`. For the example above, this produces `annotation/comparisons.tab`: + +``` +comparison compound_samples DMSO_samples +test1_control sample3_test1,sample4_test1 sample1_control,sample2_control +test2_control sample5_test2,sample6_test2 sample1_control,sample2_control +``` + +In addition, `splicekit annotation` creates the processing shell scripts and (on `"cluster"`/`"SLURM"` platforms) cluster job files under `jobs/*`. An example cluster job file: + +```bash +#!/bin/bash +#BSUB -J edgeR_junctions_sample1 # job name +#BSUB -n 4 # number of tasks +#BSUB -R "span[hosts=1]" # 1 host +#BSUB -q short # select queue +#BSUB -o logs_edgeR_junctions/sample1_control.out # output file +#BSUB -e logs_edgeR_junctions/sample1_control.err # error file + +ml R +R --no-save --args splicekit comparison_junctions_data junctions control test ... < comps_edgeR.R +``` + +When running under Snakemake instead, job submission and resourcing come from `config.yaml` and `run_snakemake_slurm.sh` — see [Configuration](configuration.md#configyaml). diff --git a/docs/pages/stylesheets/extra.css b/docs/pages/stylesheets/extra.css new file mode 100644 index 0000000..b2233e4 --- /dev/null +++ b/docs/pages/stylesheets/extra.css @@ -0,0 +1,9 @@ +/* Custom brand blue, applied via theme.palette.primary: custom in mkdocs.yml -- + mkdocs-material's built-in "blue" swatch is a fixed named color and can't be + pointed at an arbitrary hex value directly, so this overrides the CSS custom + properties it reads the header/accent color from instead. */ +[data-md-color-primary="custom"] { + --md-primary-fg-color: #3c78d7; + --md-primary-fg-color--light: #3c78d7; + --md-primary-fg-color--dark: #3c78d7; +} diff --git a/docs/splicekit_docs.pdf b/docs/splicekit_docs.pdf deleted file mode 100755 index b7c189c..0000000 Binary files a/docs/splicekit_docs.pdf and /dev/null differ diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 0000000..752cbf6 --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,62 @@ +site_name: splicekit +site_description: splicekit - splicing analysis from short-read RNA-seq +site_url: https://bedapub.github.io/splicekit/ +repo_url: https://github.com/bedapub/splicekit +repo_name: bedapub/splicekit +docs_dir: docs/pages +site_dir: docs/site + +theme: + name: material + logo: assets/splicekit_logo.png + favicon: assets/splicekit_logo.png + palette: + # "custom" + extra_css below, rather than a built-in named color, so the + # header/accent can be the exact brand blue (#3c78d7) instead of whichever + # named material swatch happens to be closest. + - scheme: default + primary: custom + toggle: + icon: material/brightness-7 + name: Switch to dark mode + - scheme: slate + primary: custom + toggle: + icon: material/brightness-4 + name: Switch to light mode + features: + - navigation.tabs + - navigation.top + - content.code.copy + - search.suggest + +extra_css: + - stylesheets/extra.css + +markdown_extensions: + - admonition + - toc: + permalink: true + - pymdownx.highlight + - pymdownx.superfences + - pymdownx.inlinehilite + +nav: + - Home: index.md + - Installation: installation.md + - Quick Start: quickstart.md + - Guides: + - Configuration: configuration.md + - Sample annotation: samples.md + - Features & count tables: features.md + - Differential splicing (edgeR): edgeR.md + - Motif & RNA-binding analysis: motifs.md + - juDGE plots: judge.md + - Additional analyses: additional-analyses.md + - Exploring results: jbrowse2.md + - Reference: + - Command-line reference: cli.md + - File formats: fileformats.md + - Genomic coordinates: coordinates.md + - Dependencies: dependencies.md + - About: about.md diff --git a/splicekit/config/__init__.py b/splicekit/config/__init__.py index f1e3cf7..04200cd 100644 --- a/splicekit/config/__init__.py +++ b/splicekit/config/__init__.py @@ -76,6 +76,27 @@ def jbrowse2_config(): except: clip = None +try: + bam_path + bam_path_defined = True +except NameError: + bam_path_defined = False + +try: + bam_column + bam_column_defined = True +except NameError: + bam_column_defined = False + +if not bam_path_defined and not bam_column_defined: + print(f"{module_desc} ERROR: neither bam_path nor bam_column is set in splicekit.config") + print(f"{module_desc} Set at least one of the following in splicekit.config:") + print(f"{module_desc} bam_path = \"/path/to/bam/files\"") + print(f"{module_desc} BAMs are then expected at {{bam_path}}/{{sample_id}}.bam") + print(f"{module_desc} bam_column = \"bam_file\"") + print(f"{module_desc} samples.tab must then have a column with this name containing the full BAM path for each sample") + sys.exit(1) + try: bam_column except: