OSCR

CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.

Code ↔ Paper

40 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 40 matches · 5 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Enhancer models › Motif analysis › Human PBMC motif analysis ↔ Figure_2/DeepBICCN2_analysis.ipynb, lines 582–594 · score 0.89 · crested.tl.modisco.create_pattern_matrix, seqlet_count_log, discard_ic_threshold, pattern_parameter, sim_threshold, crested.tl.modisco.process_patterns
  2. [2] § Methods › Enhancer models › Motif analysis › Mouse cortex motif analysis ↔ Figure_2/DeepBICCN2_analysis.ipynb, lines 582–594 · score 0.89 · crested.tl.modisco.create_pattern_matrix, seqlet_count_log, discard_ic_threshold, pattern_parameter, sim_threshold, crested.tl.modisco.process_patterns
  3. [3] § Methods › Enhancer models › Motif analysis › Mouse cortex motif analysis ↔ Figure_3/DeepPBMC_analysis.ipynb, lines 498–513 · score 0.89 · crested.tl.modisco.create_pattern_matrix, seqlet_count_log, discard_ic_threshold, pattern_parameter, sim_threshold, crested.tl.modisco.process_patterns
  4. [4] § Methods › Enhancer models › Motif analysis › Human PBMC motif analysis ↔ Figure_3/DeepPBMC_analysis.ipynb, lines 498–513 · score 0.89 · crested.tl.modisco.create_pattern_matrix, seqlet_count_log, discard_ic_threshold, pattern_parameter, sim_threshold, crested.tl.modisco.process_patterns
  5. [5] § Methods › Enhancer models › Motif analysis › Human PBMC motif enrichment analysis using pycisTarget ↔ create_fasta_with_padded_bg_from_bed.sh, lines 1–78 · score 0.86 · padded bg, create fasta, create_cistarget_motif_databases.py, bed, background, arguments
  6. [6] § A CREsted human PBMC model captures validated TFBSs ↔ Figure_3/chip_seq_analysis.ipynb, lines 1204–1343 · score 0.81 · UniBind, ChIP seq, GATA3, RUNX1, CEBPA, SPI1
  7. [7] § Methods › Enhancer models › Motif analysis › Zebrafish motif analysis ↔ src/crested/tl/modisco/_tfmodisco.py, lines 1048–1179 · score 0.79 · crested tl modisco, discard_ic_threshold, sim_threshold, trim_ic_threshold, lite, tfmodisco
  8. [8] § Methods › Enhancer models › Motif analysis › Human PBMC UniBind comparison ↔ Figure_3/chip_seq_analysis.ipynb, lines 1204–1343 · score 0.78 · recursive_seqlets, ChIP seq, Unibind, GATA3, RUNX1, CEBPA
  9. [9] § Methods › CREsted workflow › Model training › Regression models ↔ src/crested/tl/losses/_poissonmultinomial.py, lines 8–124 · score 0.77 · Poisson loss, multinomial loss, loss function, log transformed, magnitude, truth
  10. [10] § Methods › Enhancer models › Model training › CREsted cancer cell line peak regression model, DeepCCL ↔ Figure_4/evaluate_deepccl_chrombpnet.ipynb, lines 427–447 · score 0.76 · HepG2, M059J, reverse complement, GM12878, MM001, MM099
  11. [11] § CREsted identifies high similarity between MES enhancer codes in cancer ↔ Figure_4/topic_scoring.ipynb, lines 316–376 · score 0.74 · melanoma MES, GBM MES, Spearman correlation, biopsy topic, DeepCCL, NPC
  12. [12] § Methods › Enhancer models › Model evaluation › CREsted predictions against topics derived from the cisTopic analysis of human Gliomas ↔ Figure_4/topic_scoring.ipynb, lines 316–376 · score 0.74 · melanoma MES, Spearman correlation, biopsy topics, DeepCCL, NPC, OPC
  13. [13] § CREsted provides detailed insights into enhancer codes of mouse cortical cell types ↔ Figure_5/compare_models.ipynb, lines 762–806 · score 0.73 · Micro PVMs, leptomeningeal cell, GABA, VLMC, glut, vascular
  14. [14] § Methods › Enhancer models › Motif analysis › Human PBMC ChIP–seq comparison ↔ Figure_3/chip_seq_analysis.ipynb, lines 58–121 · score 0.69 · POU2F2, ChIP seq, GATA3, RUNX1, CEBPA, SPI1
  15. [15] § CREsted identifies high similarity between MES enhancer codes in cancer ↔ Figure_4/chrombpnet_porting.ipynb, lines 26–103 · score 0.67 · ChromBPnet, m059j, ensemble, HepG2, GM12878, MM001
  16. [16] § Methods › Enhancer models › Model training › CREsted Borzoi transfer learning ↔ src/crested/tl/zoo/_dilated_cnn_decoupled.py, the whole file · a weak match · score 0.66 · Conv1D, Dense layer, CNN, cropping, Softplus, head
  17. [17] § Methods › Human cancer cell lines › Data processing ↔ Figure_4/topic_scoring.ipynb, lines 258–314 · score 0.66 · HepG2, DeepCCL, MEL, GM12878, MM001, MM099
  18. [18] § Methods › Enhancer models › Motif analysis › Human PBMC motif enrichment analysis using pycisTarget ↔ create_cistarget_motif_databases.py, lines 154–217 · score 0.66 · cisTarget motif database, create cistarget motif, background padding, bg, command, fasta
  19. [19] § CREsted identifies high similarity between MES enhancer codes in cancer ↔ Figure_4/enhancer_code_deepccl.ipynb, lines 127–272 · score 0.65 · melanoma cell, RUNX, ATF, CREB, TEAD, MM029
  20. [20] § Methods › Enhancer models › Model training › CREsted Borzoi transfer learning ↔ src/crested/tl/zoo/_dilated_cnn.py, the whole file · a weak match · score 0.64 · Conv1D, Dense layer, CNN, cropping, Softplus, convolutional
  21. [21] § Methods › Enhancer models › Model training › CREsted cancer cell line peak regression model, DeepCCL ↔ Figure_4/topic_scoring.ipynb, lines 58–91 · score 0.64 · HepG2, M059J, GM12878, MM001, MM099, MM029
  22. [22] § CREsted is a software package for efficient enhancer modeling and design ↔ src/crested/tl/losses/_cosinemse_log.py, the whole file · a weak match · score 0.63 · squared error, loss function, dynamically, MSE, cosine, sum
  23. [23] § Methods › Enhancer models › Model evaluation › Gene locus predictions ↔ src/crested/tl/_tools.py, lines 256–383 · score 0.63 · transcription start site, prediction score, upstream, downstream, window, locus
  24. [24] § CREsted is a software package for efficient enhancer modeling and design ↔ src/crested/tl/losses/_cosinemse.py, the whole file · a weak match · score 0.63 · squared error, loss function, dynamically, MSE, cosine, Optionally
  25. [25] § CREsted provides detailed insights into enhancer codes of mouse cortical cell types ↔ Figure_5/compare_models.ipynb, lines 762–806 · score 0.62 · excitatory neuron, SstChodl, oligodendrocytes, microglia, model, cell
  26. [26] § A CREsted human PBMC model captures validated TFBSs ↔ Figure_3/chip_seq_analysis.ipynb, lines 449–516 · score 0.62 · Precision recall curve, ChIP seq, AR, AP, Figure 3, Thresholding
  27. [27] § Methods › CREsted workflow › Data preprocessing › Peak-height preprocessing ↔ src/crested/_datasets.py, lines 109–174 · score 0.62 · cut sites, bigWig, consensus peaks, pseudobulk, coverage, tracks
  28. [28] § Methods › Enhancer models › Model training › gReLU mouse cortex peak regression model ↔ src/crested/_datasets.py, lines 109–174 · score 0.60 · cut site, peak regression, consensus peaks, BigWigs, cortex, mouse
  29. [29] § CREsted is a software package for efficient enhancer modeling and design ↔ src/crested/pp/_normalization.py, the whole file · a weak match · score 0.60 · low variability, standard deviation, Gini, retrieved, max, peak
  30. [30] § Methods › Enhancer models › Model training › ChromBPNet cancer cell line model ↔ src/crested/pl/corr/_heatmap.py, lines 182–330 · score 0.59 · Pearson correlation coefficient, ground truth, log transformed, split, heights, trained
  31. [31] § CREsted-trained models outperform large, pretrained models on cell-type-specific chromatin accessibility predictions ↔ src/crested/_datasets.py, lines 202–336 · score 0.58 · DeepBICCN2, transfer learning, mouse cortex, brain, Borzoi, maps
  32. [32] § Methods › Enhancer models › Input size benchmark › Comparison of seqlet locations and identified patterns ↔ src/tfmindi/pp/seqlets.py, lines 203–229 · score 0.57 · TF MInDi, recursive_seqlets, tangermeme, pp, scores
  33. [33] § Methods › Enhancer models › Model evaluation › EBF1 degradation analysis in B cells ↔ src/crested/tl/_tools.py, lines 256–383 · score 0.56 · score gene locus, crested tl, upstream, downstream
  34. [34] § Methods › CREsted workflow › Model training ↔ src/crested/tl/data/utils/_datawrapper.py, lines 429–472 · score 0.55 · TensorFlow, PyTorch, backends, Keras, CREsted, accessibility
  35. [35] § Methods › Enhancer models › Model training › ChromBPNet cancer cell line model ↔ src/crested/pl/corr/_violin.py, lines 15–89 · score 0.55 · ground truth, Pearson correlation, log transformed, fine tuned, split, heights
  36. [36] § CREsted-trained models outperform large, pretrained models on cell-type-specific chromatin accessibility predictions ↔ Figure_5/compare_models.ipynb, lines 69–122 · score 0.54 · double fine tuning, fine tuned Borzoi, frozen, consensus, peaks, model
  37. [37] § Methods › CREsted workflow › Enhancer code analysis › Pattern to TF matching ↔ src/crested/tl/modisco/_tfmodisco.py, lines 1969–2019 · score 0.53 · TF candidates, correlation threshold, Pearson, kept, cell
  38. [38] § Methods › CREsted workflow › Data preprocessing › Peak-height preprocessing ↔ src/crested/_io.py, lines 567–613 · score 0.52 · bigWig, boundaries, cut, coverage, chromosome, Optionally
  39. [39] § Methods › CREsted workflow › Model training ↔ src/crested/tl/_explainer.py, lines 83–214 · score 0.52 · TensorFlow, PyTorch, backends, Keras, sequence, predicting
  40. [40] § CREsted identifies high similarity between MES enhancer codes in cancer ↔ src/crested/_datasets.py, lines 202–336 · score 0.51 · DeepGlioma, DeepCCL, gliomas, accessibility, human, CREsted

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,382 lines · 48 KB · MIT · 4 matches

  1. # %%
  2. import pandas as pd
  3. import numpy as np
  4. import keras
  5. import crested
  6. import matplotlib.pyplot as plt
  7. import matplotlib
  8. %matplotlib inline
  9. matplotlib.rcParams["pdf.fonttype"] = 42
  10. matplotlib.rcParams["ps.fonttype"] = 42
  11. from pathlib import Path
  12. import anndata
  13. import h5py
  15. # %% [markdown]
  16. # Download the data for the notebooks from the dedicated Zenodo link of the CREsted paper. Then use it below.
  18. # %%
  19. DATA_DIR ="../../../crested_data/Figure_3/pbmc/"
  21. # %% [markdown]
  22. # # Load hg38 genome
  24. # %% [markdown]
  25. # This notebook requires a hg38 fasta file and a hg38 chromosome sizes file. You can download that here for example: https://hgdownload.soe.ucsc.edu/downloads.html
  26. # Once downloaded, we load them to the notebook.
  28. # %%
  29. genome_dir = "../../../human/genome/"
  30. genome_fasta = f"{genome_dir}hg38.fa"
  31. genome_chrom_sizes = f"{genome_dir}hg38.chrom.sizes"
  33. # %%
  34. genome = crested.Genome(genome_fasta, genome_chrom_sizes)
  35. crested.register_genome(genome)
  37. # %% [markdown]
  38. # # Load DeepPBMC
  40. # %%
  41. bigwigs_folder = f"{DATA_DIR}bw/"
  42. regions_file = f"{DATA_DIR}consensus_regions.bed"
  44. # %%
  45. adata = anndata.read_h5ad(f"{DATA_DIR}pbmc_filtered.h5ad")
  46. adata.obs_names = pd.Index([
  47. 'Bcell', 'CD14_monocyte', 'CD16_monocyte',
  48. 'CD4_Tcell', 'Cytotoxic_T_cell',
  49. 'Dendritic_cell', 'Natural_killer_cell'
  50. ])
  52. # %%
  53. model_path =f"{DATA_DIR}DeepPBMC.keras"
  54. model = keras.models.load_model(
  55. model_path, compile=False
  56. )
  58. # %% [markdown]
  59. # # Analysis on modisco top regions
  61. # %% [markdown]
  62. # First we define all relevant files, patterns and directories and add them to the tf_dict dictionary.
  64. # %%
  65. tf_dict = {}
  66. tf_list = ['PAX5','EBF1','POU2F2', 'RUNX1', 'GATA3', 'ETS1', 'CEBPA', 'SPI1']
  67. ct_list= ['Bcell','Bcell','Bcell','CD4_Tcell','CD4_Tcell','CD4_Tcell', 'CD14_monocyte', 'CD14_monocyte']
  69. chip_bed_files = [
  70. f'{DATA_DIR}/chip/pax5/pax5_CHIP.bed',
  71. f"{DATA_DIR}/chip/ebf1/ebf1_CHIP.bed",
  72. f"{DATA_DIR}/chip/pou2f2/pou2f2_chip.bed",
  73. f"{DATA_DIR}/chip/runx1/runx1_CHIP.bed",
  74. f"{DATA_DIR}/chip/gata3_2/gata3_CHIP.bed",
  75. f"{DATA_DIR}/chip/ets1/ets1_CHIP.bed",
  76. f"{DATA_DIR}/chip/cebpa/cebpa_CHIP.bed",
  77. f"{DATA_DIR}/chip/spi1_2/spi1_CHIP.bed",
  78. ]
  79. unibind_bed_files = [
  80. f"{DATA_DIR}/chip/pax5/pax5.bed",
  81. f"{DATA_DIR}/chip/ebf1/ebf1_unibind.bed",
  82. None,
  83. f"{DATA_DIR}/chip/runx1/runx1.bed",
  84. f"{DATA_DIR}/chip/gata3_2/gata3.bed",
  85. f"{DATA_DIR}/chip/ets1/ets1_unibind.bed",
  86. f"{DATA_DIR}/chip/cebpa/cebpa.bed",
  87. f"{DATA_DIR}/chip/spi1_2/spi1.bed",
  88. ]
  89. chip_bw_files = [
  90. f"{DATA_DIR}chip/pax5/PAX5.bigWig",
  91. f"{DATA_DIR}chip/ebf1/EBF1.bigWig",
  92. f"{DATA_DIR}chip/pou2f2/pou2f2.bigWig",
  93. f"{DATA_DIR}chip/runx1/runx1.bw",
  94. f"{DATA_DIR}chip/gata3_2/gata3.bw",
  95. f"{DATA_DIR}chip/ets1/ets1.bw",
  96. f"{DATA_DIR}chip/cebpa/cebpa.bw",
  97. f"{DATA_DIR}chip/spi1_2/spi1.bw",
  98. ]
  100. # Modisco pattern numbers
  101. pattern_nrs = [[5, 6, 9, 12, 13, 19, 23],
  102. [3,24],
  103. [2],
  104. [2, 4, 19, 22, 25],
  105. [7], #[3,7,27]
  106. [0],
  107. [0, 4],
  108. [2]]
  110. for i, tf in enumerate(tf_list):
  111. tf_dict[tf]={}
  112. tf_dict[tf]['cell_type']=ct_list[i]
  113. tf_dict[tf]['chip_bed']=chip_bed_files[i]
  114. tf_dict[tf]['chip_bw']=chip_bw_files[i]
  115. tf_dict[tf]['unibind_bed']=unibind_bed_files[i]
  116. tf_dict[tf]['pattern_indices']=pattern_nrs[i]
  119. top_n = 1000
  120. contribution_dir=f"{DATA_DIR}modisco/"
  121. modisco_regions=f"{DATA_DIR}modisco_regions.csv"
  123. # %% [markdown]
  124. # ## Load the patterns from the modisco files
  126. # %%
  127. def recursively_load_h5_data(h5_group):
  128. """
  129. Recursively load all data from an HDF5 group or dataset into Python objects.
  130. """
  131. if isinstance(h5_group, h5py.Group):
  132. # If it's a group, recurse into its items
  133. return {key: recursively_load_h5_data(h5_group[key]) for key in h5_group.keys()}
  134. elif isinstance(h5_group, h5py.Dataset):
  135. # If it's a dataset, load its value into memory
  136. return h5_group[()]
  137. else:
  138. # Unknown type, return as-is
  139. return h5_group
  141. for tf in tf_dict:
  142. ppms = []
  143. patterns = []
  144. h5_file = contribution_dir+'/'+tf_dict[tf]['cell_type']+'_modisco_results.h5'
  145. with h5py.File(h5_file) as hdf5_results:
  146. for metacluster_name in ['pos_patterns']:
  147. pattern_idx = 0
  148. for i in range(len(list(hdf5_results[metacluster_name]))):
  149. p = "pattern_" + str(i)
  150. pattern = hdf5_results[metacluster_name][p]
  151. pattern_data = recursively_load_h5_data(pattern)
  152. patterns.append(pattern_data)
  153. tf_dict[tf]['patterns']=patterns
  155. # %% [markdown]
  156. # ## Load the regions used for modisco
  158. # %%
  159. for tf in tf_dict:
  160. file_path = modisco_regions
  161. region_df = pd.read_csv(file_path)
  162. region_df['Class name'] = region_df['Class name'].str.lstrip('Merged__')
  163. region_df['Class name'] = region_df['Class name'].replace('B_cell','Bcell')
  164. region_df['Class name'] = region_df['Class name'].replace('CD4_T_cell','CD4_Tcell')
  166. region_df = region_df.loc[region_df['Class name']==tf_dict[tf]['cell_type']]
  167. region_df = region_df.head(top_n)
  168. region_df
  169. tf_dict[tf]['region_df']=region_df
  171. # %%
  172. tf_dict['EBF1']['region_df']
  174. # %% [markdown]
  175. # ## Load the ChIP Bed files
  177. # %%
  178. # Define the base column names for a standard BED file
  179. base_column_names = ["chrom", "start", "end", "name", "score", "strand", "signal_value", "p-val", "-log10(q-value)", "peak-summit"]
  180. chip_df_list=[]
  181. for tf in tf_dict:
  182. bed_file_path = tf_dict[tf]['chip_bed']
  184. bed_preview = pd.read_csv(bed_file_path, sep="\t", header=None, nrows=1)
  186. # Extend the column names list to match the number of columns in the file
  187. extra_columns = [f"extra_col_{i}" for i in range(len(bed_preview.columns) - len(base_column_names))]
  188. column_names = base_column_names + extra_columns
  190. # Read the BED file into a DataFrame with the dynamically generated column names
  191. chip_all_df = pd.read_csv(bed_file_path, sep="\t", header=None, names=column_names)
  193. chip_df_list.append(chip_all_df)
  195. # %% [markdown]
  196. # ## Load the UniBind Bed files
  198. # %%
  199. import pandas as pd
  201. # Define the base column names for a standard BED file
  202. base_column_names = ["chrom", "start", "end", "name", "score", "strand", "signal_value", "p-val", "-log10(q-value)", "peak-summit"]
  204. unibind_df_list=[]
  205. for tf in tf_dict:
  206. # Read the BED file to determine the number of columns
  207. bed_file_path = tf_dict[tf]['unibind_bed']
  208. if bed_file_path:
  209. bed_preview = pd.read_csv(bed_file_path, sep="\t", header=None, nrows=1)
  211. # Extend the column names list to match the number of columns in the file
  212. extra_columns = [f"extra_col_{i}" for i in range(len(bed_preview.columns) - len(base_column_names))]
  213. column_names = base_column_names + extra_columns
  215. # Read the BED file into a DataFrame with the dynamically generated column names
  216. unibind_df = pd.read_csv(bed_file_path, sep="\t", header=None, names=column_names)
  218. unibind_df_list.append(unibind_df)
  219. else:
  220. unibind_df_list.append(None)
  222. # %% [markdown]
  223. # ## Calculate overlapping Unibind hits in the top regions used in modisco
  225. # %%
  226. for i,tf in enumerate(tf_dict):
  227. possible_hits = []
  228. regions_to_process = tf_dict[tf]['region_df'].head(top_n) if top_n else tf_dict[tf]['region_df']
  229. unibind_df = unibind_df_list[i]
  230. if unibind_df is not None:
  231. for _, region in regions_to_process.iterrows():
  232. overlaps = unibind_df[
  233. (unibind_df["chrom"] == region["chr"]) &
  234. (unibind_df["start"] < region["end"] - 557) &
  235. (unibind_df["end"] > region["start"] + 557)
  236. ]
  237. for _, overlap in overlaps.iterrows():
  238. possible_hits.append({
  239. "region": region['region'],
  240. "region_chr": region["chr"],
  241. "region_start": region["start"],
  242. "region_end": region["end"],
  243. "chip_chr": overlap["chrom"],
  244. "chip_start": int(overlap["start"]),
  245. "chip_end": int(overlap["end"]),
  246. "chip_name": overlap["name"],
  247. "chip_score": overlap["score"],
  248. "chip_strand": overlap["strand"],
  249. "status": "Possible Hit"
  250. })
  252. possible_hits_df_unibind = pd.DataFrame(possible_hits)
  253. tf_dict[tf]['unibind_hits_df'] = possible_hits_df_unibind
  254. else:
  255. tf_dict[tf]['unibind_hits_df'] = None
  257. # %%
  258. tf_dict['RUNX1']['unibind_hits_df']
  260. # %% [markdown]
  261. # ## Retrieve seqlet positions inside top modisco regions
  263. # %%
  264. for tf in tf_dict:
  265. seqlet_starts=[]
  266. seqlet_chroms=[]
  267. seqlet_ends=[]
  268. seqlet_regions=[]
  269. seqlet_attributions=[]
  271. pattern_nrs= tf_dict[tf]['pattern_indices']
  272. patterns = tf_dict[tf]['patterns']
  273. region_df = tf_dict[tf]['region_df']
  275. for pattern_nr in pattern_nrs:
  276. pattern = patterns[pattern_nr]
  277. for i in range(pattern['seqlets']['n_seqlets'][0]):
  278. region = pattern['seqlets']['example_idx'][i]
  279. if region<top_n:
  280. region_id = region_df.iloc[region]['region']
  281. region_chr = region_df.iloc[region]['chr']
  282. region_start = region_df.iloc[region]['start']
  283. region_end = region_df.iloc[region]['end']
  285. seqlet_start = region_start+pattern['seqlets']['start'][i]+557
  286. seqlet_end = seqlet_start+30
  287. seqlet_chrom = region_chr
  289. average_contrib=pattern['seqlets']['contrib_scores'][i]
  290. average_contrib = average_contrib[average_contrib!=0]
  291. average_contrib = np.mean(average_contrib)
  293. seqlet_starts.append(seqlet_start)
  294. seqlet_ends.append(seqlet_end)
  295. seqlet_chroms.append(seqlet_chrom)
  296. seqlet_regions.append(region_id)
  297. seqlet_attributions.append(average_contrib)
  300. seqlet_df= pd.DataFrame({
  301. 'region': seqlet_regions,
  302. 'chr':seqlet_chroms,
  303. 'start':seqlet_starts,
  304. 'end':seqlet_ends,
  305. 'average_contrib':seqlet_attributions
  306. })
  307. tf_dict[tf]['seqlet_df']=seqlet_df
  309. # %%
  310. tf='POU2F2'
  311. tf_dict[tf]['seqlet_df']
  313. # %% [markdown]
  314. # ## Calculate precision and recall of Seqlets inside ChIP-seq peaks
  316. # %%
  317. def calculate_precision_recall(seqlet_df, chip_all_df, region_df, top_n=None):
  318. """
  319. Calculate precision and recall for overlap between seqlet_df and chip_all_df, with region_df for context.
  321. Parameters:
  322. - seqlet_df (pd.DataFrame): DataFrame with seqlets.
  323. - chip_all_df (pd.DataFrame): DataFrame with ChIP-seq data.
  324. - region_df (pd.DataFrame): DataFrame with regions for possible hits.
  325. - top_n (int, optional): Limit the number of regions processed from region_df.
  327. Returns:
  328. - results_df (pd.DataFrame): DataFrame with detailed overlap information.
  329. - precision (float): Precision value.
  330. - recall (float): Recall value.
  331. """
  332. # Initialize results and counts
  333. results = []
  334. true_positives = 0
  335. false_positives = 0
  336. false_negatives = 0
  338. # Step 1: True Positives and False Positives
  339. for _, seqlet in seqlet_df.iterrows():
  340. overlaps = chip_all_df[
  341. (chip_all_df["chrom"] == seqlet["chr"]) &
  342. (chip_all_df["start"] < seqlet["end"] + 5) &
  343. (chip_all_df["end"] > seqlet["start"] - 5)
  344. ]
  345. if not overlaps.empty:
  346. for _, overlap in overlaps.iterrows():
  347. results.append({
  348. "region": seqlet['region'],
  349. "seqlet_chr": seqlet["chr"],
  350. "seqlet_start": seqlet["start"],
  351. "seqlet_end": seqlet["end"],
  352. "seqlet_average_contrib": seqlet["average_contrib"],
  353. "chip_chr": overlap["chrom"],
  354. "chip_start": int(overlap["start"]),
  355. "chip_end": int(overlap["end"]),
  356. "chip_name": overlap["name"],
  357. "chip_score": overlap["score"],
  358. "chip_strand": overlap["strand"],
  359. "status": "True Positive"
  360. })
  361. true_positives += 1
  362. else:
  363. results.append({
  364. "region": seqlet['region'],
  365. "seqlet_chr": seqlet["chr"],
  366. "seqlet_start": seqlet["start"],
  367. "seqlet_end": seqlet["end"],
  368. "seqlet_average_contrib": seqlet["average_contrib"],
  369. "chip_chr": np.nan,
  370. "chip_start": np.nan,
  371. "chip_end": np.nan,
  372. "chip_name": np.nan,
  373. "chip_score": np.nan,
  374. "chip_strand": np.nan,
  375. "status": "False Positive"
  376. })
  377. false_positives += 1
  379. # Step 2: Possible Hits (region_df vs chip_all_df)
  380. possible_hits = []
  381. regions_to_process = region_df.head(top_n) if top_n else region_df
  382. for _, region in regions_to_process.iterrows():
  383. overlaps = chip_all_df[
  384. (chip_all_df["chrom"] == region["chr"]) &
  385. (chip_all_df["start"] < region["end"] - 557) &
  386. (chip_all_df["end"] > region["start"] + 557)
  387. ]
  388. for _, overlap in overlaps.iterrows():
  389. possible_hits.append({
  390. "region": region['region'],
  391. "region_chr": region["chr"],
  392. "region_start": region["start"],
  393. "region_end": region["end"],
  394. "chip_chr": overlap["chrom"],
  395. "chip_start": int(overlap["start"]),
  396. "chip_end": int(overlap["end"]),
  397. "chip_name": overlap["name"],
  398. "chip_score": overlap["score"],
  399. "chip_strand": overlap["strand"],
  400. "status": "Possible Hit"
  401. })
  403. possible_hits_df = pd.DataFrame(possible_hits)
  405. # Step 3: False Negatives (possible hits not in seqlet_df)
  406. for _, hit in possible_hits_df.iterrows():
  407. matches = seqlet_df[
  408. (seqlet_df["chr"] == hit["chip_chr"]) &
  409. (seqlet_df["start"] < hit["chip_end"] + 5) &
  410. (seqlet_df["end"] > hit["chip_start"] - 5)
  411. ]
  412. if matches.empty:
  413. results.append({
  414. "region": hit['region'],
  415. "seqlet_chr": np.nan,
  416. "seqlet_start": np.nan,
  417. "seqlet_end": np.nan,
  418. "seqlet_average_contrib": np.nan,
  419. "chip_chr": hit["chip_chr"],
  420. "chip_start": hit["chip_start"],
  421. "chip_end": hit["chip_end"],
  422. "chip_name": hit["chip_name"],
  423. "chip_score": hit["chip_score"],
  424. "chip_strand": hit["chip_strand"],
  425. "status": "False Negative"
  426. })
  427. false_negatives += 1
  429. # Convert results to DataFrame
  430. results_df = pd.DataFrame(results)
  432. # Calculate precision and recall
  433. precision = true_positives / (true_positives + false_positives) if true_positives + false_positives > 0 else 0
  434. recall = true_positives / (true_positives + false_negatives) if true_positives + false_negatives > 0 else 0
  436. return results_df, precision, recall, true_positives, false_positives, false_negatives, possible_hits_df
  438. # %%
  439. # Assuming seqlet_df, chip_all_df, and region_df are defined
  440. for i, tf in enumerate(tf_dict):
  441. print(tf)
  442. results_df, precision, recall, tp, fp, fn, possible_hits_df = calculate_precision_recall(tf_dict[tf]['seqlet_df'], chip_df_list[i], tf_dict[tf]['region_df'], top_n)
  443. tf_dict[tf]['results_df']=results_df
  444. tf_dict[tf]['chip_hits_df']=possible_hits_df
  445. # Display results
  446. print(f"Precision: {precision:.3f}")
  447. print(f"Recall: {recall:.3f}")
  449. # %% [markdown]
  450. # ### Fig. 3e
  452. # %%
  453. from tqdm import tqdm
  454. def calculate_precision_recall_by_threshold(seqlet_df, chip_all_df, region_df, thresholds, top_n=None):
  455. """
  456. Calculate precision and recall for varying thresholds on 'average_contrib' in seqlet_df.
  458. Parameters:
  459. - seqlet_df (pd.DataFrame): DataFrame with seqlets.
  460. - chip_all_df (pd.DataFrame): DataFrame with ChIP-seq data.
  461. - region_df (pd.DataFrame): DataFrame with regions for possible hits.
  462. - thresholds (list or np.array): Thresholds for filtering seqlet_df by 'average_contrib'.
  463. - top_n (int, optional): Limit the number of regions processed from region_df.
  465. Returns:
  466. - precision_list (list): List of precision values for each threshold.
  467. - recall_list (list): List of recall values for each threshold.
  468. - thresholds (list): List of thresholds used.
  469. - avg_precision (float): Average precision across thresholds.
  470. - avg_recall (float): Average recall across thresholds.
  471. """
  472. precision_list = []
  473. recall_list = []
  475. for threshold in thresholds:
  476. # Filter seqlet_df by current threshold
  477. filtered_seqlet_df = seqlet_df.loc[seqlet_df['average_contrib'] > threshold].reset_index(drop=True)
  479. # Calculate precision and recall using the previous function
  480. _, precision, recall, _, _, _, _ = calculate_precision_recall(filtered_seqlet_df, chip_all_df, region_df, top_n)
  482. # Store results
  483. precision_list.append(precision)
  484. recall_list.append(recall)
  486. # Calculate average precision and recall
  487. avg_precision = np.mean(precision_list)
  488. avg_recall = np.mean(recall_list)
  490. return precision_list, recall_list, thresholds, avg_precision, avg_recall
  492. plt.figure(figsize=(6, 6))
  494. for i, tf in tqdm(enumerate(tf_list), total=len(tf_dict)):
  495. # Define thresholds
  496. thresholds = np.linspace(np.min(tf_dict[tf]['seqlet_df']['average_contrib'].values),
  497. np.max(tf_dict[tf]['seqlet_df']['average_contrib'].values) * 0.95, 5)
  499. # Calculate precision, recall, and averages
  500. precision_list, recall_list, thresholds, avg_precision, avg_recall = calculate_precision_recall_by_threshold(
  501. tf_dict[tf]['seqlet_df'], chip_df_list[i], tf_dict[tf]['region_df'], thresholds
  502. )
  504. # Plot all PR curves on the same figure
  505. plt.plot(recall_list, precision_list, marker='o', label=f'{tf} (AP={avg_precision:.3f}, AR={avg_recall:.3f})')
  507. # Formatting
  508. plt.title("Precision-Recall Curves by Thresholding on 'average_contrib'")
  509. plt.xlabel("Recall")
  510. plt.ylabel("Precision")
  511. plt.grid(True, linestyle='--', alpha=0.7)
  512. plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
  513. plt.xlim([0, 1])
  514. plt.ylim([0, 1])
  515. #plt.savefig('paperfigs/all_TFs_PR.pdf', bbox_inches='tight')
  516. plt.show()
  519. # %%
  520. tf_dict[tf]['chip_hits_df'] # Dataframe for ChIP-seq peaks inside top modisco regions
  522. # %% [markdown]
  523. # ## Closer look at example regions
  525. # %%
  526. tf='POU2F2'
  527. region_index = 3
  528. region_df = tf_dict[tf]['region_df']
  529. results_df = tf_dict[tf]['results_df']
  530. region_id = region_df.iloc[region_index]['region']
  531. region_data = results_df[results_df['region'] == region_id]
  533. chrom=region_id.split(':')[0]
  534. start=int(region_id.split(':')[1].split('-')[0])
  535. end=int(region_id.split(':')[1].split('-')[1])
  538. # Extract non-NaN seqlet coordinates
  539. seqlet_coords = region_data.dropna(subset=["seqlet_chr"]).apply(
  540. lambda row: (int(row["seqlet_start"]-start), int(row["seqlet_end"]-start)), axis=1
  541. ).tolist() if len(region_data.dropna(subset=["seqlet_chr"]))>0 else None
  543. # Extract non-NaN ChIP coordinates
  544. chip_coords = region_data.dropna(subset=["chip_chr"]).apply(
  545. lambda row: (int(row["chip_start"]-start), int(row["chip_end"]-start)), axis=1
  546. ).tolist() if len(region_data.dropna(subset=["chip_chr"]))>0 else None
  548. seqlet_coords
  549. chip_coords
  550. region_data
  552. # %%
  553. adata.obs_names
  555. # %%
  556. classes_of_interest = [tf_dict[tf]['cell_type']]
  557. sequence = genome.fetch(chrom, start, end).upper()
  559. class_idx = list(adata.obs_names.get_indexer(classes_of_interest))
  561. scores, one_hot_encoded_sequences = crested.tl.contribution_scores(
  562. sequence,
  563. target_idx=class_idx,
  564. model=model,
  565. )
  567. # %%
  568. # Highlight identified seqlets
  569. crested.pl.patterns.contribution_scores(
  570. scores,
  571. one_hot_encoded_sequences,
  572. sequence_labels=[''],
  573. class_labels=classes_of_interest,
  574. zoom_n_bases=500,
  575. title="Test set region",
  576. height=3.5,
  577. highlight_positions=seqlet_coords
  578. )
  580. # %%
  581. # Highlight ChIP peaks
  582. crested.pl.patterns.contribution_scores(
  583. scores,
  584. one_hot_encoded_sequences,
  585. sequence_labels=[''],
  586. class_labels=classes_of_interest,
  587. zoom_n_bases=500,
  588. title="Test set region",
  589. height=3.5,
  590. highlight_positions=chip_coords
  591. )
  593. # %%
  594. # Plot ChIP BigWigs
  595. for chip_bw in chip_bw_files:
  596. filename = chip_bw.split('/')[-1].split('.')[0].upper() # Extract and format filename
  598. # Read BigWig region
  599. vals, _ = crested.utils.read_bigwig_region(chip_bw, (chrom, start + 807, end - 807))
  601. # Plot
  602. plt.figure(figsize=(50, 2))
  603. plt.plot(np.arange(0, 500), vals)
  605. # Add filename as text in the top-left corner
  606. plt.text(0.02, 0.95, filename, transform=plt.gca().transAxes, fontsize=18, fontweight='bold',
  607. verticalalignment='top',horizontalalignment='left', bbox=dict(facecolor='white', alpha=0.7, edgecolor='black'))
  609. plt.show()
  611. # %%
  612. # Plot all possible seqlets identified by tangermeme
  614. from tangermeme.seqlet import recursive_seqlets
  615. from crested.utils import one_hot_encode_sequence
  617. oh_seq = one_hot_encode_sequence(sequence)[0].T
  618. X_attr = scores[0,0].T
  619. X_attr=X_attr*oh_seq
  620. X_attr = np.expand_dims(np.sum(X_attr, axis=0),axis=0)
  622. seqlets = recursive_seqlets(X_attr, threshold=0.05)
  624. highlight_positions=[]
  625. for i in range(len(seqlets)):
  626. start_ = seqlets.iloc[i]['start']
  627. end_ = seqlets.iloc[i]['end']
  628. highlight_positions.append((start_,end_))
  630. crested.pl.patterns.contribution_scores(
  631. scores,
  632. one_hot_encoded_sequences,
  633. sequence_labels=[''],
  634. class_labels=classes_of_interest,
  635. zoom_n_bases=500,
  636. title="Test set region",
  637. height=3.5,
  638. highlight_positions=highlight_positions
  639. )
  640. seqlets
  642. # %%
  643. # Plot region predictions
  644. sequence = genome.fetch(chrom, start, end).upper()
  646. prediction = crested.tl.predict(sequence, model)
  649. crested.pl.bar.prediction(prediction, classes=list(adata.obs_names))
  652. # %%
  653. # Plot scATAC profile
  655. bigwigs_folder = f"{DATA_DIR}/bw/"
  657. vals = np.zeros(len(list(adata.obs_names)))
  658. for i, ct in enumerate(list(adata.obs_names)):
  659. print(ct)
  660. values = crested.utils.read_bigwig_region(bigwigs_folder+ct+'.bw', (chrom,start+557,end-557))
  661. bw_values=values[0]
  662. midpoints=values[1]
  663. vals[i] = np.sum(bw_values)*adata.obsm['weights'][i]
  665. plt.figure(figsize=(20,3))
  666. plt.bar(list(adata.obs_names), vals)
  667. plt.title('ATAC profile')
  668. plt.show()
  670. # %% [markdown]
  671. # ## Recall of UniBind hits with seqlet calling
  673. # %%
  674. def get_row_with_overlap(df, input_start, input_end, overlap=0.5):
  675. """
  676. Retrieves the row with at least 30% overlap with the given input range.
  678. Parameters:
  679. df (pd.DataFrame): Input DataFrame with 'start' and 'end' columns.
  680. input_start (int): Start of the input range.
  681. input_end (int): End of the input range.
  682. overlap (float): Minimum overlap fraction required.
  684. Returns:
  685. pd.DataFrame: Row(s) with at least the specified overlap, sorted by overlap.
  686. """
  687. def overlap_fraction(row):
  688. # Calculate overlap
  689. overlap_start = max(row['start'], input_start)
  690. overlap_end = min(row['end'], input_end)
  691. overlap_length = max(0, overlap_end - overlap_start + 1) # Ensure non-negative
  692. row_length = row['end'] - row['start'] + 1
  693. return overlap_length / row_length
  695. # Calculate overlap for each row and add as a new column
  696. df['overlap'] = df.apply(overlap_fraction, axis=1)
  698. # Filter rows with at least the specified overlap
  699. overlap_rows = df[df['overlap'] >= overlap]
  701. # Sort the filtered rows by the overlap column
  702. overlap_rows = overlap_rows.sort_values(by='overlap', ascending=False)
  704. return overlap_rows
  706. # %%
  707. from tqdm import tqdm
  709. for tf in tf_dict:
  710. print(tf)
  711. filename=contribution_dir+'/'+tf_dict[tf]['cell_type']+'_contrib.npz'
  712. scores = np.load(filename)['arr_0']
  713. scores = np.transpose(scores,(0,2,1))
  714. oh_seqs=np.load(contribution_dir+'/'+tf_dict[tf]['cell_type']+'_oh.npz')['arr_0']
  715. oh_seqs = np.transpose(oh_seqs,(0,2,1))
  717. seqlet_starts=[]
  718. seqlet_ends = []
  719. p_values=[]
  720. attributions=[]
  721. attributions_exact=[]
  722. attributions_exact_padded=[]
  724. if tf_dict[tf]['unibind_hits_df'] is None:
  725. tf_dict[tf]['unibind_recall_df'] = None
  726. continue
  729. for i in tqdm(range(top_n)):
  730. oh_seq = oh_seqs[i].T
  731. X_attr = scores[i].T
  732. X_attr=X_attr*oh_seq
  733. X_attr = np.expand_dims(np.sum(X_attr, axis=0),axis=0)
  735. seqlets = recursive_seqlets(X_attr, threshold=0.05)
  736. region_id=tf_dict[tf]['region_df'].iloc[i]['region']
  738. if region_id in tf_dict[tf]['unibind_hits_df']['region'].values:
  739. rows = tf_dict[tf]['unibind_hits_df'].loc[tf_dict[tf]['unibind_hits_df']['region']==region_id]
  740. if len(rows)==1:
  741. row=rows
  742. chrom= row['region_chr']
  743. start=int(row['region_start'].iloc[0])
  744. end=int(row['region_end'].iloc[0])
  746. tf_start=int(row['chip_start'].iloc[0]-start)
  747. tf_end = int(row['chip_end'].iloc[0]-start)
  749. overlap_rows = get_row_with_overlap(seqlets, tf_start, tf_end)
  750. if len(overlap_rows)>0:
  751. r = overlap_rows.iloc[0]
  753. seq_start=start+int(r['start'])
  754. seq_end = start+int(r['end'])
  756. seqlet_starts.append(seq_start)
  757. seqlet_ends.append(seq_end)
  758. p_values.append(r['p-value'])
  759. attributions.append(r['attribution'])
  761. attribution = (np.mean(X_attr[0,int(r['start']):int(r['end'])]))
  762. padding=(30-(int(r['end'])-int(r['start'])))//2
  763. attribution_padded = (np.mean(X_attr[0,int(r['start'])-padding:int(r['end'])+padding]))
  764. attributions_exact.append(attribution)
  765. attributions_exact_padded.append(attribution_padded)
  767. else:
  768. seqlet_starts.append(np.nan)
  769. seqlet_ends.append(np.nan)
  770. p_values.append(np.nan)
  771. attributions.append(np.nan)
  772. attributions_exact.append(np.nan)
  773. attributions_exact_padded.append(np.nan)
  775. else:
  776. for i in range(len(rows)):
  777. chrom= rows['region_chr']
  778. start=int(rows['region_start'].iloc[i])
  779. end=int(rows['region_end'].iloc[i])
  781. tf_start=int(rows['chip_start'].iloc[i]-start)
  782. tf_end = int(rows['chip_end'].iloc[i]-start)
  784. overlap_rows = get_row_with_overlap(seqlets, tf_start, tf_end)
  785. if len(overlap_rows)>0:
  786. r = overlap_rows.iloc[0]
  787. seq_start=start+int(r['start'])
  788. seq_end = start+int(r['end'])
  790. seqlet_starts.append(seq_start)
  791. seqlet_ends.append(seq_end)
  792. p_values.append(r['p-value'])
  793. attributions.append(r['attribution'])
  795. attribution = (np.mean(X_attr[0,int(r['start']):int(r['end'])]))
  796. padding=(30-(int(r['end'])-int(r['start'])))//2
  797. attribution_padded = (np.mean(X_attr[0,int(r['start'])-padding:int(r['end'])+padding]))
  798. attributions_exact.append(attribution)
  799. attributions_exact_padded.append(attribution_padded)
  801. else:
  802. seqlet_starts.append(np.nan)
  803. seqlet_ends.append(np.nan)
  804. p_values.append(np.nan)
  805. attributions.append(np.nan)
  806. attributions_exact.append(np.nan)
  807. attributions_exact_padded.append(np.nan)
  808. df_final = tf_dict[tf]['unibind_hits_df'].head(top_n).copy()
  809. df_final.loc[:,'seqlet_start'] = seqlet_starts
  810. df_final.loc[:,'seqlet_end'] = seqlet_ends
  811. df_final.loc[:,'seqlet_p_val']= p_values
  812. df_final.loc[:,'seqlet_attribution'] = attributions
  813. df_final.loc[:,'chip_attribution'] = attributions_exact
  814. df_final.loc[:,'chip_attribution_padded'] = attributions_exact_padded
  816. tf_dict[tf]['unibind_recall_df']=df_final
  818. tf_dict['PAX5']['unibind_recall_df']
  821. # %%
  822. recalls ={}
  824. for tf in tf_dict:
  825. if tf_dict[tf]['unibind_recall_df'] is not None:
  826. print(tf)
  827. # Determine true positives (non-NaN in seqlet_start)
  828. true_positives = tf_dict[tf]['unibind_recall_df']['seqlet_start'].notna().sum()
  830. # Determine false negatives (NaN in seqlet_start)
  831. false_negatives = tf_dict[tf]['unibind_recall_df']['seqlet_start'].isna().sum()
  833. # Calculate recall
  834. recall = true_positives / (true_positives + false_negatives)
  836. recalls[tf] = {}
  837. recalls[tf]['recall']=recall
  838. recalls[tf]['total']=true_positives+false_negatives
  840. # Output the results
  841. print(f"True Positives (TP): {true_positives}")
  842. print(f"False Negatives (FN): {false_negatives}")
  843. print(f"Recall: {recall:.2f}")
  844. else:
  845. print('No unibind bed file provided for '+tf)
  847. # %% [markdown]
  848. # ### Fig. 3f
  850. # %%
  851. # Sorting by recall values
  852. sorted_factors = sorted(recalls, key=lambda x: recalls[x]['recall'], reverse=True)
  853. sorted_recalls = [recalls[factor]['recall'] for factor in sorted_factors]
  854. sorted_totals = [recalls[factor]['total'] for factor in sorted_factors]
  856. # Generate new x-axis labels with (n = total)
  857. x_labels = [f"{factor}\n(n = {total})" for factor, total in zip(sorted_factors, sorted_totals)]
  859. # Plotting
  860. fig, ax = plt.subplots(figsize=(8, 5))
  861. bars = ax.bar(x_labels, sorted_recalls, color='royalblue', edgecolor='black', linewidth=1.2)
  863. # Aesthetics
  864. ax.set_ylabel('Recall', fontsize=14)
  865. ax.set_xlabel('Transcription Factors', fontsize=14)
  866. ax.set_title('Recall Scores of Identified Seqlets in Unibind Sites', fontsize=16)
  867. ax.set_ylim(0, 1.05)
  868. ax.grid(axis='y', linestyle='--', alpha=0.7)
  870. # Adding value labels
  871. for bar in bars:
  872. yval = bar.get_height()
  873. ax.text(bar.get_x() + bar.get_width()/2, yval + 0.02, f'{yval:.3f}',
  874. ha='center', va='bottom', fontsize=12, color='black')
  876. # Remove top and right borders
  877. ax.spines['top'].set_visible(False)
  878. ax.spines['right'].set_visible(False)
  880. # Show the plot
  881. plt.xticks(rotation=45, fontsize=12, ha='right', va='top')
  882. plt.yticks(fontsize=12)
  883. #plt.savefig('paperfigs/unibind_recall.pdf', bbox_inches='tight')
  884. plt.show()
  886. # %% [markdown]
  887. # ## Check ChIP and ATAC peak heights in identified seqlets
  889. # %%
  890. for tf in tf_dict:
  891. tf_dict[tf]['atac_bw_file']="/home/VIB.LOCAL/niklas.kempynck/nkemp/pbmc/bw/"+tf_dict[tf]['cell_type']+".bw"
  893. # %%
  894. import pandas as pd
  895. import matplotlib.pyplot as plt
  896. from crested.utils import read_bigwig_region # Ensure you have the required module
  898. def plot_average_peak_height(bigwig_file, region_df, padding=0, seqlet_len=30, cols=None,max_carrot=None, title='', save_path=None):
  899. """
  900. Plot a histogram of average peak heights and a heatmap carrot plot for all peak values.
  902. Parameters:
  903. - bigwig_file (str): Path to the BigWig file.
  904. - region_df (pd.DataFrame): DataFrame with columns ['region', 'chr', 'start', 'end'].
  905. - padding (int): Number of bases to add as flanking around the regions.
  906. - seqlet_len (int): Expected sequence length for normalization.
  907. - cols (list): Column names for chrom, start, end. Defaults to ['chr', 'start', 'end'].
  908. """
  909. # Setup columns
  910. if cols:
  911. chrom_col = cols[0]
  912. start_col = cols[1]
  913. end_col = cols[2]
  914. else:
  915. chrom_col = 'chr'
  916. start_col = 'start'
  917. end_col = 'end'
  919. # Initialize peak values
  920. window_length = seqlet_len + 2 * padding
  921. peaks = np.zeros((len(region_df), window_length))
  922. avg_peak_heights = []
  924. for i, row in region_df.iterrows():
  925. chrom = row[chrom_col]
  926. start = int(row[start_col])
  927. end = int(row[end_col])
  928. seq_len = end - start
  930. # Normalize to seqlet_len
  931. if seq_len > seqlet_len:
  932. diff = (seq_len - seqlet_len) // 2
  933. start = start + diff
  934. end = end - diff
  935. elif seq_len < seqlet_len:
  936. diff = (seqlet_len - seq_len) // 2
  937. start = start - diff
  938. end = end + diff
  940. seq_len = end-start
  942. if seq_len<seqlet_len:
  943. end=end+1
  945. # Read the BigWig region with flanking
  946. values, _ = read_bigwig_region(bigwig_file, (chrom, start - padding, end + padding))
  948. # Handle peak values
  949. #if len(values) >= window_length:
  950. peaks[i] = values[:window_length]
  951. avg_peak_heights.append(np.mean(values))
  952. #else:
  953. # peaks[i] = np.pad(values, (0, window_length - len(values)), 'constant', constant_values=0)
  954. # avg_peak_heights.append(np.mean(values) if len(values) > 0 else 0)
  955. avg_peaks = np.mean(peaks,axis=0)
  957. num_figs= 3 if max_carrot else 2
  958. # Plot histogram of average peak heights
  959. fig, axes = plt.subplots(1, num_figs , figsize=(num_figs*4, 6))
  960. x_start = -window_length/2
  961. x_end = window_length/2
  962. # Histogram
  963. # Plot
  964. axes[0].plot(np.arange(x_start,x_end),avg_peaks)
  965. axes[0].set_title("Average Peak Height in Regions "+title)
  966. axes[0].set_xlabel("Position")
  967. axes[0].set_ylabel("Average Peak Height")
  968. axes[0].grid(axis='y', linestyle='--', alpha=0.7)
  971. peak_max_values = np.max(peaks, axis=1) # Find max value for each row
  972. sorted_indices = np.argsort(peak_max_values)[::-1] # Sort indices in descending order
  973. sorted_peaks = peaks[sorted_indices] # Reorder rows based on sorted indices
  975. # Heatmap carrot plot
  976. im = axes[1].imshow(
  977. sorted_peaks,
  978. aspect='auto',
  979. cmap='coolwarm',
  980. extent=[-window_length // 2, window_length // 2, 0, len(region_df)]
  981. )
  982. axes[1].set_title("Carrot Plot of Peak Values "+(title))
  983. axes[1].set_xlabel("Position (relative to center)")
  984. axes[1].set_ylabel("Region Index")
  985. fig.colorbar(im, ax=axes[1], orientation='vertical', label='Peak Value')
  987. if max_carrot:
  988. peaks_carrot = np.clip(peaks, 0, max_carrot)
  989. peak_max_values = np.max(peaks_carrot, axis=1) # Find max value for each row
  990. sorted_indices = np.argsort(peak_max_values)[::-1] # Sort indices in descending order
  991. sorted_peaks = peaks_carrot[sorted_indices] # Reorder rows based on sorted indices
  993. # Heatmap carrot plot
  994. im = axes[2].imshow(
  995. sorted_peaks,
  996. aspect='auto',
  997. cmap='coolwarm',
  998. extent=[-window_length // 2, window_length // 2, 0, len(region_df)]
  999. )
  1000. axes[2].set_title("Carrot Plot of Peak Values - Clipped "+(title))
  1001. axes[2].set_xlabel("Position (relative to center)")
  1002. axes[2].set_ylabel("Region Index")
  1003. fig.colorbar(im, ax=axes[2], orientation='vertical', label='Peak Value')
  1005. plt.tight_layout()
  1006. if save_path:
  1007. plt.savefig(save_path, bbox_inches='tight')
  1008. plt.show()
  1009. return peaks
  1012. # %%
  1013. tf = 'CEBPA'
  1015. bw_file =tf_dict[tf]['chip_bw']
  1016. seqlet_df = tf_dict[tf]['seqlet_df']
  1017. results_df = tf_dict[tf]['results_df']
  1019. possible_hits_df = tf_dict[tf]['chip_hits_df']
  1020. possible_hits_df_unibind = tf_dict[tf]['unibind_hits_df']
  1022. peaks=plot_average_peak_height(bw_file, seqlet_df, padding=50, max_carrot=5, title='(All seqlets)')#, save_path='paperfigs/'+tf+'_all_seqlets.pdf')
  1023. results_df_tp = results_df.loc[results_df['status']=='True Positive'].reset_index()
  1024. _=plot_average_peak_height(bw_file, results_df_tp, padding=50, cols=['seqlet_chr','seqlet_start','seqlet_end'],max_carrot=5,title='(True positives)')
  1025. results_df_fp = results_df.loc[results_df['status']=='False Positive'].reset_index()
  1026. _=plot_average_peak_height(bw_file, results_df_fp, padding=50, cols=['seqlet_chr','seqlet_start','seqlet_end'],max_carrot=0.12, title='(False positives)')#, save_path='paperfigs/'+tf+'_FP.pdf')
  1027. results_df_fn = results_df.loc[results_df['status']=='False Negative'].reset_index()
  1028. _=plot_average_peak_height(bw_file, results_df_fn, padding=50, cols=['chip_chr','chip_start','chip_end'],max_carrot=5,title='(False negatives)')
  1029. _=plot_average_peak_height(bw_file, possible_hits_df, padding=50, cols=['chip_chr','chip_start','chip_end'],max_carrot=5,title='(Chip sites)')
  1030. _=plot_average_peak_height(bw_file, possible_hits_df_unibind, padding=50, cols=['chip_chr','chip_start','chip_end'],max_carrot=5,title='(Unibind sites)')
  1034. # %%
  1035. results_df
  1037. # %%
  1038. import numpy as np
  1039. import matplotlib.pyplot as plt
  1041. def plot_combined_average_peak_heights(bigwig_file, subset_dict, padding=0, seqlet_len=30, cols_dict=None, title='', save_path=None):
  1042. """
  1043. Plot a combined average peak height plot for multiple subsets.
  1045. Parameters:
  1046. - bigwig_file (str): Path to the BigWig file.
  1047. - subset_dict (dict): Dictionary where keys are subset names and values are DataFrames with ['chr', 'start', 'end'].
  1048. - padding (int): Number of bases to add as flanking around the regions.
  1049. - seqlet_len (int): Expected sequence length for normalization.
  1050. - cols (list): Column names for chrom, start, end. Defaults to ['chr', 'start', 'end'].
  1051. - title (str): Title of the plot.
  1052. """
  1054. window_length = seqlet_len + 2 * padding
  1055. x_positions = np.arange(-window_length / 2, window_length / 2)
  1057. plt.figure(figsize=(2, 4))
  1059. max_peak = 0
  1061. for subset_name, region_df in subset_dict.items():
  1062. if region_df is None:
  1063. continue
  1064. if cols_dict:
  1065. cols = cols_dict[subset_name]
  1067. num_instances = len(region_df) # Get the number of instances
  1068. label = f"{subset_name} (n={num_instances})" # Modify the label
  1069. # Extract peak values
  1070. peaks = np.zeros((len(region_df), window_length))
  1072. for i, row in region_df.iterrows():
  1073. chrom = row[cols[0]] if cols else row['chr']
  1074. start = int(row[cols[1]]) if cols else int(row['start'])
  1075. end = int(row[cols[2]]) if cols else int(row['end'])
  1076. seq_len = end - start
  1078. diff=seq_len
  1079. # Normalize to seqlet_len
  1080. if seq_len > seqlet_len:
  1081. diff = (seq_len - seqlet_len) // 2
  1082. start += diff
  1083. end -= diff
  1084. elif seq_len < seqlet_len:
  1085. diff = (seqlet_len - seq_len) // 2
  1086. start -= diff
  1087. end += diff
  1089. if diff<seq_len:
  1090. end=end+1
  1092. # Read BigWig region
  1093. values, _ = read_bigwig_region(bigwig_file, (chrom, start - padding, end + padding))
  1094. peaks[i] = values[:window_length]
  1096. avg_peaks = np.mean(peaks, axis=0)
  1097. if np.max(avg_peaks)>max_peak:
  1098. max_peak = np.max(avg_peaks)
  1100. # Plot
  1101. plt.plot(x_positions, avg_peaks, label=label)
  1103. plt.xlabel("Position (relative to center)")
  1104. plt.ylabel("Average Peak Height")
  1105. plt.title(title)
  1106. plt.ylim([0,max_peak+0.05*max_peak])
  1107. plt.legend(loc='center left', bbox_to_anchor=(1, 0.7))
  1108. if save_path:
  1109. plt.savefig(save_path, bbox_inches='tight')
  1110. plt.show()
  1112. # %% [markdown]
  1113. # ### Fig. 3g
  1115. # %%
  1116. for tf in tf_list:
  1117. print(tf)
  1118. bw_file =tf_dict[tf]['chip_bw']
  1119. seqlet_df = tf_dict[tf]['seqlet_df']
  1120. results_df = tf_dict[tf]['results_df']
  1122. possible_hits_df = tf_dict[tf]['chip_hits_df']
  1123. possible_hits_df_unibind = tf_dict[tf]['unibind_hits_df']
  1125. subset_dict = {
  1126. "All Seqlets": seqlet_df,
  1127. "True Positives": results_df.loc[results_df['status'] == 'True Positive'].reset_index(),
  1128. "False Positives": results_df.loc[results_df['status'] == 'False Positive'].reset_index(),
  1129. "False Negatives": results_df.loc[results_df['status'] == 'False Negative'].reset_index(),
  1130. "Chip Sites": possible_hits_df,
  1131. "Unibind": possible_hits_df_unibind
  1132. }
  1134. cols_dict = {
  1135. "All Seqlets": ['chr', 'start', 'end'],
  1136. "True Positives": ['seqlet_chr','seqlet_start','seqlet_end'],
  1137. "False Positives": ['seqlet_chr','seqlet_start','seqlet_end'],
  1138. "False Negatives": ['chip_chr','chip_start','chip_end'],
  1139. "Chip Sites": ['chip_chr','chip_start','chip_end'],
  1140. "Unibind": ['chip_chr','chip_start','chip_end']
  1141. }
  1143. plot_combined_average_peak_heights(bw_file, subset_dict, padding=200, cols_dict=cols_dict)# title="Average Peak for "+tf, save_path='paperfigs/'+tf+'_peaks.pdf')
  1145. # %% [markdown]
  1146. # # Recall of seqlets on UniBind hits inside all consensus peaks
  1148. # %%
  1149. peak_bed=f"{DATA_DIR}consensus_regions.bed"
  1150. peak_df = pd.read_csv(peak_bed, sep="\t", header=None, names=column_names)
  1151. peak_df
  1153. # %% [markdown]
  1154. # ## Calculate overlapping ChIP peaks and consensus regions
  1156. # %%
  1157. from tqdm import tqdm
  1158. ct_list= ['Bcell',
  1159. 'Bcell',
  1160. 'CD4_Tcell',
  1161. 'CD4_Tcell',
  1162. 'CD4_Tcell',
  1163. 'CD14_monocyte',
  1164. 'CD14_monocyte'
  1165. ]
  1167. chip_bed_files = [
  1168. f"{DATA_DIR}/chip/pax5/pax5.bed",
  1169. f"{DATA_DIR}/chip/ebf1/ebf1_unibind.bed",
  1170. f"{DATA_DIR}/chip/runx1/runx1.bed",
  1171. f"{DATA_DIR}/chip/gata3_2/gata3.bed",
  1172. f"{DATA_DIR}/chip/ets1/ets1_unibind.bed",
  1173. f"{DATA_DIR}/chip/cebpa/cebpa.bed",
  1174. f"{DATA_DIR}/chip/spi1_2/spi1.bed",
  1175. ]
  1177. unibind_dfs = []
  1179. for chip_bed_file in tqdm(chip_bed_files):
  1180. print(chip_bed_file)
  1181. column_names = ["chrom", "start", "end", "name", "score", "strand"]
  1183. # Read the unibind bed files
  1184. chip_df = pd.read_csv(chip_bed_file, sep="\t", header=None, names=column_names)
  1186. # Merge the two DataFrames on overlapping intervals
  1187. merged_df = pd.merge(
  1188. chip_df,
  1189. peak_df,
  1190. on="chrom",
  1191. suffixes=("_chip", "_peak")
  1192. )
  1194. # Filter for overlaps
  1195. subset_chip_df = merged_df[
  1196. (merged_df["start_chip"] >= merged_df["start_peak"]) &
  1197. (merged_df["end_chip"] <= merged_df["end_peak"])
  1198. ]
  1200. sorted_chip_df = subset_chip_df.sort_values(by="score_chip", ascending=False)
  1201. unibind_dfs.append(sorted_chip_df)
  1204. # %% [markdown]
  1205. # ## Calculate contribution scores for all overlapping regions
  1206. # WARNING: TAKES A LONG TIME
  1207. #
  1208. # just load the results
  1210. # %%
  1211. #out_dirs = [
  1212. # f'{DATA_DIR}chip/pax5/',
  1213. # f"{DATA_DIR}chip/ebf1/",
  1214. # f"{DATA_DIR}chip/runx1/",
  1215. # f"{DATA_DIR}chip/gata3_2/",
  1216. # f"{DATA_DIR}chip/ets1/",
  1217. # f"{DATA_DIR}chip/cebpa/",
  1218. # f"{DATA_DIR}chip/spi1_2/",
  1219. #]
  1220. #for a, sorted_chip_df in enumerate(unibind_dfs):
  1221. # top_n=len(sorted_chip_df)
  1222. #
  1223. # sequences=[]
  1224. # for i in range(top_n):
  1225. # row= sorted_chip_df.iloc[i]
  1226. #
  1227. # chrom= row['chrom']
  1228. # start=int(row['start_peak']-807)
  1229. # end=int(row['end_peak']+807)
  1230. # sequence = genome.fetch(chrom, start, end).upper()
  1231. # sequences.append(sequence)
  1232. # scores, one_hot_encoded_sequences = evaluators[2].calculate_contribution_scores_sequence(sequences, ct_list[a], method='expected_integrated_grad', disable_tqdm=False)
  1233. # scores = np.squeeze(scores, axis=1)
  1234. # scores = np.transpose(scores,(0,2,1))
  1235. # print(scores.shape)
  1236. # one_hot_encoded_sequences = np.transpose(one_hot_encoded_sequences,(0,2,1))
  1237. # np.savez(out_dirs[a]+"/"+ct_list[a]+"_contrib.npz",scores)
  1238. # np.savez(out_dirs[a]+"/"+ct_list[a]+"_oh.npz", one_hot_encoded_sequences)
  1241. # %%
  1242. out_dirs = [
  1243. f'{DATA_DIR}/chip/pax5/',
  1244. f"{DATA_DIR}/chip/ebf1/",
  1245. f"{DATA_DIR}/chip/runx1/",
  1246. f"{DATA_DIR}/chip/gata3_2/",
  1247. f"{DATA_DIR}/chip/ets1/",
  1248. f"{DATA_DIR}/chip/cebpa/",
  1249. f"{DATA_DIR}/chip/spi1_2/",
  1250. ]
  1251. ct_list= ['Bcell',
  1252. 'Bcell',
  1253. 'CD4_Tcell',
  1254. 'CD4_Tcell',
  1255. 'CD4_Tcell',
  1256. 'CD14_monocyte',
  1257. 'CD14_monocyte'
  1258. ]
  1260. recalls={}
  1261. tfs=['PAX5','EBF1','RUNX1','GATA3','ETS1','CEBPA','SPI']
  1262. a=0
  1263. for ct, outdir, sorted_chip_df, tf in zip(ct_list, out_dirs, unibind_dfs, tfs):
  1264. scores = np.load(outdir+ct+"_contrib.npz")['arr_0']
  1265. scores = np.transpose(scores, (0,2,1))
  1266. one_hot_encoded_sequences = np.load(outdir+ct+"_oh.npz")['arr_0']
  1267. one_hot_encoded_sequences = np.transpose(one_hot_encoded_sequences, (0,2,1))
  1269. print(scores.shape)
  1270. print(sorted_chip_df.shape)
  1271. from tqdm import tqdm
  1273. seqlet_starts=[]
  1274. seqlet_ends = []
  1275. p_values=[]
  1276. attributions=[]
  1277. attributions_exact=[]
  1278. overlaps=[]
  1281. for i in tqdm(range(len(sorted_chip_df))):
  1282. oh_seq = one_hot_encoded_sequences[i].T
  1283. X_attr = scores[i].T
  1284. X_attr=X_attr*oh_seq
  1285. X_attr = np.expand_dims(np.sum(X_attr, axis=0),axis=0)
  1287. seqlets = recursive_seqlets(X_attr, threshold=0.05)
  1289. row=sorted_chip_df.iloc[i]
  1291. chrom= row['chrom']
  1292. start=int(row['start_peak']-807)
  1293. end=int(row['end_peak']+807)
  1295. tf_start=int(row['start_chip']-start)
  1296. tf_end = int(row['end_chip']-start)
  1298. attribution = (np.mean(X_attr[0,tf_start:tf_end]))
  1300. overlap_rows = get_row_with_overlap(seqlets, tf_start, tf_end)
  1301. if len(overlap_rows)>0:
  1302. r = overlap_rows.iloc[0]
  1303. seqlet_starts.append(start+int(r['start']))
  1304. seqlet_ends.append(start+int(r['end']))
  1305. p_values.append(r['p-value'])
  1306. attributions.append(r['attribution'])
  1307. overlaps.append(r['overlap'])
  1308. attributions_exact.append(attribution)
  1310. else:
  1311. seqlet_starts.append(np.nan)
  1312. seqlet_ends.append(np.nan)
  1313. p_values.append(np.nan)
  1314. attributions.append(np.nan)
  1315. overlaps.append(np.nan)
  1316. attributions_exact.append(attribution)
  1318. df_final = sorted_chip_df.copy()
  1319. df_final.loc[:,'seqlet_start'] = seqlet_starts
  1320. df_final.loc[:,'seqlet_end'] = seqlet_ends
  1321. df_final.loc[:,'seqlet_p_val']= p_values
  1322. df_final.loc[:,'seqlet_attribution'] = attributions
  1323. df_final.loc[:,'chip_attribution'] = attributions_exact
  1324. df_final.loc[:,'UniBind-Seqlet overlap fraction'] = overlaps
  1327. true_positives = df_final['seqlet_start'].notna().sum()
  1329. # Determine false negatives (NaN in seqlet_start)
  1330. false_negatives = df_final['seqlet_start'].isna().sum()
  1332. # Calculate recall
  1333. recall = true_positives / (true_positives + false_negatives)
  1335. # Output the results
  1336. print(f"True Positives (TP): {true_positives}")
  1337. print(f"False Negatives (FN): {false_negatives}")
  1338. print(f"Recall: {recall:.2f}")
  1339. recalls[tf]={}
  1340. recalls[tf]['recall']=recall
  1341. recalls[tf]['total']=true_positives+false_negatives
  1342. unibind_dfs[a]=df_final
  1343. a+=1
  1345. # %%
  1346. # Sorting by recall values
  1347. sorted_factors = sorted(recalls, key=lambda x: recalls[x]['recall'], reverse=True)
  1348. sorted_recalls = [recalls[factor]['recall'] for factor in sorted_factors]
  1349. sorted_totals = [recalls[factor]['total'] for factor in sorted_factors]
  1351. # Generate new x-axis labels with (n = total)
  1352. x_labels = [f"{factor}\n(n = {total})" for factor, total in zip(sorted_factors, sorted_totals)]
  1354. # Plotting
  1355. fig, ax = plt.subplots(figsize=(8, 5))
  1356. bars = ax.bar(x_labels, sorted_recalls, color='royalblue', edgecolor='black', linewidth=1.2)
  1358. # Aesthetics
  1359. ax.set_ylabel('Recall', fontsize=14)
  1360. ax.set_xlabel('Transcription Factors', fontsize=14)
  1361. ax.set_title('Recall Scores of Identified Seqlets in Unibind Sites', fontsize=16)
  1362. ax.set_ylim(0, 1.05)
  1363. ax.grid(axis='y', linestyle='--', alpha=0.7)
  1365. # Adding value labels
  1366. for bar in bars:
  1367. yval = bar.get_height()
  1368. ax.text(bar.get_x() + bar.get_width()/2, yval + 0.02, f'{yval:.3f}',
  1369. ha='center', va='bottom', fontsize=12, color='black')
  1371. # Remove top and right borders
  1372. ax.spines['top'].set_visible(False)
  1373. ax.spines['right'].set_visible(False)
  1375. # Show the plot
  1376. plt.xticks(rotation=45, fontsize=12, ha='right', va='top')
  1377. plt.yticks(fontsize=12)
  1378. #plt.savefig('paperfigs/unibind_recall_ALL.pdf', bbox_inches='tight')
  1379. plt.show()
  1381. # %%

chip_seq_analysis.ipynb at commit 916dee8, under MIT · at the source

Overview

Authors: Niklas Kempynck1,2,3, Seppe De Winter1,2,3, Casper H Blaauw1,2,3,4, Vasileios Konstantakos1,2,3, Eren Can Ekşi1,2,3, Sam Dieltiens1,2,3, Darina Abaffyová1,2,3, Valérie Bercier2,5, Ibrahim I Taskiran1,2,3,6, Gert Hulselmans1,2,3,7, Katina Spanier1,2,3, Valerie Christiaens1,2,3,7, Ludo Van Den Bosch2,5, Lukas Mahieu1,2,3,7, Stein Aerts1,2,3,7
ORCID iDs: Niklas Kempynck, Seppe De Winter, Casper H Blaauw, Vasileios Konstantakos, Eren Can Ekşi, Sam Dieltiens, Darina Abaffyová, Valérie Bercier, Ibrahim I Taskiran, Katina Spanier, Ludo Van Den Bosch, Lukas Mahieu, Stein Aerts
  1. Laboratory of Computational Biology, VIB Center for AI and Computational Biology (VIB.AI), Leuven, Belgium
  2. VIB-KU Leuven Center for Brain and Disease Research, Leuven, Belgium
  3. Department of Human Genetics, KU Leuven, Leuven, Belgium
  4. Oncode Institute, Hubrecht Institute-KNAW (Royal Netherlands Academy of Arts and Sciences) and University Medical Center Utrecht, Utrecht, the Netherlands
  5. Department of Neurosciences, KU Leuven, Leuven, Belgium
  6. Present Address: Illumina Artificial Intelligence Laboratory, Illumina, Foster City, CA USA
  7. Aligning Science Across Parkinson’s (ASAP) Collaborative Research Network, Chevy Chase, MD USA
Institutions: VIB-KU Leuven Center for Brain & Disease Research (Belgium); VIB.AI (Belgium); KU Leuven (Belgium); Royal Netherlands Academy of Arts and Sciences (Netherlands); University Medical Center Utrecht (Netherlands); Oncode Institute (Netherlands); Illumina (United States) (United States); Aligning Science Across Parkinson's (United States)
Journal: Nature methods, volume 23, issue 5, pages 946-959
Dates: received 2 April 2025; accepted 4 March 2026; published online 2 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41592-026-03057-2 · PMID 41927920 · PMCID PMC13167471 · OpenAlex W7147749522
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), mouse (organism), zebrafish (organism), methods / tools (subfield)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Evoked potentials, Connectivity, Machine learning
Keywords: Software, Epigenomics, Machine learning, Embryogenesis
MeSH: Enhancer Elements, Genetic*, Genomics*, Software*, Animals, Chromatin, Deep Learning, Humans, Mice, Organ Specificity, Species Specificity, Zebrafish (* major topic)
Topic: Genomics and Chromatin Dynamics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: European Research Council (101054387); Fonds Wetenschappelijk Onderzoek (Research Foundation Flanders) (1267625N, 1191323N, 1SH6J24N, G0C1620N, G088523N and G026924N, S005024N, G0I2722N EOS ID 40007513, G094121N, G044124N); Stichting Tegen Kanker (Belgian Foundation Against Cancer) (2024-140, 2020-1396); KU Leuven (Katholieke Universiteit Leuven) (C14/22/132, IDN/22/012 and "Opening the Future" Fund); Fondation Thierry Latran (Thierry Latran Foundation); Association Belge contre les Maladies Neuro-Musculaires (Belgian Association against Neuromuscular Disorders); Muscular Dystrophy Association (Muscular Dystrophy Association Inc.)
Citations: cited by 8 papers (Europe PMC); 101 references in the paper
Research resources: FASTQ files for HepG2 RRID:CVCL_0027, A172 cells RRID:CVCL_0131, LN229 cells RRID:CVCL_0393, MO59J cells RRID:CVCL_0400, GM12878 RRID:CVCL_7526, WiggleTools86 RRID:SCR_001170, RRID:SCR_003030, RRID:SCR_005553, duplicates were removed using Picard RRID:SCR_006525, 95 RRID:SCR_006793, cut-site BigWigs RRID:SCR_007708, peaks were called using MACS87 RRID:SCR_013291, NIS-element RRID:SCR_014329, RRID:SCR_014583, allowing both TensorFlow RRID:SCR_016345, M059J and LN229 using deepTools85 RRID:SCR_016366, RRID:SCR_016368, RRID:SCR_018139, 80 RRID:SCR_018209, PyTorch RRID:SCR_018536, RRID:SCR_022798, by transfer learning from Enformer10 RRID:SCR_024805, 39 RRID:SCR_024811, Keras 3.0 RRID:SCR_026159, epiAneufinder90 v.1.1.3 RRID:SCR_026269, CREsted RRID:SCR_026617, 36 RRID:SCR_026618, Borzoi15 RRID:SCR_026619, RRID:SCR_026620, RRID:SCR_026621, RRID:SCR_026622, pycisTarget34 RRID:SCR_026626, RRID:SCR_026627, RRID:SCR_026702, RRID:SCR_027381, RRID:SCR_027429, RRID:SCR_027436, RRID:SCR_027456, PyTorch Lightning RRID:SCR_027468, HyenaDNA74 RRID:SCR_027471, and Nucleotide Transformer75 RRID:SCR_027472, RRID:SCR_027473

Abstract

Sequence-based deep learning models have become the state of the art for analyzing the genomic regulatory code. Particularly for enhancers, these models excel at deciphering sequence grammar that underlies their activity. To enable end-to-end enhancer modeling and design, we developed a software package called CREsted (cis-regulatory element sequence training, explanation and design). It combines preprocessing and analysis of single-cell assay for transposase-accessible chromatin using sequencing data, modeling chromatin accessibility from sequence, sequence design and downstream analysis to decipher enhancer grammar. We demonstrate CREsted’s functionality on a mouse cortex and a human peripheral blood mononuclear cell dataset. Additionally, we use CREsted to compare mesenchymal-like cancer cell states between tumor types, and we investigate different fine-tuning strategies of genomic foundation models within CREsted. Finally, we train a model on a zebrafish development atlas and use this to design and in vivo validate cell-type-specific enhancers. For varying datasets, we demonstrate that CREsted facilitates efficient training and analyses, enabling scrutinization of the enhancer logic and design of synthetic enhancers across tissues and species.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 40 matches between paragraphs and lines of code.

aertslab/CREsted

License: other
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 3940ea94c8a56344c5478119c5e13755e4eb698e, 26 September 2026
Languages: Python (125), Jupyter (7)
Size: 251 files, 132 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (pyproject.toml), tests, continuous integration, documentation, 7 notebooks
Not found: CITATION.cff
Tools: NumPy (53 files), Keras (43 files), anndata (25 files), Matplotlib (24 files), pandas (14 files), SciPy (12 files), PyTorch (9 files), TensorFlow (8 files), seaborn (5 files), h5py (4 files), Numba (1 file), Pillow (1 file), pysam (1 file), Scanpy (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
134 files
At the source: github.com/aertslab/CREsted

aertslab/create_cisTarget_databases

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 304d5dc1b15e5c923908a50a1ec291c3faaccf9c, 18 June 2025
Languages: Python (27), Shell (1)
Size: 30 files, 28 scripts
Software Heritage: not archived
Found in: the text, “Human PBMC motif enrichment analysis using pycis”
Holds: README, environment (pyproject.toml)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (8 files), pandas (6 files), BEDTools (1 file), Numba (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
29 files
At the source: github.com/aertslab/create_cisTarget_databases.git/

Zenodo 17791463

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file, 0 scripts
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
At the source: zenodo.org/records/17791463

aertslab/CREsted-paper

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 916dee8dc9bcd9908049ce4aa48a0dbfefe5b15a, 3 December 2025
Languages: Jupyter (20), Python (19), Shell (12)
Size: 63 files, 51 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 21 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (16 files), Matplotlib (15 files), pandas (15 files), anndata (12 files), Keras (12 files), seaborn (10 files), SciPy (9 files), scikit-learn (5 files), pysam (3 files), h5py (2 files), PyTorch (2 files), Biopython (1 file), Pillow (1 file), statsmodels (1 file), TensorFlow (1 file), UMAP (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
16 files
At the source: github.com/aertslab/CREsted-paper

Zenodo 17791384

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
At the source:

Zenodo 17202107

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (22 files), pandas (22 files), anndata (14 files), Matplotlib (9 files), SciPy (6 files), seaborn (5 files), Scanpy (4 files), Keras (2 files), CuPy (1 file), h5py (1 file), Numba (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
38 files
At the source:

aertslab/tf-mindi

License: other
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: ca4467f7e020a08fcbbeedffa6f12dc870d78c40, 22 July 2026
Languages: Python (36), Jupyter (3)
Size: 95 files, 39 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, license file, environment (pyproject.toml), tests, continuous integration, documentation, 3 notebooks
Not found: CITATION.cff
Tools: NumPy (25 files), pandas (25 files), anndata (16 files), Matplotlib (9 files), SciPy (6 files), seaborn (5 files), Scanpy (4 files), h5py (2 files), Keras (2 files), CuPy (1 file), Numba (1 file), PyWavelets (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
41 files
At the source: github.com/aertslab/TF-MInDi/tree/1.1.0

Code availability

The CREsted package is available at https://github.com/aertslab/CREsted and https://crested.readthedocs.io and is stored100 at https://zenodo.org/records/15045960. All computational analyses for the main figures can be found in https://github.com/aertslab/CREsted-paper and are stored101 at https://zenodo.org/records/17791384. The following packages were used for data analyses; this information is also available in the key resource table99 (Supplementary Table 2) accompanying this paper (https://zenodo.org/records/17791463): CREsted v.1.4.0, pybigtools v.0.2.0, ChromBPNet v.1.0, anndata v.0.11.3, Enformer, Borzoi, Keras v.3.0, TensorFlow v.2.19.0, PyTorch v.2.6.0, tfmodisco-lite v.2.2.1, tangermeme v.0.4.0, gReLU v.1.0.3, SnapATAC2 v.2.6.4, Harmony v.0.0.10, statsmodels v.0.14.4, Scipy v.1.16, Python Programming Language v.3, NIS-Elements, SCENIC+ v.1.0a, scatac_fragment_tools v.0.1.4, TF-MInDi v.1.0.0, pycisTarget v.1.1, memesuite-lite v.0.2, scanpy v.1.11.4, pyChromVAR v.0.0.4, epiAneufinder v.1.1.3, HyenaDNA, Nucleotide Transformer, Transformers v.4.54.1, PyTorch Lightning v.2.5.2, create_cisTarget_databases, WiggleTools v.1.2.11, MACS v.2 and v.3, pycisTopic v.2 and NGSCheckMate v.1.0.1.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 7 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 251 scripts, each with its path and the digest of its content;
  • 40 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

Analysis data required for reproducing the findings in this paper is available at https://resources.aertslab.org/CREsted/manuscript_data. All CREsted models developed in this paper, and other legacy models, can be loaded through ‘crested.get_model’ or directly from https://resources.aertslab.org/CREsted. All raw and processed sequencing data generated in this study have been deposited in the GEO and are accessible through accession number GSE292617 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE292617). This includes the OmniATAC-seq data of three human GBM cell lines, namely A172, M059J and LN229. Raw imaging data of the zebrafish enhancer reporter assays are available on EBI Biostudies via 10.6019/S-BIAD1962. The key resources table99 (Supplementary Table 2) listing all resources needed to reproduce the results of this paper is available on Zenodo at 10.5281/zenodo.17791463 (ref. 99). The following publicly available datasets were used: Mouse cortex dataset (Zemke et al.42) downloaded from the GEO (GSE229169 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE229169)); Human PBMC dataset (De Rop et al.84) downloaded from the GEO (GSE194028 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE194028)); Zebrafish developmental dataset (Sun et al.76) downloaded from the GEO (GSE243256 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE243256)); Dataset of 100 enhancers tested in zebrafish (supplementary table of Sun et al.76; 10.1038/s41556-024-01449-0); Mouse brain pseudobulk BigWig (Zu et al.19) downloaded from the GEO (GSE246791 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE246791)); PAX5 ChIP–seq peaks downloaded from ENCODE (ENCFF827VVQ (http://www.encodeproject.org/files/ENCFF827VVQ/)); PAX5 ChIP–seq B cells—BigWig track downloaded from ENCODE (ENCFF914QGY (http://www.encodeproject.org/files/ENCFF914QGY/)); EBF1 ChIP–seq B cells—peaks downloaded from ENCODE (ENCFF895MHN (http://www.encodeproject.org/files/ENCFF895MHN/)); EBF1 ChIP–seq B cells—BigWig track downloaded from ENCODE (EBF1 ChIP–seq B cells—BigWig track); POU2F2 ChIP–seq B cells—peaks downloaded from ENCODE (ENCFF934JFA (http://www.encodeproject.org/files/ENCFF934JFA/)); POU2F2 ChIP–seq B cells—BigWig downloaded from ENCODE (ENCFF803HIP (http://www.encodeproject.org/files/ENCFF803HIP/)); GATA3 ChIP–seq CD4+ T cells—peaks and BigWig downloaded from ChIP-Atlas (SRX4705120 (http://www.ncbi.nlm.nih.gov/sra/?term=SRX4705120)); RUNX1 ChIP–seq CD4+ T cells—peaks and BigWig downloaded from ChIP-Atlas (SRX1492212 (http://www.ncbi.nlm.nih.gov/sra/?term=SRX1492212)); ETS1 ChIP–seq CD4+ T cells—peaks and BigWig downloaded from ChIP-Atlas (SRX015825 (https://www.ncbi.nlm.nih.gov/sra/?term=SRX015825)); CEBPA ChIP–seq CD14+ monocytes—peaks and BigWig downloaded from ChIP-Atlas (SRX097095 (https://www.ncbi.nlm.nih.gov/sra/?term=SRX097095)); SPI1 ChIP–seq CD14+ monocytes—peaks and BigWig downloaded from ChIP-Atlas (SRX4001818 (http://www.ncbi.nlm.nih.gov/sra/?term=SRX4001818)); PAX5 ChIP–seq direct targets B cells (https://unibind.uio.no/factor/ENCSR000BHD.GM12878_female_B-cells_lymphoblastoid_cell_line.PAX5/); EBF1 ChIP–seq direct targets B cells (https://unibind.uio.no/factor/ENCSR000DZQ.GM12878_female_B-cells_lymphoblastoid_cell_line.EBF1/); GATA3 ChIP–seq direct targets CD4+ T cells (https://unibind.uio.no/factor/GSE76181.Jurkat_T-cells.GATA3/); RUNX1 ChIP–seq direct targets CD4+ T cells (https://unibind.uio.no/factor/GSE76181.Jurkat_T-cells.RUNX1/); ETS1 ChIP–seq direct targets CD4+ T cells (https://unibind.uio.no/factor/EXP000299.Jurkat_E6_1_T-cells.ETS1); CEBPA ChIP–seq direct targets CD14+ monocytes (https://unibind.uio.no/factor/EXP000946.U937_adult_acute_monocytic_leukemia.CEBPA); SPI1 ChIP–seq direct targets CD14+ monocytes (https://unibind.uio.no/factor/EXP047756.MDMmonocyte_derived_macrophages.SPI1/); Human genome (hg38; https://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips); Mouse genome (mm10; https://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/); Zebrafish genome (danRer11; https://hgdownload.cse.ucsc.edu/goldenPath/danRer11/bigZips).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 15 authors, 4 keywords, 11 MeSH terms, 7 funders, 93 references, 42 RRIDs.

Cite

This paper

Kempynck, N., De Winter, S., Blaauw, C. H., Konstantakos, V., Ekşi, E. C., Dieltiens, S., Abaffyová, D., Bercier, V., Taskiran, I. I., Hulselmans, G., Spanier, K., Christiaens, V., Van Den Bosch, L., Mahieu, L., & Aerts, S. (2026). CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species. Nature methods, 23(5), 946-959. https://doi.org/10.1038/s41592-026-03057-2

BibTeX

@article{kempynck2026crested,
author = {Kempynck, Niklas and De Winter, Seppe and Blaauw, Casper H and Konstantakos, Vasileios and Ekşi, Eren Can and Dieltiens, Sam and Abaffyová, Darina and Bercier, Valérie and Taskiran, Ibrahim I and Hulselmans, Gert and Spanier, Katina and Christiaens, Valerie and Van Den Bosch, Ludo and Mahieu, Lukas and Aerts, Stein},
title = {{CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species}},
journal = {Nature methods},
year = {2026},
month = apr,
volume = {23},
number = {5},
pages = {946--959},
publisher = {Nature Portfolio},
issn = {1548-7091},
doi = {10.1038/s41592-026-03057-2},
url = {https://doi.org/10.1038/s41592-026-03057-2},
pmid = {41927920},
pmcid = {PMC13167471}
}

RIS

TY - JOUR
AU - Kempynck, Niklas
AU - De Winter, Seppe
AU - Blaauw, Casper H
AU - Konstantakos, Vasileios
AU - Ekşi, Eren Can
AU - Dieltiens, Sam
AU - Abaffyová, Darina
AU - Bercier, Valérie
AU - Taskiran, Ibrahim I
AU - Hulselmans, Gert
AU - Spanier, Katina
AU - Christiaens, Valerie
AU - Van Den Bosch, Ludo
AU - Mahieu, Lukas
AU - Aerts, Stein
TI - CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species
T2 - Nature methods
J2 - Nat Methods
PY - 2026
DA - 2026/04/02
VL - 23
IS - 5
SP - 946
EP - 959
SN - 1548-7091
PB - Nature Portfolio
DO - 10.1038/s41592-026-03057-2
UR - https://doi.org/10.1038/s41592-026-03057-2
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41592-026-03057-2",
"type": "article-journal",
"title": "CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species",
"container-title": "Nature methods",
"author": [
{
"family": "Kempynck",
"given": "Niklas"
},
{
"family": "De Winter",
"given": "Seppe"
},
{
"family": "Blaauw",
"given": "Casper H"
},
{
"family": "Konstantakos",
"given": "Vasileios"
},
{
"family": "Ekşi",
"given": "Eren Can"
},
{
"family": "Dieltiens",
"given": "Sam"
},
{
"family": "Abaffyová",
"given": "Darina"
},
{
"family": "Bercier",
"given": "Valérie"
},
{
"family": "Taskiran",
"given": "Ibrahim I"
},
{
"family": "Hulselmans",
"given": "Gert"
},
{
"family": "Spanier",
"given": "Katina"
},
{
"family": "Christiaens",
"given": "Valerie"
},
{
"family": "Van Den Bosch",
"given": "Ludo"
},
{
"family": "Mahieu",
"given": "Lukas"
},
{
"family": "Aerts",
"given": "Stein"
}
],
"container-title-short": "Nat Methods",
"volume": "23",
"issue": "5",
"page": "946-959",
"DOI": "10.1038/s41592-026-03057-2",
"PMID": "41927920",
"PMCID": "PMC13167471",
"ISSN": "1548-7091",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41592-026-03057-2",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
2
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10462-y [code]
SpaDC enables sequence-based integrative analysis and regulatory inference of spatial chromatin accessibility data.
Journal: Communications biology
In common: pysam, Biopython, BEDTools, 8 other tools, genetics / omics, mouse, 12 references
[2] doi:10.1038/s41467-026-75700-7 [code]
Gene regulatory innovations from transposable elements in primate cerebellum development.
Journal: Nature communications
In common: pysam, Biopython, BEDTools, 9 other tools, genetics / omics, mouse, 8 references
[3] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: pysam, BEDTools, UMAP, 10 other tools, genetics / omics, mouse, 6 references
[4] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pysam, Keras, UMAP, 12 other tools, genetics / omics, 3 references
[5] doi:10.1093/bioinformatics/btag652 [code]
mmVelo: a deep generative model for estimating cell state-dependent dynamics across multiple modalities.
Journal: Bioinformatics (Oxford, England)
In common: pysam, BEDTools, UMAP, 11 other tools, genetics / omics, mouse, 3 references
[6] doi:10.1016/j.celrep.2026.117110 [code]
Single-nucleus multiome analysis in the human prefrontal cortex identifies gene expression and cis-regulatory elements associated with aging.
Journal: Cell reports
In common: pysam, BEDTools, anndata, 8 other tools, genetics / omics, 6 references
[7] doi:10.64898/2026.03.30.714220 [code]
An integrated single cell and spatial omics atlas of human prenatal development
Journal: bioRxiv (preprint)
In common: CuPy, Biopython, UMAP, 12 other tools, 2 references
[8] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: CuPy, Biopython, Keras, 12 other tools, zebrafish, mouse
[9] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: pysam, BEDTools, anndata, 6 other tools, genetics / omics, mouse, 7 references
[10] doi:10.1038/s44320-026-00208-7 [code]
Interpretable deep generative ensemble learning for single-cell omics with Hydra.
Journal: Molecular systems biology
In common: pysam, Keras, UMAP, 12 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.