Workflow Steps and Code Snippets

import sys
sys.stderr = open(snakemake.log[0], "w")

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

calls = pd.read_table(snakemake.input[0], header=[0, 1])
samples = [name for name in calls.columns.levels[0] if name != "VARIANT"]
sample_info = calls.loc[:, samples].stack([0, 1]).unstack().reset_index(1, drop=False)
sample_info = sample_info.rename({"level_1": "sample"}, axis=1)

sample_info = sample_info[sample_info["DP"] > 0]
sample_info["freq"] = sample_info["AD"] / sample_info["DP"]
sample_info.index = np.arange(sample_info.shape[0])

plt.figure()

sns.stripplot(x="sample", y="freq", data=sample_info, jitter=True)
plt.ylabel("allele frequency")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.freqs)

plt.figure()

sns.stripplot(x="sample", y="DP", data=sample_info, jitter=True)
plt.ylabel("read depth")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.depths)

import sys
sys.stderr = open(snakemake.log[0], "w")

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

calls = pd.read_table(snakemake.input[0], header=[0, 1])
samples = [name for name in calls.columns.levels[0] if name != "VARIANT"]
sample_info = calls.loc[:, samples].stack([0, 1]).unstack().reset_index(1, drop=False)
sample_info = sample_info.rename({"level_1": "sample"}, axis=1)

sample_info = sample_info[sample_info["DP"] > 0]
sample_info["freq"] = sample_info["AD"] / sample_info["DP"]
sample_info.index = np.arange(sample_info.shape[0])

plt.figure()

sns.stripplot(x="sample", y="freq", data=sample_info, jitter=True)
plt.ylabel("allele frequency")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.freqs)

plt.figure()

sns.stripplot(x="sample", y="DP", data=sample_info, jitter=True)
plt.ylabel("read depth")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.depths)

import sys
sys.stderr = open(snakemake.log[0], "w")

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

calls = pd.read_table(snakemake.input[0], header=[0, 1])
samples = [name for name in calls.columns.levels[0] if name != "VARIANT"]
sample_info = calls.loc[:, samples].stack([0, 1]).unstack().reset_index(1, drop=False)
sample_info = sample_info.rename({"level_1": "sample"}, axis=1)

sample_info = sample_info[sample_info["DP"] > 0]
sample_info["freq"] = sample_info["AD"] / sample_info["DP"]
sample_info.index = np.arange(sample_info.shape[0])

plt.figure()

sns.stripplot(x="sample", y="freq", data=sample_info, jitter=True)
plt.ylabel("allele frequency")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.freqs)

plt.figure()

sns.stripplot(x="sample", y="DP", data=sample_info, jitter=True)
plt.ylabel("read depth")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.depths)

import sys
sys.stderr = open(snakemake.log[0], "w")

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

calls = pd.read_table(snakemake.input[0], header=[0, 1])
samples = [name for name in calls.columns.levels[0] if name != "VARIANT"]
sample_info = calls.loc[:, samples].stack([0, 1]).unstack().reset_index(1, drop=False)
sample_info = sample_info.rename({"level_1": "sample"}, axis=1)

sample_info = sample_info[sample_info["DP"] > 0]
sample_info["freq"] = sample_info["AD"] / sample_info["DP"]
sample_info.index = np.arange(sample_info.shape[0])

plt.figure()

sns.stripplot(x="sample", y="freq", data=sample_info, jitter=True)
plt.ylabel("allele frequency")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.freqs)

plt.figure()

sns.stripplot(x="sample", y="DP", data=sample_info, jitter=True)
plt.ylabel("read depth")
plt.xticks(rotation="vertical")

plt.savefig(snakemake.output.depths)

import bacchus.directional as bcd
import bacchus.insulation as bci
import bacchus.io as bcio
import cooler
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.collections as mc
import numpy as np
import os
import pandas as pd
import scipy.sparse as sp
import scipy.stats as st
import seaborn as sns
from random import choices

mat_file = snakemake.input.mat
rna_file = snakemake.input.rna
annotation_file = snakemake.input.annotation
cmap = snakemake.params.cmap
dpi = snakemake.params.dpi
text_file = str(snakemake.output.text_file)
full_mat_file = str(snakemake.output.full_mat)
zoom_1_file = str(snakemake.output.zoom1)
zoom_2_file = str(snakemake.output.zoom2)
full_mat_di_file = str(snakemake.output.full_mat_di)
zoom_1_di_file = str(snakemake.output.zoom1_di)
zoom_2_di_file = str(snakemake.output.zoom2_di)
out_files_mat = [full_mat_file, zoom_1_file, zoom_2_file]
out_files_mat_di = [full_mat_di_file, zoom_1_di_file, zoom_2_di_file]
bor_trans = str(snakemake.output.bor_trans)
trans_bor = str(snakemake.output.trans_bor)

# Create outdir if necessary.
os.makedirs(str(snakemake.params.outdir), exist_ok=True)


def write_stats(text_file, text, writing_type="a"):
    with open(text_file, writing_type) as out:
        out.write(f"{text}\n")


# Compute CIDs at 5kb resolution with 100kb window size.
M5000 = cooler.Cooler(f"{mat_file}::/resolutions/{5000}").matrix(
    balance=True, sparse=True
)[:]
di5000 = bcd.directional_index(M5000, 20)
l5000 = bcd.di_borders(di5000)
write_stats(
    text_file,
    f"Numbers of CIDs at 5kb with 100kb window size: {len(l5000)}",
    "w",
)

# Compute CIDs at 2kb resolution with 50kb window size.
M2000 = cooler.Cooler(f"{mat_file}::/resolutions/{2000}").matrix(
    balance=True, sparse=True
)[:]
di2000 = bcd.directional_index(M2000, 25)
l2000 = bcd.di_borders(di2000)
write_stats(
    text_file, f"Numbers of CIDs at 2kb with 50kb window size: {len(l2000)}"
)

# Compute CIDs at 1kb resolution using insulation score.
M1000 = cooler.Cooler(f"{mat_file}::/resolutions/{1000}").matrix(
    balance=True, sparse=True
)[:]
di1000 = bcd.directional_index(M1000, 20)
l1000 = bcd.di_borders(di1000)
write_stats(
    text_file, f"Numbers of CIDs at 1kb with 50kb window size: {len(l1000)}"
)

mat = cooler.Cooler(f"{mat_file}::/resolutions/{1000}").matrix(
    balance=True, sparse=False
)[:]
mat[np.isnan(mat)] = 0
final_borders_1000, lri_1000 = bci.get_insulation_score(
    mat, [10, 15, 20, 25, 30]
)
write_stats(
    text_file,
    f"{len(final_borders_1000)} CIDs have been detected with insulation score at 1kb resolution.",
)


rna, _ = bcio.extract_big_wig(rna_file, ztransform=False)

annotation = pd.DataFrame(
    columns=["type", "start", "end", "strand", "gene_name", "rpkm"]
)

n = 0
with open(annotation_file, "r") as file:
    for line in file:
        if line.startswith("#"):
            continue
        elif line.startswith(">"):
            break
        else:
            line = line.split("\t")
            if line[2] in ["gene", "tRNA", "rRNA"]:
                if line[2] == "gene":
                    name = line[8].split("Name=")[-1].split(";")[0]
                    annot = {
                        "type": line[2],
                        "start": int(line[3]),
                        "end": int(line[4]),
                        "strand": line[6],
                        "gene_name": name,
                        "rpkm": np.mean(rna[int(line[3]) : int(line[4])]),
                    }
                    annotation = annotation.append(annot, ignore_index=True)


###
# Print matrices with CIDs borders (Fig 1e).
###
starts = [0, 2500, 3000]
ends = [len(mat), 3000, 3500]
col = ["blue", "lime", "cyan"]

for i in range(3):
    out_file, start, end = out_files_mat[i], starts[i], ends[i]

    vmax = 99
    title = None
    fig, ax = plt.subplots(1, 1, figsize=(7, 7), dpi=dpi)

    # Axis values
    scaling_factor = 1.0
    axis = "kb"

    # No end values given.
    if end == 0:
        end = len(mat)

    # Display plots
    im = ax.imshow(
        mat[start:end, start:end],
        cmap=cmap,
        vmin=0,
        vmax=np.nanpercentile(mat, vmax),
        extent=(
            start * scaling_factor,
            end * scaling_factor,
            end * scaling_factor,
            start * scaling_factor,
        ),
    )

    # Legend
    ax.set_xlabel(f"Genomic coordinates ({axis:s})", fontsize=16)
    ax.set_ylabel(f"Genomic coordinates ({axis:s})", fontsize=16)
    ax.tick_params(axis="both", labelsize=16)

    # Title
    if title is not None:
        ax.set_title(title, size=18)

    # Colorbar
    cbar = plt.colorbar(im, shrink=0.33, anchor=(0, 0.5))
    cbar.ax.tick_params(labelsize=16)

    first = True
    for i in [x * 5 for x in l5000]:
        if i in np.arange(start, end):
            #         plt.axvline(i, lw=1, ls='dashed', c='blue')
            if first:
                s = i
                first = False
                ax.axvline(
                    s,
                    ymin=(end - s) / (end - start),
                    ymax=1,
                    c=col[0],
                    lw=1,
                    zorder=0.5,
                )
                ax.axhline(
                    s,
                    xmin=0,
                    xmax=1 - ((end - s) / (end - start)),
                    c=col[0],
                    lw=1,
                    zorder=0.5,
                )
            else:
                e = i
                ax.add_patch(
                    patches.Rectangle(
                        (s, s), e - s, e - s, edgecolor=col[0], fill=False
                    )
                )
                s = i
    ax.axvline(
        s, ymin=(end - s) / (end - start), ymax=0, c=col[0], lw=1, zorder=0.5
    )
    ax.axhline(
        s,
        xmin=1,
        xmax=1 - ((end - s) / (end - start)),
        c=col[0],
        lw=1,
        zorder=0.5,
    )

    first = True
    for i in [x * 5 for x in l2000]:
        if i in np.arange(start, end):
            #         plt.axvline(i, lw=1, ls='dashed', c='cyan')
            if first:
                s = i
                first = False
                ax.axvline(
                    s,
                    ymin=(end - s) / (end - start),
                    ymax=1,
                    c=col[1],
                    lw=1,
                    zorder=0.5,
                )
                ax.axhline(
                    s,
                    xmin=0,
                    xmax=1 - ((end - s) / (end - start)),
                    c=col[1],
                    lw=1,
                    zorder=0.5,
                )
            else:
                e = i
                ax.add_patch(
                    patches.Rectangle(
                        (s, s), e - s, e - s, edgecolor=col[1], fill=False
                    )
                )
                s = i
    ax.axvline(
        s, ymin=(end - s) / (end - start), ymax=0, c=col[1], lw=1, zorder=0.5
    )
    ax.axhline(
        s,
        xmin=1,
        xmax=1 - ((end - s) / (end - start)),
        c=col[1],
        lw=1,
        zorder=0.5,
    )

    first = True
    for i in final_borders_1000:
        if i in np.arange(start, end):
            #         plt.axvline(i, lw=1, ls='dashed', c='magenta')
            if first:
                s = i
                first = False
                ax.axvline(
                    s,
                    ymin=(end - s) / (end - start),
                    ymax=1,
                    c=col[2],
                    lw=1,
                    zorder=0.5,
                )
                ax.axhline(
                    s,
                    xmin=0,
                    xmax=1 - ((end - s) / (end - start)),
                    c=col[2],
                    lw=1,
                    zorder=0.5,
                )
            else:
                e = i
                ax.add_patch(
                    patches.Rectangle(
                        (s, s), e - s, e - s, edgecolor=col[2], fill=False
                    )
                )
                s = i
    ax.axvline(
        s, ymin=(end - s) / (end - start), ymax=0, c=col[2], lw=1, zorder=0.5
    )
    ax.axhline(
        s,
        xmin=1,
        xmax=1 - ((end - s) / (end - start)),
        c=col[2],
        lw=1,
        zorder=0.5,
    )

    # Savefig
    plt.savefig(out_file, dpi=dpi, bbox_inches="tight")
    plt.close()


###
# Print full matrices with CIDs borders using DI only (Supp Fig 1f).
###

for i in range(3):
    out_file, start, end = out_files_mat_di[i], starts[i], ends[i]

    vmax = 99
    title = None
    fig, ax = plt.subplots(1, 1, figsize=(7, 7), dpi=dpi)

    # Axis values
    scaling_factor = 1.0
    axis = "kb"

    # No end values given.
    if end == 0:
        end = len(mat)

    # Display plots
    im = ax.imshow(
        mat[start:end, start:end],
        cmap=cmap,
        vmin=0,
        vmax=np.percentile(mat, vmax),
        extent=(
            start * scaling_factor,
            end * scaling_factor,
            end * scaling_factor,
            start * scaling_factor,
        ),
    )

    # Legend
    ax.set_xlabel(f"Genomic coordinates ({axis:s})", fontsize=16)
    ax.set_ylabel(f"Genomic coordinates ({axis:s})", fontsize=16)
    ax.tick_params(axis="both", labelsize=16)

    # Title
    if title is not None:
        ax.set_title(title, size=18)

    # Colorbar
    cbar = plt.colorbar(im, shrink=0.33, anchor=(0, 0.5))
    cbar.ax.tick_params(labelsize=16)

    first = True
    for i in [x * 5 for x in l5000]:
        if i in np.arange(start, end):
            #         plt.axvline(i, lw=1, ls='dashed', c='blue')
            if first:
                s = i
                first = False
                ax.axvline(
                    s,
                    ymin=(end - s) / (end - start),
                    ymax=1,
                    c=col[0],
                    lw=1,
                    zorder=0.5,
                )
                ax.axhline(
                    s,
                    xmin=0,
                    xmax=1 - ((end - s) / (end - start)),
                    c=col[0],
                    lw=1,
                    zorder=0.5,
                )
            else:
                e = i
                ax.add_patch(
                    patches.Rectangle(
                        (s, s), e - s, e - s, edgecolor=col[0], fill=False
                    )
                )
                s = i
    ax.axvline(
        s, ymin=(end - s) / (end - start), ymax=0, c=col[0], lw=1, zorder=0.5
    )
    ax.axhline(
        s,
        xmin=1,
        xmax=1 - ((end - s) / (end - start)),
        c=col[0],
        lw=1,
        zorder=0.5,
    )

    first = True
    for i in [x * 5 for x in l2000]:
        if i in np.arange(start, end):
            # plt.axvline(i, lw=1, ls='dashed', c='cyan')
            if first:
                s = i
                first = False
                ax.axvline(
                    s,
                    ymin=(end - s) / (end - start),
                    ymax=1,
                    c=col[1],
                    lw=1,
                    zorder=0.5,
                )
                ax.axhline(
                    s,
                    xmin=0,
                    xmax=1 - ((end - s) / (end - start)),
                    c=col[1],
                    lw=1,
                    zorder=0.5,
                )
            else:
                e = i
                ax.add_patch(
                    patches.Rectangle(
                        (s, s), e - s, e - s, edgecolor=col[1], fill=False
                    )
                )
                s = i
    ax.axvline(
        s, ymin=(end - s) / (end - start), ymax=0, c=col[1], lw=1, zorder=0.5
    )
    ax.axhline(
        s,
        xmin=1,
        xmax=1 - ((end - s) / (end - start)),
        c=col[1],
        lw=1,
        zorder=0.5,
    )

    first = True
    for i in [x * 5 for x in l1000]:
        if i in np.arange(start, end):
            # plt.axvline(i, lw=1, ls='dashed', c='cyan')
            if first:
                s = i
                first = False
                ax.axvline(
                    s,
                    ymin=(end - s) / (end - start),
                    ymax=1,
                    c=col[2],
                    lw=1,
                    zorder=0.5,
                )
                ax.axhline(
                    s,
                    xmin=0,
                    xmax=1 - ((end - s) / (end - start)),
                    c=col[2],
                    lw=1,
                    zorder=0.5,
                )
            else:
                e = i
                ax.add_patch(
                    patches.Rectangle(
                        (s, s), e - s, e - s, edgecolor=col[2], fill=False
                    )
                )
                s = i
    ax.axvline(
        s, ymin=(end - s) / (end - start), ymax=0, c=col[2], lw=1, zorder=0.5
    )
    ax.axhline(
        s,
        xmin=1,
        xmax=1 - ((end - s) / (end - start)),
        c=col[2],
        lw=1,
        zorder=0.5,
    )

    # Savefig
    plt.savefig(out_file, dpi=dpi, bbox_inches="tight")
    plt.close()

###
# Transcription at borders.
###


def get_borders_transcriptions(annotation, borders, size, step):
    borders_transcriptions = []
    for i in annotation.index:
        if annotation.loc[i, "strand"] == "-":
            pos = annotation.loc[i, "end"]
        else:
            pos = annotation.loc[i, "start"]
        for j in borders:
            a = j * size + step
            if abs(pos - a) < 2500:
                borders_transcriptions.append(annotation.loc[i, "rpkm"])
    return borders_transcriptions


borders_trans_5000 = get_borders_transcriptions(annotation, l5000, 5000, 2500)
borders_trans_2000 = get_borders_transcriptions(annotation, l2000, 5000, 1000)
borders_trans_1000 = get_borders_transcriptions(
    annotation, final_borders_1000, 1000, 500
)
data = {
    "RPKM (log)": np.log(
        list(annotation.rpkm)
        + list(np.log(borders_trans_5000))
        + list(np.log(borders_trans_2000))
        + list(np.log(borders_trans_1000))
    ),
    "Genes": list(np.repeat(f"All genes\nn={len(annotation)}", len(annotation)))
    + list(
        np.repeat(f"5kb\nn={len(borders_trans_5000)}", len(borders_trans_5000))
    )
    + list(
        np.repeat(f"2kb\nn={len(borders_trans_2000)}", len(borders_trans_2000))
    )
    + list(
        np.repeat(f"1kb\nn={len(borders_trans_1000)}", len(borders_trans_1000))
    ),
}
data = pd.DataFrame(data)
data.replace([np.inf, -np.inf, np.nan], 0, inplace=True)

sns.violinplot(x="Genes", y="RPKM (log)", data=data, palette="tab10")
plt.xlabel("Resolution of the matrix to detect the borders", size=14)
plt.ylabel("RPKM (log)", size=14)
plt.title(
    "Genes transcription of genes at less\nthan 5kb of a detected border",
    size=16,
)
plt.xticks(size=12)
plt.yticks(size=12)
plt.savefig(bor_trans, bbox_inches="tight")
plt.close()

write_stats(
    text_file,
    f"p-value between 5kb borders and whole genome: {st.mannwhitneyu(np.log(list(annotation.rpkm)), np.log(borders_trans_5000))[1]}",
)
write_stats(
    text_file,
    f"p-value between 2kb borders and whole genome: {st.mannwhitneyu(np.log(list(annotation.rpkm)), np.log(borders_trans_2000))[1]}",
)
write_stats(
    text_file,
    f"p-value between 1kb borders and whole genome: {st.mannwhitneyu(np.log(list(annotation.rpkm)), np.log(borders_trans_1000))[1]}",
)

##
# Transcription depending on distance borders
###

annotation["dist_5kb"] = 0
annotation["dist_2kb"] = 0
annotation["dist_1kb"] = 0
annotation["log1p_rpkm"] = np.log1p(annotation["rpkm"])
for i in annotation.index:
    value = np.inf
    if annotation.loc[i, "strand"] == "+":
        for x in l5000:
            pos = annotation.loc[i, "start"]
            if pos < x * 5000:
                a = x * 5000 - pos
            elif pos > (x + 1) * 5000:
                a = pos - (x + 1) * 5000
            else:
                a = -1
            value = min(a, value)
    else:
        for x in l5000:
            pos = annotation.loc[i, "end"]
            if pos < x * 5000:
                a = x * 5000 - pos
            elif pos > (x + 1) * 5000:
                a = pos - (x + 1) * 5000
            else:
                a = -1
            value = min(a, value)
    annotation.loc[i, "dist_5kb"] = value

for i in annotation.index:
    value = np.inf
    if annotation.loc[i, "strand"] == "+":
        for x in l2000:
            pos = annotation.loc[i, "start"]
            if pos < x * 5000:
                a = x * 5000 - pos
            elif pos > (x + 0.4) * 5000:
                a = pos - (x + 0.4) * 5000
            else:
                a = -1
            value = min(a, value)
    else:
        for x in l2000:
            pos = annotation.loc[i, "end"]
            if pos < x * 5000:
                a = x * 5000 - pos
            elif pos > (x + 0.4) * 5000:
                a = pos - (x + 0.4) * 5000
            else:
                a = -1
            value = min(a, value)
    annotation.loc[i, "dist_2kb"] = value

for i in annotation.index:
    value = np.inf
    if annotation.loc[i, "strand"] == "+":
        for x in final_borders_1000:
            pos = annotation.loc[i, "start"]
            if pos < x * 1000:
                a = x * 1000 - pos
            elif pos > (x + 1) * 1000:
                a = pos - (x + 1) * 1000
            else:
                a = -1
            value = min(a, value)
    else:
        for x in final_borders_1000:
            pos = annotation.loc[i, "end"]
            if pos < x * 1000:
                a = x * 1000 - pos
            elif pos > (x + 1) * 1000:
                a = pos - (x + 1) * 1000
            else:
                a = -1
            value = min(a, value)
    annotation.loc[i, "dist_1kb"] = value

bin_means_l3, bin_edges_l3, binnumber_l3 = st.binned_statistic(
    annotation.dist_1kb,
    annotation.rpkm,
    statistic="mean",
    bins=41,
    range=(-5000, 200000),
)
bin_means_l2, bin_edges_l2, binnumber_l2 = st.binned_statistic(
    annotation.dist_2kb,
    annotation.rpkm,
    statistic="mean",
    bins=41,
    range=(-5000, 200000),
)
bin_means_l1, bin_edges_l1, binnumber_l1 = st.binned_statistic(
    annotation.dist_5kb,
    annotation.rpkm,
    statistic="mean",
    bins=41,
    range=(-5000, 200000),
)


def bootstrap_down(data):
    l = np.zeros((10000))
    for i in range(10000):
        l[i] = np.mean(choices(data, k=len(data)))
    return np.percentile(l, 5)


def bootstrap_up(data):
    l = np.zeros((10000))
    for i in range(10000):
        l[i] = np.mean(choices(data, k=len(data)))
    return np.percentile(l, 95)


bin_perc5_l1, bin_edges_l1, binnumber_l1 = st.binned_statistic(
    annotation.dist_5kb,
    annotation.rpkm,
    statistic=bootstrap_down,
    bins=41,
    range=(-5000, 200000),
)
bin_perc95_l1, bin_edges_l1, binnumber_l1 = st.binned_statistic(
    annotation.dist_5kb,
    annotation.rpkm,
    statistic=bootstrap_up,
    bins=41,
    range=(-5000, 200000),
)
bin_perc5_l2, bin_edges_l2, binnumber_l2 = st.binned_statistic(
    annotation.dist_2kb,
    annotation.rpkm,
    statistic=bootstrap_down,
    bins=41,
    range=(-5000, 200000),
)
bin_perc95_l2, bin_edges_l2, binnumber_l2 = st.binned_statistic(
    annotation.dist_2kb,
    annotation.rpkm,
    statistic=bootstrap_up,
    bins=41,
    range=(-5000, 200000),
)
bin_perc5_l3, bin_edges_l3, binnumber_l3 = st.binned_statistic(
    annotation.dist_1kb,
    annotation.rpkm,
    statistic=bootstrap_down,
    bins=41,
    range=(-5000, 200000),
)
bin_perc95_l3, bin_edges_l3, binnumber_l3 = st.binned_statistic(
    annotation.dist_1kb,
    annotation.rpkm,
    statistic=bootstrap_up,
    bins=41,
    range=(-5000, 200000),
)

plt.plot(bin_edges_l1[1:] / 1000, bin_means_l1, label="5kb")
plt.fill_between(
    bin_edges_l1[1:] / 1000, bin_perc5_l1, bin_perc95_l1, alpha=0.5
)
plt.plot(bin_edges_l2[1:] / 1000, bin_means_l2, label="2kb")
plt.fill_between(
    bin_edges_l2[1:] / 1000, bin_perc5_l2, bin_perc95_l2, alpha=0.5
)
plt.plot(bin_edges_l3[1:] / 1000, bin_means_l3, label="1kb")
plt.fill_between(
    bin_edges_l3[1:] / 1000, bin_perc5_l3, bin_perc95_l3, alpha=0.5
)
plt.xlabel("Distance from the closest border (kb)", size=14)
plt.ylabel("Mean RPKM", size=14)
plt.xlim(0, 50)
plt.title(
    "Genes transcription depending on\nthe distance to the closest border",
    size=16,
)
plt.xticks(size=12)
plt.yticks(size=12)
plt.legend()
plt.savefig(trans_bor, bbox_inches="tight")
plt.close()

import bacchus.hic as bch
import bacchus.io as bcio
import bacchus.plot as bcp
import cooler
import copy
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import os
from os.path import join
import pandas as pd
from scipy import ndimage
import seaborn as sns

# Import snakemake parameters.
annotation_file = snakemake.input.annotation
mat_wt_file = str(snakemake.input.cool_wt)
mat_rf_file = str(snakemake.input.cool_rf)
rna_wt_file = snakemake.input.rna_wt
rna_rf_file = snakemake.input.rna_rf
cov_wt_file = snakemake.input.cov_wt
cov_rf_file = snakemake.input.cov_rf
gc_file = snakemake.input.gc
epod_file = snakemake.input.EPODs
res = int(snakemake.params.res)
cmap = snakemake.params.cmap
outdir = str(snakemake.params.outdir)
width = snakemake.params.width
pos = snakemake.params.positions

# Create outdir if necessary.
os.makedirs(outdir, exist_ok=True)


def import_annotation_gff(annotation_file):
    """Function to create a table of the gene positions from the gff file."""
    annotation = pd.DataFrame(
        columns=["type", "start", "end", "strand", "name"]
    )
    with open(annotation_file, "r") as file:
        for line in file:
            # Header.
            if line.startswith("#"):
                continue
            # Stop at the fasta sequences.
            elif line.startswith(">"):
                break
            else:
                line = line.split("\t")
                if line[2] in ["gene", "tRNA", "rRNA"]:
                    if line[2] == "gene":
                        name = line[8].split("Name=")[-1].split(";")[0]
                        # Extract gene position.
                        annot = {
                            "type": line[2],
                            "start": int(line[3]),
                            "end": int(line[4]),
                            "strand": line[6],
                            "gene_name": name,
                        }
                        annotation = annotation.append(annot, ignore_index=True)
    return annotation


# Import files.
annotation = import_annotation_gff(annotation_file)
mat_wt = cooler.Cooler(f"{mat_wt_file}::/resolutions/{res}").matrix(
    balance=True, sparse=False
)[:]
mat_wt[np.isnan(mat_wt)] = 0
mat_rf = cooler.Cooler(f"{mat_rf_file}::/resolutions/{res}").matrix(
    balance=True, sparse=False
)[:]
mat_rf[np.isnan(mat_rf)] = 0
rna_wt, _ = bcio.extract_big_wig(rna_wt_file, binning=100)
rna_rf, _ = bcio.extract_big_wig(rna_rf_file, binning=100)
cov_wt, _ = bcio.extract_big_wig(cov_wt_file, binning=1000)
cov_rf, _ = bcio.extract_big_wig(cov_wt_file, binning=1000)
gc_content = pd.read_csv(gc_file, sep="\t", header=None).iloc[:, 1]
epod = pd.read_csv(epod_file, sep="\t", header=None)

# Make the rotation of the matrix to plot them
mat_wt_rot = copy.copy(mat_wt)
mat_wt_rot = bch.interpolate_white_lines(mat_wt_rot)
mat_wt_rot[np.isnan(mat_wt_rot)] = 0
mat_wt_rot = ndimage.rotate(mat_wt_rot, 45, reshape=True)
mat_rf_rot = copy.copy(mat_rf)
mat_rf_rot = bch.interpolate_white_lines(mat_rf_rot)
mat_rf_rot[np.isnan(mat_rf_rot)] = 0
mat_rf_rot = ndimage.rotate(mat_rf_rot, 45, reshape=True)


def plot_region(
    M1,
    M1_rot,
    M2,
    M2_rot,
    rna1,
    rna2,
    cov1,
    cov2,
    annotation,
    binning,
    width,
    zoom_ini,
    outfile,
    split=False,
):
    rna_binning = 100
    pal = sns.color_palette("Paired")

    # Defined values to specify the borders of the matrices and tables
    zoom = [zoom // binning for zoom in zoom_ini]
    width = width / 1000
    row1 = len(M1_rot) // 2 - int(np.sqrt(2) * width)
    row2 = len(M1_rot) // 2 + int(np.sqrt(2) * width)
    col1 = int(zoom[0] * np.sqrt(2))
    col2 = int(zoom[1] * np.sqrt(2))
    annotation_zoom = annotation[
        np.logical_and(
            annotation.start > zoom_ini[0], annotation.end < zoom_ini[1]
        )
    ]
    ymax = np.nanmax(
        np.concatenate(
            (
                rna1[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
                rna2[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
            )
        )
    )
    max_cov = np.nanmax(cov1) * 1.02
    min_cov = np.nanmin(cov1) * 0.98

    # Define panels depending on split or not.
    if split and ymax > 1500:
        a = 4
        fig, ax = plt.subplots(
            7,
            2,
            figsize=(16, 12),
            gridspec_kw={"height_ratios": [1, 30, 1, 6, 12, 3, 3]},
            sharex=False,
        )

        # Parameters to merge panel 2 and 3.
        d = 0.02
        D = 0.04
        for j in range(2):
            ax[a - 1, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
            ax[a - 1, j].tick_params(axis="both", labelsize=15)
            ax[a - 1, j].spines["bottom"].set_visible(False)
            ax[a - 1, j].spines["top"].set_visible(False)
            ax[a - 1, j].spines["right"].set_visible(False)
            ax[a - 1, j].set_ylim(ymax - 525, ymax)
            ax[a - 1, j].get_xaxis().set_visible(False)

            # Add the small cut between the two panels
            kwargs = dict(
                transform=ax[a - 1, j].transAxes, color="k", clip_on=False
            )
            ax[a - 1, j].plot((-d, +d), (-D, +D), **kwargs)  # top-left diagonal
            kwargs.update(
                transform=ax[a, j].transAxes
            )  # switch to the bottom axes
            ax[a, j].plot(
                (-d, +d), (1 - d, 1 + d), **kwargs
            )  # bottom-left diagonal

        # Plot the RNA on the split panel
        ax[a - 1, 0].fill_between(
            np.arange(
                zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000
            ),
            rna1[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
            color="black",
        )
        ax[a - 1, 1].fill_between(
            np.arange(
                zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000
            ),
            rna2[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
            color="black",
        )

    else:
        a = 3
        fig, ax = plt.subplots(
            6,
            2,
            figsize=(16, 12),
            gridspec_kw={"height_ratios": [1, 30, 1, 15, 3, 3]},
            sharex=False,
        )

    # Plot expressed genes - blue are forward - red are reversed
    for j in range(2):
        ax[0, j].set_xlim(zoom_ini[0], zoom_ini[1])
        ax[0, j].get_xaxis().set_visible(False)
        ax[0, j].get_yaxis().set_visible(False)
        ax[2, j].set_xlim(zoom_ini[0], zoom_ini[1])
        ax[2, j].get_xaxis().set_visible(False)
        ax[2, j].get_yaxis().set_visible(False)
    pos = 1
    for i in range(len(annotation_zoom)):
        # Extract annotation information
        annot = annotation_zoom.iloc[
            i,
        ]
        strand = annot.strand
        start = annot.start // rna_binning
        end = annot.end // rna_binning
        name = annot.gene_name
        # Defined color depending on the strand
        if strand == "+":
            color = pal.as_hex()[1]
        else:
            color = pal.as_hex()[5]
        # Print it only if it transcribed (10% most transcribed genes)
        if np.mean(rna1[start:end]) >= 120.7628998139441:
            pos = pos * -1
            for j in range(2):
                ax[0, j].add_patch(
                    patches.Rectangle(
                        (start * 100, 0),
                        (end - start) * 100,
                        1,
                        edgecolor=color,
                        facecolor=color,
                        fill=True,
                    )
                )
                ax[0, j].text(
                    x=(start + ((end - start) / 2)) * 100,
                    y=pos + 0.5,
                    s=name,
                    rotation=90 * pos,
                    wrap=True,
                )
    color = pal.as_hex()[3]
    for i in epod.index:
        start = epod.loc[i, 3]
        end = epod.loc[i, 4]
        if (
            (start > zoom_ini[0])
            and (start < zoom_ini[1])
            or (end > zoom_ini[0])
            and (end < zoom_ini[1])
        ):
            start = max(start, zoom_ini[0]) // rna_binning
            end = min(end, zoom_ini[1]) // rna_binning
            for j in range(2):
                ax[2, j].add_patch(
                    patches.Rectangle(
                        (start * 100, 0),
                        (end - start) * 100,
                        1,
                        edgecolor=color,
                        facecolor=color,
                        fill=True,
                    )
                )
    #                 ax[2, j].text(
    #                     x=(start + ((end - start) / 2)) * 100,
    #                     y=pos + 0.5,
    #                     s="EPOD",
    #                     rotation=90 * pos,
    #                     wrap=True,
    #                 )

    # Plot the matrices
    ax[1, 0].set_ylabel("Genomic distance (kb)", fontsize=15)
    for j in range(2):
        ax[1, j].get_xaxis().set_visible(False)
        ax[1, j].tick_params(axis="both", labelsize=15)

    # Plot the matrice 1
    im = ax[1, 0].imshow(
        M1_rot[row1:row2, col1:col2],
        cmap=cmap,
        interpolation="none",
        vmin=0,
        vmax=np.nanpercentile(M1[zoom[0] : zoom[1], zoom[0] : zoom[1]], 97),
        extent=(col1 / np.sqrt(2), col2 / np.sqrt(2), width, -width),
    )

    # Plot the matrice 2
    ax[1, 1].imshow(
        M2_rot[row1:row2, col1:col2],
        cmap=cmap,
        interpolation="none",
        vmin=0,
        vmax=np.nanpercentile(M1[zoom[0] : zoom[1], zoom[0] : zoom[1]], 97),
        extent=(col1 / np.sqrt(2), col2 / np.sqrt(2), width, -width),
    )

    # RNAseq legends and plot at the bottom panel.
    ax[a, 0].set_ylabel("RNA count (CPM)", fontsize=15)
    for j in range(2):
        #         ax[a, j].set_xlabel("Genomic coordinates (kb)", size=15)
        ax[a, j].get_xaxis().set_visible(False)
        ax[a, j].tick_params(axis="both", labelsize=15)
        ax[a, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
        ax[a, j].spines["top"].set_visible(False)
        ax[a, j].spines["right"].set_visible(False)
        if split:
            if ymax > 1500:
                ax[a, j].set_ylim(0, 1050)
            else:
                ax[a, j].set_ylim(0, 1500)
        else:
            ax[a, j].set_ylim(0, ymax)

    # Plot RNA at the bottom panel
    ax[a, 0].fill_between(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000),
        rna1[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
        color="black",
    )
    ax[a, 1].fill_between(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000),
        rna2[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
        color="black",
    )

    # GC content
    gc_binning = 100
    for j in range(2):
        ax[a + 1, j].plot(
            np.arange(
                zoom_ini[0] // 1000, zoom_ini[1] // 1000, gc_binning / 1000
            ),
            gc_content[
                (zoom_ini[0] + 250)
                // gc_binning : (zoom_ini[1] + 250)
                // gc_binning
            ],
            color="black",
        )
        ax[a + 1, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
        ax[a + 1, j].tick_params(axis="both", labelsize=15)
        ax[a + 1, j].get_xaxis().set_visible(False)
    ax[a + 1, 0].set_ylabel("GC", fontsize=15)

    # Coverage
    cov_binning = 1000
    ax[a + 2, 0].plot(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, cov_binning / 1000),
        cov1[
            (zoom_ini[0] + 250)
            // cov_binning : (zoom_ini[1] + 250)
            // cov_binning
        ],
        color="black",
    )
    ax[a + 2, 1].plot(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, cov_binning / 1000),
        cov2[
            (zoom_ini[0] + 250)
            // cov_binning : (zoom_ini[1] + 250)
            // cov_binning
        ],
        color="black",
    )
    for j in range(2):
        ax[a + 2, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
        ax[a + 2, j].set_ylim(min_cov, max_cov)
        ax[a + 2, j].set_xlabel("Genomic coordinates (kb)", size=15)
        ax[a + 2, j].tick_params(axis="both", labelsize=15)
    ax[a + 2, 0].set_ylabel("HiC\ncoverage\n(CPM)", fontsize=15)

    # Colorbar
    cbar = fig.colorbar(
        im, ax=ax.ravel().tolist(), shrink=0.3, anchor=(1.2, 0.75)
    )
    cbar.ax.tick_params(labelsize=15)

    # Save the fig and adjust margins
    if split and ymax > 1500:
        plt.subplots_adjust(wspace=0.2, hspace=0.1)
    else:
        plt.subplots_adjust(wspace=0.2, hspace=0.1)
    plt.savefig(outfile, bbox_inches="tight", dpi=200)


for split in [True, False]:
    if split:
        outfile = join(
            outdir,
            f"region_{pos[0]}_{pos[1]}_split.pdf",
        )
    else:
        outfile = join(
            outdir,
            f"region_{pos[0]}_{pos[1]}.pdf",
        )
    position = [int(p) * 1000 for p in pos]
    plot_region(
        mat_wt,
        mat_wt_rot,
        mat_rf,
        mat_rf_rot,
        rna_wt,
        rna_rf,
        cov_wt,
        cov_rf,
        annotation,
        res,
        width,
        position,
        outfile,
        split=split,
    )