dLEM on one chromosome using Python¶

Get input data¶

In [ ]:

%%bash
mkdir -p ../data
cd ../data
wget --quiet https://4dn-open-data-public.s3.amazonaws.com/fourfront-webprod/wfoutput/d13aa1ea-053c-4113-94fe-a1e7ab1dbbab/4DNFI9GMP2J8.mcool

%%bash mkdir -p ../data cd ../data wget --quiet https://4dn-open-data-public.s3.amazonaws.com/fourfront-webprod/wfoutput/d13aa1ea-053c-4113-94fe-a1e7ab1dbbab/4DNFI9GMP2J8.mcool

Train and assess dLEM model¶

In [1]:

import sys
sys.path.append('../')

import sys sys.path.append('../')

In [2]:

from dlem.api import (fetch_band, 
                        train_dlem,
                        generate_dlem_prediction,
                        fit_band_row_profile_sliding,
                        plot_map,
                        compute_contact_heatmap,
                        band_generate_pixels)

from dlem.api import (fetch_band, train_dlem, generate_dlem_prediction, fit_band_row_profile_sliding, plot_map, compute_contact_heatmap, band_generate_pixels)

In [3]:

import re
import hictkpy as htk
import numpy as np

resolution = 10000
cool_filename = f"../data/4DNFI9GMP2J8.mcool::{resolution}"
cool_handle = htk.File(cool_filename, resolution)
cur_region = "chr10:1-133797422"
split_region = re.split(r"[-:]", cur_region)
cur_chrom = split_region[0]
cur_end = cool_handle.chromosomes()[cur_chrom]
width = 700

import re import hictkpy as htk import numpy as np resolution = 10000 cool_filename = f"../data/4DNFI9GMP2J8.mcool::{resolution}" cool_handle = htk.File(cool_filename, resolution) cur_region = "chr10:1-133797422" split_region = re.split(r"[-:]", cur_region) cur_chrom = split_region[0] cur_end = cool_handle.chromosomes()[cur_chrom] width = 700

In [4]:

  
band = fetch_band(cool_filename, 
                  resolution, 
                  cur_region, 
                  width)

band = fetch_band(cool_filename, resolution, cur_region, width)

WARNING:2026-02-05 16:31:39,132:jax._src.xla_bridge:876: An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.

In [6]:

decay_fit = fit_band_row_profile_sliding(band, 
                                       start_row=5, 
                                       extent=500, 
                                       window_size=300, 
                                       window_step=200, 
                                       log_y=True)

d_vec = np.array([p['d'] for p in decay_fit['window_params']])
slowdown = -np.median(d_vec)

decay_fit = fit_band_row_profile_sliding(band, start_row=5, extent=500, window_size=300, window_step=200, log_y=True) d_vec = np.array([p['d'] for p in decay_fit['window_params']]) slowdown = -np.median(d_vec)

In [7]:

#train with MSE
resMSE = train_dlem(band,
                   steps=10, 
                   start_row=5, 
                   slowdown=slowdown,
                   train_steps=300, 
                   verbose=False, 
                   loss_type="mse")

#train with MSE resMSE = train_dlem(band, steps=10, start_row=5, slowdown=slowdown, train_steps=300, verbose=False, loss_type="mse")

Output()

In [8]:

# Generate prediction across whole chromosome
pred_this, left_this, right_this = generate_dlem_prediction(resMSE,
                                                            start=0,
                                                            span=cur_end//resolution,
                                                            mode="mse",
                                                            slowdown=slowdown,
                                                            rows=200)

# Generate prediction across whole chromosome pred_this, left_this, right_this = generate_dlem_prediction(resMSE, start=0, span=cur_end//resolution, mode="mse", slowdown=slowdown, rows=200)

Subset to a specific region and plot prediction vs input¶

In [9]:

#Subset prediction and generate normalized, matched vs input values
span = 800
start = 2050
plot_cur_start = resolution*start
plot_cur_end = resolution*(start+span)
plot_region_name = f"{cur_chrom}:{plot_cur_start}-{plot_cur_end}"
band_height = min(span, 
                  pred_this.shape[0],
                  band.shape[0])
true_vals = band[0:band_height,start:(start+span)]
pred_vals = pred_this[0:band_height,start:(start+span)]
sub_l = left_this[start:start+span]
sub_r = right_this[start:start+span]

heatmap = compute_contact_heatmap(
    pred_vals,
    lower_vals=true_vals,
    normalize=True,
    log=True,
    match_lower_range=True
)

#Subset prediction and generate normalized, matched vs input values span = 800 start = 2050 plot_cur_start = resolution*start plot_cur_end = resolution*(start+span) plot_region_name = f"{cur_chrom}:{plot_cur_start}-{plot_cur_end}" band_height = min(span, pred_this.shape[0], band.shape[0]) true_vals = band[0:band_height,start:(start+span)] pred_vals = pred_this[0:band_height,start:(start+span)] sub_l = left_this[start:start+span] sub_r = right_this[start:start+span] heatmap = compute_contact_heatmap( pred_vals, lower_vals=true_vals, normalize=True, log=True, match_lower_range=True )

In [10]:

fig = plot_map(heatmap,
                sub_l, 
                sub_r, 
                plot_region_name, 
                resolution)

fig = plot_map(heatmap, sub_l, sub_r, plot_region_name, resolution)

In [11]:

plot_width = 700
fig.update_layout(
    width=plot_width,
    height=plot_width*1.61803,
)
file_path = f"plotly_heatmap.html"
fig.write_html(file_path)

plot_width = 700 fig.update_layout( width=plot_width, height=plot_width*1.61803, ) file_path = f"plotly_heatmap.html" fig.write_html(file_path)

In [12]:

from IPython.display import HTML, display
display(HTML("<style>.container { width:100% !important; }</style>"))
display(HTML(filename='../docs/plotly_heatmap.html'))

from IPython.display import HTML, display display(HTML("")) display(HTML(filename='../docs/plotly_heatmap.html'))

Output predictions as matrix and tracks¶

In [13]:

import os
import re
import hictkpy as htk
import matplotlib.pyplot as plt
import numpy as np

#Save predicted output as cool using input as reference

out_cool_name = "test_data.cool"
ref_cool_name = "../data/4DNFI9GMP2J8.mcool"

if os.path.exists(out_cool_name):
    os.remove(out_cool_name)
cool_handle = htk.File(ref_cool_name,resolution)
bins = cool_handle.bins().to_df()[:]
resolution = cool_handle.resolution()
bins = bins.copy()
bins["weight"] = 1.0
columns_inc = ['bin1_id', 'bin2_id', 'count']
custom_dtypes = {'count': np.float32}
chrom_sizes = cool_handle.chromosomes()
chr_sub = ["chr10", "chr12", "chr13"]
pred_cool_filehandle = htk.cooler.FileWriter(out_cool_name,
                                              bins=cool_handle.bins()
                                              )
cur_pred = {chrom_name:None for chrom_name in chr_sub}
cur_pred[split_region[0]] = pred_this

for out_chrom in chr_sub:
    pred_cool_filehandle.add_pixels(sorted=False,
                                    pixels= band_generate_pixels(
                                            band=cur_pred[out_chrom],
                                            ref_cool_filename=ref_cool_name,
                                            resolution=resolution,
                                            chrom_name=out_chrom,
                                            ))
pred_cool_filehandle.finalize()

import os import re import hictkpy as htk import matplotlib.pyplot as plt import numpy as np #Save predicted output as cool using input as reference out_cool_name = "test_data.cool" ref_cool_name = "../data/4DNFI9GMP2J8.mcool" if os.path.exists(out_cool_name): os.remove(out_cool_name) cool_handle = htk.File(ref_cool_name,resolution) bins = cool_handle.bins().to_df()[:] resolution = cool_handle.resolution() bins = bins.copy() bins["weight"] = 1.0 columns_inc = ['bin1_id', 'bin2_id', 'count'] custom_dtypes = {'count': np.float32} chrom_sizes = cool_handle.chromosomes() chr_sub = ["chr10", "chr12", "chr13"] pred_cool_filehandle = htk.cooler.FileWriter(out_cool_name, bins=cool_handle.bins() ) cur_pred = {chrom_name:None for chrom_name in chr_sub} cur_pred[split_region[0]] = pred_this for out_chrom in chr_sub: pred_cool_filehandle.add_pixels(sorted=False, pixels= band_generate_pixels( band=cur_pred[out_chrom], ref_cool_filename=ref_cool_name, resolution=resolution, chrom_name=out_chrom, )) pred_cool_filehandle.finalize()

Massaging band data for chr10 into pixels

Passing through data for chr12 into pixels

Passing through data for chr13 into pixels

Out[13]:

File(test_data.cool)

Use cooler to visualize region¶

In [14]:

%%bash

cooler show \
    ./test_data.cool \
    chr10:20500000-22500000 \
    -o test_img.png

%%bash cooler show \ ./test_data.cool \ chr10:20500000-22500000 \ -o test_img.png

In [15]:

from IPython.display import Image
Image("test_img.png")

from IPython.display import Image Image("test_img.png")

Out[15]:

In [17]:

%%bash

# Zoomify dLEM cool output
hictk zoomify \
    test_data.cool \
    test_data.mcool

%%bash # Zoomify dLEM cool output hictk zoomify \ test_data.cool \ test_data.mcool

[2026-02-05 21:47:56.248] [info]: Running hictk v2.2.0-bioconda
[2026-02-05 21:47:56.278] [info]: coarsening cooler at test_data.cool 10 times (10000 -> 10000 -> 20000 -> 50000 -> 100000 -> 200000 -> 500000 -> 1000000 -> 2000000 -> 5000000 -> 10000000)
[2026-02-05 21:47:56.278] [info]: copying 10000 resolution from test_data.cool
[2026-02-05 21:47:56.726] [info]: generating 20000 resolution from 10000 (2x)
[2026-02-05 21:47:57.342] [info]: [10000 -> 20000] processing chr12:29780000-29800000 at 1655629 pixels/s...
[2026-02-05 21:47:57.898] [info]: [10000 -> 20000] processing chr13:25100000-25120000 at 1769912 pixels/s...
[2026-02-05 21:47:58.672] [info]: generating 50000 resolution from 10000 (5x)
[2026-02-05 21:47:59.544] [info]: [10000 -> 50000] processing chr13:46950000-47000000 at 1141553 pixels/s...
[2026-02-05 21:48:00.059] [info]: generating 100000 resolution from 50000 (2x)
[2026-02-05 21:48:00.544] [info]: generating 200000 resolution from 100000 (2x)
[2026-02-05 21:48:00.772] [info]: generating 500000 resolution from 100000 (5x)
[2026-02-05 21:48:00.901] [info]: generating 1000000 resolution from 500000 (2x)
[2026-02-05 21:48:00.930] [info]: generating 2000000 resolution from 1000000 (2x)
[2026-02-05 21:48:00.940] [info]: generating 5000000 resolution from 1000000 (5x)
[2026-02-05 21:48:00.960] [info]: generating 10000000 resolution from 5000000 (2x)
[2026-02-05 21:48:00.970] [info]: DONE! Processed 10 resolution(s) in 4.72s!

Use Higlass to visualize¶

In [18]:

import higlass as hg
from IPython.display import display
import os

import higlass as hg from IPython.display import display import os

In [ ]:

out_mcool_name = "test_data.mcool"
# Input data
tileset1 = hg.cooler(ref_cool_name)
# dLEM output
tileset2 = hg.cooler(out_mcool_name)

# Create a `hg.HeatmapTrack` for each tileset
track1 = tileset1.track("heatmap")
track2 = tileset2.track("heatmap")

# Create two independent `hg.View`s, one for each heatmap
view1 = hg.view(track1, width=6)
view2 = hg.view(track2, width=6)

# Lock zoom & location for each `View`
view_lock = hg.lock(view1, view2)

# Concatenate views horizontally and apply synchronization lock
(view1 | view2).locks(view_lock)

out_mcool_name = "test_data.mcool" # Input data tileset1 = hg.cooler(ref_cool_name) # dLEM output tileset2 = hg.cooler(out_mcool_name) # Create a `hg.HeatmapTrack` for each tileset track1 = tileset1.track("heatmap") track2 = tileset2.track("heatmap") # Create two independent `hg.View`s, one for each heatmap view1 = hg.view(track1, width=6) view2 = hg.view(track2, width=6) # Lock zoom & location for each `View` view_lock = hg.lock(view1, view2) # Concatenate views horizontally and apply synchronization lock (view1 | view2).locks(view_lock)

Output L and R predictions as bigwig tracks¶

In [20]:

import pyranges as pr

target_regions_gr = pr.PyRanges(
    chromosomes=[cur_chrom],
    starts=[0],
    ends=[cool_handle.chromosomes()[cur_chrom]]
)

import pyranges as pr target_regions_gr = pr.PyRanges( chromosomes=[cur_chrom], starts=[0], ends=[cool_handle.chromosomes()[cur_chrom]] )

In [21]:

cur_chr_gr = target_regions_gr[cur_chrom].as_df().iloc[0]
cur_chr_name = f"{cur_chr_gr.Chromosome}:{cur_chr_gr.Start}-{cur_chr_gr.End}"

out_gr_l = []
bw_res = 50
cur_l_gr = target_regions_gr.tile(resolution)
cur_l_gr.Score = np.append(left_this, 0)
cur_l_gr = cur_l_gr.tile(bw_res)
cur_l_gr.End = cur_l_gr.End-1

cur_chr_gr = target_regions_gr[cur_chrom].as_df().iloc[0] cur_chr_name = f"{cur_chr_gr.Chromosome}:{cur_chr_gr.Start}-{cur_chr_gr.End}" out_gr_l = [] bw_res = 50 cur_l_gr = target_regions_gr.tile(resolution) cur_l_gr.Score = np.append(left_this, 0) cur_l_gr = cur_l_gr.tile(bw_res) cur_l_gr.End = cur_l_gr.End-1

In [ ]:

cur_l_gr.to_bigwig("test.bw", 
                   dict(cool_handle.chromosomes()), 
                   rpm=False, 
                   value_col="Score")

cur_l_gr.to_bigwig("test.bw", dict(cool_handle.chromosomes()), rpm=False, value_col="Score")