### Setup Conda Environment for UniversalEPI

Source: https://context7.com/boevalab/universalepi/llms.txt

Creates and activates a conda environment with all necessary dependencies for UniversalEPI, including PyTorch and PyTorch Lightning. Requires a `environment.yml` file and data/model checkpoints to be downloaded.

```bash
conda env create -f environment.yml

conda activate universalepi

# Download required data directory from Zenodo and place in root
# https://zenodo.org/records/15079298
# Expected structure: ./data/

# Download model checkpoints from Zenodo
# https://doi.org/10.5281/zenodo.14622040
# Place in ./checkpoints/ directory
```

--------------------------------

### Create Conda Environment for UniversalEPI

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Installs all necessary packages for UniversalEPI by creating a conda environment from the provided YAML file. This ensures a reproducible and isolated environment for the project.

```bash
conda env create -f environment.yml
```

--------------------------------

### Define Genomic Region of Interest (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

Specifies the target genomic region for analysis. This is a string representing a chromosome and its start and end coordinates.

```python
# Define the region of interest
region = "chr12:43600000-44300000"
```

--------------------------------

### Train Stage 1 CNN Model

Source: https://context7.com/boevalab/universalepi/llms.txt

Trains the Stage 1 Convolutional Neural Network (CNN) model for genomic feature encoding. Configuration is managed via Hydra, allowing for flexible setup of training cell lines and model parameters. Key parameters include kernel sizes, channel dimensions, and pooling sizes.

```bash
# Configure training cell lines in:
# ./Stage1/configs/datamodule/validation/cross-cell.yaml
# Set: train: [gm12878, k562]

# Run training
python ./Stage1/train.py

# The model is configured via ./Stage1/configs/demo.yaml
# Key parameters in Stage2/configs/configs.yaml:
#   ksizes: [11, 11, 11, 5, 5]        # CNN kernel sizes
#   channels: [5, 30, 60, 60, 90, 90, 90, 4]  # Channel dimensions
#   poolings: [4, 5, 5, 4, 2]         # Pooling sizes
```

--------------------------------

### Process GENCODE GTF for Protein-Coding Genes (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

This Python function reads a GENCODE GTF file, filters for 'gene' features, extracts 'gene_name' and 'gene_type' from attributes, keeps only 'protein_coding' genes, selects relevant columns, calculates the transcription start site (TSS), removes duplicate gene names, drops rows with missing values, and saves the result to a CSV file. It depends on the pandas and os libraries.

```python
import pandas as pd
import os

def process_gencode():
    gene_df = pd.read_csv("../data/genome/gencode.v47.annotation.gtf.gz", sep="\t", comment="#", 
                            names=["chr", "source", "feature", "start", "end", "score", "strand", "frame", "attributes"])
    gene_df = gene_df[gene_df.feature == "gene"]
    gene_df["gene_name"] = gene_df.attributes.str.extract(r'gene_name \"(.+?)\";')
    gene_df['gene_type'] = gene_df.attributes.str.extract(r'gene_type \"(.+?)\";')
    gene_df = gene_df[gene_df.gene_type == "protein_coding"]
    gene_df = gene_df[['chr', 'start', 'end', 'gene_name', 'strand']]
    gene_df["tss"] = gene_df.apply(lambda x: x["start"] if x["strand"] == "+" else x["end"], axis=1)
    gene_df = gene_df.drop_duplicates(subset="gene_name")
    gene_df = gene_df.dropna()
    gene_df.to_csv("../data/genome/gencode_processed.csv", index=False)
    return gene_df

if os.path.isfile("../data/genome/gencode_processed.csv"):
    gene_df = pd.read_csv("../data/genome/gencode_processed.csv")
else:
    gene_df = process_gencode()
gene_df
```

--------------------------------

### Download and Unzip Macrophage Predictions (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

Downloads prediction data for macrophage differentiation from Zenodo and extracts the zip files. It creates a 'results' directory if it doesn't exist and cleans up the zip files after extraction.

```python
# Download Macrophage predictions from Zenodo
def download_predictions(cell_line):
    download_path = f"https://zenodo.org/records/15079440/files/{cell_line}.zip?download=1"
    os.makedirs('../results', exist_ok=True)
    os.system(f"wget {download_path} -O ../results/{cell_line}.zip")
    os.system(f"unzip ../results/{cell_line}.zip -d ../results/")
    os.system(f"rm ../results/{cell_line}.zip")
download_predictions("Reed_0000_merged")
download_predictions("Reed_1440_merged")
```

--------------------------------

### Stage 2 Configuration File Reference (YAML)

Source: https://context7.com/boevalab/universalepi/llms.txt

This is a reference for the main Stage 2 configuration file, detailing runtime settings, data paths, model architecture parameters for both Stage 1 CNN and Stage 2 transformer, feature flags, uncertainty estimation settings, and positional encoding configurations.

```yaml
# ./Stage2/configs/configs.yaml

# Runtime settings
print_freq: 10000
seed: 1337
wandb_log: False
run_name: paper-hg38-map-concat-stage1024-rf-lrelu-eval-stg-newsplit-newdata-atac-var-beta-neg-s1337

# Data paths
hic_res: 5000                              # Hi-C resolution in bp
hic_path: data/hic/
atac_path: data/stage1_outputs/
cell_lines_train: ["gm12878", "k562"]
cell_lines_val: ["gm12878", "k562"]
cell_lines_test: hepg2
data_save_dir: null                        # Auto-generate NPZ files
train_dir: null                            # Or path to pre-made train NPZ
val_dir: null
test_dir: null
logs_dir: checkpoints
stage_1_model: checkpoints/stage1_model_gm12878_k562.pth

# Stage 1 CNN architecture
ksizes: [11, 11, 11, 5, 5]
channels: [5, 30, 60, 60, 90, 90, 90, 4]
poolings: [4, 5, 5, 4, 2]

# Stage 2 transformer architecture
n_epochs: 20
batch_size: 32
n_feat: 180
seq_len: 401                               # Number of ATAC peaks in window
binning: 1000                              # Peak region length in bp
dropout: 0.1
num_heads: 4
num_layers: 4
embed_dim: 32
hidden_dim: 96
lr: 0.001
log_scale: True

# Feature flags
MAP: True                                  # Include mappability track
ATAC: True                                 # Include ATAC signal

# Uncertainty estimation
var_flg: True                              # Enable aleatoric uncertainty
beta: 0.5                                  # Beta-NLL loss hyperparameter

# Feature selection (STG)
STG: True
stg_reg: 0.01
input_feat_dims: [7500, 3000, 600, 180, 90]

# Positional encoding
pe_res: 500                                # PE resolution in bp
max_len: 12001
PEU: True                                  # Enable PE upscaling
```

--------------------------------

### Create Stage 2 Datasets

Source: https://context7.com/boevalab/universalepi/llms.txt

Generates NPZ datasets required for Stage 2 training. This process combines Stage 1 parquet outputs with Hi-C target data. Datasets can be created for training, validation, and testing, with output files named according to cell line and data mode.

```bash
# Create training datasets for each cell line
python ./Stage2/create_dataset.py \
    -g ./data/stage1_outputs/predict_gm12878 \
    -s ./data/processed_data/ \
    --hic_data_dir ./data/hic/ \
    --mode train

python ./Stage2/create_dataset.py \
    -g ./data/stage1_outputs/predict_k562 \
    -s ./data/processed_data/ \
    --hic_data_dir ./data/hic/ \
    --mode train

# Create validation datasets
python ./Stage2/create_dataset.py \
    -g ./data/stage1_outputs/predict_gm12878 \
    -s ./data/processed_data/ \
    --hic_data_dir ./data/hic/ \
    --mode val

# Output: ./data/processed_data/{cell_line}_{train|val|test}.npz
# Contains: dnase, sequence, mappability, meta, indexing, target
```

--------------------------------

### Remove Cell-Type Specific Blacklisted Regions - Python

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

This Python script removes cell-type-specific blacklisted genomic regions. It reads a blacklist file specific to a cell line, identifies bins to be dropped based on start positions, and filters the main DataFrame. This is useful for fine-tuning data based on specific cell characteristics.

```python
import pandas as pd
from tqdm import tqdm

# Load cell-type specific blacklist
blacklist = pd.read_csv(f"../data/blacklist/{cell_line.lower()}_blacklist.bed", sep="\t", names=["chr", "start", "end"])
blacklist = blacklist[blacklist.chr.isin(test_chr)]

drop_info = {}
for chr in test_chr:
    blacklist_chr = blacklist[blacklist.chr == chr]
    chr = int(chr.split("chr")[-1])
    blacklist_chr['bin1'] = blacklist_chr.start//5000
    
    drop_list = []
    for _,row in tqdm(blacklist_chr.iterrows()):
        start = row['bin1']
        skip_regions = [start, start+1, start-1]
        drop_list.extend(skip_regions)
    
    drop_info[chr] = drop_list

df_list = []
for chr, drop_list in drop_info.items():
    df_filter = final_df[final_df.chrom==chr]
    df_filter = df_filter[~df_filter.bin1.isin(drop_list)]
    df_filter = df_filter[~df_filter.bin2.isin(drop_list)]
    df_list.append(df_filter)
final_df = pd.concat(df_list, ignore_index=True)
final_df.head()
```

--------------------------------

### Navigate to ATAC-seq Preprocessing Directory

Source: https://github.com/boevalab/universalepi/wiki/Generating-predictions-on-new-data

Change the current directory to the ATAC-seq preprocessing scripts. This is a prerequisite for running the normalization script.

```bash
cd preprocessing/atac
```

--------------------------------

### Load and Analyze Predictions with NumPy

Source: https://context7.com/boevalab/universalepi/llms.txt

Demonstrates how to load prediction results from an NPZ file using NumPy and perform basic analysis. It shows how to access prediction arrays, filter data by chromosome, and convert log-scale predictions to a linear scale.

```python
import numpy as np

# Load predictions
results = np.load('./results/imr90/paper-hg38-map-concat-stage1024-rf-lrelu-eval-stg-newsplit-newdata-atac-var-beta-neg-s1337/results.npz')

# Access prediction arrays
chromosomes = results['chr']      # Shape: (N,) - chromosome numbers
positions_1 = results['pos1']     # Shape: (N,) - peak 1 positions in bp
positions_2 = results['pos2']     # Shape: (N,) - peak 2 positions in bp
log_hic = results['predictions']  # Shape: (N,) - predicted log Hi-C
uncertainty = results['variance'] # Shape: (N,) - aleatoric uncertainty

# Filter for specific chromosome
chr2_mask = chromosomes == 2
chr2_predictions = log_hic[chr2_mask]
chr2_positions = list(zip(positions_1[chr2_mask], positions_2[chr2_mask]))

# Convert log predictions to linear scale
linear_hic = np.exp(log_hic)
```

--------------------------------

### Load and Process Prediction Results (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

Loads prediction results from .npz files, calculates mean predictions and uncertainties (aleatoric and epistemic), and structures them into a pandas DataFrame. It filters results for a specific chromosome and distance between bins.

```python
# Load results and filter for the region of interest
chrom = int((region.split(":")[0]).split('chr')[-1])
seeds = [1337, 2714, 5040, 5452, 5895, 5994, 6286, 8737, 9354, 9597]
pred_dim = 200

def get_results(cell_line):
    variance_list = []
    predictions_list = []
    for i, seed in enumerate(seeds):
        results = np.load(f'../results/{cell_line}/paper-hg38-map-concat-stage1024-rf-lrelu-eval-stg-newsplit-newdata-atac-var-beta-neg-s{seed}/results.npz')
        if i == 0:
            targets = np.reshape(results['targets'], [-1, pred_dim])
        predictions_seed = np.reshape(results['predictions'], [-1, pred_dim])
        variance_seed = np.reshape(results['variance'], [-1, pred_dim])
        variance_list.append(variance_seed)
        predictions_list.append(predictions_seed)
    aleatoric_uncertainty = np.mean(variance_list, axis=0)
    predictions = np.mean(predictions_list, axis=0)
    predictions_std = np.std(predictions_list, axis=0, ddof=1)
    epistemic_uncertainty = predictions_std**2
    variance = aleatoric_uncertainty + epistemic_uncertainty
    predictions = np.clip(predictions, 0, None)

    res = 5000
    chr_data = np.reshape(results['chr'], [-1, pred_dim])
    pos1 = np.reshape(results['pos1'], [-1, pred_dim])//res
    pos2 = np.reshape(results['pos2'], [-1, pred_dim])//res
    print(targets.shape, predictions.shape, chr_data.shape, pos1.shape, pos2.shape)
    # Create result dataframe
    result_df = pd.DataFrame({
        'chrom': chr_data.flatten(),
        'bin1': pos1.flatten(),
        'bin2': pos2.flatten(),
        'target': targets.flatten(),
        'preds': predictions.flatten(),
        'variance': variance.flatten(),
    })
    result_df['dist'] = abs(result_df['bin2']-result_df['bin1'])
    result_df = result_df[result_df.dist < 203]
    result_df = result_df[result_df.chrom == chrom]
    return result_df
```

--------------------------------

### Import Libraries for UniversalEPI Analysis (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

Imports essential Python libraries for data manipulation, numerical operations, plotting, and specialized bioinformatics tools like pyBigWig. It also configures warnings and sets a random seed for reproducibility.

```python
# Import all the necessary libraries
import sys
if 'src' not in sys.path:
    sys.path.append('src')
import os
import warnings
warnings.filterwarnings("ignore")
    
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors
from matplotlib.patches import Polygon

import pyBigWig

np.random.seed(42)
```

--------------------------------

### Configure Cell Lines for UniversalEPI

Source: https://context7.com/boevalab/universalepi/llms.txt

Creates configuration files for new cell lines by defining ATAC-seq data paths. This involves creating a new YAML file in the `dataset/cells` directory and updating the `multicell.yaml` and `cross-cell.yaml` configuration files to include the new cell line.

```yaml
# Create ./Stage1/configs/datamodule/dataset/cells/new_cell.yaml
# Template structure:
species: human
cell_line: new_cell
atac:
  bigwig_path: data/atac/raw/new_cell_normalized.bw
  peaks_path: data/atac/raw/new_cell_dedup.bed
# Then add "new_cell" to ./Stage1/configs/datamodule/multicell.yaml
# And set as prediction cell in ./Stage1/configs/datamodule/validation/cross-cell.yaml
```

--------------------------------

### Prepare and Filter DataFrames (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

Loads and processes prediction results for untreated and treated conditions, sorts them, removes duplicates, and prints their shapes. This prepares the data for merging and further analysis.

```python
untreated_df = get_results('Reed_0000_merged')
treated_df = get_results('Reed_1440_merged')
treated_df.sort_values('preds', ascending=True, inplace=True)
untreated_df.sort_values('preds', ascending=True, inplace=True)
treated_df.drop_duplicates(['chrom', 'bin1', 'bin2'], keep='last', inplace=True)
untreated_df.drop_duplicates(['chrom', 'bin1', 'bin2'], keep='last', inplace=True)
print(untreated_df.shape, treated_df.shape)
```

--------------------------------

### Create Target Dataset (Python)

Source: https://github.com/boevalab/universalepi/wiki/Training-UniversalEPI

Combines ATAC-seq and Hi-C data to extract targets for training. Adds pseudopeaks and saves updated ATAC-seq peaks. Assumes 5Kb Hi-C resolution and hg38 genome by default.

```python
python ./preprocessing/prepare_target_data.py --cell_line <cell_line> --atac_bed_path ./data/atac/raw/<cell_line>_dedup.bed --hic_data_dir ./data/hic/
```

--------------------------------

### Download and Unzip HepG2 Predictions (Bash)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

Downloads prediction data for the HepG2 cell line from a Zenodo repository and unzips the downloaded file into a specified results directory. This step prepares the necessary data for subsequent analysis.

```bash
# Download HepG2 predictions from Zenodo
download_path = f"https://zenodo.org/records/15079440/files/{cell_line}.zip?download=1"
os.makedirs('../results', exist_ok=True)
os.system(f"wget {download_path} -O ../results/{cell_line}.zip")
os.system(f"unzip ../results/{cell_line}.zip -d ../results/")
```

--------------------------------

### Convert Hi-C Data to Sparse Pairwise Interactions

Source: https://context7.com/boevalab/universalepi/llms.txt

Converts Hi-C contact matrices from either `.hic` or `.cool` format into sparse pairwise interaction files. This script supports ICE normalization and outputs files in a format suitable for downstream processing.

```bash
cd preprocessing/hic

# Convert from .hic format with ICE normalization
./hic2sparse.sh ../../data/hic/K562.hic ../../data/hic/k562 5000 --ice

# Convert from .cool format with ICE normalization
./cool2sparse.sh ../../data/hic/K562.cool ../../data/hic/k562 --ice

# Output: ../../data/hic/k562/raw_iced/chr{1-22}_raw.bed
# Format: pos1<tab>pos2<tab>hic_score
```

--------------------------------

### Test Stage 2 Model

Source: https://context7.com/boevalab/universalepi/llms.txt

Evaluates the trained Stage 2 model on test cell lines using ground truth Hi-C data. Configuration requires setting the test cell lines and optionally specifying a directory for pre-generated test NPZ files. The output includes predictions, variance, and ground truth values.

```bash
# Configure in ./Stage2/configs/configs.yaml:
#   cell_lines_test: hepg2
#   test_dir: null  # or path to pre-generated test NPZ

python ./Stage2/main.py \
    --config_dir ./Stage2/configs/configs.yaml \
    --mode test

# Output: ./results/hepg2/{run_name}/results.npz
# Contains:
#   - chr: Chromosome numbers
#   - pos1: Position of ATAC peak 1
#   - pos2: Position of ATAC peak 2
#   - predictions: Log Hi-C predictions
#   - variance: Aleatoric uncertainty
#   - targets: Ground truth log Hi-C values
```

--------------------------------

### Create and Display Color Bar for Confidence Fold Change

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

This snippet demonstrates how to create a color bar axis and display it using Matplotlib. It's used to visually represent the 'Max. Confidence Fold Change (treated / untreated)' based on a specified colormap and normalization.

```python
cbar_ax = fig.add_axes([0.9, 0.15, 0.025, 0.7])
sm = cm.ScalarMappable(cmap=cmap, norm=norm)
plt.colorbar(sm, cax=cbar_ax, label="Max. Confidence Fold Change (treated / untreated)")
plt.show()
```

--------------------------------

### Prepare Target Data for Training and Testing

Source: https://context7.com/boevalab/universalepi/llms.txt

Combines ATAC-seq peaks with Hi-C interaction scores to generate training targets. For training, it includes negative sampling by adding pseudopeaks. For testing, it applies Hi-C smoothing without pseudopeaks. Requires specifying cell line, ATAC-seq peak paths, Hi-C data directory, and genomic parameters.

```bash
# For training cell lines (adds 10% pseudopeaks for negative sampling)
python ./preprocessing/prepare_target_data.py \
    --cell_line gm12878 \
    --atac_bed_path ./data/atac/raw/GM12878_dedup.bed \
    --hic_data_dir ./data/hic/ \
    --resolution 5000 \
    --max_len 2000000 \
    --genome hg38

# For test cell lines (applies Hi-C smoothing, no pseudopeaks)
python ./preprocessing/prepare_target_data.py \
    --cell_line hepg2 \
    --atac_bed_path ./data/atac/raw/HEPG2_dedup.bed \
    --hic_data_dir ./data/hic/ \
    --test

# Output:
# - data/atac/raw/{cell_line}_dedup_neg.bed (peaks with pseudopeaks)
# - data/hic/{cell_line}/target_data/chr{N}_target.pkl (interaction matrices)
```

--------------------------------

### Estimate Epistemic Uncertainty with Model Checkpoints (Python)

Source: https://context7.com/boevalab/universalepi/llms.txt

This snippet demonstrates how to compute epistemic uncertainty by running predictions across multiple model checkpoints. It loads predictions from various seeds, stacks them, and calculates the variance to represent epistemic uncertainty, combining it with aleatoric uncertainty.

```python
import numpy as np
import os

# Run predictions with each of 10 model checkpoints
cell_line = 'imr90'
seeds = [1337, 42, 123, 456, 789, 1234, 5678, 9012, 3456, 7890]
all_predictions = []

for seed in seeds:
    result_path = f'./results/{cell_line}/paper-hg38-map-concat-stage1024-rf-lrelu-eval-stg-newsplit-newdata-atac-var-beta-neg-s{seed}/results.npz'
    results = np.load(result_path)
    all_predictions.append(results['predictions'])

# Stack predictions and compute epistemic uncertainty
predictions_stack = np.stack(all_predictions, axis=0)  # Shape: (10, N)
mean_prediction = np.mean(predictions_stack, axis=0)
epistemic_uncertainty = np.var(predictions_stack, axis=0)

# Total uncertainty = aleatoric + epistemic
aleatoric_uncertainty = results['variance']  # From any single run
total_uncertainty = aleatoric_uncertainty + epistemic_uncertainty
```

--------------------------------

### Train Stage 2 Model

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Trains the Stage 2 model of UniversalEPI. Requires correct configuration of genomic data and HiC paths in `configs.yaml`. Supports training with pre-generated npz files if `create_dataset.py` and `merge_dataset.py` have been executed.

```python
python ./Stage2/main.py --config_dir ./Stage2/configs/configs.yaml --mode train
```

--------------------------------

### Test Stage 1 Model

Source: https://context7.com/boevalab/universalepi/llms.txt

Tests the trained Stage 1 model on held-out cell lines. The cell line for prediction is configured in the cross-cell validation YAML file. This step evaluates the model's performance on unseen data.

```bash
# Configure test cell line in:
# ./Stage1/configs/datamodule/validation/cross-cell.yaml
# Set: predict: hepg2

python ./Stage1/test.py
```

--------------------------------

### Test Stage 1 Model

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Tests the Stage 1 model of UniversalEPI. It uses test cell lines defined in the `cross_cell.yaml` configuration file.

```python
python ./Stage1/test.py
```

--------------------------------

### Prepare Target Data for Training and Testing

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Processes Hi-C and ATAC-seq data to extract target interactions for training and testing UniversalEPI. It handles normalization and combines data for specific cell lines, adding pseudopeaks for negative sampling.

```python
python ./preprocessing/prepare_target_data.py --cell_line gm12878 --atac_bed_path ./data/atac/raw/GM12878_dedup.bed --hic_data_dir ./data/hic/
python ./preprocessing/prepare_target_data.py --cell_line hepg2 --atac_bed_path ./data/atac/raw/HEPG2_dedup.bed --hic_data_dir ./data/hic/ --test
```

--------------------------------

### Convert .cool to Pairwise Interactions (Shell)

Source: https://github.com/boevalab/universalepi/blob/main/preprocessing/hic/README.md

Converts a .cool file to pairwise interaction files (tab-separated: pos1, pos2, hic_score) for each chromosome. Supports ICE normalization. The output is stored in a specified directory.

```shell
./cool2sparse.sh ../../data/hic/HEPG2.cool ../../data/hic/hepg2 --ice

```

--------------------------------

### Train Stage 1 Model

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Trains the Stage 1 model of UniversalEPI. It utilizes training cell lines specified in the `cross_cell.yaml` configuration file.

```python
python ./Stage1/train.py
```

--------------------------------

### Train Stage 2 Transformer Model

Source: https://context7.com/boevalab/universalepi/llms.txt

Trains the Stage 2 transformer encoder model for Hi-C prediction, incorporating uncertainty estimation. Configuration involves setting paths to Hi-C and ATAC data, specifying cell lines for training and validation, and pointing to the Stage 1 model checkpoint. Key training parameters include epochs, batch size, learning rate, and model architecture details.

```bash
# Configure paths in ./Stage2/configs/configs.yaml:
#   hic_path: data/hic/
#   atac_path: data/stage1_outputs/
#   cell_lines_train: ["gm12878", "k562"]
#   cell_lines_val: ["gm12878", "k562"]
#   stage_1_model: checkpoints/stage1_model_gm12878_k562.pth

python ./Stage2/main.py \
    --config_dir ./Stage2/configs/configs.yaml \
    --mode train

# Key training parameters (configs.yaml):
#   n_epochs: 20
#   batch_size: 32
#   lr: 0.001
#   num_heads: 4
#   num_layers: 4
#   embed_dim: 32
#   hidden_dim: 96
#   var_flg: True        # Enable aleatoric uncertainty
#   STG: True            # Enable feature selection
```

--------------------------------

### Test Stage 2 Model

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Tests the Stage 2 model of UniversalEPI. Requires correct configuration of genomic data and test directory paths in `configs.yaml`. Generates a `.npz` file with predictions, variance, and targets. Evaluation plots can be generated using `plot_scores.ipynb`.

```python
python ./Stage2/main.py --config_dir ./Stage2/configs/configs.yaml --mode test
```

--------------------------------

### Merge and Clean DataFrames (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

Merges the untreated and treated dataframes based on genomic coordinates and removes rows with any missing values. It then prints the shape of the merged dataframe.

```python
# Merge the two dataframes
merged_df = pd.merge(untreated_df, treated_df, on=['chrom', 'bin1', 'bin2'], suffixes=('_untreated', '_treated'))
merged_df.dropna(inplace=True)
print(merged_df.shape)
```

--------------------------------

### Store Genomic Inputs for Stage 1

Source: https://github.com/boevalab/universalepi/wiki/Generating-predictions-on-new-data

Execute the script to store genomic inputs for a specified cell line. This script processes DNA-sequence, ATAC-seq, and mappability data, saving them as parquet files. It uses all chromosomes by default, but a subset can be specified in the configuration file.

```python
python ./Stage1/store_inputs.py --cell_line <cell_line>
```

--------------------------------

### Process Predictions and Calculate Uncertainties (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

Loads prediction results from multiple seeds, calculates mean predictions, aleatoric uncertainty, and epistemic uncertainty. It also extracts chromosome and position information, reshaping arrays for further processing.

```python
pred_dim = 200
variance_list = []
predictions_list = []
for i, seed in enumerate(seeds):
    results = np.load(f'../results/{cell_line}/paper-hg38-map-concat-stage1024-rf-lrelu-eval-stg-newsplit-newdata-atac-var-beta-neg-s{seed}/results.npz')
    if i == 0:
        targets = np.reshape(results['targets'], [-1, pred_dim])
    predictions_seed = np.reshape(results['predictions'], [-1, pred_dim])
    variance_seed = np.reshape(results['variance'], [-1, pred_dim])
    variance_list.append(variance_seed)
    predictions_list.append(predictions_seed)
aleatoric_uncertainty = np.mean(variance_list, axis=0)
predictions = np.mean(predictions_list, axis=0)
predictions_std = np.std(predictions_list, axis=0, ddof=1)
epistemic_uncertainty = predictions_std**2
variance = aleatoric_uncertainty + epistemic_uncertainty
predictions = np.clip(predictions, 0, None)

res = 5000
chr_data = np.reshape(results['chr'], [-1, pred_dim])
pos1 = np.reshape(results['pos1'], [-1, pred_dim])//res
pos2 = np.reshape(results['pos2'], [-1, pred_dim])//res
print(targets.shape, predictions.shape, chr_data.shape, pos1.shape, pos2.shape)
```

--------------------------------

### Plot Differential Arcs for a Given Region

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

This function, `plot_arcs_differential`, is designed to visualize differential arcs for a specified region. It likely takes a region object as input and generates corresponding arc plots, crucial for analyzing spatial or temporal relationships in epidemiological data.

```python
plot_arcs_differential(region)
```

--------------------------------

### Generate Stage 2 Hi-C Predictions

Source: https://context7.com/boevalab/universalepi/llms.txt

Generates Hi-C predictions for new cell types using pretrained Stage 2 models, without requiring ground truth data. Predictions can be made for all chromosomes or specific ones. The output NPZ file contains chromosome numbers, peak positions, predicted log Hi-C scores, and uncertainty estimates.

```bash
# Predict for all chromosomes
python ./Stage2/predict.py \
    --config_dir ./Stage2/configs/configs.yaml \
    --cell_line_predict imr90

# Predict for specific chromosomes
python ./Stage2/predict.py \
    --config_dir ./Stage2/configs/configs.yaml \
    --cell_line_predict imr90 \
    --chroms_predict 2 6 19

# Output: ./results/imr90/{run_name}/results.npz
# Contains:
#   - chr: Chromosome numbers (1-22)
#   - pos1: Genomic position of first ATAC peak (bp)
#   - pos2: Genomic position of second ATAC peak (bp)
#   - predictions: Log Hi-C interaction scores
#   - variance: Aleatoric uncertainty estimates
```

--------------------------------

### Generate Scatter Plot on Test Data - Python

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

This Python script generates a scatter plot to visualize the correlation between predictions and target values on test data. It uses a subset of the data for clarity, employs Gaussian KDE for density estimation, and displays the plot with Pearson (R) and Spearman (ρ) correlation coefficients. This provides a visual assessment of the model's predictive accuracy.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import pearsonr, spearmanr
from scipy.stats import gaussian_kde

t_b = np.array(final_df['target'].values.tolist())
t_hat_b = np.array(final_df['preds'].values.tolist())
r = np.round(pearsonr(t_b, t_hat_b)[0],3)
rho = np.round(spearmanr(t_b, t_hat_b)[0],3)
idx_subset = np.random.choice(np.arange(len(t_b)), size=50000, replace=False)
t_b_subset = t_b[idx_subset]
t_hat_b_subset = t_hat_b[idx_subset]
xy = np.vstack([t_hat_b_subset, t_b_subset])
z = gaussian_kde(xy)(xy)
idx = z.argsort()
t_hat_b_subset, t_b_subset, z = t_hat_b_subset[idx], t_b_subset[idx], z[idx]
plt.scatter(t_hat_b_subset, t_b_subset, c=z, s=2)
plt.title(f"UniversalEPI ({cell_line}): R={r}; ρ:{rho}")
plt.xlabel("Predictions")
plt.ylabel("Targets")
plt.colorbar()
plt.show()
```

--------------------------------

### Plot Differential Genomic Arcs and Signals (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/differential_interactions.ipynb

This Python function visualizes differential Hi-C interactions as arcs, overlaid with ATAC-seq signals and gene annotations for a given genomic region. It requires pre-processed dataframes for interactions and genes, and ATAC-seq data from normalized bigWig files. The output is a matplotlib figure with four subplots.

```python
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.cm as cm
import pyBigWig
import numpy as np
from matplotlib.patches import Polygon


def plot_arcs_differential(region, gene_df=gene_df, result_df=merged_df_filtered):
    
    res = 5000
    chr, start, end = int((region.split(":")[0]).split('chr')[-1]), \
        int(region.split(":")[1].split("-")[0]), int(region.split(":")[1].split("-")[1])
    bin_start = start//res
    bin_end = end//res
    df_region = result_df[(result_df.chrom == chr) & (result_df.bin1 >= bin_start) & (result_df.bin2 <= bin_end)]
    df_region['pos1'] = df_region['bin1']*res
    df_region['pos2'] = df_region['bin2']*res
    

    # Normalize preds values for shading
    norm = mcolors.Normalize(vmin=-0.5, vmax=0.5)    
    cmap = cm.RdBu
    
    # Define genomic range
    chrom = "chr"+str(int(df_region["chrom"].iloc[0]))

    # Load ATAC-seq bigWig
    bw1 = pyBigWig.open(f"../data/atac/raw/0000_normalized.bw") 
    bw2 = pyBigWig.open(f"../data/atac/raw/1440_normalized.bw")
    
    positions = np.arange(start, end, 100, dtype=int)  # Sample every 100bp
    atac_signal1 = np.array([bw1.values(chrom, pos, pos + 100) for pos in positions])
    atac_signal1 = np.nan_to_num(atac_signal1)  # Replace NaNs with zeros
    bw1.close()

    atac_signal2 = np.array([bw2.values(chrom, pos, pos + 100) for pos in positions])
    atac_signal2 = np.nan_to_num(atac_signal2)  # Replace NaNs with zeros
    bw2.close()

    atac_signal1 = atac_signal1.max(axis=1)
    atac_signal2 = atac_signal2.max(axis=1)

    # Normalize genomic positions for plotting
    positions_kb = positions / 1000

    # Create figure with two subplots
    fig, axs = plt.subplots(4, 1, figsize=(12, 10), sharex=True, gridspec_kw={"height_ratios": [3, 1, 1, 0.5]}, dpi=100)

    # Plot 1: Hi-C Arcs
    axs[0].hlines(y=0, xmin=start/1000, xmax=end/1000, color='black', linewidth=1)
    for _, row in df_region.iterrows():
        x1, x2 = row["pos1"] / 1000, row["pos2"] / 1000
        arc_height = (x2 - x1) / 2  
        arc = np.linspace(0, np.pi, 100)
        arc_x = np.linspace(x1, x2, 100)
        arc_y = np.sin(arc) * arc_height
        color = cmap(norm(row["max_conf_fc"]))
        axs[0].plot(arc_x, arc_y, color=color, alpha=0.9)
    axs[0].set_title("{}: {:,}—{:,}".format(chrom, start, end))
    axs[0].set_ylabel("Predicted Differential Interactions")
    axs[0].set_yticks([])
    axs[0].spines["top"].set_visible(False)
    axs[0].spines["right"].set_visible(False)
    axs[0].spines["bottom"].set_visible(False)
    axs[0].spines["left"].set_visible(False)

    # Plot 2: ATAC-seq Signal
    axs[1].fill_between(positions_kb, atac_signal1, color="red", alpha=0.7)
    axs[1].set_ylabel("ATAC-seq\n(before treatment)")
    axs[1].set_xlim(start / 1000, end / 1000)  # Keep X-axis consistent
    axs[1].spines["top"].set_visible(False)
    axs[1].spines["right"].set_visible(False)
    axs[1].spines["bottom"].set_visible(False)

    # Plot 3: ATAC-seq Signal
    axs[2].fill_between(positions_kb, atac_signal2, color="blue", alpha=0.7)
    axs[2].set_ylabel("ATAC-seq\n(after treatment)")
    axs[2].set_xlim(start / 1000, end / 1000)  # Keep X-axis consistent
    axs[2].spines["top"].set_visible(False)
    axs[2].spines["right"].set_visible(False)
    axs[2].spines["bottom"].set_visible(False)

    # Plot 4: Gene Annotations
    gene_df = gene_df[gene_df.chr == chrom]
    gene_df = gene_df[(gene_df.end >= start) & (gene_df.start <= end)]
    for i, row in gene_df.iterrows():
        y_pos = (i+1) % 2 * 0.3  # Alternates between 0 and 0.3 to stagger labels
        j=0.25
        axs[3].hlines(y=j+y_pos, xmin=row["start"] / 1000, xmax=row["end"] / 1000, color="black", linewidth=3)
        if row['tss'] < start or row['tss'] > end:
            pass
        else:
            axs[3].text(row['tss']/1000, j-0.15-y_pos, row["gene_name"], ha="center", va="top", fontsize=8)
        tss_kb = row["tss"] / 1000
        if row["strand"] == '+':  
            triangle_coords = [(tss_kb, j+y_pos), (tss_kb, j+y_pos + 0.4), (tss_kb + 0.3, j+y_pos)]  # Right-facing
        else:
            triangle_coords = [(tss_kb, j+y_pos), (tss_kb, j+y_pos + 0.4), (tss_kb - 0.3, j+y_pos)]  # Left-facing

        # Add triangle to plot
        triangle = Polygon(triangle_coords, closed=True, facecolor="black", edgecolor="black")
        axs[3].add_patch(triangle)

    axs[3].set_yticks([])
    axs[3].set_xlabel("Genomic Position (Kb)")
    axs[3].set_xlim(start / 1000, end / 1000)  # Keep X-axis consistent
    axs[3].set_ylim(-1, 1)
    axs[3].set_ylabel("Genes")
    axs[3].spines["top"].set_visible(False)
    axs[3].spines["right"].set_visible(False)
    axs[3].spines["left"].set_visible(False)


    # Remove top subplot's X-axis ticks for cleaner look
    axs[0].tick_params(axis="x", which="both", bottom=False, labelbottom=False)
    axs[1].tick_params(axis="x", which="both", bottom=False, labelbottom=False)
    axs[2].tick_params(axis="x", which="both", bottom=False, labelbottom=False)

    # Add colorbar to whole figure
    fig.subplots_adjust(right=0.89)

```

--------------------------------

### Normalize ATAC-seq Data with Python

Source: https://github.com/boevalab/universalepi/wiki/Generating-predictions-on-new-data

Run the ATAC-seq normalization script. It supports input from BAM files or BigWig and peak files. The script normalizes ATAC-seq data using GM12878 as a reference and generates normalized bigwig and deduplicated peak files.

```python
python normalize_atac.py -p ../../data/atac/raw/ --input_bam ../../data/atac/raw/<cell_line1>.bam ../../data/atac/raw/<cell_line2>.bam
```

```python
python normalize_atac.py -p ../../data/atac/raw/ --input_bw ../../data/atac/raw/<cell_line1>.bigWig ../../data/atac/raw/<cell_line2>.bigWig --input_bed ../../data/atac/raw/<cell_line1>.bed ../../data/atac/raw/<cell_line2>.bed
```

--------------------------------

### Import Libraries for UniversalEPI Evaluation (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

Imports necessary Python libraries for data manipulation, scientific computing, plotting, and specialized bioinformatics tools. It also configures warning filters and sets a random seed for reproducibility.

```python
# Import all the necessary libraries
import sys
if 'src' not in sys.path:
    sys.path.append('src')
import os
import warnings
warnings.filterwarnings("ignore")
    
import numpy as np
import pandas as pd
from scipy.stats import kendalltau,spearmanr,pearsonr,gaussian_kde

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors
from matplotlib.patches import Polygon

import pyBigWig
from tqdm import tqdm

np.random.seed(42)
```

--------------------------------

### Generate Hi-C Predictions with Python

Source: https://github.com/boevalab/universalepi/wiki/Generating-predictions-on-new-data

Run the prediction script for Hi-C data generation. This script requires the ATAC path to be correctly set in its configuration file. It can optionally predict for a subset of chromosomes. The output includes predicted Hi-C values and their associated uncertainties.

```python
python ./Stage2/predict.py --config_dir ./Stage2/configs/configs.yaml --cell_line_predict <cell_line>
```

```python
python ./Stage2/predict.py --config_dir ./Stage2/configs/configs.yaml --cell_line_predict <cell_line> --chroms_predict 2 6 19
```

--------------------------------

### Generate Hi-C Predictions (Stage 2)

Source: https://github.com/boevalab/universalepi/blob/main/README.md

Generates Hi-C predictions from Stage 2 outputs. Requires correct configuration of `atac_path` in `configs.yaml`. Can optionally select a subset of chromosomes for prediction. Outputs a `.npz` file containing chromosome, peak positions, predictions, and variance.

```python
python ./Stage2/predict.py --config_dir ./Stage2/configs/configs.yaml --cell_line_predict imr90
```

```python
python ./Stage2/predict.py --config_dir ./Stage2/configs/configs.yaml --cell_line_predict imr90 --chroms_predict 2 6 19
```

--------------------------------

### Distance Stratified Correlation Analysis (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

This script calculates and visualizes Pearson and Spearman correlations between predicted and actual values, stratified by genomic distance. It iterates through distance bins, computes correlation coefficients and sample counts, and then plots these metrics against genomic distance in kilobases.

```python
r_d_200 = []
rho_d_200 = []
n_samples = []

for index in range(np.uint(1000000/res)):
    t_b_i = np.array(final_df[final_df.dist == index]['target'].values.tolist())
    t_hat_b_i = np.array(final_df[final_df.dist == index]['preds'].values.tolist())

    n_samples.append(len(np.array(t_b_i).flatten()))
    
    r_d_200.append(pearsonr(np.array(t_b_i).flatten(), np.array(t_hat_b_i).flatten())[0])
    rho, _ = spearmanr(np.array(t_b_i).flatten(), np.array(t_hat_b_i).flatten())
    rho_d_200.append(rho)

x_label = np.arange(0,len(r_d_200)) * res / 1000

f, ax = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))
ax[0].plot(x_label, np.array(r_d_200), label="UniversalEPI")
ax[0].grid(axis='y')
ax[0].set_ylim(-0.1,1)
ax[0].set_xlabel("Genomic Distance [kb]")
ax[0].set_ylabel("Pearson Correlation")

ax[1].plot(x_label, np.array(rho_d_200), label="UniversalEPI")
ax[1].grid(axis='y')
ax[1].set_ylim(-0.1,1)
ax[1].set_xlabel("Genomic Distance [kb]")
ax[1].set_ylabel("Spearman Correlation")

ax[2].plot(x_label, np.array(n_samples), label="UniversalEPI")
ax[2].grid(axis='y')
ax[2].set_xlabel("Genomic Distance [kb]")
ax[2].set_ylabel("Counts")

lines, labels = ax[1].get_legend_handles_labels()
f.legend(lines, labels, loc="upper right", title="method", fontsize='large', title_fontsize='large')
f.suptitle("Distance Stratified Correlation")
```

--------------------------------

### Define Cell Line and Seeds for Prediction (Python)

Source: https://github.com/boevalab/universalepi/blob/main/Stage2/plot_scores.ipynb

Sets the target cell line to 'hepg2' and defines a list of seeds to be used for generating predictions. These parameters are crucial for reproducible analysis.

```python
# Cell line of interest
cell_line = "hepg2"
seeds = [1337, 2714, 5040, 5452, 5895, 5994, 6286, 8737, 9354, 9597]
```