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