OSCR

Gene regulatory innovations from transposable elements in primate cerebellum development.

Code ↔ Paper

18 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 18 matches
  1. [1] § Methods › Regression analysis of sequence conservation at high attribution sites ↔ 03_HERVL/08_nucleotide_preservation.ipynb, lines 417–510 · score 0.97 · hot encoded nucleotide, absolute attribution score, low attribution positions, high attribution, outcome variable, linear regression
  2. [2] § Results › Lineage-specific co-option of TE subfamilies as regulatory elements ↔ 02_CRE_potential_screening/03_visualize_screening_candidates.ipynb, lines 159–172 · score 0.76 · LTR9B, MLT1M, ERVL B4, LTR1A2, MER72, MER130
  3. [3] § Results › Complex cis-regulatory sequences in ancestral TE sequences facilitate co-option of TEs ↔ 02_CRE_potential_screening/03_visualize_screening_candidates.ipynb, lines 159–172 · score 0.76 · LTR9B, MLT1M, ERVL B4, LTR1A2, MER72, MER130
  4. [4] § Methods › Logistic regression analysis of TE–cCRE overlap ↔ 01_CRE_TE_overlap/06a_regression_analysis_human.ipynb, lines 334–368 · score 0.74 · log_dist_TSS, GC_content, TE_overlap, regression model, mappability
  5. [5] § Methods › Logistic regression analysis of TE–cCRE overlap ↔ 01_CRE_TE_overlap/06b_regression_analysis_marmoset.ipynb, lines 334–368 · score 0.74 · log_dist_TSS, GC_content, TE_overlap, regression model, mappability
  6. [6] § Methods › Overlap between TEs and cCREs in mammalian cerebellum development ↔ 01_CRE_TE_overlap/01_make_repeat_beds.ipynb, lines 142–149 · score 0.70 · nf LO, liftOver, chain, transposable element, UCSC, Dfam
  7. [7] § Results › Complex cis-regulatory sequences in ancestral TE sequences facilitate co-option of TEs ↔ 03_HERVL/01_find_accessible_copies.py, lines 393–431 · score 0.67 · MamRTE1, MamTip2b, MLT1M, ancestral sequence, MER130, copies
  8. [8] § Results › Lineage-specific co-option of TE subfamilies as regulatory elements ↔ 03_HERVL/01_find_accessible_copies.py, lines 433–482 · score 0.65 · MamRTE1, MamTip2b, MLT1M, MER130, TEs, mice
  9. [9] § Methods › Logistic regression analysis of TE–cCRE overlap ↔ 01_CRE_TE_overlap/07_regression_analysis_catlas.ipynb, lines 233–254 · score 0.64 · accessibility_class, phyloP, GC_content, TE_overlap, regression, model
  10. [10] § Results › Determinants of copy-level TE co-option ↔ 03_HERVL/08_nucleotide_preservation.ipynb, lines 417–510 · score 0.58 · high attribution positions, linear regression, HERVL copies, ancestral sequences, motif, scores
  11. [11] § Methods › Logistic regression analysis of TE copy accessibility ↔ 01_CRE_TE_overlap/06a_regression_analysis_human.ipynb, lines 334–368 · score 0.58 · Odds ratios, GC content, confidence intervals, coefficients, TSS, fitted
  12. [12] § Methods › Logistic regression analysis of TE copy accessibility ↔ 01_CRE_TE_overlap/06b_regression_analysis_marmoset.ipynb, lines 334–368 · score 0.58 · Odds ratios, GC content, confidence intervals, coefficients, TSS, fitted
  13. [13] § Methods › Logistic regression analysis of TE–cCRE overlap ↔ 01_CRE_TE_overlap/03a_celltype_specific_peaks_human.ipynb, lines 113–138 · score 0.57 · peakType, n_NMF, intronic, exonic, distal, TE overlap
  14. [14] § Methods › Logistic regression analysis of TE–cCRE overlap ↔ 01_CRE_TE_overlap/06a_regression_analysis_human.ipynb, lines 523–558 · score 0.57 · n_NMF, postnatal, embryonic, adult, fetal, cons
  15. [15] § Methods › Evaluation of the contributions of species-specific TE copies to cross-species gene expression divergence ↔ 03_HERVL/14_HERVL_gene_expression.ipynb, lines 140–154 · score 0.55 · cross species comparisons, gene expression, macaque, marmoset, mouse, cells
  16. [16] § Methods › Mapping chromatin accessibility to the consensus sequences and identifying accessible copies of the screened TE subfamilies ↔ 03_HERVL/01_find_accessible_copies.py, lines 152–163 · score 0.54 · pyBigWig, tracks, copies, accessibility, cell
  17. [17] § Methods › Single-cell multi-omics atlas of cerebellum development ↔ 03_HERVL/11_distance_analysis.ipynb, lines 31–34 · score 0.52 · cross species comparison, highly variable, distance, NMF, peaks, human
  18. [18] § Results › Lineage-specific co-option of TE subfamilies as regulatory elements ↔ 03_HERVL/03_HERVL_accessibilitiy.ipynb, lines 424–438 · score 0.51 · UBC_diff, GC_diff_1, GCP, accessibility, HERVL

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 · 582 lines · 18 KB · CC-BY-NC-SA-4.0 · 3 matches

  1. # %%
  2. import pandas as pd
  3. import numpy as np
  4. import pyranges as pr
  5. import pyBigWig
  6. from Bio import SeqIO
  7. from Bio.SeqUtils import gc_fraction as calc_GC
  8. from scipy import stats
  9. import statsmodels.api as sm
  10. import statsmodels.formula.api as smf
  11. from statsmodels.stats.outliers_influence import variance_inflation_factor
  13. # %%
  14. path_to_lsdf2 = '/work/tetsuya/sds/sd17d003/'
  16. # %% [markdown]
  17. # ## 1. Load data
  19. # %%
  20. # Load CREs
  21. cCREs = pr.read_bed(f"{path_to_lsdf2}/Ioannis/Comparative_Cereb/peak_annotations/Human/human_peaks.bed")
  22. print(f"Loaded {len(cCREs)} CREs")
  24. # Load TEs
  25. TEs = pr.read_bed(f"{path_to_lsdf2}/Tetsuya/data/transposable_elements/dfam/3.8/annotations/hg38/hg38.nrph.hits.te.bed")
  26. print(f"Loaded {len(TEs)} TEs")
  28. # %%
  29. TEs
  31. # %% [markdown]
  32. # ## 2. Create random background regions
  34. # %%
  35. def create_random_regions(n_regions, width, genome_file, exclude_regions=None):
  36. """
  37. Create random genomic regions
  39. Parameters:
  40. n_regions: number of regions to create
  41. width: width of each region (500 for your case)
  42. genome_file: path to chromosome sizes file
  43. exclude_regions: PyRanges object of regions to avoid
  44. """
  45. # Load chromosome sizes
  46. chrom_sizes = pd.read_csv(genome_file, sep='\t', header=None,
  47. names=['Chromosome', 'End'])
  48. chrom_sizes['Start'] = 0
  50. # Filter to main chromosomes
  51. chrom_sizes = chrom_sizes[chrom_sizes['Chromosome'].str.contains('^chr[0-9XY]+$')]
  53. regions = []
  54. attempts = 0
  55. max_attempts = n_regions * 10
  57. while len(regions) < n_regions and attempts < max_attempts:
  58. # Randomly select chromosome weighted by size
  59. weights = chrom_sizes['End'].values / chrom_sizes['End'].sum()
  60. chrom = np.random.choice(chrom_sizes['Chromosome'], p=weights)
  61. chrom_size = chrom_sizes[chrom_sizes['Chromosome'] == chrom]['End'].values[0]
  63. # Random position
  64. start = np.random.randint(0, max(1, chrom_size - width))
  65. end = start + width
  67. # Create region
  68. region = pr.PyRanges(chromosomes=[chrom], starts=[start], ends=[end])
  70. # Check if it overlaps excluded regions
  71. if exclude_regions is not None:
  72. if len(region.overlap(exclude_regions)) == 0:
  73. regions.append({'Chromosome': chrom, 'Start': start, 'End': end})
  74. else:
  75. regions.append({'Chromosome': chrom, 'Start': start, 'End': end})
  77. attempts += 1
  79. background = pr.PyRanges(pd.DataFrame(regions))
  80. return background
  82. # %%
  83. # Create background (5x CREs)
  84. np.random.seed(123)
  85. n_background = len(cCREs)
  87. background = create_random_regions(
  88. n_regions=n_background,
  89. width=500,
  90. genome_file=f'{path_to_lsdf2}/Tetsuya/data/genomes/hg38/hg38.chrom.sizes',
  91. exclude_regions=cCREs
  92. )
  94. print(f"Created {len(background)} background regions")
  96. # %% [markdown]
  97. # ## 3. Calculate TE overlap
  99. # %%
  100. def calculate_overlap(regions, features, min_overlap=250):
  101. """
  102. Calculate binary overlap: does region overlap any feature by at least min_overlap bp?
  104. Parameters:
  105. regions: PyRanges object of genomic regions
  106. features: PyRanges object of features (TEs)
  107. min_overlap: minimum overlap width in bp (default 250)
  109. Returns:
  110. PyRanges object with TE_overlap column added
  111. """
  112. # Get all overlaps
  113. overlaps = regions.join(features, how='left', suffix='_TE')
  115. if len(overlaps) == 0:
  116. # No overlaps at all
  117. regions_df = regions.df.copy()
  118. regions_df['TE_overlap'] = 0
  119. return pr.PyRanges(regions_df)
  121. # Calculate intersection width for each overlap
  122. overlaps_df = overlaps.df.copy()
  124. # Intersection start = max of the two starts
  125. # Intersection end = min of the two ends
  126. overlaps_df['intersection_start'] = overlaps_df[['Start', 'Start_TE']].max(axis=1)
  127. overlaps_df['intersection_end'] = overlaps_df[['End', 'End_TE']].min(axis=1)
  128. overlaps_df['overlap_width'] = overlaps_df['intersection_end'] - overlaps_df['intersection_start']
  130. # Filter to overlaps >= min_overlap
  131. significant_overlaps = overlaps_df[overlaps_df['overlap_width'] > min_overlap]
  133. # Get unique region indices with significant overlaps
  134. if len(significant_overlaps) > 0:
  135. # Need to match back to original regions
  136. # Create unique identifier for regions
  137. overlaps_df['region_id'] = (overlaps_df['Chromosome'].astype(str) + ':' +
  138. overlaps_df['Start'].astype(str) + '-' +
  139. overlaps_df['End'].astype(str))
  140. significant_overlaps['region_id'] = (significant_overlaps['Chromosome'].astype(str) + ':' +
  141. significant_overlaps['Start'].astype(str) + '-' +
  142. significant_overlaps['End'].astype(str))
  144. overlapped_ids = set(significant_overlaps['region_id'].unique())
  145. else:
  146. overlapped_ids = set()
  148. # Add binary overlap to original regions
  149. regions_df = regions.df.copy()
  150. regions_df['region_id'] = (regions_df['Chromosome'].astype(str) + ':' +
  151. regions_df['Start'].astype(str) + '-' +
  152. regions_df['End'].astype(str))
  153. regions_df['TE_overlap'] = regions_df['region_id'].isin(overlapped_ids).astype(int)
  154. regions_df = regions_df.drop('region_id', axis=1)
  156. return pr.PyRanges(regions_df)
  158. # %%
  159. # Calculate for both
  160. cCREs_with_TE = calculate_overlap(cCREs, TEs)
  161. background_with_TE = calculate_overlap(background, TEs)
  163. # Add group labels
  164. cCREs_df = cCREs_with_TE.df
  165. cCREs_df['is_CRE'] = 1
  167. background_df = background_with_TE.df
  168. background_df['is_CRE'] = 0
  170. # Combine
  171. all_regions_df = pd.concat([cCREs_df, background_df], ignore_index=True)
  172. print(f"Total regions: {len(all_regions_df)}")
  173. print(f"TE overlap in CREs: {all_regions_df[all_regions_df['is_CRE']==1]['TE_overlap'].mean():.3f}")
  174. print(f"TE overlap in background: {all_regions_df[all_regions_df['is_CRE']==0]['TE_overlap'].mean():.3f}")
  176. # %%
  177. len(cCREs_df[cCREs_df['TE_overlap']==1])
  179. # %% [markdown]
  180. # ## 4. Calculate GC content
  182. # %%
  183. from pyfaidx import Fasta
  185. def calculate_gc_content(regions_df, genome_fasta):
  186. """Calculate GC content for genomic regions"""
  187. genome = Fasta(genome_fasta)
  188. gc_values = []
  190. for idx, row in regions_df.iterrows():
  191. chrom = row['Chromosome']
  192. start = row['Start']
  193. end = row['End']
  195. try:
  196. seq = genome[chrom][start:end].seq.upper()
  197. gc = (seq.count('G') + seq.count('C')) / len(seq) if len(seq) > 0 else 0
  198. gc_values.append(gc)
  199. except:
  200. gc_values.append(np.nan)
  202. return gc_values
  204. # %%
  205. all_regions_df['GC_content'] = calculate_gc_content(all_regions_df, f'{path_to_lsdf2}/Tetsuya/data/genomes/hg38/hg38.fa')
  206. print(f"GC content calculated, mean: {all_regions_df['GC_content'].mean():.3f}")
  208. # %% [markdown]
  209. # ## 5. Calculate mappability
  211. # %%
  212. def calculate_mappability(regions_df, bigwig_file):
  213. """Calculate mean mappability score for each region"""
  214. bw = pyBigWig.open(bigwig_file)
  215. mappability_scores = []
  217. for idx, row in regions_df.iterrows():
  218. chrom = row['Chromosome']
  219. start = int(row['Start'])
  220. end = int(row['End'])
  222. try:
  223. # Get values for this region
  224. values = bw.values(chrom, start, end)
  225. # Filter out None values and calculate mean
  226. valid_values = [v for v in values if v is not None]
  227. if len(valid_values) > 0:
  228. mappability_scores.append(np.mean(valid_values))
  229. else:
  230. mappability_scores.append(0)
  231. except:
  232. mappability_scores.append(0)
  234. bw.close()
  235. return mappability_scores
  238. # %%
  239. all_regions_df['mappability'] = calculate_mappability(
  240. all_regions_df,
  241. f'{path_to_lsdf2}/Tetsuya/data/mappability/hg38.k50.Umap.MultiTrackMappability.bw'
  242. )
  243. print(f"Mappability calculated, mean: {all_regions_df['mappability'].mean():.3f}")
  245. # %% [markdown]
  246. # ## 6. Calculate distance to nearest TSS
  248. # %%
  249. def extract_tss_from_gtf(gtf_file, feature_type='gene'):
  250. """
  251. Extract TSS positions from GTF file
  253. Parameters:
  254. gtf_file: path to GTF file (e.g., gencode.v38.annotation.gtf.gz)
  255. feature_type: 'transcript' or 'gene' (transcript gives more TSSs)
  257. Returns:
  258. PyRanges object with TSS positions (single bp)
  259. """
  260. # Read GTF
  261. print(f"Reading {gtf_file}...")
  262. gtf = pr.read_gtf(gtf_file)
  264. # Filter to desired feature type
  265. features = gtf[gtf.Feature == feature_type]
  267. print(f"Found {len(features)} {feature_type} features")
  269. # Extract TSS based on strand
  270. # TSS = 5' end of transcript
  271. # For + strand: TSS = Start
  272. # For - strand: TSS = End
  274. df = features.df.copy()
  276. # Create TSS positions
  277. df['TSS'] = np.where(df['Strand'] == '+', df['Start'], df['End'])
  279. # Create single-bp regions for TSS
  280. # PyRanges uses 0-based half-open coordinates [start, end)
  281. # So a single position is represented as [pos, pos+1)
  282. tss_df = pd.DataFrame({
  283. 'Chromosome': df['Chromosome'],
  284. 'Start': df['TSS'],
  285. 'End': df['TSS'] + 1,
  286. 'Strand': df['Strand'],
  287. 'gene_id': df.get('gene_id', ''),
  288. 'gene_name': df.get('gene_name', ''),
  289. 'transcript_id': df.get('transcript_id', '')
  290. })
  292. # Remove duplicates (same TSS position)
  293. tss_df = tss_df.drop_duplicates(subset=['Chromosome', 'Start', 'End'])
  295. print(f"Extracted {len(tss_df)} unique TSS positions")
  297. return pr.PyRanges(tss_df)
  299. def calculate_distance_to_tss(regions, gtf_file):
  300. """
  301. Calculate distance from regions to nearest TSS
  303. Parameters:
  304. regions: PyRanges object of genomic regions
  305. gtf_file: path to GTF annotation file
  307. Returns:
  308. Array of distances to nearest TSS
  309. """
  310. # Extract TSS positions
  311. tss = extract_tss_from_gtf(gtf_file, feature_type='transcript')
  313. # Find nearest TSS for each region
  314. print("Calculating distances to nearest TSS...")
  315. nearest = regions.nearest(tss, overlap=True)
  317. # Extract distances
  318. # PyRanges nearest() returns Distance column
  319. # Negative distances mean overlap (upstream/downstream)
  320. # We want absolute distance
  321. distances = nearest.df['Distance'].abs().values
  323. return distances
  326. # %%
  327. # Usage
  328. gtf_file = f"{path_to_lsdf2}/Ioannis/Comparative_Cereb/genome_annotations/Human/hg38_ens92_filtered.gtf"
  329. all_regions_pr = pr.PyRanges(all_regions_df[['Chromosome', 'Start', 'End']])
  330. all_regions_df['dist_to_TSS'] = calculate_distance_to_tss(all_regions_pr, gtf_file)
  331. print(f"Distance to TSS - median: {np.median(all_regions_df['dist_to_TSS']):.0f} bp")
  332. print(f"Distance to TSS - mean: {np.mean(all_regions_df['dist_to_TSS']):.0f} bp")
  334. # %% [markdown]
  335. # ## 7. Main regression model
  337. # %%
  338. import statsmodels.api as sm
  340. # Remove any rows with missing data
  341. df_clean = all_regions_df.dropna(subset=['TE_overlap', 'is_CRE', 'GC_content', 'mappability', 'dist_to_TSS'])
  342. df_clean['log_dist_TSS'] = np.log(df_clean['dist_to_TSS'] + 1)
  344. # Main model: GC + mappability
  345. formula_main = 'TE_overlap ~ is_CRE + GC_content + mappability + log_dist_TSS'
  346. model_main = smf.logit(formula_main, data=df_clean).fit()
  348. print("\n" + "="*60)
  349. print("MAIN MODEL: TE_overlap ~ is_CRE + GC_content + mappability + log_dist_TSS")
  350. print("="*60)
  351. print(model_main.summary())
  353. # Get odds ratios
  354. odds_ratios = np.exp(model_main.params)
  355. conf_int = np.exp(model_main.conf_int())
  356. conf_int['OR'] = odds_ratios
  357. conf_int.columns = ['2.5%', '97.5%', 'OR']
  359. print("\nOdds Ratios and 95% Confidence Intervals:")
  360. print(conf_int)
  362. # Focus on is_CRE coefficient
  363. or_cre = odds_ratios['is_CRE']
  364. ci_low = conf_int.loc['is_CRE', '2.5%']
  365. ci_high = conf_int.loc['is_CRE', '97.5%']
  366. pval = model_main.pvalues['is_CRE']
  368. print(f"\nCRE effect: OR = {or_cre:.3f}, 95% CI: [{ci_low:.3f}, {ci_high:.3f}], P = {pval:.2e}")
  370. # %%
  371. print(model_main.pvalues)
  373. # %% [markdown]
  374. # ---
  376. # %% [markdown]
  377. # # Within cCRE comparisons:
  379. # %%
  380. human_peaks = pd.read_csv(f'{path_to_lsdf2}/Tetsuya/project/cerebellum/ioannis/Comparative_Cereb/peak_annotations/Human/human_peaks_repeat_overlap.tsv',
  381. sep='\t')
  382. human_peaks['Name'] = human_peaks['peak'].str.replace('hg38_', '')
  384. # %%
  385. human_peaks = human_peaks.merge(all_regions_df[all_regions_df['is_CRE']==1][['Name', 'TE_overlap', 'GC_content', 'mappability', 'dist_to_TSS']],
  386. how='left', on='Name')
  388. # %%
  389. mapping = {
  390. 1: "human_specific",
  391. 3: "primate_specific",
  392. 7: "conserved",
  393. }
  395. human_peaks["Cons_group_label"] = human_peaks["Cons_group"].map(mapping).fillna("others")
  397. # Set "Others" as reference
  398. human_peaks['Cons_group_label'] = pd.Categorical(
  399. human_peaks['Cons_group_label'],
  400. categories=['others', 'human_specific', 'primate_specific', 'conserved'],
  401. ordered=False
  402. )
  404. human_peaks.to_csv(f'{path_to_lsdf2}/Tetsuya/project/cerebellum/ioannis/Comparative_Cereb/peak_annotations/Human/human_peaks_regression.tsv',
  405. sep='\t', index=False)
  407. # %%
  408. human_peaks
  410. # %%
  411. print("\n" + "="*70)
  412. print("MAIN MODEL: Peak Type + Conservation + Covariates")
  413. print("="*70)
  415. formula = 'TE_overlap ~ C(peakType) + C(Cons_group_label) + GC_content'
  416. model_main = smf.logit(formula, data=human_peaks).fit()
  418. print(model_main.summary())
  420. # Get odds ratios
  421. odds_ratios = np.exp(model_main.params)
  422. conf_int = np.exp(model_main.conf_int())
  423. conf_int['OR'] = odds_ratios
  424. conf_int.columns = ['2.5%', '97.5%', 'OR']
  426. print("\n" + "="*70)
  427. print("ODDS RATIOS AND 95% CONFIDENCE INTERVALS")
  428. print("="*70)
  429. print(conf_int)
  431. model_main.pvalues
  433. # %% [markdown]
  434. # ### n_NMF
  436. # %%
  437. human_peaks = pd.read_csv(f'{path_to_lsdf2}/Tetsuya/project/cerebellum/ioannis/Comparative_Cereb/peak_annotations/Human/human_peaks_regression.tsv',
  438. sep='\t')
  440. # %%
  441. nmf_list = [f'NMF_{i}' for i in range(1, 19)]
  442. nmf_size_list = []
  443. for nmf in nmf_list:
  444. nmf_peak_df = pd.read_csv(f'{path_to_lsdf2}/Ioannis/Comparative_Cereb/crossSpecies_comparisons/atac_nmf/allFilteredPeaks_nnls/human/{nmf}.bed',
  445. sep='\t', names=['chrom', 'start', 'end', 'peak'])
  446. print(f'{nmf}: {len(nmf_peak_df)}')
  447. nmf_size_list.append(len(nmf_peak_df))
  448. nmf_peak_df[nmf] = 1
  449. human_peaks = human_peaks.merge(nmf_peak_df[['peak', nmf]], on='peak', how='left')
  450. human_peaks[nmf] = human_peaks[nmf].fillna(0).astype(int)
  451. human_peaks['n_NMF'] = human_peaks[nmf_list].sum(axis=1)
  452. human_peaks
  454. # %%
  455. # Add n_NMF_label column with classification
  456. human_peaks_subset = human_peaks[human_peaks['n_NMF'] > 0].copy()
  458. human_peaks_subset['n_NMF_label'] = human_peaks_subset['n_NMF'].apply(
  459. lambda x: str(int(x)) if x < 10 else '>9'
  460. )
  462. # Set "Others" as reference
  463. human_peaks_subset['Cons_group_label'] = pd.Categorical(
  464. human_peaks_subset['Cons_group_label'],
  465. categories=['others', 'human_specific', 'primate_specific', 'conserved'],
  466. ordered=False
  467. )
  469. # Convert to categorical with '0' as the reference level
  470. human_peaks_subset['n_NMF_label'] = pd.Categorical(
  471. human_peaks_subset['n_NMF_label'],
  472. categories=['1', '2', '3', '4', '5', '6', '7', '8', '9', '>9'],
  473. ordered=True
  474. )
  476. # %%
  477. print("\n" + "="*70)
  478. print("MAIN MODEL: Peak Type + Conservation + Covariates")
  479. print("="*70)
  481. formula = 'TE_overlap ~ C(peakType) + C(Cons_group_label) + C(n_NMF_label) + GC_content'
  482. model_main = smf.logit(formula, data=human_peaks_subset).fit()
  484. print(model_main.summary())
  486. # Get odds ratios
  487. odds_ratios = np.exp(model_main.params)
  488. conf_int = np.exp(model_main.conf_int())
  489. conf_int['OR'] = odds_ratios
  490. conf_int.columns = ['2.5%', '97.5%', 'OR']
  492. print("\n" + "="*70)
  493. print("ODDS RATIOS AND 95% CONFIDENCE INTERVALS")
  494. print("="*70)
  495. print(conf_int)
  497. model_main.pvalues
  499. # %% [markdown]
  500. # ### Highly variable NMF
  502. # %%
  503. nmf_list = [f'NMF_{i}' for i in range(1, 19)]
  504. nmf_size_list = []
  505. for nmf in nmf_list:
  506. nmf_peak_df = pd.read_csv(f'{path_to_lsdf2}/Ioannis/Comparative_Cereb/crossSpecies_comparisons/atac_nmf/highly_variable_peaks/human/{nmf}.bed',
  507. sep='\t', names=['chrom', 'start', 'end', 'peak'])
  508. print(f'{nmf}: {len(nmf_peak_df)}')
  509. nmf_size_list.append(len(nmf_peak_df))
  510. nmf_peak_df[nmf] = 1
  511. human_peaks = human_peaks.merge(nmf_peak_df[['peak', nmf]], on='peak', how='left')
  512. human_peaks[nmf] = human_peaks[nmf].fillna(0).astype(int)
  513. human_peaks['n_NMF'] = human_peaks[nmf_list].sum(axis=1)
  514. human_peaks
  516. # %%
  517. peak_info = pd.read_csv(f'{path_to_lsdf2}/Ioannis/Comparative_Cereb/peak_annotations/Human/human_peaks_info.txt',
  518. sep='\t')
  519. peak_info['maxAcc_devState'] = peak_info['maxAcc_Sample'].str.split(':', expand=True)[1]
  520. peak_info_subset = peak_info[['peak', 'maxAcc_devState']].copy()
  521. peak_info_subset
  523. # %%
  524. human_peaks_subset = human_peaks[human_peaks['n_NMF'] > 0]
  526. dev_state_mapping = {
  527. 'CS18-19': 'embryonic',
  528. 'CS20': 'embryonic',
  529. 'CS22-23': 'embryonic',
  530. '11wpc': 'fetal',
  531. '15-17wpc': 'fetal',
  532. 'newborn': 'postnatal',
  533. 'infant': 'postnatal',
  534. 'toddler': 'postnatal',
  535. 'adult': 'adult',
  536. 'adult_DN': 'adult'
  537. }
  539. # Map developmental states to numeric values
  540. peak_info_subset['dev_state_label'] = peak_info_subset['maxAcc_devState'].map(dev_state_mapping)
  542. human_peaks_subset = human_peaks_subset.merge(peak_info_subset[['peak', 'dev_state_label']], on='peak', how='left')
  544. # Set "Others" as reference
  545. human_peaks_subset['Cons_group_label'] = pd.Categorical(
  546. human_peaks_subset['Cons_group_label'],
  547. categories=['others', 'human_specific', 'primate_specific', 'conserved'],
  548. ordered=False
  549. )
  551. # Convert to categorical with '0' as the reference level
  552. human_peaks_subset['dev_state_label'] = pd.Categorical(
  553. human_peaks_subset['dev_state_label'],
  554. categories=['adult', 'postnatal', 'fetal', 'embryonic'],
  555. ordered=True
  556. )
  558. human_peaks_subset
  560. # %%
  561. print("\n" + "="*70)
  562. print("MAIN MODEL: Peak Type + Conservation + Covariates")
  563. print("="*70)
  565. formula = 'TE_overlap ~ C(peakType) + C(Cons_group_label) + C(dev_state_label) + GC_content'
  566. model_main = smf.logit(formula, data=human_peaks_subset).fit()
  568. print(model_main.summary())
  570. # Get odds ratios
  571. odds_ratios = np.exp(model_main.params)
  572. conf_int = np.exp(model_main.conf_int())
  573. conf_int['OR'] = odds_ratios
  574. conf_int.columns = ['2.5%', '97.5%', 'OR']
  576. print("\n" + "="*70)
  577. print("ODDS RATIOS AND 95% CONFIDENCE INTERVALS")
  578. print("="*70)
  579. print(conf_int)
  581. # %%
  582. model_main.pvalues

06a_regression_analysis_human.ipynb at commit d89c357, under CC-BY-NC-SA-4.0 · at the source

Overview

Authors: Tetsuya Yamada1, Mari Sepp1,2, Ioannis Sarropoulos1,3,4, Henrik Kaessmann1
ORCID iDs: Tetsuya Yamada, Mari Sepp, Ioannis Sarropoulos, Henrik Kaessmann
  1. Center for Molecular Biology of Heidelberg University (ZMBH), DKFZ-ZMBH Alliance, Heidelberg, Germany
  2. Present Address: Centre of Genomics, Evolution and Medicine (cGEM), Institute of Genomics, University of Tartu, Tartu, Estonia
  3. Present Address: Wellcome Sanger Institute, Cambridge, UK
  4. Present Address: Cambridge Stem Cell Institute and Department of Medicine, University of Cambridge, Cambridge, UK
Institutions: Heidelberg University (Germany); DKFZ-ZMBH Alliance (Germany); University of Tartu (Estonia); Wellcome/MRC Cambridge Stem Cell Institute (United Kingdom); University of Cambridge (United Kingdom); Wellcome Sanger Institute (United Kingdom)
Journal: Nature communications, volume 17, issue 1, article 7598
Dates: received 12 September 2025; accepted 8 July 2026; published online 30 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-75700-7 · PMID 42527403 · PMCID PMC13421694 · OpenAlex W7171675596
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), mouse (organism), non-human primate (organism), cellular / molecular (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, fMRI & imaging
Keywords: Gene regulation, Evolutionary genetics, Transcriptomics, Machine learning, Genetics of the nervous system
MeSH: Cerebellum*, DNA Transposable Elements*, Gene Regulatory Networks*, Primates*, Animals, Callithrix, Chromatin, Evolution, Molecular, Gene Expression Regulation, Developmental, Humans, Macaca, Mice (* major topic)
Topic: Chromosomal and Genetic Variations (Plant Science, Agricultural and Biological Sciences), according to OpenAlex
Funding: European Research Council (101019268); European Molecular Biology Organization (EMBO) (Scientific Exchange Grant (9231), Postdoctoral Fellowship (ALTF 769-2022))
Citations: cited by 1 paper (Europe PMC); 124 references in the paper

Abstract

Transposable elements are hypothesized to have driven gene regulatory innovation, yet their contributions to primate brain development at the cell type level remain underexplored. Here, we use single-cell multiomics data from human, macaque, marmoset, and mouse cerebella to show that transposable element contributions to different cell types are shaped by varying degrees of constraints across cell types, as well as the preferential co-option of certain transposable elements in specific cell states. Using a sequence-based deep-learning model that predicts cell-type-specific chromatin accessibility, we systematically assess the co-option potential of transposable elements into cerebellar gene regulatory networks, identifying twelve transposable element subfamilies with complex regulatory sequences in their ancestral states that facilitate their co-option as cell-type-specific cis-regulatory elements. Preservation of these ancestral regulatory sequences, as well as the active chromatin environment surrounding the insertion site, is the major determinant of the accessibility of extant copies. Lineage-specific accessible copies contribute to human-specific gene expression. Broadly, we demonstrate how transposable elements can be flexibly co-opted into cell-type-specific gene regulatory networks, and introduce a generalizable analytical framework for dissecting their contribution to mammalian regulatory evolution.

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 18 matches between paragraphs and lines of code.

kaessmannlab/cerebellum_te

License: CC-BY-NC-SA-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d89c3571724b09827196d19ed0ca22d19a6cd6aa, 10 June 2026
Languages: Jupyter (30), Python (2)
Size: 35 files, 32 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 30 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: pandas (32 files), NumPy (30 files), Matplotlib (27 files), seaborn (27 files), BEDTools (21 files), Biopython (19 files), SciPy (19 files), pysam (13 files), statsmodels (8 files), TensorFlow (2 files), Keras (1 file), scikit-learn (1 file), SHAP (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
34 files
At the source: gitlab.com/kaessmannlab/cerebellum_te

Zenodo 19354896

License: CC-BY-4.0
State: the link answers, verified on 27 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: pandas (22 files), NumPy (20 files), Matplotlib (17 files), seaborn (17 files), BEDTools (16 files), Biopython (13 files), SciPy (10 files), pysam (8 files), statsmodels (6 files), SHAP (1 file), TensorFlow (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
22 files
At the source:

Code availability

All original code is available on GitLab [https://gitlab.com/kaessmannlab/cerebellum_te]. An archived snapshot corresponding to the manuscript is deposited on Zenodo [10.5281/zenodo.19354896]124.

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:

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

Previously published datasets used in this study are available in the ArrayExpress database under accession codes E-MTAB-9765 (https://www.ebi.ac.uk/biostudies/ArrayExpress/studies/E-MTAB-9765) and E-MTAB-10533 (https://www.ebi.ac.uk/biostudies/ArrayExpress/studies/E-MTAB-10533)37, and in the heiData repository under accession codes QDOC4E (https://doi.org/10.11588/data/QDOC4E)38,122 and GDSZG9 (https://doi.org/10.11588/DATA/GDSZG9)39,123. Source data are provided with this paper.

All original code is available on GitLab [https://gitlab.com/kaessmannlab/cerebellum_te]. An archived snapshot corresponding to the manuscript is deposited on Zenodo [10.5281/zenodo.19354896]124.

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

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 5 keywords, 12 MeSH terms, 2 funders, 116 references.

Cite

This paper

Yamada, T., Sepp, M., Sarropoulos, I., & Kaessmann, H. (2026). Gene regulatory innovations from transposable elements in primate cerebellum development. Nature communications, 17(1), 7598. https://doi.org/10.1038/s41467-026-75700-7

BibTeX

@article{yamada2026gene,
author = {Yamada, Tetsuya and Sepp, Mari and Sarropoulos, Ioannis and Kaessmann, Henrik},
title = {{Gene regulatory innovations from transposable elements in primate cerebellum development}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {7598},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75700-7},
url = {https://doi.org/10.1038/s41467-026-75700-7},
pmid = {42527403},
pmcid = {PMC13421694}
}

RIS

TY - JOUR
AU - Yamada, Tetsuya
AU - Sepp, Mari
AU - Sarropoulos, Ioannis
AU - Kaessmann, Henrik
TI - Gene regulatory innovations from transposable elements in primate cerebellum development
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/30
VL - 17
IS - 1
SP - 7598
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75700-7
UR - https://doi.org/10.1038/s41467-026-75700-7
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75700-7",
"type": "article-journal",
"title": "Gene regulatory innovations from transposable elements in primate cerebellum development",
"container-title": "Nature communications",
"author": [
{
"family": "Yamada",
"given": "Tetsuya"
},
{
"family": "Sepp",
"given": "Mari"
},
{
"family": "Sarropoulos",
"given": "Ioannis"
},
{
"family": "Kaessmann",
"given": "Henrik"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "7598",
"DOI": "10.1038/s41467-026-75700-7",
"PMID": "42527403",
"PMCID": "PMC13421694",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75700-7",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
30
]
]
}
}

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/s41592-026-03057-2 [code]
CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: pysam, Biopython, BEDTools, 9 other tools, genetics / omics, mouse, 8 references
[2] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: pysam, BEDTools, TensorFlow, 6 other tools, non-human primate, genetics / omics, mouse, 1 other category, 6 references
[3] doi:10.1186/s13059-026-04050-w
Transposable element-mediated evolutionary expansion of Sox2- and Brn2-binding regulatory modules for mammalian neural-cell differentiation.
Journal: Genome biology
In common: 12 references
[4] 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, SHAP, 6 other tools, genetics / omics, mouse, 3 references
[5] doi:10.1126/sciadv.aed2952 [code]
Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: pysam, Biopython, BEDTools, 5 other tools, cellular / molecular, 4 references
[6] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Biopython, SHAP, Keras, 8 other tools, mouse, cellular / molecular
[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, statsmodels, 5 other tools, genetics / omics, cellular / molecular, 3 references
[8] doi:10.1038/s44400-026-00094-8 [code]
Haplotype-resolved DNA methylation at the &lt;i&gt;APOE&lt;/i&gt; locus identifies allele-specific epigenetic signatures relevant to Alzheimer's disease risk.
Journal: NPJ dementia
In common: pysam, BEDTools, statsmodels, 6 other tools, genetics / omics, cellular / molecular, 2 references
[9] doi:10.1016/j.xhgg.2026.100629 [code]
Positive selection on brain cis-regulatory elements in the human lineage drives changes in gene expression and susceptibility to neuropsychiatric disorders.
Journal: HGG advances
In common: Biopython, statsmodels, seaborn, 5 other tools, cellular / molecular, 3 references
[10] doi:10.1038/s41592-026-03211-w [code]
Spatial isoform sequencing at single-cell resolution reveals cell-type-specific spatial isoform variability in multiple brain cell types.
Journal: Nature methods
In common: pysam, Biopython, BEDTools, 7 other tools, genetics / omics, mouse

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.