simulate_data.py

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df = pd.read_csv("annotations.csv", header=7)

df_sorted = df[df["CHR"] == 1].sort_values(by="MAPINFO")

plt.figure(figsize=(30, 3))

LENGTH = 10000

regions = []

dmrs = []
non_DMRS = []
outside_tss_site = []

distances = np.empty(LENGTH - 1)
true_states = np.empty(LENGTH, dtype=int)

last = None

biggest_loc = None

prev_loc = None

for idx, loc in enumerate(df_sorted["MAPINFO"][:LENGTH]):
    if prev_loc is not None:
        distances[idx - 1] = loc - prev_loc
    prev_loc = loc

np.save("distances.npy", distances)

for idx, (loc, highlight) in enumerate(
    zip(
        df_sorted["MAPINFO"][:LENGTH],
        df_sorted["UCSC_RefGene_Group"][:LENGTH] == "TSS1500",
    )
):
    if highlight and (last is None or last + 1500 < loc):
        # plt.axvspan(loc, loc+1500, color='red', alpha=0.7)
        last = loc
        regions.append([])

    if highlight:
        regions[-1].append(idx)

    else:
        outside_tss_site.append(idx)

    if biggest_loc is None or loc > biggest_loc:
        biggest_loc = loc

    # plt.axvline(loc, color="green", linewidth = 0.1)


for i in regions:
    if np.random.random() > 0.5:
        non_DMRS.append(i)
    else:
        dmrs.append(i)

colors = np.full((LENGTH), None)

sample_1_m = np.empty(LENGTH)
sample_2_m = np.empty(LENGTH)


sample1_data = np.empty((2, LENGTH))


def beta2M(b):
    delta = 0.01
    b = (b + delta) / (1 + 2 * delta)
    return np.log2(b / (1 - b))


def get_outside_tss_site():
    # if going to be a methylated region
    if np.random.random() > 0.5:
        # Numbers from peters et al supp section
        return beta2M(np.random.beta(2.4, 20, 2)), "tab:purple"
    else:
        return beta2M(np.random.beta(14, 3, 2)), "tab:blue"


for idx in outside_tss_site:
    idx = int(idx)
    (m1, m2), colors[idx] = get_outside_tss_site()

    true_states[idx] = 1 if colors[idx] == "tab:purple" else 2
    sample_1_m[idx] = m1
    sample_2_m[idx] = m2

for r in non_DMRS:
    for idx in r:
        idx = int(idx)

        (m1, m2), colors[idx] = get_outside_tss_site()

        true_states[idx] = 1 if colors[idx] == "tab:purple" else 2
        sample_1_m[idx] = m1
        sample_2_m[idx] = m2


def draw_beta(mode):
    K = 100

    mu = mode
    r = mu / (1 - mu)
    B = K / (1 + r)
    A = r * B
    a = A + 1
    b = B + 1
    return np.random.beta(a, b)


for r in dmrs:
    # Hyper half the time, hypo the other half
    hyper = np.random.random() > 0.5

    m1 = np.random.uniform(0, 0.8)
    m2 = m1 + 0.2

    if hyper:
        temp = m1
        m1 = m2
        m2 = temp

    for idx in r:
        idx = int(idx)

        sample_1_m[idx] = beta2M(draw_beta(m1))
        sample_2_m[idx] = beta2M(draw_beta(m2))

        colors[idx] = "tab:red" if hyper else "tab:green"
        true_states[idx] = 3 if hyper else 4

sample1_data = np.empty((2, len(sample_1_m)))

for idx in range(len(sample_1_m)):
    sample1_data[0, idx] = sample_1_m[idx]
    sample1_data[1, idx] = sample_2_m[idx]
np.save("simulated_data.npy", sample1_data)

np.save("true_states.npy", true_states)

plt.figure()
plt.scatter(
    sample_1_m,
    sample_2_m,
    color=colors,
    s=0.6,
)
plt.xlabel("Sample 1 M Value")
plt.ylabel("Sample 2 M Value")
plt.title("Distribution of Simulated Paired Samples")

plt.xticks((-6, -3, 0, 3, 6))
plt.yticks((-6, -3, 0, 3, 6))

import matplotlib.patches as mpatches

handles = []
patch = mpatches.Patch(color="tab:purple", label="Low methylation")
patch2 = mpatches.Patch(color="tab:blue", label="High methylation")
patch3 = mpatches.Patch(color="tab:green", label="DMR (hypermethylation)")
patch4 = mpatches.Patch(color="tab:red", label="DMR (hypomethylation)")

# handles is a list, so append manual patch
handles = [patch, patch2, patch3, patch4]

# plot the legend
plt.legend(handles=handles, loc="upper left")

plt.tight_layout()

plt.savefig("sim_distr.pdf")

plt.figure()
plt.hist(np.concatenate((sample_1_m, sample_2_m)), bins=50)
plt.title("Distribution of Simulated Methylation Amounts")
plt.xlabel("M Value")
plt.ylabel("Count")

plt.savefig("distr_methy_sim.pdf")

plt.figure()

real_world = np.genfromtxt("sample.txt", skip_header=1, delimiter=",")
plt.hist(beta2M(real_world[:, range(1, 1501, 2)].flatten()), bins=50, color="tab:green")
plt.title("Distribution of Biological Methylation Amounts")
plt.xlabel("M Value")
plt.ylabel("Count")
plt.savefig("distr_methly_real.pdf")

plt.figure()

plt.xscale("log")
logbins = np.geomspace(distances.min(), distances.max(), 40)
plt.hist(distances, bins=logbins)
plt.title("Distribution of Probe Distances")
plt.xlabel("Probe Distance (in base pairs)")
plt.ylabel("Count")
plt.axvline(10000, linestyle="dashed", color="grey")
plt.savefig("probe_dist_distr.pdf")