RBP Edit Denoiser (RED) is a tool for separating RBP-driven editing events from deaminase background in editing-based RBP profiling experiments (STAMP, INSCRIBE, and related assays). Given per-site MARINE profiled edit counts for matched treatment (RBP-deaminase fusion) and control (deaminase-only) replicates, it produces a per-site score that quantifies how likely each candidate site is to be a real RBP-binding-driven edit rather than baseline editing.

Scoring

train.py consumes matched treatment/control MARINE inputs and writes treatment_scored.tsv / control_scored.tsv, ranked by a combined_score column (higher = more likely a real RBP-driven edit). The score is ccs4 — a weighted sum of four control-calibrated z-axes, where each axis is the mean control-z over a family of differential (treatment-vs-control) features:

combined_score = W1·z_density + W2·z_edit_rate + W3·z_repro + W4·z_motif

• A1 z_density (clustering) — mean-z over multi-scale pooled_density_l2fc_{W} + local_enrichment_l2fc_{W}.
• A2 z_edit_rate (editing level) — mean-z over beta_edit_rate_l2fc + window_exceed_prob_{W}.
• A3 z_repro (reproducibility) — mean-z over beta_treat_reproducibility + n_treatment_reps + Beta confidence.
• A4 z_motif (sequence context) — control-background k-mer motif enrichment: enriched k-mers are called from log2(treat_freq / ctrl_freq) over the ±200 bp context, so the deaminase-only control is the background and genuinely treatment-specific sequence context is enriched.

Installation

Create a fresh conda/mamba environment named RED with Python 3.10.18 and the required packages:

mamba create -n RED python=3.10.18 -c conda-forge -c bioconda \
    numpy=2.2.6 pandas=2.3.3 scikit-learn=1.7.2 scipy=1.15.2 \
    matplotlib=3.10.9 seaborn=0.13.2 pysam=0.24.0 pyfaidx=0.9.0.4
mamba activate RED

External resources

scripts/annotate_sites.py (called by train.py) uses three external files for site filtering: ANNOVAR for gene-region annotation (drops intergenic sites), a dbSNP common-variants BED (drops known SNPs), and a RepeatMasker BED (drops Alu overlaps for A→I editing). The bundled scripts look for them under resources/ by default — drop each file into the path below or pass a custom path with the matching CLI flag.

Resource	Default path	CLI override	Source
ANNOVAR perl scripts	resources/table_annovar.pl (+ sibling *.pl)	—	annovar.openbioinformatics.org (academic registration required)
ANNOVAR hg38 refGene database	resources/humandb/	--annovar-db <dir>	downloaded with ANNOVAR (see below)
dbSNP common-variants BED3	resources/b151_GRCh38p7_common.bed3	--snp-bed <path>	NCBI dbSNP FTP — flatten to chrom\tstart\tend
RepeatMasker BED (Alu)	resources/alu_elements.bed	--alu-bed <path>	UCSC Table Browser → group Repeats → track RepeatMasker, output BED with columns chrom, start, end, family, subfamily, strand

After downloading ANNOVAR, there are also *.pl scripts, drop all *.pl files into resources/ (so table_annovar.pl can find its sibling helper scripts), we will need them to do annotation later. Now we can fetch the hg38 refGene database:

cd resources/
perl annotate_variation.pl -downdb -webfrom annovar refGene -buildver hg38 humandb/

This produces humandb/hg38_refGene.txt, humandb/hg38_refGeneMrna.fa, and humandb/hg38_refGeneVersion.txt (≈ 800 MB total).

Quick test (chr22 example)

example/ ships chr22-only MARINE outputs from the paper's INSCRIBE APOBEC1 TDP43 dataset (2 control + 2 treatment replicates, ≈ 13 MB total). It runs end-to-end in under a minute on 4 cores and is a good sanity check that your install + reference genome are set up correctly.

PYTHON=$(which python)            # from the activated `RED` env
REF=/path/to/hg38/genome.fa

$PYTHON scripts/train.py --edit-type 'C>T' \
    --control-files   example/control_rep1.tsv example/control_rep2.tsv \
    --treatment-files example/treatment_rep1.tsv example/treatment_rep2.tsv \
    --reference $REF \
    --output-dir example_out/ \
    --no-annotate \
    --min-coverage 7 \
    --density-windows 25,50,100,200 --n-jobs 4

This writes example_out/treatment_scored.tsv (and control_scored.tsv), ranked by combined_score. Add --axis-weights W1 W2 W3 W4 to reweight the four ccs4 axes (default 1 1 1 1).

--no-annotate skips the entire annotation step (ANNOVAR + intergenic removal + SNP/Alu filtering) so the example runs without any of the external resources above — only the reference FASTA is needed. Look for example_out/treatment_scored.tsv — the scored chr22 sites.

Sanity-check the scoring with a motif density plot

TDP43 binds UGUGUG on RNA, so ccs4-preferred edits should sit close to a TGTGTG motif on the genome, while control edits should not. scripts/motif_density.py plots this directly:

$PYTHON scripts/motif_density.py \
    --treatment example_out/treatment_scored.tsv \
    --control   example_out/control_scored.tsv \
    --reference $REF \
    --motif TGTGTG --top-n 3000 \
    --output example_out/motif_density.pdf

Takes the top 3,000 sites by combined_score from each table, computes the strand-aware distance from each edit to the nearest TGTGTG within ±400 bp, and saves a KDE comparing the two distributions. A working ccs4 score produces a sharp treatment peak around 0 bp and a flat / much weaker control curve.

Step 1 — Score

PYTHON=/tscc/nfs/home/s5xu/miniforge3/envs/ML/bin/python
REF=/tscc/projects/ps-yeolab3/ekofman/ReferenceData/hg38/cellranger-GRCh38-3.0.0/fasta/genome.fa
ANNOVAR=/tscc/nfs/home/s5xu/projects/resource/annovar/humandb 
SNP=/tscc/projects/ps-yeolab3/s5xu/deprecated/edit_classifier/resources/b151_GRCh38p7_common.bed3 
ALU=/tscc/projects/ps-yeolab3/s5xu/deprecated/edit_classifier/resources/alu_elements.bed

$PYTHON scripts/train.py --edit-type 'C>T' \
    --control-files   <ctrl_rep1.tsv> <ctrl_rep2.tsv> \
    --treatment-files <treat_rep1.tsv> <treat_rep2.tsv> \
    --manual-filtered-data <output_folder/marine/analysis/treat_rep1/filtered_site_info.tsv> <output_folder/marine/analysis/treat_rep2/filtered_site_info.tsv> \
    --reference $REF \
    --output-dir <out>/ \
    --min-coverage 7 \
    --annovar-db $ANNOVAR \
    --snp-bed $SNP \
    --alu-bed $ALU \
    --density-windows 25,50,100,200 --n-jobs 16

Outputs in <out>/:

• treatment_scored.tsv — the scored sites table, ranked by combined_score
• control_scored.tsv — the same scoring applied to the control sites
• annotated/ — gene/region annotations (ANNOVAR cache)

Step 2 — Cluster benchmark (optional, only with manual-filtered data)

If the user has manual-filtered (rule-based) call sets to compare against, cluster_edits_benchmark.py computes the five "yield setpoint" comparisons on top of treatment_scored.tsv:

$PYTHON scripts/cluster_edits_benchmark.py \
    --denoiser-output      <out>/treatment_scored.tsv \
    --manual-filtered-data <out>/annotated/<treat_rep1>_manual_filtered.tsv \
                           <out>/annotated/<treat_rep2>_manual_filtered.tsv \
    --control-data         <out>/annotated/<ctrl_rep1>_control_filtered.tsv \
                           <out>/annotated/<ctrl_rep2>_control_filtered.tsv \
    --sample-names         <treat_rep1>_annotated <treat_rep2>_annotated \
    --reference $REF --motif <MOTIF> --score-feature combined_score \
    --edit-type 'C>T' --annovar-db <annovar_humandb> \
    --annovar-reuse-dir <out>/annotated \
    --keep-pct 50 \
    --output-dir <out>/cluster_benchmark_keep50

The --keep-pct / --match-count / --threshold flag picks one of five modes (keep25, keep50, keep75, matchcount, threshold 0.7). Each call writes one mode's directory; the paper used five parallel calls.

Step 3 — Run FlARE Edit Peak Caller

flag	default
--flare-regions	/tscc/projects/ps-yeolab3/ekofman/ReferenceData/peakcalling_regions/hg38/hg38_cellranger
--flare-fdr	0.1
--flare-max-merge-dist	100
count_coverage_on_both_strands	true

For analysis on FLARE results, please refers to RED_Analysis repo at https://github.com/howardxu520/RED_Analysis