scraps

Snakemake Check Github Clones

scraps extracts mRNA polyadenylation sites from “TVN”-primed single-cell RNA-seq libraries at near-nucleotide resolution.

scraps (Single Cell RNA PolyA Site Discovery) is currently implemented as a Snakemake pipeline for 10X Genomics 3’ end v2/v3/v4 libraries (and other platforms with similar library structure, including Drop-seq, Microwell-seq, and BD Rhapsody), sequenced on Illumina, Element or Ultima instruments, as well as pre-aligned, barcode-tagged BAMs (e.g. ONT / PacBio long reads). If long Read1 is available (estimated ~6% of SRA-deposited data, or now planning new experiments), positional information will be calculated from paired realignment; otherwise, the less optimal anchored Read2 approach is used. scraps will eventually be expanded for analyzing a range of RNA processing changes in single-cell RNA-seq data.

For additional discussions and usage cases, please see bioRxiv preprint.

Example usage

scraps requires the following as input (defined in config.yaml):

Quick Start

  1. Set up conda environment: 
    conda env create -f scraps_conda.yml
conda activate scraps_conda

  2. Configure your samples in config.yaml under the SAMPLES section

  3. Run the pipeline: 
    snakemake -j 24 --resources total_impact=5 --keep-going

Detailed Usage

To run the bundled test samples (after downloading the STAR index and whitelists from Setup; the dry-run-only fixtures chromiumv3_ultima_fastq_test and the bam_*_test samples in config.yaml should be removed from SAMPLES first, since they have no FASTQs / only a synthetic BAM):

snakemake --snakefile Snakefile \
  -j 24 \
  --resources total_impact=5 \
  --keep-going

config.yaml is read automatically (configfile: in the Snakefile). To pass a different one, use --configfile other.yaml; put any explicit target before that flag, because --configfile accepts multiple files.

DAG steps illustration (rule graph of all four alignment modes with discovery enabled)

Cluster execution (SLURM): either run the whole pipeline inside one batch allocation with the command above, or submit each rule as its own SLURM job with the bundled profile (set your account and partition in profiles/slurm/config.v8+.yaml first):

snakemake --profile profiles/slurm --configfile config.yaml

See profiles/slurm/README.md and the SLURM executor plugin docs. Other schedulers use their own executor plugins (--executor <name>); Snakemake 8+ no longer has --cluster/--drmaa.

Notes: each STAR alignment job (one per sample × alignment mode) uses total_impact=5; --resources total_impact=5*N allows N concurrent STAR jobs.

STAR memory: each STARsolo job requests the genome index size (read from STAR_INDEX) plus 8 GB, at least 48 GB. STAR frees the index after mapping, before STARsolo counting and BAM sorting, and --limitBAMsortRAM is set to half the request, the other half kept for STARsolo’s per-read data, so sorting stays inside the allocation. For large libraries set STAR_mem_mb (MB) in DEFAULTS, a platform block, or a sample block, or override mem_mb with set-resources/--set-resources; the sort buffer follows either way. The SLURM profile retries a failed job once, and STAR jobs get twice the memory on the retry.

Configuration

scraps uses two main configuration files for flexible pipeline setup:

config.yaml

Main pipeline configuration file containing:

Sample Configuration Example:

SAMPLES:
  sample_name:
    basename: sample-           # FASTQ file prefix
    platform: illumina          # Sequencing platform
    chemistry: chromiumV3       # Platform chemistry
    alignments:                 # Optional: override default alignments
      - R2
      - paired

chemistry.yaml

Platform and chemistry-specific parameters organized hierarchically. Each chemistry type (chromiumV3, chromiumV4, chromiumV2, ddseq, dropseq, microwellseq, bd, indrop) contains:

Configuration Hierarchy

The pipeline uses hierarchical configuration lookup to determine parameters for each sample:

┌─────────────────────────────────────────────────────────┐
│  1. Sample-specific settings (config.yaml SAMPLES)     │
│     Highest priority - overrides everything             │
└────────────────────┬────────────────────────────────────┘
                     │ If not found ↓
┌─────────────────────────────────────────────────────────┐
│  2. Chemistry + Platform (chemistry.yaml)               │
│     e.g., chromiumV3 → illumina → STAR_R1               │
└────────────────────┬────────────────────────────────────┘
                     │ If not found ↓
┌─────────────────────────────────────────────────────────┐
│  3. Chemistry defaults (chemistry.yaml)                 │
│     e.g., chromiumV3 → bc_whitelist                     │
└────────────────────┬────────────────────────────────────┘
                     │ If not found ↓
┌─────────────────────────────────────────────────────────┐
│  4. Global defaults (config.yaml DEFAULTS)              │
│     Lowest priority - fallback values                   │
└─────────────────────────────────────────────────────────┘

This allows platform-specific customization (e.g., Illumina vs Ultima Genomics) while maintaining chemistry-specific defaults.

Supported scRNA-seq platforms

| Platform | Library (BC+UMI+A) | Setting | Test data | | :——–|:————| :————| :———| | 10x Chromium V3 | 16 + 12 + 30 | chromiumV3 | (via the Ultima and synthetic BAM fixtures) | | 10x Chromium V4 (GEM-X) | 16 + 12 + 30 | chromiumV4 | | | 10x V3/V4 - Ultima Genomics, CRAM-derived reads | adapter + 16 + 9 + 3 ignored + 8 | chromiumV3 / chromiumV4, platform: ultima_cram (R1 mode) | ✓ | | 10x V3/V4 - Ultima Genomics, delivered FASTQs | R1 = 16 + 9, R2 = cDNA (pseudo-paired); priming site = R2 3’ end + 5 bp | chromiumV3 / V4, platform: ultima_fastq (R2 mode) | | | Pre-aligned BAM (ONT / PacBio long reads, or any compatible single- or paired-end BAM) | CB/UB-tagged, full-length molecules; see Pre-aligned BAM input | any 10x chemistry, platform: generic, alignments: [bam] | ✓ (synthetic) | | 10x Chromium V2 | 16 + 10 + 30 | chromiumV2 | ✓ | | ddSEQ (Bio-Rad / Illumina SureCell) | 21 + 8 (after bc_cut linker removal) | ddseq | | | Drop-seq | 12 + 8 + 30 | dropseq | ✓ | | Microwell-seq | 6x3 + 6 + 30 | microwellseq | ✓ | | BD Rhapsody | 9x3 + 8 + 18 | bd | | | inDrop | 8 + 6 + 18 | indrop | |

Custom chemistry supported, by editing chemistry.yaml. Also see synthetic FASTQ tool.

Ultima Genomics notes. Reads converted from Ultima CRAM files are one long read (adapter + CB + UMI + cDNA) and use platform: ultima_cram with alignments: [R1] (no cutadapt step). FASTQs delivered by Ultima for 10x libraries are instead pseudo-paired (R1 = cell barcode + UMI, R2 = cDNA) and use platform: ultima_fastq with alignments: [R2]. Because Ultima trims the 5 bp adjacent to the poly(A) homopolymer, ultima_fastq disables the soft-clipped poly(A) filter and shifts the priming site 5 bp downstream of the R2 alignment end (R2_polyA_filter: false, R2_priming_offset: 5). No detectable poly(A) is required for these reads at any step. The FASTQ file names must follow the usual *_R1* / *_R2* pattern.

Pre-aligned BAM input

The bam alignment mode counts a BAM that was aligned and barcode-tagged elsewhere (ONT or PacBio long reads, or any compatible short-read BAM) without a STAR step. Use platform: generic, alignments: [bam], and a per-sample bam: path (absolute, or relative to DATA); basename is not needed.

SAMPLES:
  my_longread:
    bam: my_longread.tagged.bam
    chemistry: chromiumV3
    platform: generic
    alignments: [bam]
    # bam_strand: same        # see table
    # bam_polyA_filter: false
    # bam_mate: 0

Which end of the read is the priming site, and whether a poly(A)/poly(T) soft clip is required there, is declared per sample (defaults in config.yaml DEFAULTS):

class bam_strand counted end bam_polyA_filter featureCounts
1 same (read is on the transcript strand) read 3’ end false -s 1 --read2pos 3
2 same read 3’ end true (3’ poly(A) soft clip required, like R2) -s 1 --read2pos 3
3 reverse (read is antisense) read 5’ end false -s 2 --read2pos 5
4 reverse read 5’ end true (5’ poly(T) soft clip required) -s 2 --read2pos 5

Paired-end BAMs: set bam_mate: 1 or 2 to name the mate whose end is the priming site. Only proper pairs are kept, the other mate is dropped, and the selected mate is rewritten as a single-end record, so every fragment contributes one position. bam_strand and bam_polyA_filter then describe that mate (a 10x-style paired BAM is class 3 on mate 1). bam_priming_offset (bp the priming site lies beyond the counted end) works as for R2_priming_offset.

Assumptions (documented, not enforced): the BAM contains full-length molecules (adapters at both ends), carries CellRanger/STARsolo-style tags (only CB and UB are read), and its strand assignment is accurate for the declared bam_strand. Reads without CB/UB, and secondary/supplementary alignments, are always dropped. featureCounts runs in long-read mode (-L, single-threaded) so spliced long reads keep their full CIGAR. The soft-clip filter only recognises M/S CIGAR operations; BAMs with =/X or hard clips (e.g. pbmm2, minimap2 --eqx) must use bam_polyA_filter: false.

Output

  1. Aligned BAMs: {RESULTS}/{sample}/{sample}_{mode}_Aligned.sortedByCoord.out.bam (STARsolo output with CB/UB tags; not produced for the bam mode)
  2. Count table: {RESULTS}/counts/{sample}_{mode}_counts.tsv.gz, directional-UMI-collapsed counts per cell barcode for each POLYA_SITES feature (bundled hg38 SAF: a transcript-oriented -10/+5 bp window around each polyAdb cleavage site). GeneID fields are ;-delimited for hg38 and _-delimited for mm10; overlapping features are comma-joined. 
    gene    cell    count
RPL26;6154;ENSG00000161970;chr17;8377421;-;3'UTR(L)     AACTCCCGTTCCTCCA        1
RPL26;6154;ENSG00000161970;chr17;8377421;-;3'UTR(L)     CCCATACGTTAAAGAC        1
ACTG1;71;ENSG00000184009;chr17;81509999;-;3'UTR(M)      ACATCAGGTGATGTCT        1
ADRM1;11047;ENSG00000130706;chr20;62308862;+;3'UTR(M)   CAGCGACTCTGCCCTA        1

  3. (Discovery only) Strand-aware single-base priming pileups {RESULTS}/discovery/beds/{sample}_{mode}.stranded.bed.gz (chrom, start, end, UMI count, transcript strand) and the annotated site tables described under De novo priming site discovery.

  4. Benchmarks: {RESULTS}/benchmarks/<rule>_<sample>_<mode>.tsv (discovery: <rule>_<unit>.tsv), Snakemake’s per-job wall time, peak memory and CPU time for every rule. Summarize them per sample with python3 inst/scripts/benchmark_summary.py {RESULTS}/benchmarks.

A MultiQC report is not produced by default (the qc_check rule in rules/qc.snake is disabled in rule all); an example report from an earlier version is kept under inst/test_output/.

R functions available for importing results into Seurat object, and finding differential PA site usage; source() inst/scripts/R/geneid.R first (it provides the delimiter-aware GeneID parser the other files use; they also source it automatically when loaded with source()). These copies are authoritative: the separate scrapR package (remotes::install_github("rnabioco/scrapR"), also in the Docker image) still carries the older _-only GeneID parsing and does not work with the ;-delimited hg38 output.

De novo priming site discovery

By default scraps quantifies poly(A) usage against validated poly(A) sites (polyAdb), which is the recommended reference for downstream quantitative analyses. Optionally, scraps can also use its single-base resolution to discover RT priming sites de novo and annotate them.

The discovery module (off by default) applies kernel density estimation to pseudobulked, UMI-deduplicated, strand-aware single-base priming positions (per sample, or per a defined group of samples), calls priming-site summits, and classifies each as:

Note: discovered sites are for hypothesis generation. For quantitative analyses, use validated poly(A) sites. Both novel classes require orthogonal validation. suspected_internal_priming sites are ambiguous rather than artifacts: interpret them with caution, report results with and without them where it matters, and confirm with an orthogonal method before treating them as real sites.

Enable it by setting GENOME_FASTA and DISCOVERY.enabled: true in config.yaml. Outputs are written to {RESULTS}/discovery/ (annotated TSV, BED, and a reusable de novo SAF per discovery unit, plus the per-(sample, mode) stranded priming beds under discovery/beds/). See docs/discovery.md for the full method, configuration, caveats, and re-quantitation workflow.

Setup

1. Clone repository

git clone https://github.com/rnabioco/scraps
cd scraps

conda env create -f scraps_conda.yml
conda activate scraps_conda

Alternatively, ensure all dependencies are installed and available in your PATH.

3. Prepare reference files

STAR genome index

Place STAR index in the ref/ directory or specify custom path in config.yaml (STAR_INDEX must point at the extracted index directory; the default is ref/cr2020A_star/). Not needed for BAM input.

Download link (extract after download):

Barcode whitelists (required for the 10x chemistries as configured)

Whitelist paths are configured per chemistry in chemistry.yaml and passed to STARsolo --soloCBwhitelist, so a missing file fails the alignment; set bc_whitelist: None to run without one. Place downloaded whitelists in the ref/ directory.

Download links (extract after download):

Update chemistry.yaml with the correct paths:

chromiumV3:
  bc_whitelist: ref/3M-february-2018.txt
chromiumV4:
  bc_whitelist: ref/3M-3pgex-may-2023_TRU.txt   # 10x GEM-X 3' whitelist (shipped with Cell Ranger)
chromiumV2:
  bc_whitelist: ref/737K-august-2016.txt

4. Configure your samples

Edit config.yaml to specify:

Example:

SAMPLES:
  my_sample:
    basename: SRR9887775_        # Matches SRR9887775_R1.fastq.gz, SRR9887775_R2.fastq.gz
    chemistry: chromiumV3
    platform: illumina

Note: scraps does not download from SRA; FASTQs must already be present under DATA (e.g. after fasterq-dump --split-files, renamed to the *R1* / *R2* pattern).

5. Run the pipeline

# Dry-run to check configuration (also validates chemistry/platform/mode and bam-mode keys)
snakemake -np

# Run pipeline (--cores/-j is required by Snakemake for local execution)
snakemake -j 24 --resources total_impact=5 --keep-going

# Dry-run a single target (target BEFORE any --configfile flag)
snakemake -np results/counts/chromiumv2_test_R2_counts.tsv.gz

Sample test results from an earlier version can be found at inst/test_output/ (legacy _-delimited GeneIDs and bed tracks).

Dependencies

scraps requires the following executables in your PATH and Python packages in the running environment (all in scraps_conda.yml):

Recommended: Use Conda to manage these dependencies:

conda env create -f scraps_conda.yml
conda activate scraps_conda

All required dependencies (including zsh) will be installed automatically.

Docker image for automated deployment can also be found at https://hub.docker.com/r/rnabioco/scraps.

Please also see the Snakemake documentation for general information on executing and manipulating snakemake pipelines.

For Developers

For detailed development guidelines including code style conventions, testing procedures, and instructions for adding new rules or chemistry configurations, see AGENTS.md.

Key resources:

Extended function

1) Measuring internal priming as indicator of apoptotic cytoplasmic poly(A) RNA decay

(Based on widespread RNA decay during apoptosis: Liu and Fu et al.) Build a region SAF marking all gene regions (5’UTR, intron, CDS, 3’UTR) from the polyAdb SAF and a GTF with expand_region_saf.R, set it as POLYA_SITES, and process the output with the helper R functions. Please see Rmarkdown notebook for more.

2) Accurate intron/exon quantification for RNA velocity

(See discussions on quantification approaches and pitfalls: Soneson et al.)

Consideration scraps
Avoid feature double-counting ✓
Take strandedness into account ✓
Avoid count subtraction ✓
Resolve spliced vs unspliced target ✓
Speed ✓
This site is open source. Improve this page.