<jupyter_start><jupyter_text>ATAC - DNase - EcoDAM Comparison Let's try to find the overlapping regions of all three methods in a Venn style way. We'll first normalize the EcoDAM data acording to the (smoothed) theoretical value.<jupyter_code>%load_ext autoreload
%autoreload 2
%load_ext nb_black

import pathlib

import seaborn as sns
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pybedtools

from ecodam_py.bedgraph import BedGraphAccessor, equalize_loci
from ecodam_py.eco_atac_normalization import *
dnase = pybedtools.BedTool(
    "/mnt/saphyr/Saphyr_Data/DAM_DLE_VHL_DLE/Michael/R/Refrences/Wu_merged_hotspots_as_windows.sort.bed"
)

dnase = (
    dnase.to_dataframe()
    .rename({"chrom": "chr", "start": "start_locus", "end": "end_locus"}, axis=1)
    .reset_index(drop=True)
    .assign(intensity=1)
)
dnase
nfr_all_chrom = pathlib.Path(
    "/mnt/saphyr/Saphyr_Data/DAM_DLE_VHL_DLE/Hagai/hmm_output/ENCFF240YRV.sorted.nfr.bedgraph"
)
nfr_all_chrom = read_bedgraph(nfr_all_chrom).query("intensity != 0")
nfr_all_chrom<jupyter_output><empty_output><jupyter_text>Before performing the intersection we want to use the data from the peak calling procedure rather than the entire dataset, as per Michal's advice.<jupyter_code>peak_called_atac = read_bedgraph(
    pathlib.Path(
        "/mnt/saphyr/Saphyr_Data/DAM_DLE_VHL_DLE/Hagai/hmm_output/atac_peaks.bedgraph"
    )
)
peak_called_atac<jupyter_output><empty_output><jupyter_text>We'lll do an unweighted intersection - we're satisfied even if a small portion of the peak touches on an NFR.<jupyter_code>nfr = []
peaks = []
for chr_, nfr_group, peak_group in iter_over_bedgraphs_chromosomes(
    nfr_all_chrom, peak_called_atac
):
    cur_peak, cur_nfr = peak_group.bg.unweighted_overlap(nfr_group)
    nfr.append(cur_nfr.bg.index_to_columns())
    peaks.append(cur_peak.bg.index_to_columns())
nfr_all_chrom = pd.concat(nfr, axis=0, ignore_index=True)<jupyter_output><empty_output><jupyter_text>Now we'll find the intersection of ATAC and DNase, and we'll finally intersect that with our own data:<jupyter_code>total_intersection = {}
for overlap_pct in [0.0001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95]:
    dnase_intersection, atac_intersection = [], []
    for chr_, dnase_group, atac_group in iter_over_bedgraphs_chromosomes(
        dnase, nfr_all_chrom
    ):
        cur_dnase, cur_atac = dnase_group.bg.weighted_overlap(
            atac_group, overlap_pct=overlap_pct
        )
        dnase_intersection.append(cur_dnase.bg.index_to_columns())
        atac_intersection.append(cur_atac.bg.index_to_columns())
    dnase_intersection = pd.concat(dnase_intersection, axis=0, ignore_index=True)
    atac_intersection = pd.concat(atac_intersection, axis=0, ignore_index=True)
    total_intersection[overlap_pct] = len(atac_intersection) / len(nfr_all_chrom)
    print(f"Finished {overlap_pct}")
total_intersection
len(dnase_intersection) / len(dnase)
total_intersection_df = pd.DataFrame.from_dict(
    total_intersection, orient="index", columns=["Overlap percentage requirement"]
)
ax0 = total_intersection_df.reset_index().plot.scatter(
    x="index", y="Overlap percentage requirement"
)
ax0.set_ylabel("ATAC NFR Coverage %")
ax0.set_xlabel("Overlap %")
ax0.set_title("ATAC NFR Coverage as a function of overlap requirement with DNase")<jupyter_output><empty_output><jupyter_text>In broad strokes, all DNase areas are covered with ATAC NFR areas, while the DNase only covers about 7% of the NFR. We can check the data on a chromosome-by-chromosome level:<jupyter_code>atac_per_chr = {}
for chr_, dnase_group, atac_group in iter_over_bedgraphs_chromosomes(
    dnase, nfr_all_chrom
):
    dnase_int_grouped = dnase_intersection.query("chr == @chr_")
    atac_int_grouped = atac_intersection.query("chr == @chr_")
    atac_per_chr[chr_] = len(atac_int_grouped) / len(atac_group)

atac_per_chr_df = pd.DataFrame.from_dict(
    atac_per_chr, orient="index", columns=["Overlap fraction"]
)
ax = atac_per_chr_df.plot.bar()
ax.set_ylabel("Overlap fraction")
ax.set_title("ATAC vs DNase overlap")<jupyter_output><empty_output><jupyter_text>Not too informative, but definitely reassuring.  Adding the EcoDAM data To do that we'll need to normalize that data with the theoretical site density first.<jupyter_code>smoothed_theo_fname = pathlib.Path(
    "/mnt/saphyr/Saphyr_Data/DAM_DLE_VHL_DLE/Theoretical_EcoDam/New folder/Hg38.EcoDam.insilico.canonical.keynames.bg_smoothed_7kb_resampled_with_700_overlapping_bp.bedgraph"
)
smoothed_theo = read_bedgraph(smoothed_theo_fname)
smoothed_theo
chromatin = read_bedgraph(
    pathlib.Path(
        "/mnt/saphyr/Saphyr_Data/DAM_DLE_VHL_DLE/Michael/Chromatin_rep1/filtered_xmap/pipeline_bedgraphs/NoBlackMask/Chromatin_WholeGenome.bedgraph"
    )
)
chromatin
chrom_intersection, nfr_intersection, dnase_chrom_intersection = [], [], []
for (
    chr_,
    chrom_group,
    theo_group,
    dnase_group,
    nfr_group,
) in iter_over_bedgraphs_chromosomes(chromatin, smoothed_theo, dnase, nfr_all_chrom):
    equalized_chrom, equalized_theo = equalize_loci(
        chrom_group.copy(), theo_group.copy()
    )
    # Normalze the chromatin data
    mask, norm_by = prepare_site_density_for_norm(equalized_theo.even)
    equalized_chrom.even.loc[~mask, "intensity"] *= norm_by
    equalized_chrom.even.loc[mask, "intensity"] = 0
    # Find open areas and filter
    chrom_top = equalized_chrom.even.intensity.quantile(0.9)
    chrom_top = equalized_chrom.even.query("intensity > @chrom_top")
    #         dnase_even, dnase_at_1bp, _, _ = equalize_loci(dnase_group.copy(), theo_group.copy())
    #         nfr_even, nfr_at_1bp, _, _ = equalize_loci(nfr_group.copy(), theo_group.copy())

    try:
        nfr_chrom, cur_nfr = chrom_top.bg.weighted_overlap(nfr_group)
        dnase_chrom, _ = chrom_top.bg.weighted_overlap(dnase_group)
        chrom_intersection.append(nfr_chrom.bg.index_to_columns())
        nfr_intersection.append(cur_nfr.bg.index_to_columns())
        dnase_chrom_intersection.append(dnase_chrom.bg.index_to_columns())
    except AssertionError:
        continue

chrom_intersection = pd.concat(chrom_intersection, axis=0, ignore_index=True)
nfr_intersection = pd.concat(nfr_intersection, axis=0, ignore_index=True)
dnase_chrom_intersection = pd.concat(
    dnase_chrom_intersection, axis=0, ignore_index=True
)
nfr_per_chr = {}
dnase_nfr = {}
for chr_, chrom_group, nfr_group, dnase_group in iter_over_bedgraphs_chromosomes(
    chromatin, nfr_all_chrom, dnase
):
    nfr_int_grouped = nfr_intersection.query("chr == @chr_")
    nfr_per_chr[chr_] = len(nfr_int_grouped) / len(nfr_group)
    dnase_chrom_int_grouped = dnase_chrom_intersection.query("chr == @chr_")
    dnase_nfr[chr_] = len(dnase_chrom_int_grouped) / len(dnase_group)
nfr_per_chr_df = pd.DataFrame.from_dict(
    nfr_per_chr, orient="index", columns=["Overlap fraction"]
)
ax = nfr_per_chr_df.plot.bar()
ax.set_ylabel("Overlap fraction (ATAC vs. Chromatin)")
ax.set_title("ATAC Overlap")

dnase_chrom_df = pd.DataFrame.from_dict(
    dnase_nfr, orient="index", columns=["Overlap fraction"]
)
axx = dnase_chrom_df.plot.bar()
axx.set_ylabel("Overlap fraction (DNase vs Chromatin)")
axx.set_title("DNase overlap")
chromatin.groupby("chr").get_group("chr10").intensity.hist(bins=100)
chromatin.groupby("chr").get_group("chr8").intensity.hist(bins=100)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10557580666936682, 'fraction_numerical': 0.009180504927771027, 'mean_word_length': 5.532627865961199, 'pattern_counts': {'":': 3, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '536543', 'license': 'MIT', 'n_tokens_mistral': 2937, 'n_tokens_neox': 2779, 'n_words': 470, 'path': 'notebooks/venn.ipynb', 'repo_name': 'ebensteinLab/EcoDAM', 'size': '7407'}