OSCR

Dissecting epigenetic heterogeneity in single-cell DNA methylomes with a unified framework.

Code ↔ Paper

12 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 12 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › Identification of hypo-DMRs ↔ scMethCraft/function/DMR.py, the whole file · a weak match · score 0.73 · log2 fold change, hypo DMRs, thresholds, rank, methylation
  2. [2] § Methods › The model architecture of scMethCraft ↔ scMethCraft/model/scmethcraft_model.py, lines 114–199 · score 0.63 · hot encodes, DNA sequence, scMethCraft models, mer, chromosome, position
  3. [3] § Methods › The model architecture of scMethCraft ↔ scMethCraft/model/scmethcraft_model.py, lines 114–199 · score 0.63 · hot encodes, DNA sequence, scMethCraft models, mer, chromosome, position
  4. [4] § Results › Design principles of scMethCraft inspire single-cell omics data modeling ↔ tutorial/tutorial_model_training.ipynb, lines 27–83 · score 0.62 · entire training process, genomic sequences, epoch, sequence feature, encoding, mer
  5. [5] § Results › Design principles of scMethCraft inspire single-cell omics data modeling ↔ scMethCraft/model/scmethcraft_model.py, lines 114–199 · score 0.58 · hot encoding, genomic position, single cell, mer, DNA, sequence
  6. [6] § Results › Design principles of scMethCraft inspire single-cell omics data modeling ↔ scMethCraft/model/scmethcraft_model.py, lines 114–199 · score 0.58 · hot encoding, genomic position, single cell, mer, DNA, sequence
  7. [7] § Methods › Implementation details of downstream analyses ↔ ldscore/regressions.py, lines 538–677 · score 0.57 · score regression, LD scores, genetic, heritability, enrichment, partitioned
  8. [8] § Methods › The model architecture of scMethCraft ↔ scMethCraft/model/scmethcraft_model.py, lines 311–426 · score 0.57 · dense layer, batch normalization, Hidden, GELU, dropout, activation
  9. [9] § Methods › The model architecture of scMethCraft ↔ scMethCraft/model/scmethcraft_model.py, lines 311–426 · score 0.57 · dense layer, batch normalization, Hidden, GELU, dropout, activation
  10. [10] § Methods › Implementation details of downstream analyses ↔ ldsc.py, lines 1–43 · score 0.54 · LD scores, score regression, command, genetic, heritability, partitioned
  11. [11] § Methods › The model architecture of scMethCraft ↔ scMethCraft/model/scmethcraft_model.py, lines 311–426 · score 0.54 · Dense layer, batch normalization, GELU, dropout, activation, model
  12. [12] § Methods › The model architecture of scMethCraft ↔ scMethCraft/model/scmethcraft_model.py, lines 311–426 · score 0.54 · Dense layer, batch normalization, GELU, dropout, activation, model

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

Python · 762 lines · 25 KB · MIT · 4 matches

  1. import math
  2. from typing import List
  3. import sys
  4. import random
  6. import torch
  7. import torch.nn as nn
  8. import numpy as np
  9. import h5py
  10. import pandas as pd
  11. import sklearn.metrics as metrics
  12. import scanpy as sc
  13. import anndata as ad
  14. from scipy.special import expit
  16. from .utils_model import *
  17. from .layers import ConvLayer, DenseLayer
  18. from .utils_model import m_round
  20. device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
  22. def sigmoid(x):
  23. s = 1 / (1 + np.exp(-x))
  24. return s
  26. def _col_round(x):
  27. frac = x - math.floor(x)
  28. if frac <= 0.5:
  29. return math.floor(x)
  30. return math.ceil(x)
  33. def m_get_filter_dim(seq_length: int, pooling_sizes: List[int]):
  34. filter_dim = seq_length
  35. for ps in pooling_sizes:
  36. filter_dim = _col_round(filter_dim / ps)
  37. return int(filter_dim/2)
  42. # =========================
  43. # Chromosome length maps
  44. # =========================
  45. # Used for genomic position normalization in positional encoding.
  46. # Key: chromosome ID
  47. # Value: chromosome length (bp)
  48. p_max_map_human = dict(pd.DataFrame([[1, 2,3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25],[248956422, 242193529, 198295559, 190214555,
  49. 181538259, 170805979, 159345973, 145138636, 138394717,
  50. 133797422, 135086622, 133275309, 114364328, 107043718,
  51. 101991189, 90338345, 83257441, 80373285, 58617616, 64444167,
  52. 46709983, 50818468, 156040895, 57227415,16569]]).T.values)
  54. p_max_map_mouse = dict(pd.DataFrame([['1', '2', '3', '4', '5', '6', '7', '8', '9', '10',
  55. '11', '12', '13', '14', '15', '16', '17', '18', '19','23', '24','25'],[195471971, 182113224, 160039680, 156508116,
  56. 151834684, 149736546, 145441459, 129401213, 124595110,
  57. 130694993, 122082543, 120129022, 120421639, 124902244,
  58. 104043685, 98207768, 94987271, 90702639, 61431566, 171031299, 91744698,10000]]).astype(int).T.values)
  61. # =========================
  62. # Sequence loading utilities
  63. # =========================
  64. def load_seq(input_path,filename,load_range,mode= "onehot"):
  65. """
  66. Load sequence data from HDF5 file.
  68. Args:
  69. input_path (str): Path to input directory
  70. filename (str): HDF5 file name
  71. load_range (tuple or None): Subset range (start, end)
  72. mode (str): "onehot", "kmer", or "both"
  74. Returns:
  75. Tensor(s) depending on mode:
  76. - onehot: (N, 4, L)
  77. - kmer: (N, K)
  78. - both: (onehot, kmer, position)
  79. """
  80. if mode == "onehot":
  81. with h5py.File(input_path+filename) as file:
  82. if load_range != False:
  83. onehot = file["X"][load_range[0]:load_range[1]]
  84. else:
  85. onehot = file["X"]
  86. onehot = torch.nn.functional.one_hot(torch.tensor(onehot).to(torch.int64)).transpose(-2, -1)
  87. return onehot.to(torch.float)
  89. if mode == "kmer":
  90. with h5py.File(input_path+filename) as file:
  91. if load_range != False:
  92. kmer = file["Kmer"][load_range[0]:load_range[1]]
  93. else:
  94. kmer = file["Kmer"]
  95. return torch.tensor(kmer).to(torch.float)
  97. if mode == "both":
  98. with h5py.File(input_path+filename) as file:
  99. if load_range != False:
  100. onehot = file["X"][load_range[0]:load_range[1]]
  101. kmer = file["Kmer"][load_range[0]:load_range[1]]
  102. else:
  103. onehot = file["X"]
  104. kmer = file["Kmer"]
  105. pos = file["Pos"]
  107. onehot = torch.nn.functional.one_hot(torch.tensor(onehot).to(torch.int64)).transpose(-2, -1)
  108. return onehot.to(torch.float),torch.tensor(kmer).to(torch.float),torch.tensor(pos)
  110. # =========================
  111. # Dataset definition
  112. # =========================
  114. class MethyDataset(torch.utils.data.Dataset):
  115. """
  116. Dataset for single-cell methylation modeling.
  118. Each sample consists of:
  119. - DNA sequence (one-hot encoded)
  120. - Cell-level methylation state
  121. - k-mer representation
  122. - Genomic positional embedding
  123. """
  124. def __init__(self, h5_path, state_path, d_emb=64, load_range=None,p_max_map = "HUMAN"):
  125. self.h5_path = h5_path
  126. self.d_emb = d_emb
  127. if p_max_map == "HUMAN":
  128. self.p_max_map = p_max_map_human
  129. elif p_max_map == "MOUSE":
  130. self.p_max_map = p_max_map_mouse
  131. else:
  132. self.p_max_map = p_max_map
  134. self.load_range = load_range
  135. self.state_data = np.load(state_path)
  136. self.cell = self.state_data.shape[1]
  137. with h5py.File(self.h5_path, 'r') as f:
  138. self.total_len = f['X'].shape[0]
  140. def __len__(self):
  141. return self.total_len if self.load_range is None else (self.load_range[1] - self.load_range[0])
  143. def get_pos_embedding(self, raw_pos):
  144. """
  145. Generate sinusoidal-like positional embedding.
  147. Args:
  148. raw_pos: [chromosome_id, start, end]
  150. Returns:
  151. Tensor of shape (d_emb,)
  152. """
  153. chr_id = int(raw_pos[0])
  154. p_max = self.p_max_map.get(chr_id, 248956422) # 防止 key 错误
  156. p = (raw_pos[1] + raw_pos[2]) / 2
  157. t = torch.tensor([p / p_max], dtype=torch.float)
  158. w = 2 * 3.14159 * t
  160. # 计算编码
  161. pe_list = []
  162. for j in range(self.d_emb):
  163. if j <= 1:
  164. pe_list.append(t)
  165. elif j % 2 == 0:
  166. f_val = 1e-4 + (j/2 - 1) * (self.d_emb/2 - 1 - 1e-4) / (self.d_emb/2 - 1)
  167. pe_list.append(torch.cos(f_val * w))
  168. else:
  169. f_val = 1e-4 + ((j-1)/2 - 1) * (self.d_emb/2 - 1 - 1e-4) / (self.d_emb/2 - 1)
  170. pe_list.append(torch.sin(f_val * w))
  172. return torch.cat(pe_list) # 返回 [d_emb] 向量
  174. def __getitem__(self, index):
  175. """
  176. Retrieve a single sample.
  178. Returns:
  179. onehot: (4, L)
  180. state: (n_cells,)
  181. kmer: (K,)
  182. pos_embedding: (d_emb,)
  183. """
  185. if self.load_range:
  186. index += self.load_range[0]
  188. with h5py.File(self.h5_path, 'r') as f:
  189. # 读取原始数据
  190. x = f['X'][index]
  191. kmer = f['Kmer'][index]
  192. pos_info = f['Pos'][index] # 获取 [chr_id, start, end]
  194. onehot = torch.nn.functional.one_hot(torch.tensor(x, dtype=torch.long), num_classes=4).transpose(-2, -1).to(torch.float)
  195. pe_emb = self.get_pos_embedding(pos_info)
  196. state = torch.tensor(self.state_data[index], dtype=torch.float)
  197. kmer_tensor = torch.tensor(kmer, dtype=torch.float)
  199. return onehot, state, kmer_tensor, pe_emb
  201. class SimilarityLayer_fast(nn.Module):
  202. """
  203. Fast approximation of similarity propagation.
  205. Instead of explicitly constructing an NxN similarity matrix,
  206. this layer factorizes it into two low-rank matrices:
  207. A ≈ A1 @ A2^T
  209. This reduces computation from O(N^2) to O(N * dim_n).
  210. """
  213. def __init__(
  214. self,
  215. n_cells: int,
  216. dim_n :int = 256,
  217. dropout_rate: float = 0.1,
  218. batch_norm = True,
  219. device = device,
  220. ):
  222. super().__init__()
  223. self.similarity_matrix=torch.nn.Parameter(torch.rand(n_cells,dim_n))
  224. self.n_cells = n_cells
  225. self.dim_n = dim_n
  226. self.dropout_rate = dropout_rate
  227. self.dropout = self.dropout_rate
  229. self.device = device
  230. self.batch_norm = (
  231. nn.BatchNorm1d(n_cells) if batch_norm else nn.Identity()
  232. )
  237. def forward(
  238. self,
  239. input_vector: torch.Tensor,
  241. ):
  243. A_fixed = torch.abs(self.similarity_matrix)
  244. m1 = (torch.rand(self.n_cells, 1, device=self.device) > self.dropout_rate).float()
  245. m2 = (torch.rand(self.n_cells, 1, device=self.device) > self.dropout_rate).float()
  247. A1 = A_fixed * m1 / (1 - self.dropout_rate)
  248. A2 = A_fixed * m2 / (1 - self.dropout_rate)
  249. XA = torch.matmul(input_vector, A1)
  250. SX = torch.matmul(XA, A2.t())
  251. diag_S = torch.sum(A1 * A2, dim=1)
  252. output_vector = SX - input_vector * diag_S
  253. return self.batch_norm(output_vector)
  255. class SimilarityLayer(nn.Module):
  256. """
  257. Full similarity matrix version (exact but O(N^2)).
  259. Learns a symmetric similarity matrix and applies dropout masking.
  260. """
  262. def __init__(
  263. self,
  264. n_cells: int,
  265. dropout_rate: float = 0.1,
  266. batch_norm = True
  267. ):
  269. super().__init__()
  270. self.similarity_matrix=torch.nn.Parameter(torch.rand(n_cells,n_cells))
  271. self.n_cells = n_cells
  272. self.alpha = 0
  273. self.dropout_rate = dropout_rate
  274. self.dropout = self.dropout_rate
  275. self.eyematrix = 1-torch.eye(self.n_cells,self.n_cells, device=device)
  276. self.batch_norm = (
  277. nn.BatchNorm1d(n_cells) if batch_norm else nn.Identity()
  278. )
  280. def my_dropout(self,input_matrix):
  281. dropout_matrix = (torch.rand(self.n_cells,self.n_cells, device=device)>self.dropout).float()
  282. dropout_matrix = torch.mul(dropout_matrix,self.eyematrix)
  283. input_matrix = torch.mul(dropout_matrix, input_matrix)
  284. input_matrix = input_matrix/(1-self.dropout)
  285. if self.alpha>0:
  286. input_matrix = input_matrix+self.alpha*self.n_cells*torch.mean(input_matrix)*torch.eye(self.n_cells,self.n_cells, device=device)
  287. return input_matrix
  290. def forward(
  291. self,
  292. input_vector: torch.Tensor,
  294. ):
  296. fixed_similarity_matrix = self.similarity_matrix
  297. fixed_similarity_matrix = torch.abs(fixed_similarity_matrix+fixed_similarity_matrix.T)/2
  298. fixed_similarity_matrix = self.my_dropout(fixed_similarity_matrix)
  299. output_vector = torch.matmul(input_vector,fixed_similarity_matrix)/(1+self.alpha)
  300. output_vector = self.batch_norm(output_vector)
  301. return output_vector
  303. class Sequence_extraction(nn.Module):
  304. """
  305. Core feature extraction module combining:
  306. - sequence (one-hot CNN)
  307. - k-mer representation
  308. - positional embedding
  309. """
  311. def __init__(
  312. self,
  313. n_cells: int,
  314. K:int = 8,
  315. n_filters_init: int = 256,
  316. n_repeat_blocks_tower: int =2,
  317. filters_mult: float = 1.41421,
  318. n_filters_pre_bottleneck: int = 256,
  319. n_bottleneck_layer: int = 25,
  320. dropout_rate_similarity: float = 0.3,
  321. batch_norm: bool = True,
  322. embedding_dim: int = 16,
  323. dropout: float = 0.0,
  324. genomic_seq_length: int = 10000, ):
  325. super().__init__()
  327. self.stem = ConvLayer(
  328. in_channels=4,
  329. out_channels=n_filters_init,
  330. kernel_size=12,
  331. pool_size=4,
  332. dropout=dropout,
  333. batch_norm=batch_norm,
  334. )
  338. tower_layers = []
  339. curr_n_filters = n_filters_init
  341. for i in range(n_repeat_blocks_tower):
  342. tower_layers.append(
  343. ConvLayer(
  344. in_channels=curr_n_filters,
  345. out_channels=m_round(curr_n_filters * filters_mult),
  346. kernel_size=5,
  347. pool_size=2,
  348. dropout=dropout,
  349. batch_norm=batch_norm,
  350. )
  351. )
  352. curr_n_filters = m_round(curr_n_filters * filters_mult)
  354. self.tower = nn.Sequential(*tower_layers)
  356. self.pre_bottleneck = ConvLayer(
  357. in_channels=curr_n_filters,
  358. out_channels=n_filters_pre_bottleneck,
  359. kernel_size=1,
  360. dropout=dropout,
  361. batch_norm=batch_norm,
  362. pool_size=2,
  363. )
  365. # get pooling sizes of the upstream conv layers
  366. pooling_sizes = [4] + [2] * n_repeat_blocks_tower + [1]
  367. # get filter dimensionality to account for variable sequence length
  368. filter_dim = m_get_filter_dim(seq_length=genomic_seq_length, pooling_sizes=pooling_sizes)
  370. self.hidden_onehot = DenseLayer(
  371. in_features=n_filters_pre_bottleneck * filter_dim,
  372. out_features=n_bottleneck_layer,
  373. use_bias=True,
  374. batch_norm=True,
  375. dropout=0.2,
  376. activation_fn=nn.Identity(),
  377. )
  379. self.hidden_pos_1 = KANLinear(
  380. in_features=64,
  381. out_features=32,
  382. )
  384. self.hidden_pos_2 = KANLinear(
  385. in_features=32,
  386. out_features=n_bottleneck_layer*2,
  387. )
  388. self.Kmer1 = DenseLayer(
  389. in_features=4 ** K,
  390. out_features=4 ** int(K/2),
  391. use_bias=True,
  392. batch_norm=True,
  393. dropout=0.2,
  394. activation_fn=nn.GELU(),
  395. )
  396. self.pos_embedding_attention = nn.Embedding(256,embedding_dim)
  397. self.hidden_Kmer = DenseLayer(
  398. in_features=4**int(K/2),
  399. out_features=n_bottleneck_layer,
  400. use_bias=True,
  401. batch_norm=True,
  402. dropout=0.2,
  403. activation_fn=nn.GELU(),
  404. )
  406. self.kmer_embedding_linear = DenseLayer(
  407. in_features=1,
  408. out_features=embedding_dim,
  409. use_bias=True,
  410. batch_norm=False,
  411. dropout=0,
  412. activation_fn=nn.GELU(),
  413. )
  415. self.transformer = torch.nn.TransformerEncoderLayer(embedding_dim,4,batch_first=True)
  417. self.kmer_embedding_linear_inverse = DenseLayer(
  418. in_features=embedding_dim,
  419. out_features=1,
  420. use_bias=True,
  421. batch_norm=False,
  422. dropout=0,
  423. activation_fn=nn.GELU(),
  424. )
  426. self.final = nn.Linear(n_bottleneck_layer*2, n_cells)
  428. def forward(
  429. self,
  430. onehot: torch.Tensor,
  431. kmer: torch.Tensor,
  432. pos:torch.Tensor
  433. ):
  436. onehot = self.stem(onehot)
  437. onehot = self.tower(onehot)
  438. onehot = self.pre_bottleneck(onehot)
  439. onehot = onehot.view(onehot.shape[0], -1)
  440. onehot = self.hidden_onehot(onehot)
  442. kmer = self.Kmer1(kmer)
  443. res_kmer = kmer
  444. kmer = self.kmer_embedding_linear(kmer.unsqueeze(2))
  445. kmer_seq_embedding = kmer + self.pos_embedding_attention.weight
  446. kmer_seq_embedding = self.transformer(kmer_seq_embedding)
  447. kmer = self.kmer_embedding_linear_inverse(kmer_seq_embedding).squeeze(2)
  448. kmer = kmer + res_kmer
  451. kmer = self.hidden_Kmer(kmer)
  453. pos = self.hidden_pos_1(pos)
  454. pos = self.hidden_pos_2(pos)
  456. latent = torch.cat((onehot,kmer),1)
  458. res_latent = latent
  459. latent = torch.mul(latent, pos)
  460. latent = latent + res_latent
  462. latent = self.final(latent)
  465. return latent
  467. class Similarity_weighting(nn.Module):
  468. def __init__(
  469. self,
  470. n_cells: int,
  471. dropout_rate: float = 0.1,
  472. batch_norm = True
  473. ):
  474. super().__init__()
  475. self.SimilarityLayer1 = SimilarityLayer(n_cells,dropout_rate,batch_norm)
  477. def train(self):
  478. self.SimilarityLayer1.train()
  479. self.SimilarityLayer1.dropout = self.SimilarityLayer1.dropout_rate
  480. self.SimilarityLayer1.alpha = 0
  482. def eval(self):
  483. self.SimilarityLayer1.eval()
  484. self.SimilarityLayer1.dropout = 0
  485. self.SimilarityLayer1.alpha = 0.2
  487. def forward(self, methy_level_vector: torch.Tensor):
  488. methy_level_vector = self.SimilarityLayer1(methy_level_vector)
  489. return methy_level_vector
  491. class Similarity_weighting_fast(nn.Module):
  492. def __init__(
  493. self,
  494. n_cells: int,
  495. dim_n :int = 256,
  496. dropout_rate: float = 0.1,
  497. batch_norm = True,
  498. device = device
  499. ):
  500. super().__init__()
  502. self.SimilarityLayer1 = SimilarityLayer_fast(n_cells,dim_n,dropout_rate,batch_norm,device = device)
  504. def train(self):
  505. self.SimilarityLayer1.dropout = self.SimilarityLayer1.dropout_rate
  507. def eval(self):
  508. self.SimilarityLayer1.dropout = 0
  511. def forward(self, methy_level_vector: torch.Tensor):
  513. methy_level_vector = self.SimilarityLayer1(methy_level_vector)
  515. return methy_level_vector
  517. class KANLinear(torch.nn.Module):
  518. def __init__(
  519. self,
  520. in_features,
  521. out_features,
  522. grid_size=5,
  523. spline_order=3,
  524. scale_noise=0.1,
  525. scale_base=1.0,
  526. scale_spline=1.0,
  527. enable_standalone_scale_spline=True,
  528. base_activation=torch.nn.SiLU,
  529. grid_eps=0.02,
  530. grid_range=[-1, 1],
  531. ):
  532. super(KANLinear, self).__init__()
  533. self.in_features = in_features
  534. self.out_features = out_features
  535. self.grid_size = grid_size
  536. self.spline_order = spline_order
  538. h = (grid_range[1] - grid_range[0]) / grid_size
  539. grid = (
  540. (
  541. torch.arange(-spline_order, grid_size + spline_order + 1) * h
  542. + grid_range[0]
  543. )
  544. .expand(in_features, -1)
  545. .contiguous()
  546. )
  547. self.register_buffer("grid", grid)
  549. self.base_weight = torch.nn.Parameter(torch.Tensor(out_features, in_features))
  550. self.spline_weight = torch.nn.Parameter(
  551. torch.Tensor(out_features, in_features, grid_size + spline_order)
  552. )
  553. if enable_standalone_scale_spline:
  554. self.spline_scaler = torch.nn.Parameter(
  555. torch.Tensor(out_features, in_features)
  556. )
  558. self.scale_noise = scale_noise
  559. self.scale_base = scale_base
  560. self.scale_spline = scale_spline
  561. self.enable_standalone_scale_spline = enable_standalone_scale_spline
  562. self.base_activation = base_activation()
  563. self.grid_eps = grid_eps
  565. self.reset_parameters()
  567. def reset_parameters(self):
  568. torch.nn.init.kaiming_uniform_(self.base_weight, a=math.sqrt(5) * self.scale_base)
  569. with torch.no_grad():
  570. noise = (
  571. (
  572. torch.rand(self.grid_size + 1, self.in_features, self.out_features)
  573. - 1 / 2
  574. )
  575. * self.scale_noise
  576. / self.grid_size
  577. )
  578. self.spline_weight.data.copy_(
  579. (self.scale_spline if not self.enable_standalone_scale_spline else 1.0)
  580. * self.curve2coeff(
  581. self.grid.T[self.spline_order : -self.spline_order],
  582. noise,
  583. )
  584. )
  585. if self.enable_standalone_scale_spline:
  586. # torch.nn.init.constant_(self.spline_scaler, self.scale_spline)
  587. torch.nn.init.kaiming_uniform_(self.spline_scaler, a=math.sqrt(5) * self.scale_spline)
  589. def b_splines(self, x: torch.Tensor):
  590. """
  591. Compute the B-spline bases for the given input tensor.
  593. Args:
  594. x (torch.Tensor): Input tensor of shape (batch_size, in_features).
  596. Returns:
  597. torch.Tensor: B-spline bases tensor of shape (batch_size, in_features, grid_size + spline_order).
  598. """
  599. assert x.dim() == 2 and x.size(1) == self.in_features
  601. grid: torch.Tensor = (
  602. self.grid
  603. ) # (in_features, grid_size + 2 * spline_order + 1)
  604. x = x.unsqueeze(-1)
  605. bases = ((x >= grid[:, :-1]) & (x < grid[:, 1:])).to(x.dtype)
  606. for k in range(1, self.spline_order + 1):
  607. bases = (
  608. (x - grid[:, : -(k + 1)])
  609. / (grid[:, k:-1] - grid[:, : -(k + 1)])
  610. * bases[:, :, :-1]
  611. ) + (
  612. (grid[:, k + 1 :] - x)
  613. / (grid[:, k + 1 :] - grid[:, 1:(-k)])
  614. * bases[:, :, 1:]
  615. )
  617. assert bases.size() == (
  618. x.size(0),
  619. self.in_features,
  620. self.grid_size + self.spline_order,
  621. )
  622. return bases.contiguous()
  624. def curve2coeff(self, x: torch.Tensor, y: torch.Tensor):
  625. """
  626. Compute the coefficients of the curve that interpolates the given points.
  628. Args:
  629. x (torch.Tensor): Input tensor of shape (batch_size, in_features).
  630. y (torch.Tensor): Output tensor of shape (batch_size, in_features, out_features).
  632. Returns:
  633. torch.Tensor: Coefficients tensor of shape (out_features, in_features, grid_size + spline_order).
  634. """
  635. assert x.dim() == 2 and x.size(1) == self.in_features
  636. assert y.size() == (x.size(0), self.in_features, self.out_features)
  638. A = self.b_splines(x).transpose(
  639. 0, 1
  640. ) # (in_features, batch_size, grid_size + spline_order)
  641. B = y.transpose(0, 1) # (in_features, batch_size, out_features)
  642. solution = torch.linalg.lstsq(
  643. A, B
  644. ).solution # (in_features, grid_size + spline_order, out_features)
  645. result = solution.permute(
  646. 2, 0, 1
  647. ) # (out_features, in_features, grid_size + spline_order)
  649. assert result.size() == (
  650. self.out_features,
  651. self.in_features,
  652. self.grid_size + self.spline_order,
  653. )
  654. return result.contiguous()
  656. @property
  657. def scaled_spline_weight(self):
  658. return self.spline_weight * (
  659. self.spline_scaler.unsqueeze(-1)
  660. if self.enable_standalone_scale_spline
  661. else 1.0
  662. )
  664. def forward(self, x: torch.Tensor):
  665. assert x.dim() == 2 and x.size(1) == self.in_features
  667. base_output = torch.nn.functional.linear(self.base_activation(x), self.base_weight)
  668. spline_output = torch.nn.functional.linear(
  669. self.b_splines(x).view(x.size(0), -1),
  670. self.scaled_spline_weight.view(self.out_features, -1),
  671. )
  672. return base_output + spline_output
  674. @torch.no_grad()
  675. def update_grid(self, x: torch.Tensor, margin=0.01):
  676. assert x.dim() == 2 and x.size(1) == self.in_features
  677. batch = x.size(0)
  679. splines = self.b_splines(x) # (batch, in, coeff)
  680. splines = splines.permute(1, 0, 2) # (in, batch, coeff)
  681. orig_coeff = self.scaled_spline_weight # (out, in, coeff)
  682. orig_coeff = orig_coeff.permute(1, 2, 0) # (in, coeff, out)
  683. unreduced_spline_output = torch.bmm(splines, orig_coeff) # (in, batch, out)
  684. unreduced_spline_output = unreduced_spline_output.permute(
  685. 1, 0, 2
  686. ) # (batch, in, out)
  688. # sort each channel individually to collect data distribution
  689. x_sorted = torch.sort(x, dim=0)[0]
  690. grid_adaptive = x_sorted[
  691. torch.linspace(
  692. 0, batch - 1, self.grid_size + 1, dtype=torch.int64, device=x.device
  693. )
  694. ]
  696. uniform_step = (x_sorted[-1] - x_sorted[0] + 2 * margin) / self.grid_size
  697. grid_uniform = (
  698. torch.arange(
  699. self.grid_size + 1, dtype=torch.float32, device=x.device
  700. ).unsqueeze(1)
  701. * uniform_step
  702. + x_sorted[0]
  703. - margin
  704. )
  706. grid = self.grid_eps * grid_uniform + (1 - self.grid_eps) * grid_adaptive
  707. grid = torch.concatenate(
  708. [
  709. grid[:1]
  710. - uniform_step
  711. * torch.arange(self.spline_order, 0, -1, device=x.device).unsqueeze(1),
  712. grid,
  713. grid[-1:]
  714. + uniform_step
  715. * torch.arange(1, self.spline_order + 1, device=x.device).unsqueeze(1),
  716. ],
  717. dim=0,
  718. )
  720. self.grid.copy_(grid.T)
  721. self.spline_weight.data.copy_(self.curve2coeff(x, unreduced_spline_output))
  723. def regularization_loss(self, regularize_activation=1.0, regularize_entropy=1.0):
  724. """
  725. Compute the regularization loss.
  727. This is a dumb simulation of the original L1 regularization as stated in the
  728. paper, since the original one requires computing absolutes and entropy from the
  729. expanded (batch, in_features, out_features) intermediate tensor, which is hidden
  730. behind the F.linear function if we want an memory efficient implementation.
  732. The L1 regularization is now computed as mean absolute value of the spline
  733. weights. The authors implementation also includes this term in addition to the
  734. sample-based regularization.
  735. """
  736. l1_fake = self.spline_weight.abs().mean(-1)
  737. regularization_loss_activation = l1_fake.sum()
  738. p = l1_fake / regularization_loss_activation
  739. regularization_loss_entropy = -torch.sum(p * p.log())
  740. return (
  741. regularize_activation * regularization_loss_activation
  742. + regularize_entropy * regularization_loss_entropy
  743. )
  746. def output_model(MethyBasset_part1,MethyBasset_part2,savepath = f"../sample_data/output/"):
  748. import os
  749. if not os.path.exists(savepath):
  750. os.makedirs(savepath)
  751. print("Folder created")
  752. else:
  753. print("Folder already exists")
  756. torch.save(MethyBasset_part1.state_dict(), f"{savepath}/scMethCraft_part1.pth")
  757. torch.save(MethyBasset_part2.state_dict(), f"{savepath}/scMethCraft_part2.pth")

scmethcraft_model.py at commit 6c2ef56, under MIT · at the source

Overview

Authors: Songming Tang1, Siyu Li1, Guangxin Zhang2, Aoran Lyu3, Han Li1,4, Shengquan Chen1,5
ORCID iDs: Han Li, Shengquan Chen
  1. School of Mathematical Sciences and LPMC, Nankai University,Tianjin, China
  2. SeekGene BioSciences, Beijing, China
  3. Department of Mechanical and Aerospace Engineering, The University of Manchester,Manchester, UK
  4. Sheffield Institute for Translational Neuroscience, University of Sheffield,Sheffield, UK
  5. Academy for Advanced Interdisciplinary Studies, Nankai University,Tianjin, China
Institutions: Nankai University (China); University of Manchester (United Kingdom); University of Sheffield (United Kingdom)
Journal: Nature communications, volume 17, issue 1, article 6469
Dates: received 12 August 2025; accepted 30 April 2026; published online 15 May 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-73171-4 · PMID 42140984 · PMCID PMC13376737 · OpenAlex W7161287819
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity, Machine learning
Keywords: Computational models, Data mining, Machine learning, Genome informatics, Epigenetics
MeSH: DNA Methylation*, Epigenesis, Genetic*, Epigenome*, Epigenomics*, Single-Cell Analysis*, Animals, Genetic Heterogeneity, Humans (* major topic)
Topic: Epigenetics and DNA Methylation (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: National Natural Science Foundation of China (National Science Foundation of China) (62473212, 62203236)
Citations: cited by 2 papers (Europe PMC); 56 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

bulik/ldsc

License: GPL-3.0
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 2fdeeb3b44379408794154993dbd6101b8946b7e, 16 January 2026
Languages: Python (19), MATLAB (4), Perl (1), R (1)
Size: 1,093 files, 25 scripts
Software Heritage: archived
Found in: the text, “Implementation details of downstream analyses”
Holds: README, license file, environment (environment.yml, requirements.txt, setup.py), tests
Not found: CITATION.cff, continuous integration, documentation
Tools: NumPy (18 files), pandas (11 files), SciPy (5 files), BEDTools (2 files)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
27 files
At the source: github.com/bulik/ldsc/wiki/LD-Score-Estimation-Tutorial

BioX-NKU/scMethCraft

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 6c2ef5613c1edd87036d1a982d2ffea118955f05, 10 September 2026
Languages: Python (21), Jupyter (16)
Size: 71 files, 37 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (setup.py), 8 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: Scanpy (17 files), NumPy (12 files), PyTorch (11 files), SciPy (9 files), pandas (8 files), anndata (5 files), scikit-learn (5 files), h5py (4 files), pysam (3 files), seaborn (2 files), Matplotlib (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
39 files
At the source: github.com/BioX-NKU/scMethCraft

Zenodo 19549979

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 (11 files), Scanpy (11 files), PyTorch (9 files), pandas (7 files), SciPy (6 files), anndata (4 files), h5py (4 files), scikit-learn (4 files), pysam (3 files), seaborn (2 files), Matplotlib (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
31 files
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-73171-4.

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:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 91 scripts, each with its path and the digest of its content;
  • 12 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 statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-73171-4.

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 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 5 keywords, 8 MeSH terms, 1 funder, 55 references.

Cite

This paper

Tang, S., Li, S., Zhang, G., Lyu, A., Li, H., & Chen, S. (2026). Dissecting epigenetic heterogeneity in single-cell DNA methylomes with a unified framework. Nature communications, 17(1), 6469. https://doi.org/10.1038/s41467-026-73171-4

BibTeX

@article{tang2026dissecting,
author = {Tang, Songming and Li, Siyu and Zhang, Guangxin and Lyu, Aoran and Li, Han and Chen, Shengquan},
title = {{Dissecting epigenetic heterogeneity in single-cell DNA methylomes with a unified framework}},
journal = {Nature communications},
year = {2026},
month = may,
volume = {17},
number = {1},
pages = {6469},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-73171-4},
url = {https://doi.org/10.1038/s41467-026-73171-4},
pmid = {42140984},
pmcid = {PMC13376737}
}

RIS

TY - JOUR
AU - Tang, Songming
AU - Li, Siyu
AU - Zhang, Guangxin
AU - Lyu, Aoran
AU - Li, Han
AU - Chen, Shengquan
TI - Dissecting epigenetic heterogeneity in single-cell DNA methylomes with a unified framework
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/05/15
VL - 17
IS - 1
SP - 6469
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-73171-4
UR - https://doi.org/10.1038/s41467-026-73171-4
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-73171-4",
"type": "article-journal",
"title": "Dissecting epigenetic heterogeneity in single-cell DNA methylomes with a unified framework",
"container-title": "Nature communications",
"author": [
{
"family": "Tang",
"given": "Songming"
},
{
"family": "Li",
"given": "Siyu"
},
{
"family": "Zhang",
"given": "Guangxin"
},
{
"family": "Lyu",
"given": "Aoran"
},
{
"family": "Li",
"given": "Han"
},
{
"family": "Chen",
"given": "Shengquan"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "6469",
"DOI": "10.1038/s41467-026-73171-4",
"PMID": "42140984",
"PMCID": "PMC13376737",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-73171-4",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
15
]
]
}
}

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.1002/advs.77524 [code]
MethyAnno: An Interpretable Automated Annotation Method Leveraging Multi-Scale Information and Metric Learning Framework for scDNAm Data.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: anndata, Scanpy, PyTorch, 6 other tools, genetics / omics, 10 references, author Shengquan Chen
[2] doi:10.1101/gr.281350.125 [code]
High-fidelity bidirectional translation between single-cell transcriptomes and DNA methylomes with scBOND.
Journal: Genome research
In common: BEDTools, anndata, Scanpy, 8 other tools, genetics / omics, cellular / molecular, 7 references
[3] 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, BEDTools, anndata, 7 other tools, genetics / omics, 7 references
[4] doi:10.1038/s41592-026-03057-2 [code]
CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: pysam, BEDTools, anndata, 9 other tools, genetics / omics, 3 references
[5] doi:10.1038/s41467-026-71803-3 [code]
Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: BEDTools, anndata, Scanpy, 8 other tools, genetics / omics, cellular / molecular, 2 references
[6] 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, anndata, 8 other tools, genetics / omics, 1 reference
[7] 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, 7 other tools, genetics / omics, cellular / molecular, 1 reference
[8] 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, anndata, Scanpy, 8 other tools, cellular / molecular, 1 reference
[9] doi:10.1093/bib/bbag485 [code]
Single-cell-level perturbation-induced and condition-related signal estimation with batch effect removal using NDreamer.
Journal: Briefings in bioinformatics
In common: anndata, Scanpy, PyTorch, 6 other tools, genetics / omics, 4 references
[10] 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, cellular / molecular, 2 references

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.