SAGAconf

SAGAconf gives a calibrated confidence score to each call in a chromatin state annotation. Segmentation and genome annotation (SAGA) methods, such as ChromHMM and Segway, label every genomic bin with a chromatin state. Many of these labels do not reproduce across replicates. SAGAconf compares a base annotation with a verification annotation and gives each bin an r-value: a score for how well its state call reproduces. You can then keep only the reproducible calls for downstream analysis. SAGAconf works with any SAGA method, because it reads only the posterior probability matrix.

• Paper: Foroozandeh Shahraki, Farahbod and Libbrecht, Robust chromatin state annotation, Genome Research 34(3):469–483, 2024.
• Blog post: Reproducibility as a measure of confidence.

Install

pip install sagaconf

SAGAconf needs Python 3.9 or later. The active-regions mode also needs bedtools on your PATH.

Quick start

Parse the posteriors of each annotation, then compare the two parsed files:

sagaconf parse --saga chmm base_posteriors/ 200 base/
sagaconf parse --saga chmm verif_posteriors/ 200 verif/
sagaconf run base/parsed_posterior.bed verif/parsed_posterior.bed results/

results/r_values.bed holds one r-value per bin. results/confident_segments.bed holds the reproducible subset of the base annotation (this file needs mnemonics; see Known issues).

To run the full example on ChromHMM's sample data, run example/run.sh.

Input

SAGAconf compares two annotations of the same genome:

• The base annotation is the one you want to score.
• The verification annotation comes from a replicate. It can come from different data and a different model (setting S1 in the paper), the same model trained on both replicates (S2), or the same data with a different random initialization (S3).

sagaconf run reads one parsed posterior file per annotation. The file has one row per genomic bin and K + 3 columns: chr, start, end, and one posterior probability for each of the K states. The file can be BED (tab-separated) or CSV. sagaconf parse writes this file from the output of ChromHMM or Segway:

• ChromHMM: run LearnModel with -printposterior. Give sagaconf parse the directory of *_posterior.txt files (one file per cell type and chromosome).
• Segway: give sagaconf parse the directory of posterior<i>.bedGraph files (one file per state).

Commands

sagaconf parse

sagaconf parse --saga {chmm,segway} [--out-format {bed,csv}] POSTERIORDIR RESOLUTION OUTDIR

sagaconf run

sagaconf run BASE VERIF OUTDIR [--mode MODE] [options]

python -m sagaconf runs the same command as sagaconf.

Mnemonics file

A mnemonics file gives each state a biological name. It is a tab-separated file with the header old and new. The old column holds the state number, in the column order of the parsed posterior file. The numbers can start at 0 or at 1. The new column holds the name. SAGAconf shortens each name to its first 4 characters (and to 3 characters after an underscore), so Enhancer_low becomes Enha_low. example/example_mnemonics.txt is an example:

old	new
1	Enhancer_low
2	Enhancer
3	Promoter_flanking

Output

sagaconf run writes these files in OUTDIR in full mode:

Each plot is written as PDF and SVG. Each .txt file next to a plot holds the plotted values.

Known issues

These issues come from the version of SAGAconf that the paper used. This release keeps them, so that its results match that version exactly. Fixes will come in later releases.

• In full mode, -k does not merge states.
• In full mode without mnemonics, SAGAconf does not write the genome-wide results (including confident_segments.bed and the UCSC track). Give -bm and -vm to get them.
• In merge mode, SAGAconf can stop with a KeyError after it writes the merge reports.

Legacy interface

The original scripts still work, with the same flags and the same output:

python SAGAconf_parser.py --saga chmm POSTERIORDIR 200 OUTDIR
python SAGAconf.py [-q | --r_only | --ct_only | ...] BASE VERIF OUTDIR

sagaconf run also accepts the original flag names (--r_only, --windowsize, --base_mnemonics, and so on). In the original scripts, the --active_regions mode reads src/biointerpret/GRCh38-cCREs.bed and src/biointerpret/Meuleman.tsv relative to the working directory.

The code as the paper used it is on the legacy branch (tag v0-legacy). That branch also keeps the scripts that produced the paper's figures.

Tests

tests/golden/ checks that this version gives byte-identical output to the legacy version. It runs both versions on the same input and compares every file they write, and everything they print. tests/golden/gate.py describes the method. The check passes on synthetic data and on the paper's GM12878 ChromHMM and Segway annotations.

Citation

@article{foroozandeh2024robust,
  title   = {Robust chromatin state annotation},
  author  = {Foroozandeh Shahraki, Mehdi and Farahbod, Marjan and Libbrecht, Maxwell W.},
  journal = {Genome Research},
  volume  = {34},
  number  = {3},
  pages   = {469--483},
  year    = {2024},
  doi     = {10.1101/gr.278343.123}
}

License

SAGAconf is available under the MIT License.