Note
The relevant data used in this tutorial is accessible at:
Warning
Make sure to download the relevant index file for each .bam (.bam.bai index), .cram (.cram.crai index), or .gz (.tbi index).
-
BH01.hg19.frag.bed.gz: DOWNLOAD -
BH01.hg19.frag.bed.gz.tbi: DOWNLOAD -
ctcf.hg19.bed: DOWNLOAD -
tss.hg19.bed: DOWNLOAD
We start by reading the relevant locations of the CTCF and the TSS sites and making the -2kb and +2kb intervals of 1bp resolution.
import os
from tqdm import tqdm
def read_ctcf_create_intervals(path):
os.makedirs('./tmp', exist_ok=True)
with open(path, 'rt') as ctcf_file:
for line in tqdm(ctcf_file, desc="Processing CTCF data"):
fields = line.strip().split('\t')
if len(fields) < 6:
raise ValueError(f"Invalid input line: {line.strip()} - expected at least 6 columns.")
chrom, start, end, _, _, strand = fields
pos = (int(start) + int(end)) // 2
if strand == "+":
interval_range = range(pos - 2000, pos + 2001)
elif strand == "-":
interval_range = range(pos + 2000, pos - 2001, -1)
else:
raise ValueError(f"Invalid strand value: {strand} - expected '+' or '-'.")
intervals = [[chrom[3:], start, start + 1] for start in interval_range]
filename = f"./tmp/{chrom}_{start}_{end}.ctcf.intervals.bed"
with open(filename, 'w') as output:
output.writelines(f"{interval[0]}\t{interval[1]}\t{interval[2]}\n" for interval in intervals)
def read_tss_create_intervals(path):
output_dir = "./tmp"
os.makedirs(output_dir, exist_ok=True)
with open(path, 'rt') as tss_file:
for line in tqdm(tss_file):
# Parse the line into relevant fields
chrom, start, end, _, _, strand = line.strip().split('\t')
start, end = int(start), int(end)
intervals = []
if strand == "+":
interval_range = range(start - 2000, start + 2001)
elif strand == "-":
interval_range = range(end + 2000, end - 2001, -1)
else:
raise ValueError(
f"Invalid strand value: '{strand}'. Expected '+' or '-'."
)
intervals = [[chrom[3:], i, i + 1] for i in interval_range]
filename = os.path.join(output_dir, f"{chrom}_{start}_{end}.tss.intervals.bed")
with open(filename, 'w') as output_file:
output_file.writelines(
f"{interval[0]}\t{interval[1]}\t{interval[2]}\n" for interval in intervals
)
read_ctcf_create_intervals("ctcf.hg19.bed")
read_tss_create_intervals("tss.hg19.bed")Generate the relevant FinaleToolkit functions:
with open("ctcf.hg19.bed", 'rt') as ctcf, open('ctcf.commands.txt', 'w') as command_file:
for line in tqdm(ctcf):
chrom, start, end, _, _, strand = line.strip().split('\t')
filename = f"./tmp/{chrom}_{start}_{end}.ctcf.intervals.bed"
out_filename = f"./tmp/{chrom}_{start}_{end}.ctcf.intervals.wcov.bed"
command = [
'finaletoolkit', 'coverage',
'-o', out_filename,
'-q', '30',
'-w', '1',
'-p', 'any',
'BH01.hg19.frag.bed.gz',
filename
]
command_str = ' '.join(command)
command_file.write(command_str + '\n')
with open("tss.hg19.bed", 'rt') as tss, open('tss.commands.txt', 'w') as command_file:
for line in tqdm(tss):
chrom, start, end, _, _, strand = line.strip().split('\t')
filename = f"./tmp/{chrom}_{start}_{end}.tss.intervals.bed"
out_filename = f"./tmp/{chrom}_{start}_{end}.tss.intervals.wcov.bed"
command = [
'finaletoolkit', 'coverage',
'-o', out_filename,
'-q', '30',
'-w', '1',
'-p', 'any',
'BH01.hg19.frag.bed.gz',
filename
]
command_str = ' '.join(command)
command_file.write(command_str + '\n')Run the commands generated.
Then, compile and combine the results.
directory = ./tmp/"
pattern = "*.ctcf.intervals.wcov.bed" # REPEAT FOR *.tss.intervals.wcov.bed
file_paths = glob.glob(os.path.join(directory, pattern))
all_arrays = []
for file_path in tqdm(file_paths):
fifth_column = []
with open(file_path, 'r') as file:
for line in file:
fields = line.strip().split('\t')
if len(fields) >= 5:
fifth_column.append(float(fields[4]))
fifth_column = np.array(fifth_column)
total_sum = fifth_column.sum()
if total_sum > 0:
normalized_column = fifth_column
all_arrays.append(normalized_column)
if len(all_arrays) == 0:
raise ValueError("No valid arrays found.")
array_lengths = [len(arr) for arr in all_arrays]
if len(set(array_lengths)) > 1:
raise ValueError(f"All arrays are not the same length. Lengths: {array_lengths}")
stacked_arrays = np.stack(all_arrays)
np.save('ctcf.normalizedaverages.npy', average_array) # CHANGE FOR *.tss.intervals.wcov.bedPlot the figure:
import numpy as np
import matplotlib.pyplot as plt
def transform_array(array):
selected_columns = np.concatenate([array[:, :2000], array[:, -2000:]], axis=1)
row_means = selected_columns.mean(axis=1, keepdims=True)
non_zero_indices = np.where(row_means != 0)[0]
normalized_array = array[non_zero_indices] / row_means[non_zero_indices]
transformed_array = normalized_array.mean(axis=0)
return transformed_array[3000:-3000]
ctcf = np.load('./ctcf.normalizedaverages.npy')#.mean(axis=0)
tss = np.load('./tss.normalizedaverages.npy')#.mean(axis=0)
ctcf = transform_array(ctcf)
tss = transform_array(tss)
x = np.linspace(-2000, 2000, 4001)
fig, axes = plt.subplots(1, 2, figsize=(8,4), dpi=300)
xticks = [-2000, -1000, 0, 1000, 2000]
xtick_labels = ['-2kb', '-1kb', '0', '1kb', '2kb']
axes[0].plot(x, ctcf, color="#C19DB9", linewidth=1.5)
axes[0].axvline(0, color="black", linewidth=0.8)
axes[0].set_xlabel("Distance to CTCF motif (bp)", fontsize=12)
axes[0].set_ylabel("Normalized fragment coverage", fontsize=12)
axes[0].set_xticks(xticks)
#axes[0].set_yticks(ctcf_y_ticks)
axes[0].set_xticklabels(xtick_labels, fontsize=12)
axes[0].spines['top'].set_visible(False)
axes[0].spines['right'].set_visible(False)
axes[1].plot(x, tss, color="#5BA872", linewidth=1.5)
axes[1].axvline(0, color="black", linewidth=0.8)
axes[1].set_xlabel("Distance to CGI TSS (bp)", fontsize=12)
axes[1].set_xticks(xticks)
axes[1].set_xticklabels(xtick_labels, fontsize=12)
#axes[1].set_yticks(tss_y_ticks)
axes[1].spines['top'].set_visible(False)
axes[1].spines['right'].set_visible(False)
plt.tight_layout()
plt.show()
FinaleToolkit