# SAGAconf User Manual

## Introduction

SAGAconf is a software that assigns calibrated confidence scores to chromatin state annotations, addressing the problem of reproducibility of SAGA annotations. It uses an integrative approach to derive a reproducibility score (r-value) at each genomic position, allowing for the identification of a confident and reliable subset from genome annotation while removing irreproducible predictions. Thus, SAGAconf allows a researcher to select only the reliable predictions from a chromatin annotation for use in downstream analyses. SAGAconf is independent of SAGA methodology, taking as input only posterior probability matrix.

## Input Data

SAGAconf reproducibility analysis pipeline gets as input two sets of posterior probability matrices (base and verification) each with **`G`** rows and **`K+3`** columns. Here, **`G`** denotes the genome size that is the number of rows and **`K`** corresponds to the number of chromatin states identified by SAGA method with 3 additional columns representing genomic coordinates of each bin (**`chr`**, **`start`**, **`end`**). These posterior files can be in BED or CSV format.

---

---

The **`SAGAconf_parser.py`** script is a command-line tool that allows you to parse posterior files generated by the SAGA model (ChromHMM and Segway). The script takes several arguments, some of which are required and others are optional.

## Required Arguments

- **`posteriordir`**: This argument specifies the directory containing all posterior files. You must provide the path to this directory when running the script. 
- **`resolution`**: This argument specifies the resolution of the SAGA model in base pairs (bp). You must provide this value when running the script.
- **`savedir`**: This argument specifies the directory where the parsed posterior file will be saved. You must provide the path to this directory when running the script.
- **`--saga`**: This argument specifies which SAGA model to use. You can choose between **`segway`** and **`chmm`**.

## Optional Arguments

- **`-h`** or **`--help`**: This argument displays a help message and exits the script. You can use this argument if you need more information about how to use the script.
- **`--out-format`**: This argument specifies the format for saving parsed posteriors. You can choose between **`bed`** and **`csv`** formats. If you do not provide this argument, the default format **`bed`**  will be used.

### Posterior directory
The **`posteriordir`** argument specifies the directory containing all posterior files. This directory can have two formats: one is the default format outputted by the ChromHMM, and the other one is the format outputted by Segway. For ChromHMM, the posterior probabilities over state assignments are printed in a file ending with **`_posterior.txt`** or **`_posterior.txt.gz`** if the **`-gzip`** flag is specified. One file is generated per cell type and chromosome. For Segway, this information is delivered in a series of numbered **`BED`** files, one for each segment label.

## Example Usage

Here is an example of how you might run the **`SAGAconf_parser.py`** script:

`python SAGAconf_parser.py --out_format bed --saga chmm /path/to/posteriordir 200 /path/to/savedir`

In this example, we are specifying that we want to use the **`bed`** format for saving parsed posteriors, and we want to use the **`chmm`** SAGA model. We are also providing the required arguments for the **`posteriordir`**, **`resolution`**, and **`savedir`**.

---

---

Given the parserd_posterior files, SAGAconf.py can provide following analysis:

1. ratio of naive overlap between chromatin states of base and verification annotations
2. Granularity of chromatin states vs. overlap
3. Spatial misalignment of chromatin states vs. overlap
4. Posterior probability vs. overlap
5. Post-clustering of chromatin states by merging into <K chromatin states
6. r-Values per genomic position
7. Robust chromatin states (confident segments) 

---

---

The **`SAGAconf.py`** script is a command-line tool that allows you to compare base and verification annotations generated by the SAGA model. The script takes several arguments, some of which are required and others are optional.

## Required Arguments

- **`base`**: This argument specifies the path of the parsed posterior CSV or BED file of the base replicate annotation. You must provide the path to this file when running the script.
- **`verif`**: This argument specifies the path of the parsed posterior CSV or BED file of the verification replicate annotation. You must provide the path to this file when running the script.
- **`savedir`**: This argument specifies the directory where the SAGAconf results will be saved. You must provide the path to this directory when running the script.

## Optional Arguments

- **`-h`** or **`--help`**: This argument displays a help message and exits the script. You can use this argument if you need more information about how to use the script.
- **`-v`** or **`--verbosity`**: This argument increases output verbosity. You can use this argument if you want more detailed output from the script.
- **`-bm`** or **`--base_mnemonics`**: This argument specifies a file to use as mnemonics (biological label interpretations) for the base replicate. If you provide this argument, you must also provide the path to the text file.
- **`-vm`** or **`--verif_mnemonics`**: This argument specifies a file to use as mnemonics (biological label interpretations) for the verification replicate. If you provide this argument, you must also provide the path to the text file.
- **`-s`** or **`--subset`**: This argument specifies that SAGA should be run on just one chromosome (chr21). If you provide this argument, only one chromosome (chr21) will be analyzed. Note that, if your annotations do not include chr21, you should not use this argument.
- **`-w`** or **`--windowsize`**: This argument specifies the window size in base pairs (bp) to account for around each genomic bin. The default value is 1000 bp.
- **`-to`** or **`--iou_threshold`**: This argument specifies the threshold on the Intersection over Union (IoU) of overlap for considering a pair of labels as corresponding. You can provide a value for this threshold when running the script. The default value is 0.75.
- **`-tr`** or **`--repr_threshold`**: This argument specifies the threshold on the reproducibility score for considering a segment as reproduced. You can provide a value for this threshold when running the script. In the original paper, this argument is referred to as **`alpha`.** The default value is 0.8.
- **`-k`** or **`--merge_clusters`**: This argument specifies a k value to merge base annotation chromatin states until k states. You can provide a value for k when running the script.
- **`-q`** or **`--quick`**: This argument specifies that only a subset of essential analysis should be performed for a quick report. If you provide this argument, only essential analysis will be performed.

## Example Usage

Here is an example of how you might run the **`SAGAconf.py`** script:

`python SAGAconf/SAGAconf.py -v -bm /path/to/base_mnemonics.txt -vm /path/to/verif_mnemonics.txt -s -w 2000 -to 0.5 -tr 0.8 -k 5 -q /path/to/base/parsed_posterior.bed /path/to/verif/parsed_posterior.bed /path/to/savedir`

In this example, we are specifying that we want increased output verbosity (**`-v`**), and we are providing paths to text files containing mnemonics for both base and verification replicates (**`-bm`** and **`-vm`**). We are also specifying that we want to run SAGA on just chr21 chromosome (**`-s`**), and we are providing values for window size (**`-w`**), IoU threshold (**`-to`**), reproducibility threshold (**`-tr`**), and k value (**`-k`**). Finally, we are specifying that only essential analysis should be performed (**`-q`**) and providing paths to base and verification replicate files and save directory.

## Mnemonics file

A mnemonics file is a tab-separated text file with two columns. The first column contains the integer label of the original chromatin states, and the second column contains the biological label assignments. Here is an example of a mnemonics file:


```
old	new
1	Enhancer_low
2	Enhancer
3	Promoter_flanking
4	Promoter
5	Promoter
6	Promoter
7	Enhancer
8	Enhancer
9	Enhancer
10	Enhancer
```


To generate your own mnemonics file, you can create a new text file and enter the integer labels and biological label assignments in the appropriate columns. Make sure to separate the columns with a tab character and to include a header row with the column names **`old`** and **`new`**.

---

---

## Real Use Case Example

In this section, we will go through a real use case of SAGAconf starting from obtaining chromatin state annotations from ChromHMM, followed by parsing the posteriors and running SAGAconf. Note that, in the following example, everything is happening within the main `SAGAconf/` directory.

1. First, download ChromHMM and unzip it by running the following commands:
    
    ```
    curl -o ChromHMM.zip http://compbio.mit.edu/ChromHMM/ChromHMM.zip
    unzip ChromHMM.zip
    ```
    
2. Next, generate test annotations from ChromHMM by running the following command:
    
    `java -mx1600M -jar ChromHMM/ChromHMM.jar LearnModel -printposterior ChromHMM/SAMPLEDATA_HG18 ChromHMM/OUTPUTSAMPLE 10 hg18`
    
    This command runs the `LearnModel` command of the `ChromHMM` program with the `-printposterior` option. This command takes a set of binarized data files, learns chromatin state models, and produces chromatin state annotations with posterior files where `k=10`. In the ChromHMM sample run, it concatenates the data from `K562_chr11` and `GM12878_chr11`. For demonstration purposes, we will treat these two as base and verification replicates.
    
3. Next, we manually move the posterior files for `GM12878_chr11` and `K562_chr11` into separate directories for base and verification replicates using the following command:

    ```
    mkdir -p ChromHMM/OUTPUTSAMPLE/base ChromHMM/OUTPUTSAMPLE/verif && \
    mv ChromHMM/OUTPUTSAMPLE/POSTERIOR/GM12878_10_chr11_posterior.txt ChromHMM/OUTPUTSAMPLE/base && \
    mv ChromHMM/OUTPUTSAMPLE/POSTERIOR/K562_10_chr11_posterior.txt ChromHMM/OUTPUTSAMPLE/verif && \
    rm -r ChromHMM/OUTPUTSAMPLE/POSTERIOR
    ```
    
    This command creates two directories, `ChromHMM/OUTPUTSAMPLE/base` and `ChromHMM/OUTPUTSAMPLE/verif`, moves the posterior files for `GM12878_chr11` and `K562_chr11` into these directories, respectively, and removes the original `POSTERIOR` directory.
    
4. Then, we use the following commands to parse posteriors from `ChromHMM/OUTPUTSAMPLE/base` and `ChromHMM/OUTPUTSAMPLE/verif` into the standard format required by SAGAconf:
    
    ```
    python SAGAconf_parser.py --saga chmm ChromHMM/OUTPUTSAMPLE/base 200 ChromHMM/OUTPUTSAMPLE/base

    python SAGAconf_parser.py --saga chmm ChromHMM/OUTPUTSAMPLE/verif 200 ChromHMM/OUTPUTSAMPLE/verif
    ```
    
    Here, using `--saga chmm`, we are specifying that the annotations are generated by ChromHMM and we are specifying resolution to be 200bp. Ultimately, we will have `parsed_posterior.bed` files in both `ChromHMM/OUTPUTSAMPLE/base` and `ChromHMM/OUTPUTSAMPLE/verif`.
    
5. Lastly, we run SAGAconf to obtain a full reproducibility report with default parameters by running the following command:
    
    `python SAGAconf.py ChromHMM/OUTPUTSAMPLE/base/parsed_posterior.bed ChromHMM/OUTPUTSAMPLE/verif/parsed_posterior.bed ChromHMM/OUTPUTSAMPLE/sagaconf_base`


You can also, execute `example/run.sh` to run the whole real use case example. Make the `example/run.sh` file executable by running the following command:

`chmod +x example/run.sh`

`./example/run.sh`

---

---

## Outputs:

### Output files:

1. `AUC_mAUC.txt` : Area under the state_merging curve (auSMC) per chromatin state
2. `barplot` : Area under the state_merging curve (auSMC) per chromatin state
3. `binned_posterior_heatmap` : IoU overlap with binned posterior of base annotation
4. `conf_progress_entropy` : Ratio_robust vs annotation entropy during state merging
5. `conf_progress` :  Ratio_robust vs number of chromatin states during state merging
6. `confident_segments_dense.bed` : Robust subset (dense list) of the annotation according to SAGAconf
7. `confident_segments.bed` : Robust subset of the annotation according to SAGAconf
8. `coverages1.txt` : Genome coverage per chromatin state (base annotation)
9. `coverages2.txt` : Genome coverage per chromatin state (verification annotation)
10. `granularity` : state_merging curve per chromatin state
11. `heatmap_w` : IoU overlap with maximum-a-posterior state with w > 0
12. `heatmap` : IoU overlap with maximum-a-posterior state with w = 0
13. `len_bound_overall` : genome-wide overlap as a function of w 
14. `len_bound` : overlap as a function of w per chromatin state
15. `MI_post_clustering_progress.txt` : Change in Mutual information during state merging
16. `naive_post_clustering_progress.txt` : Change in naive overlap during state merging
17. `NMI.txt` : Mutual information with and without taking posterior into account
18. `overlap_ratio.txt` : Per chromatin state naive overlap
19. `post_clustering_progress.txt` : Change in ratio_robust during state merging
20. `r_values_report.txt` : average r-value genome-wide and per chromatin state
21. `r_values.bed` : r-values per genomic position
22. `ratio_robust.txt` : ratio_robust genome-wide and per chromatin state

### Output sub-directories:
1. `calib/` : posterior calibration curves
2. `Dist_vs_Corresp/` : overlap as a function of w
3. `Dist_vs_Corresp3/` : correspondence as a function of w
