# README

Generated on: `2023-03-29`
MACS2 version: `macs2 2.2.7.1`

Generated narrowPeaks (`macs2 bdgpeakcall -c 2 --no-trackline`) and broadPeaks (`macs2 bdgbroadcall -c 2 -C 1`) from imputed signal tracks for the IHEC dataset.

Any peaks overlapping the ENCODE blacklist (hg38, accession: `ENCFF356LFX`) are filtered out (`bedtools subtract -A`). 

Peaks are ranked by the `score` field (BED column 5). Ties are broken by genomic coordinate sort order.

Questions: dincer@ucla.edu

## Downloading data and setting up environment
- Download ENCODE blacklist (ENCFF356LFX)
- Download bigWigToBedGraph
- Set up different version of MACS2 (optional)

```bash
# Download encode blacklist hg38 (ENCFF356LFX)
wget https://www.encodeproject.org/files/ENCFF356LFX/@@download/ENCFF356LFX.bed.gz
gunzip ENCFF356LFX.bed
mv ENCFF356LFX.bed encode_blacklist_hg38.ENCFF356LFX.bed

# Download bigWigToBedGraph
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bigWigToBedGraph
chmod +x ./bigWigToBedGraph

# Set up different version of MACS2 (optional)
# mamba create -c conda-forge -c bioconda --prefix ./conda_envs/ihec_macs210 macs2=2.1.0
```

Setting up the conda environment: `conda env create -f environment.yml`

```yaml
# environment.yml
name: ihec_imputed_peaks
channels:
    - conda-forge
    - bioconda
dependencies:
    - python=3.10
    - ucsc-bigwigtobedgraph
    - numpy
    - pandas
    - bedtools
    - pyBedTools
    - seaborn
    - tqdm
    - macs2=2.2.7.1
```

## Calling peaks from imputed signal tracks
Calling peaks script:
```py
# python ihec_imputed_peak_calling.py
import subprocess
import tempfile
from pathlib import Path

import pandas as pd

blacklist_fn = "encode_blacklist_hg38.ENCFF356LFX.bed"
path_bigWigToBedGraph = "./bigWigToBedGraph"
macs2_version = "2.2.7.1"


def generate_peaks_imputed_bw(
    bw_fn,
    out_prefix_fn,
    blacklist_fn,
    path_bigWigToBedGraph="bigWigToBedGraph",
    conda_run_cmd=None,
):
    if conda_run_cmd is None:
        conda_run_cmd = ""

    with tempfile.TemporaryDirectory() as tmpdir:
        bw_in_name = Path(bw_fn).name

        out_narrowpeak_fn = f"{out_prefix_fn}/narrowPeak/{bw_in_name}.narrowPeak"
        out_broadpeak_fn = f"{out_prefix_fn}/broadPeak/{bw_in_name}.broadPeak"

        Path(out_narrowpeak_fn).parent.mkdir(exist_ok=True, parents=True)
        Path(out_broadpeak_fn).parent.mkdir(exist_ok=True, parents=True)

        # generate BEDGRAPH
        bedg_out = f"{tmpdir}/{bw_in_name}.bedgraph"
        cmd_ = f"{path_bigWigToBedGraph} {bw_fn} {bedg_out}"
        print(cmd_)
        subprocess.run(cmd_, shell=True, check=True)

        # generate NARROWPEAKS
        narrowpeak_tmp_fn = f"{tmpdir}/{bw_in_name}.tmp.narrowpeak"
        cmd_ = f"{conda_run_cmd} macs2 bdgpeakcall -i {bedg_out} -c 2 --no-trackline --ofile {narrowpeak_tmp_fn}".strip()
        print(cmd_)
        subprocess.run(cmd_, shell=True, check=True)

        # filter out blacklist from narrowpeaks
        cmd_ = f"bedtools subtract -A -a {narrowpeak_tmp_fn} -b {blacklist_fn} > {out_narrowpeak_fn}"
        print(cmd_)
        subprocess.run(cmd_, shell=True, check=True)

        # rename narrowpeaks by order
        df = pd.read_csv(out_narrowpeak_fn, sep="\t", header=None)
        i_ = df.sort_values([4, 0, 1, 2], ascending=[False, True, True, True]).index
        df.loc[i_, 3] = range(df.shape[0])
        df.loc[:, 3] = "nPk_" + df.loc[:, 3].map(str)
        df.to_csv(out_narrowpeak_fn, sep="\t", header=None, index=None)

        # generate BROADPEAKS
        broadpeak_tmp_fn = f"{tmpdir}/{bw_in_name}.tmp.broadpeak"
        cmd_ = f"{conda_run_cmd} macs2 bdgbroadcall -i {bedg_out} -c 2 -C 1 --ofile {broadpeak_tmp_fn}".strip()
        print(cmd_)
        subprocess.run(cmd_, shell=True, check=True)

        # filter out blacklist from broadpeaks
        cmd_ = f"bedtools subtract -A -a {broadpeak_tmp_fn} -b {blacklist_fn} > {out_broadpeak_fn}"
        print(cmd_)
        subprocess.run(cmd_, shell=True, check=True)

        # rename peaks by order
        df = pd.read_csv(out_broadpeak_fn, sep="\t", header=None)
        i_ = df.sort_values([4, 0, 1, 2], ascending=[False, True, True, True]).index
        df.loc[i_, 3] = range(df.shape[0])
        df.loc[:, 3] = "bPk_" + df.loc[:, 3].map(str)
        df.to_csv(out_broadpeak_fn, sep="\t", header=None, index=None)


if __name__ == "__main__":
    import sys

    bw_fn = sys.argv[1]
    peak_out_dir = sys.argv[2]

    # conda environment containing desired macs2 version, default: 2.2.7.1
    macs2_versions_prefix_mapping = {
        "2.2.7.1": None,
        "2.1.0.20150731": "conda_envs/ihec_macs210",
        "2.1.1.20160309": "conda_envs/ihec_macs211",
    }

    macs2_conda_prefix = macs2_versions_prefix_mapping[macs2_version]
    conda_run_cmd = (
        f"conda run --prefix {macs2_conda_prefix} " if macs2_conda_prefix else ""
    )

    # print macs2 version
    subprocess.run(
        f"{conda_run_cmd} macs2 --version",
        check=True,
        text=True,
        shell=True,
    )

    generate_peaks_imputed_bw(
        bw_fn,
        peak_out_dir,
        blacklist_fn=blacklist_fn,
        path_bigWigToBedGraph=path_bigWigToBedGraph,
        conda_run_cmd=conda_run_cmd,
    )
```

Generating checksums:
```sh
md5sum broadPeak/*.broadPeak > broadPeak/checksum.md5
md5sum narrowPeak/*.narrowPeak > narrowPeak/checksum.md5
```