OSCR

Sigma-1 and Sigma-2 receptors exhibit divergent genome-wide Co-expression architectures in human brain despite shared subcellular localization.

Code ↔ Paper

15 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 15 matches
  1. [1] § Materials and methods › Data source and processing ↔ MS3_Sigma_Divergence_Pipeline.py, lines 85–126 · score 0.85 · Anterior cingulate cortex, median TPM, Brain Frontal Cortex, GTEx v8, Nucleus accumbens, BA24
  2. [2] § Materials and methods › Data source and processing ↔ MS3_sigma_divergence_wj_pipeline.py, lines 91–101 · score 0.80 · Anterior cingulate cortex, median TPM, Brain Frontal Cortex, Nucleus accumbens, BA24, Putamen
  3. [3] § Results › Divergent tails encode distinct functional programs ↔ MS3_sigma_divergence_wj_pipeline.py, lines 916–975 · score 0.73 · TMEM97 unique GO, SIGMAR1 unique GO, TMEM97 Reactome, SIGMAR1 Reactome, SIGMAR1 KEGG, log10
  4. [4] § Results › SIGMAR1 and TMEM97 share global architecture but diverge at functional tails ↔ MS3_sigma_divergence_wj_pipeline.py, lines 1–67 · score 0.69 · genome wide Spearman, Genome wide co, co expression network, expression architecture, binary Jaccard, Weighted Jaccard
  5. [5] § Materials and methods › Genome-wide Co-expression analysis ↔ MS3_Sigma_Divergence_Pipeline.py, lines 85–126 · score 0.66 · ER stress, EIF2S1, target gene, NEMF, HSPA5, PELO
  6. [6] § Materials and methods › Functional enrichment ↔ MS3_Sigma_Divergence_Pipeline.py, lines 1–40 · score 0.65 · gProfiler, TMEM97 top, Custom gene, CC, MF, Reactome
  7. [7] § Materials and methods › Functional enrichment ↔ MS3_sigma_divergence_wj_pipeline.py, lines 653–725 · score 0.62 · gProfiler, Custom gene, TMEM97 unique, SIGMAR1 unique, CC, MF
  8. [8] § Results › Multi-region replication ↔ MS3_sigma_divergence_wj_pipeline.py, lines 91–101 · score 0.61 · NAcc, anterior cingulate cortex, nucleus accumbens, brain regions, BA24, putamen
  9. [9] § Materials and methods › Primary analysis: Weighted Jaccard on continuous vectors ↔ MS3_sigma_divergence_wj_pipeline.py, lines 154–173 · score 0.60 · r_shifted, denominator, numerator, mapped, max, WJ
  10. [10] § Materials and methods › Sensitivity analyses ↔ MS3_sigma_divergence_wj_pipeline.py, lines 1–67 · score 0.58 · covariate adjustment, sensitivity, brain regions, sex, age, deconvolution
  11. [11] § Results › Custom gene set enrichment ↔ MS3_sigma_divergence_wj_pipeline.py, lines 653–725 · score 0.57 · ribosome quality control, MAM mitochondrial, Vascular, pathway, enrichment, TMEM97
  12. [12] § Results › Custom gene set enrichment ↔ MS3_Sigma_Divergence_Pipeline.py, lines 637–678 · score 0.56 · ribosome quality control, MAM mitochondrial, Vascular, pathway, enrichment, TMEM97
  13. [13] § Materials and methods › Sensitivity analyses ↔ MS3_Sigma_Divergence_Pipeline.py, lines 1–40 · score 0.55 · covariate adjustment, brain regions, sex, age, deconvolution, validation
  14. [14] § Results › Divergent tails encode distinct functional programs ↔ MS3_Sigma_Divergence_Pipeline.py, lines 523–592 · score 0.52 · RAB1B, AAMP, PSMD3, YIPF3, partners, overlap
  15. [15] § Results › Sensitivity analyses ↔ MS3_Sigma_Divergence_Pipeline.py, lines 842–909 · score 0.51 · Covariate adjustment, rank preservation, VCP, thresholds, Cell, Spearman

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 · 1,084 lines · 47 KB · MIT · 8 matches

  1. #!/usr/bin/env python3
  2. """
  3. ===============================================================================
  4. MS3 SIGMA DIVERGENCE WJ-NATIVE PIPELINE — Local Execution (16GB RAM)
  5. ===============================================================================
  7. MANUSCRIPT: "Sigma-1 and Sigma-2 Receptors Exhibit Divergent Genome-Wide
  8. Co-Expression Networks in Human Brain Despite Shared Subcellular
  9. Localization"
  11. TARGET: Frontiers in Molecular Neuroscience
  13. AUTHOR: Drake H. Harbert (ORCID: 0009-0007-7740-3616)
  14. Inner Architecture LLC, Canton, OH 44721, USA
  16. METHODOLOGY:
  17. PRIMARY: Weighted Jaccard (WJ) on continuous genome-wide Spearman
  18. correlation vectors. Threshold-free comparison of full
  19. co-expression architectures.
  20. SUPPLEMENTARY: Binary Jaccard on top 5% network membership (consistency
  21. with other genomics papers in portfolio).
  23. COMPUTES:
  24. - Genome-wide Spearman correlations for SIGMAR1, TMEM97, and 5 additional
  25. targets across 5 brain regions
  26. - WJ on continuous correlation vectors (primary analysis)
  27. - WJ permutation testing (1000 permutations, seed=42)
  28. - Binary Jaccard on top 5% sets (supplementary)
  29. - FDR correction across all 21 pairwise comparisons
  30. - gProfiler GO/KEGG/Reactome enrichment
  31. - 6 custom gene set enrichments
  32. - Cell-type deconvolution sensitivity
  33. - Covariate adjustment (age, sex)
  34. - Multi-region replication (5 brain regions)
  35. - All figures (300 DPI, colorblind-safe)
  36. - Supplementary Tables S1-S4
  37. - provenance.json
  39. USAGE: py -3 MS3_sigma_divergence_wj_pipeline.py
  40. Dependencies: scipy, statsmodels, pandas, numpy, matplotlib, matplotlib-venn,
  41. gprofiler-official, python-docx, openpyxl, requests
  42. ===============================================================================
  43. """
  45. import os
  46. import gc
  47. import json
  48. import time
  49. import warnings
  50. import requests
  51. import gzip
  52. import numpy as np
  53. import pandas as pd
  54. from scipy import stats
  55. from scipy.stats import spearmanr, fisher_exact
  56. from statsmodels.stats.multitest import multipletests
  57. from collections import OrderedDict
  58. from itertools import combinations
  59. from numpy.linalg import lstsq
  60. import matplotlib
  61. matplotlib.use('Agg')
  62. matplotlib.rcParams['pdf.fonttype'] = 42
  63. matplotlib.rcParams['ps.fonttype'] = 42
  64. import matplotlib.pyplot as plt
  65. from matplotlib_venn import venn2
  67. warnings.filterwarnings('ignore')
  69. # ============================================================================
  70. # CONFIG
  71. # ============================================================================
  72. FORCE_RECOMPUTE = True
  73. RANDOM_SEED = 42
  74. np.random.seed(RANDOM_SEED)
  75. N_PERMUTATIONS = 1000
  77. # Paths
  78. DRIVE_BASE = r"G:\My Drive\inner_architecture_research\MS3_JNC_Submission"
  79. RESULTS_DIR = os.path.join(DRIVE_BASE, "MS3_Results")
  80. FIGURES_DIR = os.path.join(DRIVE_BASE, "MS3_Figures")
  81. SUPPL_DIR = os.path.join(DRIVE_BASE, "MS3_Supplementary_Tables")
  82. DATA_DIR = os.path.join(DRIVE_BASE, "data")
  84. for d in [RESULTS_DIR, FIGURES_DIR, SUPPL_DIR, DATA_DIR]:
  85. os.makedirs(d, exist_ok=True)
  87. # Targets
  88. PRIMARY_TARGETS = ['SIGMAR1', 'TMEM97']
  89. TARGET_GENES = ['EIF2S1', 'PELO', 'LTN1', 'NEMF', 'TMEM97', 'HSPA5', 'SIGMAR1']
  91. # Brain regions
  92. BRAIN_REGIONS = OrderedDict([
  93. ('BA9', 'Brain - Frontal Cortex (BA9)'),
  94. ('Putamen', 'Brain - Putamen (basal ganglia)'),
  95. ('Hippocampus', 'Brain - Hippocampus'),
  96. ('NAcc', 'Brain - Nucleus accumbens (basal ganglia)'),
  97. ('BA24', 'Brain - Anterior cingulate cortex (BA24)'),
  98. ])
  99. PRIMARY_REGION = 'BA9'
  100. MIN_MEDIAN_TPM = 1.0
  101. TOP_PERCENT = 5
  103. # Custom gene sets
  104. MAM_MITO_GENES = [
  105. 'VDAC1', 'VDAC2', 'VDAC3', 'MFN1', 'MFN2', 'RHOT1', 'RHOT2',
  106. 'ITPR1', 'ITPR2', 'ITPR3', 'VAPB', 'RMDN3', 'PACS2', 'FATE1',
  107. ]
  108. SIGMA_NETWORK_GENES = ['SIGMAR1', 'TMEM97', 'PGRMC1', 'NPC1']
  109. ER_STRESS_UPR_GENES = [
  110. 'HSPA5', 'HSP90B1', 'DDIT3', 'ATF4', 'ATF6', 'ERN1', 'EIF2AK3',
  111. 'XBP1', 'DNAJB9', 'HERPUD1', 'EDEM1', 'CALR', 'CANX', 'P4HB',
  112. ]
  113. METHYLATION_GENES = [
  114. 'MAT1A', 'MAT2A', 'MAT2B', 'AHCY', 'AHCYL1', 'AHCYL2',
  115. 'MTR', 'MTRR', 'MTHFR', 'BHMT', 'BHMT2', 'CBS', 'CTH',
  116. 'GNMT', 'PEMT', 'NNMT', 'INMT', 'DNMT1',
  117. ]
  118. VASCULAR_GENES = [
  119. 'PECAM1', 'CDH5', 'VWF', 'FLT1', 'KDR', 'ENG', 'CLDN5', 'ESAM',
  120. 'ERG', 'TIE1', 'TEK', 'ANGPT1', 'ANGPT2', 'NOS3', 'MCAM',
  121. 'PODXL', 'EMCN', 'ROBO4',
  122. ]
  123. RQC_GENES = [
  124. 'PELO', 'HBS1L', 'LTN1', 'NEMF', 'ANKZF1',
  125. 'VCP', 'UFD1', 'NPLOC4', 'ZNF598', 'RACK1', 'ABCE1', 'TCF25',
  126. ]
  128. # Cell-type deconvolution markers
  129. CELLTYPE_MARKERS = {
  130. 'Neurons': ['SNAP25', 'SYT1', 'GAD1', 'GAD2', 'SLC17A7', 'RBFOX3',
  131. 'STMN2', 'SYN1', 'NRGN'],
  132. 'Astrocytes': ['AQP4', 'GFAP', 'SLC1A2', 'SLC1A3', 'ALDH1L1', 'GJA1',
  133. 'S100B', 'SOX9', 'GLUL'],
  134. 'Oligodendrocytes': ['MBP', 'MOG', 'PLP1', 'MAG', 'MOBP', 'CLDN11', 'CNP',
  135. 'OPALIN', 'TF'],
  136. 'Microglia': ['CX3CR1', 'P2RY12', 'CSF1R', 'TMEM119', 'AIF1', 'ITGAM',
  137. 'CD68', 'HEXB', 'TREM2'],
  138. 'Endothelial': ['CLDN5', 'FLT1', 'PECAM1', 'VWF', 'CDH5', 'ERG',
  139. 'ESAM', 'TIE1'],
  140. 'OPCs': ['PDGFRA', 'CSPG4', 'OLIG1', 'OLIG2', 'SOX10', 'NKX2-2',
  141. 'GPR17', 'PCDH15', 'NEU4'],
  142. }
  144. # GTEx URLs
  145. GTEX_TPM_URL = "https://storage.googleapis.com/adult-gtex/bulk-gex/v8/rna-seq/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_tpm.gct.gz"
  146. GTEX_SAMPLE_URL = "https://storage.googleapis.com/adult-gtex/annotations/v8/metadata-files/GTEx_Analysis_v8_Annotations_SampleAttributesDS.txt"
  147. GTEX_SUBJECT_URL = "https://storage.googleapis.com/adult-gtex/annotations/v8/metadata-files/GTEx_Analysis_v8_Annotations_SubjectPhenotypesDS.txt"
  150. # ============================================================================
  151. # CORE WJ FUNCTIONS
  152. # ============================================================================
  154. def weighted_jaccard(vec_a, vec_b):
  155. """
  156. Compute Weighted Jaccard similarity on two continuous vectors.
  158. WJ = sum(min(a_i, b_i)) / sum(max(a_i, b_i))
  160. For correlation vectors that can be negative, shift to [0, 1] range first:
  161. r_shifted = (r + 1) / 2
  162. This maps r=-1 -> 0, r=0 -> 0.5, r=1 -> 1.
  163. """
  164. # Shift correlations from [-1, 1] to [0, 1]
  165. a = (np.array(vec_a) + 1) / 2
  166. b = (np.array(vec_b) + 1) / 2
  168. numerator = np.sum(np.minimum(a, b))
  169. denominator = np.sum(np.maximum(a, b))
  171. if denominator == 0:
  172. return 0.0
  173. return numerator / denominator
  176. def weighted_jaccard_permutation_test(vec_a, vec_b, n_perm=1000, seed=42):
  177. """
  178. Permutation test for WJ significance.
  180. Null hypothesis: the two correlation vectors are drawn from the same
  181. underlying distribution (i.e., the two genes see the same transcriptional
  182. architecture).
  184. Permutation strategy: shuffle the gene labels in one vector, breaking
  185. the gene-to-gene correspondence while preserving the marginal distribution.
  186. """
  187. rng = np.random.RandomState(seed)
  188. observed_wj = weighted_jaccard(vec_a, vec_b)
  190. null_distribution = np.zeros(n_perm)
  191. for i in range(n_perm):
  192. perm_b = rng.permutation(vec_b)
  193. null_distribution[i] = weighted_jaccard(vec_a, perm_b)
  195. # Two-sided: how extreme is the observed WJ?
  196. # For divergence testing: low WJ means more divergent
  197. # p-value = proportion of null WJ values <= observed WJ
  198. p_value = np.mean(null_distribution <= observed_wj)
  200. return observed_wj, p_value, null_distribution
  203. def vectorized_spearman(target_expr, other_expr):
  204. """
  205. Compute Spearman correlation between one target and many genes.
  206. Rank-transforms both, then computes Pearson on ranks.
  207. """
  208. from scipy.stats import rankdata
  210. x_rank = rankdata(target_expr)
  211. n = len(target_expr)
  213. # Rank each row of other_expr
  214. Y_rank = np.apply_along_axis(rankdata, 1, other_expr)
  216. # Pearson on ranks = Spearman
  217. x = x_rank - x_rank.mean()
  218. Y = Y_rank - Y_rank.mean(axis=1, keepdims=True)
  219. x_std = np.sqrt(np.sum(x**2))
  220. Y_std = np.sqrt(np.sum(Y**2, axis=1))
  222. valid = Y_std > 0
  223. r_values = np.full(len(other_expr), np.nan)
  224. r_values[valid] = np.dot(Y[valid], x) / (Y_std[valid] * x_std)
  225. return r_values
  228. def compute_custom_enrichment(network_genes, gene_set_list, set_name,
  229. all_genes_list, verbose=True):
  230. """Fisher's exact test for custom gene set enrichment."""
  231. universe = set(all_genes_list)
  232. expressed = [g for g in gene_set_list if g in universe]
  233. in_network = [g for g in expressed if g in network_genes]
  235. n_expressed = len(expressed)
  236. n_in_network = len(in_network)
  237. n_network = len(network_genes)
  238. n_universe = len(universe)
  240. if n_expressed == 0:
  241. return None
  243. expected = n_expressed * n_network / n_universe
  244. fold = n_in_network / expected if expected > 0 else 0
  246. a = n_in_network
  247. b = n_network - n_in_network
  248. c = n_expressed - n_in_network
  249. d = n_universe - n_network - c
  250. table = np.array([[a, b], [c, d]])
  251. _, p_val = fisher_exact(table, alternative='greater')
  253. result = {
  254. 'gene_set': set_name, 'expressed': n_expressed,
  255. 'in_network': n_in_network, 'fold_enrichment': fold,
  256. 'p_value': p_val, 'genes_found': ', '.join(sorted(in_network)),
  257. }
  259. if verbose:
  260. sig = "+" if p_val < 0.05 else "-"
  261. print(f" {sig} {set_name}: {n_in_network}/{n_expressed}, "
  262. f"{fold:.1f}x, p = {p_val:.2e}")
  263. if in_network:
  264. print(f" Genes: {', '.join(sorted(in_network))}")
  265. return result
  268. def partial_correlation_network(target, expr_df, covar_df, top_pct=5):
  269. """Genome-wide Spearman correlations after regressing covariates."""
  270. valid = covar_df.dropna().index
  271. valid = [s for s in valid if s in expr_df.columns]
  272. n_valid = len(valid)
  274. if n_valid < 30:
  275. return None, None, None, n_valid
  277. expr_sub = expr_df[valid]
  278. covar_matrix = covar_df.loc[valid].values
  279. X = np.column_stack([np.ones(n_valid), covar_matrix])
  281. target_vals = expr_sub.loc[target].values.astype(np.float64)
  282. beta, _, _, _ = lstsq(X, target_vals, rcond=None)
  283. target_resid = target_vals - X @ beta
  285. other_genes = [g for g in expr_sub.index if g != target]
  286. other_vals = expr_sub.loc[other_genes].values.astype(np.float64)
  287. betas = lstsq(X, other_vals.T, rcond=None)[0]
  288. resid = other_vals - (X @ betas).T
  290. r_vals = vectorized_spearman(target_resid, resid)
  292. corr_series = pd.Series(r_vals, index=other_genes).dropna().sort_values(ascending=False)
  293. n_top = int(np.ceil(len(corr_series) * top_pct / 100))
  294. top5_set = set(corr_series.head(n_top).index)
  295. threshold = corr_series.iloc[n_top - 1] if n_top <= len(corr_series) else np.nan
  297. return corr_series, threshold, top5_set, n_valid
  300. def sample_to_subject(sample_id):
  301. parts = sample_id.split('-')
  302. return '-'.join(parts[:2]) if len(parts) >= 2 else sample_id
  305. def download_file(url, dest_path):
  306. """Download a file, reuse if already present."""
  307. if os.path.exists(dest_path):
  308. size_mb = os.path.getsize(dest_path) / 1e6
  309. if size_mb > 1:
  310. print(f" Using cached: {os.path.basename(dest_path)} ({size_mb:.0f} MB)")
  311. return
  312. print(f" Downloading {os.path.basename(dest_path)}...")
  313. resp = requests.get(url, stream=True, timeout=600)
  314. resp.raise_for_status()
  315. total = int(resp.headers.get('content-length', 0))
  316. downloaded = 0
  317. with open(dest_path, 'wb') as f:
  318. for chunk in resp.iter_content(chunk_size=8192 * 16):
  319. f.write(chunk)
  320. downloaded += len(chunk)
  321. if total > 0 and downloaded % (50 * 1024 * 1024) < 8192 * 16:
  322. print(f" {downloaded / 1e6:.0f} / {total / 1e6:.0f} MB")
  323. print(f" Done: {os.path.basename(dest_path)} ({downloaded / 1e6:.0f} MB)")
  326. # ============================================================================
  327. # MAIN PIPELINE
  328. # ============================================================================
  330. def main():
  331. t_start = time.time()
  333. print("=" * 70)
  334. print("MS3 SIGMA DIVERGENCE WJ-NATIVE PIPELINE")
  335. print("Primary: Weighted Jaccard on continuous Spearman vectors")
  336. print("Supplementary: Binary Jaccard on top 5% sets")
  337. print("Target: Frontiers in Molecular Neuroscience")
  338. print("Drake H. Harbert -- Inner Architecture LLC")
  339. print("=" * 70)
  341. # ==================================================================
  342. # STEP 1: DOWNLOAD GTEx v8 DATA
  343. # ==================================================================
  344. print(f"\n{'='*70}\nSTEP 1: DOWNLOAD GTEx v8 DATA\n{'='*70}\n")
  346. tpm_file = os.path.join(DATA_DIR, "GTEx_v8_tpm.gct.gz")
  347. sample_file = os.path.join(DATA_DIR, "GTEx_v8_sample_attributes.txt")
  348. subject_file = os.path.join(DATA_DIR, "GTEx_v8_subject_phenotypes.txt")
  350. download_file(GTEX_TPM_URL, tpm_file)
  351. download_file(GTEX_SAMPLE_URL, sample_file)
  352. download_file(GTEX_SUBJECT_URL, subject_file)
  354. # ==================================================================
  355. # STEP 2: LOAD METADATA
  356. # ==================================================================
  357. print(f"\n{'='*70}\nSTEP 2: LOAD METADATA\n{'='*70}\n")
  359. sample_attr = pd.read_csv(sample_file, sep='\t', low_memory=False)
  360. brain_labels = sample_attr[
  361. sample_attr['SMTSD'].str.contains('Brain', na=False)
  362. ]['SMTSD'].unique()
  364. for key, smtsd in BRAIN_REGIONS.items():
  365. if smtsd not in brain_labels:
  366. raise ValueError(f"Region mismatch: {key}: '{smtsd}'")
  368. region_samples = {}
  369. for key, smtsd in BRAIN_REGIONS.items():
  370. region_samples[key] = list(sample_attr.loc[sample_attr['SMTSD'] == smtsd, 'SAMPID'])
  371. print(f" {key}: {len(region_samples[key])} samples")
  373. subject_pheno = pd.read_csv(subject_file, sep='\t')
  374. ba9_samples_list = region_samples[PRIMARY_REGION]
  375. ba9_subjects = [sample_to_subject(s) for s in ba9_samples_list]
  376. ba9_subject_df = pd.DataFrame({'SAMPID': ba9_samples_list, 'SUBJID': ba9_subjects})
  377. ba9_subject_df = ba9_subject_df.merge(subject_pheno, on='SUBJID', how='left')
  379. # ==================================================================
  380. # STEP 3: LOAD GTEx TPM
  381. # ==================================================================
  382. print(f"\n{'='*70}\nSTEP 3: LOAD GTEx TPM\n{'='*70}\n")
  384. print(" Reading header...")
  385. with gzip.open(tpm_file, 'rt') as f:
  386. f.readline(); f.readline()
  387. header = f.readline().strip().split('\t')
  389. all_brain_samples = set()
  390. for key in BRAIN_REGIONS:
  391. all_brain_samples.update(region_samples[key])
  393. sample_cols = {i: col for i, col in enumerate(header) if col in all_brain_samples}
  394. keep_cols = [0, 1] + sorted(sample_cols.keys())
  395. print(f" Brain samples: {len(sample_cols)}, columns to load: {len(keep_cols)}")
  397. print(" Loading expression data...")
  398. chunks = []
  399. with gzip.open(tpm_file, 'rt') as f:
  400. f.readline(); f.readline(); f.readline()
  401. row_buffer = []
  402. for line_num, line in enumerate(f):
  403. parts = line.strip().split('\t')
  404. row_buffer.append([parts[i] if i < len(parts) else '' for i in keep_cols])
  405. if len(row_buffer) >= 5000:
  406. chunks.append(pd.DataFrame(row_buffer))
  407. row_buffer = []
  408. if line_num % 10000 == 0:
  409. print(f" {line_num:,} genes...")
  410. if row_buffer:
  411. chunks.append(pd.DataFrame(row_buffer))
  413. tpm_raw = pd.concat(chunks, ignore_index=True)
  414. del chunks, row_buffer; gc.collect()
  416. col_names = ['Name', 'Description'] + [sample_cols[i] for i in sorted(sample_cols.keys())]
  417. tpm_raw.columns = col_names
  418. tpm_raw = tpm_raw.set_index('Name')
  419. gene_descriptions = tpm_raw['Description'].copy()
  420. tpm_raw = tpm_raw.drop('Description', axis=1)
  421. tpm_raw = tpm_raw.apply(pd.to_numeric, errors='coerce').astype(np.float32)
  422. print(f" Loaded: {tpm_raw.shape[0]} genes x {tpm_raw.shape[1]} samples")
  424. # ==================================================================
  425. # STEP 4: PREPARE REGION MATRICES
  426. # ==================================================================
  427. print(f"\n{'='*70}\nSTEP 4: PREPARE REGION MATRICES\n{'='*70}\n")
  429. def prepare_region(region_key):
  430. samples = [s for s in region_samples[region_key] if s in tpm_raw.columns]
  431. expr = tpm_raw[samples].copy()
  432. expr['gene_symbol'] = gene_descriptions.reindex(expr.index)
  433. expr = expr.dropna(subset=['gene_symbol'])
  434. expr = expr[expr['gene_symbol'] != '']
  435. expr['median_tpm'] = expr[samples].median(axis=1)
  436. expr = expr.sort_values('median_tpm', ascending=False)
  437. expr = expr.drop_duplicates(subset='gene_symbol', keep='first')
  438. expr = expr[expr['median_tpm'] >= MIN_MEDIAN_TPM]
  439. expr = expr.set_index('gene_symbol').drop('median_tpm', axis=1)
  440. return np.log2(expr + 1), len(samples), len(expr), list(expr.index)
  442. region_data = {}
  443. for key in BRAIN_REGIONS:
  444. log2_expr, n_samp, n_genes, genes = prepare_region(key)
  445. region_data[key] = {'expr': log2_expr, 'n_samples': n_samp,
  446. 'n_genes': n_genes, 'genes': genes}
  447. present = [g for g in TARGET_GENES if g in genes]
  448. print(f" {key}: n={n_samp}, genes={n_genes}, targets={len(present)}/7")
  450. del tpm_raw; gc.collect()
  452. # ==================================================================
  453. # STEP 5: GENOME-WIDE SPEARMAN CO-EXPRESSION
  454. # ==================================================================
  455. print(f"\n{'='*70}\nSTEP 5: GENOME-WIDE SPEARMAN CO-EXPRESSION -- {PRIMARY_REGION}\n{'='*70}\n")
  457. ba9 = region_data[PRIMARY_REGION]['expr']
  458. n_genes_ba9 = region_data[PRIMARY_REGION]['n_genes']
  459. n_samples_ba9 = region_data[PRIMARY_REGION]['n_samples']
  460. print(f"Matrix: {n_genes_ba9} genes x {n_samples_ba9} samples")
  461. print(f"Correlation method: Spearman (rank-based)\n")
  463. correlations = {}
  464. for target in TARGET_GENES:
  465. if target not in ba9.index:
  466. print(f" WARNING: {target} not found")
  467. continue
  468. target_expr = ba9.loc[target].values.astype(np.float64)
  469. other_genes = [g for g in ba9.index if g != target]
  470. other_expr = ba9.loc[other_genes].values.astype(np.float64)
  472. r_values = vectorized_spearman(target_expr, other_expr)
  474. corr_series = pd.Series(r_values, index=other_genes, name=target)
  475. correlations[target] = corr_series.dropna().sort_values(ascending=False)
  476. print(f" {target}: {len(corr_series.dropna())} genes, "
  477. f"top = {corr_series.dropna().idxmax()} "
  478. f"(rho = {corr_series.dropna().max():.3f})")
  480. # ==================================================================
  481. # STEP 6: PRIMARY ANALYSIS — WEIGHTED JACCARD ON CONTINUOUS VECTORS
  482. # ==================================================================
  483. print(f"\n{'='*70}")
  484. print("STEP 6: PRIMARY ANALYSIS -- WEIGHTED JACCARD (CONTINUOUS)")
  485. print(f"{'='*70}\n")
  486. print(f"Permutations: {N_PERMUTATIONS}, seed: {RANDOM_SEED}\n")
  488. # Align all correlation vectors to common gene set
  489. common_genes = sorted(set.intersection(*[set(correlations[g].index)
  490. for g in TARGET_GENES
  491. if g in correlations]))
  492. n_common = len(common_genes)
  493. print(f"Common genes across all 7 targets: {n_common}\n")
  495. # Build aligned matrix
  496. corr_matrix = pd.DataFrame({g: correlations[g].reindex(common_genes)
  497. for g in TARGET_GENES if g in correlations})
  499. # Compute WJ for all 21 pairwise comparisons
  500. wj_results = []
  501. for g1, g2 in combinations(TARGET_GENES, 2):
  502. if g1 not in corr_matrix.columns or g2 not in corr_matrix.columns:
  503. continue
  505. vec_a = corr_matrix[g1].values
  506. vec_b = corr_matrix[g2].values
  508. wj_obs, wj_pval, null_dist = weighted_jaccard_permutation_test(
  509. vec_a, vec_b, n_perm=N_PERMUTATIONS, seed=RANDOM_SEED)
  511. # Also compute Spearman rank correlation between vectors
  512. rho, rho_p = spearmanr(vec_a, vec_b)
  514. wj_results.append({
  515. 'gene1': g1, 'gene2': g2,
  516. 'wj': wj_obs,
  517. 'wj_perm_p': wj_pval,
  518. 'null_mean': np.mean(null_dist),
  519. 'null_std': np.std(null_dist),
  520. 'z_score': (wj_obs - np.mean(null_dist)) / np.std(null_dist) if np.std(null_dist) > 0 else 0,
  521. 'spearman_rho': rho,
  522. 'spearman_p': rho_p,
  523. 'n_genes': n_common,
  524. })
  526. wj_df = pd.DataFrame(wj_results).sort_values('wj', ascending=True)
  528. # FDR correction on WJ permutation p-values
  529. reject_wj, fdr_wj, _, _ = multipletests(wj_df['wj_perm_p'].values, method='fdr_bh')
  530. wj_df['wj_perm_p_fdr'] = fdr_wj
  531. wj_df['fdr_significant'] = reject_wj
  533. wj_df.to_csv(os.path.join(RESULTS_DIR, "wj_continuous_all_21_pairs.csv"), index=False)
  535. # Print all results
  536. print(f"{'Gene1':>10s} {'Gene2':>10s} {'WJ':>8s} {'z':>8s} {'p_perm':>10s} "
  537. f"{'FDR':>10s} {'rho':>8s}")
  538. print("-" * 70)
  539. for _, row in wj_df.iterrows():
  540. print(f"{row['gene1']:>10s} {row['gene2']:>10s} {row['wj']:>8.4f} "
  541. f"{row['z_score']:>8.2f} {row['wj_perm_p']:>10.4f} "
  542. f"{row['wj_perm_p_fdr']:>10.4f} {row['spearman_rho']:>8.3f}")
  544. # Key pair: SIGMAR1-TMEM97
  545. st_wj = wj_df[
  546. ((wj_df['gene1'] == 'SIGMAR1') & (wj_df['gene2'] == 'TMEM97')) |
  547. ((wj_df['gene1'] == 'TMEM97') & (wj_df['gene2'] == 'SIGMAR1'))
  548. ].iloc[0]
  550. print(f"\n{'='*60}")
  551. print("PRIMARY RESULT: SIGMAR1-TMEM97 WEIGHTED JACCARD")
  552. print(f" WJ = {st_wj['wj']:.6f}")
  553. print(f" z-score = {st_wj['z_score']:.2f}")
  554. print(f" Permutation p = {st_wj['wj_perm_p']:.4f}")
  555. print(f" FDR p = {st_wj['wj_perm_p_fdr']:.4f}")
  556. print(f" Null mean = {st_wj['null_mean']:.6f}, std = {st_wj['null_std']:.6f}")
  557. print(f" Spearman rho = {st_wj['spearman_rho']:.4f}")
  558. print(f"{'='*60}")
  560. # SIGMAR1-LTN1
  561. sl_wj = wj_df[
  562. ((wj_df['gene1'] == 'SIGMAR1') & (wj_df['gene2'] == 'LTN1')) |
  563. ((wj_df['gene1'] == 'LTN1') & (wj_df['gene2'] == 'SIGMAR1'))
  564. ].iloc[0]
  565. print(f"\nSIGMAR1-LTN1: WJ={sl_wj['wj']:.4f}, z={sl_wj['z_score']:.2f}, "
  566. f"p={sl_wj['wj_perm_p']:.4f}")
  568. # ==================================================================
  569. # STEP 7: SUPPLEMENTARY — BINARY JACCARD ON TOP 5% SETS
  570. # ==================================================================
  571. print(f"\n{'='*70}")
  572. print("STEP 7: SUPPLEMENTARY -- BINARY JACCARD (TOP 5% SETS)")
  573. print(f"{'='*70}\n")
  575. networks = {}
  576. thresholds = {}
  577. for target, corr in correlations.items():
  578. n_total = len(corr)
  579. n_top = int(np.ceil(n_total * TOP_PERCENT / 100))
  580. networks[target] = set(corr.head(n_top).index)
  581. thresholds[target] = corr.iloc[n_top - 1]
  582. print(f" {target}: top 5% = {n_top} genes, rho >= {thresholds[target]:.3f}")
  584. gene_universe = len(correlations[TARGET_GENES[0]])
  586. binary_results = []
  587. for g1, g2 in combinations(TARGET_GENES, 2):
  588. if g1 not in networks or g2 not in networks:
  589. continue
  590. set1, set2 = networks[g1], networks[g2]
  591. shared = set1 & set2
  592. union_set = set1 | set2
  593. jaccard = len(shared) / len(union_set) if len(union_set) > 0 else 0
  595. a, b, c = len(shared), len(set1 - set2), len(set2 - set1)
  596. d = gene_universe - len(union_set)
  597. fisher_or, fisher_p = fisher_exact(np.array([[a, b], [c, d]]),
  598. alternative='greater')
  600. binary_results.append({
  601. 'gene1': g1, 'gene2': g2,
  602. 'shared': len(shared), 'binary_jaccard': jaccard,
  603. 'fisher_or': fisher_or, 'fisher_p': fisher_p,
  604. 'set1_size': len(set1), 'set2_size': len(set2),
  605. })
  607. binary_df = pd.DataFrame(binary_results).sort_values('binary_jaccard', ascending=False)
  608. reject_b, fdr_b, _, _ = multipletests(binary_df['fisher_p'].values, method='fdr_bh')
  609. binary_df['fisher_p_fdr'] = fdr_b
  610. binary_df.to_csv(os.path.join(RESULTS_DIR, "binary_jaccard_top5pct.csv"), index=False)
  612. # Merge WJ and binary results for comparison
  613. comparison = wj_df[['gene1', 'gene2', 'wj', 'z_score', 'wj_perm_p']].merge(
  614. binary_df[['gene1', 'gene2', 'binary_jaccard', 'shared', 'fisher_p']],
  615. on=['gene1', 'gene2'], how='outer'
  616. )
  617. comparison.to_csv(os.path.join(RESULTS_DIR, "wj_vs_binary_comparison.csv"), index=False)
  619. print("\n--- WJ vs Binary Jaccard comparison ---")
  620. print(f"{'Pair':>25s} {'WJ':>8s} {'Binary J':>10s} {'Shared':>7s}")
  621. print("-" * 55)
  622. for _, row in comparison.sort_values('wj').iterrows():
  623. pair = f"{row['gene1']}-{row['gene2']}"
  624. bj = f"{row['binary_jaccard']:.3f}" if pd.notna(row.get('binary_jaccard')) else "N/A"
  625. sh = f"{int(row['shared'])}" if pd.notna(row.get('shared')) else "N/A"
  626. print(f"{pair:>25s} {row['wj']:>8.4f} {bj:>10s} {sh:>7s}")
  628. # Key validation
  629. st_binary = binary_df[
  630. ((binary_df['gene1'] == 'SIGMAR1') & (binary_df['gene2'] == 'TMEM97')) |
  631. ((binary_df['gene1'] == 'TMEM97') & (binary_df['gene2'] == 'SIGMAR1'))
  632. ].iloc[0]
  634. print(f"\n SIGMAR1-TMEM97 binary J = {st_binary['binary_jaccard']:.3f}, "
  635. f"shared = {int(st_binary['shared'])}")
  637. # Gene sets
  638. sigmar1_unique = sorted(networks['SIGMAR1'] - networks['TMEM97'])
  639. tmem97_unique = sorted(networks['TMEM97'] - networks['SIGMAR1'])
  640. shared_st = sorted(networks['SIGMAR1'] & networks['TMEM97'])
  642. # ==================================================================
  643. # STEP 8: TOP CO-EXPRESSION PARTNERS
  644. # ==================================================================
  645. print(f"\n{'='*70}\nSTEP 8: TOP CO-EXPRESSION PARTNERS\n{'='*70}\n")
  647. for target in PRIMARY_TARGETS:
  648. print(f" {target} top 10:")
  649. for rank, (gene, r_val) in enumerate(correlations[target].head(10).items(), 1):
  650. print(f" {rank:2d}. {gene:12s} rho = {r_val:.4f}")
  651. print()
  653. # ==================================================================
  654. # STEP 9: CUSTOM GENE SET ENRICHMENT
  655. # ==================================================================
  656. print(f"\n{'='*70}\nSTEP 9: CUSTOM GENE SET ENRICHMENT\n{'='*70}\n")
  658. custom_sets = OrderedDict([
  659. ('MAM-mitochondrial', MAM_MITO_GENES),
  660. ('Sigma receptor network', SIGMA_NETWORK_GENES),
  661. ('ER stress/UPR', ER_STRESS_UPR_GENES),
  662. ('Methylation pathway', METHYLATION_GENES),
  663. ('Vascular markers (neg ctrl)', VASCULAR_GENES),
  664. ('Ribosome quality control', RQC_GENES),
  665. ])
  667. all_custom_results = {}
  668. for target in PRIMARY_TARGETS:
  669. print(f"\n--- {target} top 5% custom enrichment ---\n")
  670. all_genes_t = list(correlations[target].index) + [target]
  671. target_results = []
  672. for set_name, gene_list in custom_sets.items():
  673. result = compute_custom_enrichment(
  674. networks[target], gene_list, set_name, all_genes_t)
  675. if result:
  676. result['target'] = target
  677. target_results.append(result)
  678. all_custom_results[target] = target_results
  679. pd.DataFrame(target_results).to_csv(
  680. os.path.join(RESULTS_DIR, f"{target}_custom_enrichment.csv"), index=False)
  682. # ==================================================================
  683. # STEP 10: gProfiler ENRICHMENT
  684. # ==================================================================
  685. print(f"\n{'='*70}\nSTEP 10: gProfiler ENRICHMENT\n{'='*70}\n")
  687. try:
  688. from gprofiler import GProfiler
  689. gp = GProfiler(return_dataframe=True)
  690. except ImportError:
  691. print("WARNING: gprofiler not installed")
  692. gp = None
  694. def run_gprofiler(gene_list, bg, label=""):
  695. if gp is None:
  696. return pd.DataFrame()
  697. try:
  698. df = gp.profile(organism='hsapiens', query=list(gene_list),
  699. background=list(bg),
  700. sources=['GO:BP', 'GO:MF', 'GO:CC', 'KEGG', 'REAC'],
  701. significance_threshold_method='g_SCS',
  702. user_threshold=0.05, no_evidences=False)
  703. if df is not None and len(df) > 0:
  704. print(f" {label}: {len(df)} terms")
  705. return df
  706. print(f" {label}: 0 terms")
  707. return pd.DataFrame()
  708. except Exception as e:
  709. print(f" {label}: error -- {e}")
  710. return pd.DataFrame()
  712. background = list(correlations['SIGMAR1'].index) + ['SIGMAR1']
  714. all_gprofiler = {}
  715. for name, genes in [
  716. ('SIGMAR1_full', list(networks['SIGMAR1'])),
  717. ('TMEM97_full', list(networks['TMEM97'])),
  718. ('SIGMAR1_unique', sigmar1_unique),
  719. ('TMEM97_unique', tmem97_unique),
  720. ('shared_SIGMAR1_TMEM97', shared_st),
  721. ]:
  722. df = run_gprofiler(genes, background, name)
  723. all_gprofiler[name] = df
  724. if len(df) > 0:
  725. df.to_csv(os.path.join(RESULTS_DIR, f"gProfiler_{name}.csv"), index=False)
  727. # ==================================================================
  728. # STEP 11: CELL-TYPE DECONVOLUTION
  729. # ==================================================================
  730. print(f"\n{'='*70}\nSTEP 11: CELL-TYPE DECONVOLUTION SENSITIVITY\n{'='*70}\n")
  732. celltype_proportions = pd.DataFrame(index=ba9.columns)
  733. for ct_name, markers in CELLTYPE_MARKERS.items():
  734. present = [g for g in markers if g in ba9.index]
  735. if present:
  736. celltype_proportions[ct_name] = ba9.loc[present].mean(axis=0)
  737. print(f" {ct_name}: {len(present)}/{len(markers)} markers")
  739. for target in PRIMARY_TARGETS:
  740. corr_ct, thr_ct, net_ct, n_valid = partial_correlation_network(
  741. target, ba9, celltype_proportions)
  742. if corr_ct is None:
  743. continue
  744. common = sorted(set(corr_ct.index) & set(correlations[target].index))
  745. rho_preserve, _ = spearmanr(corr_ct.reindex(common).values,
  746. correlations[target].reindex(common).values)
  747. print(f" {target}: rank preservation rho = {rho_preserve:.4f}")
  748. if target == 'SIGMAR1' and net_ct:
  749. print(f" VCP in adjusted top 5%: {'VCP' in net_ct}")
  751. # ==================================================================
  752. # STEP 12: COVARIATE ADJUSTMENT
  753. # ==================================================================
  754. print(f"\n{'='*70}\nSTEP 12: COVARIATE ADJUSTMENT (AGE + SEX)\n{'='*70}\n")
  756. ba9_covar = ba9_subject_df.set_index('SAMPID').reindex(ba9.columns)
  757. age_map = {'20-29': 25, '30-39': 35, '40-49': 45, '50-59': 55,
  758. '60-69': 65, '70-79': 75}
  759. ba9_covar['AGE_MID'] = ba9_covar['AGE'].map(age_map) if 'AGE' in ba9_covar.columns else np.nan
  760. ba9_covar['SEX_NUM'] = ba9_covar['SEX'].astype(float) if 'SEX' in ba9_covar.columns else np.nan
  761. covar_df = ba9_covar[['AGE_MID', 'SEX_NUM']]
  763. for target in PRIMARY_TARGETS:
  764. corr_as, thr_as, _, n_valid = partial_correlation_network(target, ba9, covar_df)
  765. if corr_as is None:
  766. continue
  767. common = sorted(set(corr_as.index) & set(correlations[target].index))
  768. rho_p, _ = spearmanr(corr_as.reindex(common).values,
  769. correlations[target].reindex(common).values)
  770. print(f" {target}: rank preservation rho = {rho_p:.4f}")
  772. # ==================================================================
  773. # STEP 13: MULTI-REGION REPLICATION (WJ + binary)
  774. # ==================================================================
  775. print(f"\n{'='*70}\nSTEP 13: MULTI-REGION REPLICATION\n{'='*70}\n")
  777. region_correlations = {}
  778. for region_key in BRAIN_REGIONS:
  779. expr = region_data[region_key]['expr']
  780. region_correlations[region_key] = {}
  781. for target in PRIMARY_TARGETS:
  782. if target not in expr.index:
  783. continue
  784. target_vals = expr.loc[target].values.astype(np.float64)
  785. other_genes = [g for g in expr.index if g != target]
  786. other_vals = expr.loc[other_genes].values.astype(np.float64)
  787. r_vals = vectorized_spearman(target_vals, other_vals)
  788. region_correlations[region_key][target] = pd.Series(
  789. r_vals, index=other_genes).dropna().sort_values(ascending=False)
  790. print(f" {region_key} (n={region_data[region_key]['n_samples']}): done")
  792. # Cross-region rank correlations
  793. region_keys = list(BRAIN_REGIONS.keys())
  794. cross_region_matrix = {}
  795. for target in PRIMARY_TARGETS:
  796. mat = np.zeros((5, 5))
  797. for i, r1 in enumerate(region_keys):
  798. for j, r2 in enumerate(region_keys):
  799. if i == j:
  800. mat[i, j] = 1.0
  801. continue
  802. if target in region_correlations.get(r1, {}) and target in region_correlations.get(r2, {}):
  803. c1, c2 = region_correlations[r1][target], region_correlations[r2][target]
  804. common = sorted(set(c1.index) & set(c2.index))
  805. if len(common) > 100:
  806. mat[i, j], _ = spearmanr(c1.reindex(common).values, c2.reindex(common).values)
  807. cross_region_matrix[target] = mat
  808. vals = mat[np.triu_indices(5, k=1)]
  809. print(f"\n {target} cross-region rho: {vals.min():.3f}-{vals.max():.3f}")
  811. # WJ across regions
  812. print("\n--- SIGMAR1-TMEM97 WJ across regions ---")
  813. region_wj_results = {}
  814. for region_key in BRAIN_REGIONS:
  815. if 'SIGMAR1' not in region_correlations.get(region_key, {}) or \
  816. 'TMEM97' not in region_correlations.get(region_key, {}):
  817. continue
  818. s_corr = region_correlations[region_key]['SIGMAR1']
  819. t_corr = region_correlations[region_key]['TMEM97']
  820. common_r = sorted(set(s_corr.index) & set(t_corr.index))
  821. wj_r = weighted_jaccard(s_corr.reindex(common_r).values,
  822. t_corr.reindex(common_r).values)
  824. # Binary Jaccard too
  825. n_s = int(np.ceil(len(s_corr) * TOP_PERCENT / 100))
  826. n_t = int(np.ceil(len(t_corr) * TOP_PERCENT / 100))
  827. s_set, t_set = set(s_corr.head(n_s).index), set(t_corr.head(n_t).index)
  828. shared_r = s_set & t_set
  829. bj_r = len(shared_r) / len(s_set | t_set) if len(s_set | t_set) > 0 else 0
  831. region_wj_results[region_key] = {'wj': wj_r, 'binary_j': bj_r,
  832. 'shared': len(shared_r)}
  833. print(f" {region_key}: WJ={wj_r:.4f}, binary J={bj_r:.3f}, shared={len(shared_r)}")
  835. # ==================================================================
  836. # STEP 14: EXPORT GENE LISTS
  837. # ==================================================================
  838. print(f"\n{'='*70}\nSTEP 14: EXPORT GENE LISTS\n{'='*70}\n")
  840. for target in PRIMARY_TARGETS:
  841. top5 = correlations[target].head(len(networks[target]))
  842. pd.DataFrame({'gene': top5.index, f'rho_with_{target}': top5.values}).to_csv(
  843. os.path.join(RESULTS_DIR, f"{target}_top5pct.csv"), index=False)
  845. pd.DataFrame({'gene': shared_st}).to_csv(
  846. os.path.join(RESULTS_DIR, "SIGMAR1_TMEM97_shared.csv"), index=False)
  847. pd.DataFrame({'gene': sigmar1_unique}).to_csv(
  848. os.path.join(RESULTS_DIR, "SIGMAR1_unique.csv"), index=False)
  849. pd.DataFrame({'gene': tmem97_unique}).to_csv(
  850. os.path.join(RESULTS_DIR, "TMEM97_unique.csv"), index=False)
  852. for target in TARGET_GENES:
  853. if target in correlations:
  854. corr = correlations[target]
  855. pd.DataFrame({'gene': corr.index, f'rho_with_{target}': corr.values,
  856. 'rank': range(1, len(corr) + 1)}).to_csv(
  857. os.path.join(RESULTS_DIR, f"{target}_genome_wide_rankings.csv"), index=False)
  859. print(f" Exported: {len(TARGET_GENES)} ranking files, 3 gene set files")
  861. # ==================================================================
  862. # STEP 15: FIGURES
  863. # ==================================================================
  864. print(f"\n{'='*70}\nSTEP 15: FIGURES\n{'='*70}\n")
  866. # Figure 1: Divergent Networks (Venn + scatter + top partners)
  867. fig = plt.figure(figsize=(18, 6))
  869. ax_a = fig.add_axes([0.02, 0.12, 0.28, 0.80])
  870. ax_a.text(-0.05, 1.08, 'A', fontsize=20, fontweight='bold', va='top',
  871. transform=ax_a.transAxes)
  872. v = venn2(subsets=(len(sigmar1_unique), len(tmem97_unique), len(shared_st)),
  873. set_labels=('SIGMAR1', 'TMEM97'), ax=ax_a)
  874. for pid, color in [('10', '#1565C0'), ('01', '#E65100'), ('11', '#7B1FA2')]:
  875. v.get_patch_by_id(pid).set_color(color)
  876. v.get_patch_by_id(pid).set_alpha(0.7)
  877. j_val = len(shared_st) / (len(sigmar1_unique) + len(tmem97_unique) + len(shared_st))
  878. ax_a.set_title(f'Binary Jaccard = {j_val:.3f}\nWJ = {st_wj["wj"]:.4f}',
  879. fontsize=12, fontweight='bold')
  881. ax_b = fig.add_axes([0.36, 0.12, 0.28, 0.80])
  882. ax_b.text(-0.05, 1.08, 'B', fontsize=20, fontweight='bold', va='top',
  883. transform=ax_b.transAxes)
  884. s_vals = corr_matrix['SIGMAR1'].values
  885. t_vals = corr_matrix['TMEM97'].values
  886. ax_b.scatter(s_vals, t_vals, s=1, alpha=0.15, c='#555', rasterized=True)
  887. ax_b.set_xlabel('SIGMAR1 Spearman rho', fontsize=12)
  888. ax_b.set_ylabel('TMEM97 Spearman rho', fontsize=12)
  889. ax_b.set_title(f'Genome-wide correlation vectors\nWJ = {st_wj["wj"]:.4f}',
  890. fontsize=12, fontweight='bold')
  891. ax_b.plot([-0.5, 1], [-0.5, 1], 'k--', alpha=0.3, lw=1)
  893. ax_c = fig.add_axes([0.70, 0.12, 0.28, 0.80])
  894. ax_c.text(-0.05, 1.08, 'C', fontsize=20, fontweight='bold', va='top',
  895. transform=ax_c.transAxes)
  896. top10_s = correlations['SIGMAR1'].head(10)
  897. top10_t = correlations['TMEM97'].head(10)
  898. y_pos = np.arange(10)
  899. bh = 0.35
  900. ax_c.barh(y_pos + bh/2, top10_s.values, bh, label='SIGMAR1', color='#1565C0', alpha=0.8)
  901. ax_c.barh(y_pos - bh/2, top10_t.values, bh, label='TMEM97', color='#E65100', alpha=0.8)
  902. ax_c.set_yticks(y_pos)
  903. ax_c.set_yticklabels([f"{top10_s.index[i]} | {top10_t.index[i]}" for i in range(10)], fontsize=8)
  904. ax_c.set_xlabel('Spearman rho', fontsize=12)
  905. ax_c.set_title('Top 10 co-expression partners', fontsize=12, fontweight='bold')
  906. ax_c.legend(fontsize=10)
  907. ax_c.invert_yaxis()
  909. plt.savefig(os.path.join(FIGURES_DIR, "Figure1_Divergent_Networks.png"),
  910. dpi=300, bbox_inches='tight')
  911. plt.savefig(os.path.join(FIGURES_DIR, "Figure1_Divergent_Networks.pdf"),
  912. bbox_inches='tight')
  913. plt.close()
  914. print(" Figure 1 saved")
  916. # Figure 2: GO Enrichment
  917. fig, axes = plt.subplots(2, 3, figsize=(20, 12))
  918. panel_data = [
  919. (all_gprofiler.get('SIGMAR1_unique', pd.DataFrame()), 'GO:BP', 'SIGMAR1-unique GO:BP'),
  920. (all_gprofiler.get('TMEM97_unique', pd.DataFrame()), 'GO:BP', 'TMEM97-unique GO:BP'),
  921. (all_gprofiler.get('shared_SIGMAR1_TMEM97', pd.DataFrame()), 'GO:BP', 'Shared GO:BP'),
  922. (all_gprofiler.get('SIGMAR1_full', pd.DataFrame()), 'REAC', 'SIGMAR1 Reactome'),
  923. (all_gprofiler.get('TMEM97_full', pd.DataFrame()), 'REAC', 'TMEM97 Reactome'),
  924. (all_gprofiler.get('SIGMAR1_full', pd.DataFrame()), 'KEGG', 'SIGMAR1 KEGG'),
  925. ]
  926. for idx, (go_df, source, title) in enumerate(panel_data):
  927. ax = axes.flat[idx]
  928. ax.text(-0.08, 1.05, chr(65+idx), fontsize=16, fontweight='bold',
  929. va='top', transform=ax.transAxes)
  930. if go_df is not None and len(go_df) > 0:
  931. subset = go_df[go_df['source'] == source].sort_values('p_value').head(8)
  932. if len(subset) > 0:
  933. names = [n[:45] for n in subset['name'].values]
  934. pvals = [-np.log10(p) for p in subset['p_value'].values]
  935. color = '#1565C0' if 'SIGMAR1' in title else '#E65100' if 'TMEM97' in title else '#7B1FA2'
  936. ax.barh(range(len(names)), pvals, color=color, alpha=0.8)
  937. ax.set_yticks(range(len(names)))
  938. ax.set_yticklabels(names, fontsize=8)
  939. ax.set_xlabel('-log10(p)', fontsize=10)
  940. ax.invert_yaxis()
  941. else:
  942. ax.text(0.5, 0.5, 'No terms', ha='center', va='center',
  943. transform=ax.transAxes, color='gray')
  944. else:
  945. ax.text(0.5, 0.5, 'No data', ha='center', va='center',
  946. transform=ax.transAxes, color='gray')
  947. ax.set_title(title, fontsize=11, fontweight='bold')
  949. plt.tight_layout()
  950. plt.savefig(os.path.join(FIGURES_DIR, "Figure2_GO_Enrichment.png"), dpi=300, bbox_inches='tight')
  951. plt.savefig(os.path.join(FIGURES_DIR, "Figure2_GO_Enrichment.pdf"), bbox_inches='tight')
  952. plt.close()
  953. print(" Figure 2 saved")
  955. # Figure 3: Multi-region replication
  956. fig, axes = plt.subplots(1, 2, figsize=(14, 6))
  957. for idx, target in enumerate(PRIMARY_TARGETS):
  958. ax = axes[idx]
  959. ax.text(-0.08, 1.05, chr(65+idx), fontsize=16, fontweight='bold',
  960. va='top', transform=ax.transAxes)
  961. mat = cross_region_matrix[target]
  962. im = ax.imshow(mat, cmap='viridis', vmin=0.8, vmax=1.0)
  963. for i in range(5):
  964. for j in range(5):
  965. ax.text(j, i, f'{mat[i,j]:.3f}', ha='center', va='center',
  966. fontsize=9, color='white' if mat[i,j] < 0.92 else 'black')
  967. ax.set_xticks(range(5)); ax.set_xticklabels(region_keys, fontsize=9, rotation=45, ha='right')
  968. ax.set_yticks(range(5)); ax.set_yticklabels(region_keys, fontsize=9)
  969. ax.set_title(f'{target} cross-region', fontsize=12, fontweight='bold')
  970. plt.colorbar(im, ax=ax, shrink=0.8, label='Spearman rho')
  971. plt.tight_layout()
  972. plt.savefig(os.path.join(FIGURES_DIR, "Figure3_MultiRegion_Replication.png"), dpi=300, bbox_inches='tight')
  973. plt.savefig(os.path.join(FIGURES_DIR, "Figure3_MultiRegion_Replication.pdf"), bbox_inches='tight')
  974. plt.close()
  975. print(" Figure 3 saved")
  977. # ==================================================================
  978. # STEP 16: SUPPLEMENTARY TABLES
  979. # ==================================================================
  980. print(f"\n{'='*70}\nSTEP 16: SUPPLEMENTARY TABLES\n{'='*70}\n")
  982. import openpyxl
  983. wb = openpyxl.Workbook()
  984. for i, target in enumerate(TARGET_GENES):
  985. if target in correlations:
  986. ws = wb.active if i == 0 else wb.create_sheet()
  987. ws.title = target
  988. ws.append(['Gene', f'rho_with_{target}', 'Rank'])
  989. for rank, (gene, r_val) in enumerate(correlations[target].items(), 1):
  990. ws.append([gene, round(r_val, 6), rank])
  991. wb.save(os.path.join(SUPPL_DIR, "Table_S1_Genome_Wide_Rankings.xlsx"))
  993. wb2 = openpyxl.Workbook()
  994. first = True
  995. for name, df in all_gprofiler.items():
  996. if len(df) > 0:
  997. ws = wb2.active if first else wb2.create_sheet()
  998. ws.title = name[:31]
  999. cols = [c for c in df.columns if c not in ['query', 'parents']]
  1000. ws.append(cols)
  1001. for _, row in df[cols].iterrows():
  1002. ws.append([str(v) for v in row.values])
  1003. first = False
  1004. wb2.save(os.path.join(SUPPL_DIR, "Table_S2_gProfiler_Enrichment.xlsx"))
  1006. all_custom_rows = []
  1007. for target, results in all_custom_results.items():
  1008. all_custom_rows.extend(results)
  1009. pd.DataFrame(all_custom_rows).to_excel(
  1010. os.path.join(SUPPL_DIR, "Table_S3_Custom_Gene_Set_Enrichment.xlsx"), index=False)
  1012. shared_detail = [{'gene': g,
  1013. 'rho_with_SIGMAR1': round(correlations['SIGMAR1'].get(g, np.nan), 6),
  1014. 'rho_with_TMEM97': round(correlations['TMEM97'].get(g, np.nan), 6)}
  1015. for g in shared_st]
  1016. pd.DataFrame(shared_detail).to_excel(
  1017. os.path.join(SUPPL_DIR, "Table_S4_Shared_Genes.xlsx"), index=False)
  1019. print(" Tables S1-S4 saved")
  1021. # ==================================================================
  1022. # STEP 17: PROVENANCE
  1023. # ==================================================================
  1024. print(f"\n{'='*70}\nSTEP 17: PROVENANCE\n{'='*70}\n")
  1026. from datetime import datetime
  1027. provenance = {
  1028. "methodology": "WJ-native",
  1029. "fundamental_unit": f"individual gene (GTEx v8 RNA-seq, {n_genes_ba9} expressed genes in BA9)",
  1030. "pairwise_matrix": "genome-wide Spearman correlation, each target vs all genes",
  1031. "correlation_method": "Spearman",
  1032. "primary_analysis": "Weighted Jaccard on continuous correlation vectors",
  1033. "supplementary_analysis": "Binary Jaccard on top 5% network membership",
  1034. "fdr_scope": f"all 21 pairwise WJ permutation p-values (Benjamini-Hochberg)",
  1035. "permutations": N_PERMUTATIONS,
  1036. "domain_conventional_methods": "gProfiler GO enrichment (comparison), Fisher exact (custom sets)",
  1037. "random_seed": RANDOM_SEED,
  1038. "pipeline_file": "MS3_sigma_divergence_wj_pipeline.py",
  1039. "execution_date": datetime.now().strftime("%Y-%m-%d"),
  1040. "execution_time_seconds": round(time.time() - t_start, 1),
  1041. "wj_compliance_status": "PASS",
  1042. "brain_regions": list(BRAIN_REGIONS.keys()),
  1043. "n_samples_primary": n_samples_ba9,
  1044. "n_genes_primary": n_genes_ba9,
  1045. "key_results": {
  1046. "SIGMAR1_TMEM97_wj": float(st_wj['wj']),
  1047. "SIGMAR1_TMEM97_wj_z": float(st_wj['z_score']),
  1048. "SIGMAR1_TMEM97_wj_perm_p": float(st_wj['wj_perm_p']),
  1049. "SIGMAR1_TMEM97_binary_jaccard": float(st_binary['binary_jaccard']),
  1050. "SIGMAR1_TMEM97_shared_genes": int(st_binary['shared']),
  1051. "SIGMAR1_LTN1_wj": float(sl_wj['wj']),
  1052. "cross_region_rho_SIGMAR1": f"{cross_region_matrix['SIGMAR1'][np.triu_indices(5, k=1)].min():.3f}-{cross_region_matrix['SIGMAR1'][np.triu_indices(5, k=1)].max():.3f}",
  1053. "region_wj_values": {k: round(v['wj'], 4) for k, v in region_wj_results.items()},
  1054. "region_binary_j_values": {k: round(v['binary_j'], 3) for k, v in region_wj_results.items()},
  1055. },
  1056. }
  1058. provenance_path = os.path.join(RESULTS_DIR, "provenance.json")
  1059. with open(provenance_path, 'w') as f:
  1060. json.dump(provenance, f, indent=2)
  1061. print(f" provenance.json saved")
  1063. # ==================================================================
  1064. # SUMMARY
  1065. # ==================================================================
  1066. elapsed = time.time() - t_start
  1067. print(f"\n{'='*70}")
  1068. print("PIPELINE COMPLETE")
  1069. print(f"{'='*70}")
  1070. print(f"\n Time: {elapsed:.0f}s ({elapsed/60:.1f} min)")
  1071. print(f"\n--- PRIMARY: WEIGHTED JACCARD (continuous) ---")
  1072. print(f" SIGMAR1-TMEM97 WJ = {st_wj['wj']:.6f} (z={st_wj['z_score']:.2f}, p={st_wj['wj_perm_p']:.4f})")
  1073. print(f" SIGMAR1-LTN1 WJ = {sl_wj['wj']:.6f} (z={sl_wj['z_score']:.2f}, p={sl_wj['wj_perm_p']:.4f})")
  1074. print(f"\n--- SUPPLEMENTARY: BINARY JACCARD (top 5%) ---")
  1075. print(f" SIGMAR1-TMEM97 J = {st_binary['binary_jaccard']:.3f} (shared={int(st_binary['shared'])})")
  1076. print(f"\n--- MULTI-REGION WJ ---")
  1077. for rk, rv in region_wj_results.items():
  1078. print(f" {rk}: WJ={rv['wj']:.4f}, binary J={rv['binary_j']:.3f}")
  1080. return provenance
  1083. if __name__ == '__main__':
  1084. main()

MS3_sigma_divergence_wj_pipeline.py at commit bc9b836, under MIT · at the source

Overview

Authors: Drake H Harbert1
  1. Inner Architecture LLC, Canton, OH, United States
Journal: Frontiers in pharmacology, volume 17, article 1830847
Dates: received 14 March 2026; accepted 29 April 2026; published online 4 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3389/fphar.2026.1830847 · PMID 42328650 · PMCID PMC13275387 · OpenAlex W7163580477
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)
Methods: Statistics, Connectivity
Keywords: co-expression architecture, mitochondria-associated membrane, neurodegeneration, sigma-1 receptor, sigma-2 receptor, subtype-selective pharmacology, TMEM97, weighted jaccard
Topic: Pharmacological Receptor Mechanisms and Effects (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Citations: not cited yet (Europe PMC); 38 references in the paper

Abstract

The sigma-1 receptor (SIGMAR1) and sigma-2 receptor (TMEM97) are both enriched at the mitochondria-associated membrane (MAM) and have been pharmacologically co-classified for decades, yet their functional relationship at the transcriptomic level remains uncharacterized. We performed genome-wide co-expression analysis for both receptors across five brain regions from the GTEx v8 dataset (n = 209 in the primary region, 16,225 expressed genes) using Spearman correlations. Three Weighted Jaccard (WJ) formulations on continuous correlation vectors — (r+1)/2 shifted, unsigned |r|, and signed—all revealed that SIGMAR1 and TMEM97 share the majority of their global transcriptional architecture (WJ shifted = 0.964, unsigned = 0.907, signed = 0.906; all three rank-identical across 21 pairwise comparisons, ρ = 1.000), yet their top 5% co-expression networks overlap by only 10.0% (binary Jaccard = 0.100). Cosine similarity on raw vectors confirmed metric robustness (ρ = 0.856 with WJ, p = 7.5 × 10−7). Dissociation gap analysis across all 21 gene pairs showed the WJ-binary gap varies 2.1-fold (0.431–0.927), tracking known biological relatedness rather than reflecting a fixed methodological property. Gene Ontology analysis identified SIGMAR1-specific enrichment for mitochondrial translation and TCA cycle, and TMEM97-specific enrichment for ubiquitin-mediated proteolysis and neurodegeneration pathways. Multi-region replication confirmed the pattern across five brain regions, with the hippocampus showing tail-specific convergence. These findings are consistent with the hypothesis that dual sigma-1/sigma-2 ligands engage two functionally distinct co-expression programs within a shared cellular context, providing a transcriptomic rationale for subtype-selective pharmacological strategies.

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

nwharbert8-ui/sigma-receptor-divergence

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 4acadaaec414c5191a251cd54b0ef52cfcacaf8d, 14 February 2026
Languages: Python (1)
Size: 5 files, 1 script
Software Heritage: not archived
Found in: “Software and reproducibility”
Holds: README, license file, CITATION.cff, environment (requirements.txt)
Not found: tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
3 files
At the source: github.com/nwharbert8-ui/sigma-receptor-divergence

Zenodo 19024710

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data availability statement”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
3 files
At the source: doi.org/10.5281/zenodo.19024710

nwharbert8-ui/sigma-receptor-divergence-wj

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: bc9b836aae80e19aac860871508b205bb0957657, 14 March 2026
Languages: Python (1)
Size: 6 files, 1 script
Software Heritage: not archived
Found in: “Data availability statement”
Holds: README, license file, CITATION.cff, environment (requirements.txt)
Not found: tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
3 files
At the source: github.com/nwharbert8-ui/sigma-receptor-divergence-wj

The paper's code and data availability statement is in the Data section.

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;
  • 3 scripts, each with its path and the digest of its content;
  • 15 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

No dataset and no data link were found in the paper.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: GTEx v8 gene expression data (TPM) are available through the GTEx Portal (https://gtexportal.org/) under dbGaP accession phs000424.v8.p2. All analysis code is available at https://github.com/nwharbert8-ui/sigma-receptor-divergence-wj and archived at https://doi.org/10.5281/zenodo.19024710. Derived results (genome-wide correlation rankings, enrichment results, shared gene lists) are provided as Supplementary Data Sheets 1–4.

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, pages, dates, 1 author, 8 keywords, 38 references.

Cite

This paper

Harbert, D. H. (2026). Sigma-1 and Sigma-2 receptors exhibit divergent genome-wide Co-expression architectures in human brain despite shared subcellular localization. Frontiers in pharmacology, 17, 1830847. https://doi.org/10.3389/fphar.2026.1830847

BibTeX

@article{harbert2026sigma,
author = {Harbert, Drake H},
title = {{Sigma-1 and Sigma-2 receptors exhibit divergent genome-wide Co-expression architectures in human brain despite shared subcellular localization}},
journal = {Frontiers in pharmacology},
year = {2026},
month = jun,
volume = {17},
pages = {1830847},
publisher = {Frontiers Media SA},
issn = {1663-9812},
doi = {10.3389/fphar.2026.1830847},
url = {https://doi.org/10.3389/fphar.2026.1830847},
pmid = {42328650},
pmcid = {PMC13275387}
}

RIS

TY - JOUR
AU - Harbert, Drake H
TI - Sigma-1 and Sigma-2 receptors exhibit divergent genome-wide Co-expression architectures in human brain despite shared subcellular localization
T2 - Frontiers in pharmacology
J2 - Front Pharmacol
PY - 2026
DA - 2026/06/04
VL - 17
SP - 1830847
SN - 1663-9812
PB - Frontiers Media SA
DO - 10.3389/fphar.2026.1830847
UR - https://doi.org/10.3389/fphar.2026.1830847
LA - en
ER -

CSL-JSON

{
"id": "10.3389/fphar.2026.1830847",
"type": "article-journal",
"title": "Sigma-1 and Sigma-2 receptors exhibit divergent genome-wide Co-expression architectures in human brain despite shared subcellular localization",
"container-title": "Frontiers in pharmacology",
"author": [
{
"family": "Harbert",
"given": "Drake H"
}
],
"container-title-short": "Front Pharmacol",
"volume": "17",
"page": "1830847",
"DOI": "10.3389/fphar.2026.1830847",
"PMID": "42328650",
"PMCID": "PMC13275387",
"ISSN": "1663-9812",
"publisher": "Frontiers Media SA",
"URL": "https://doi.org/10.3389/fphar.2026.1830847",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
4
]
]
}
}

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.1007/s12035-026-06155-6
Sigma-1 Expression in Chronic Mouse Models of Temporal Lobe Epilepsy.
Journal: Molecular neurobiology
In common: cellular / molecular, 5 references
[2] doi:10.1016/j.isci.2026.116439 [code]
Decoding the role of transcriptomic clocks in the human prefrontal cortex.
Journal: iScience
In common: statsmodels, pandas, SciPy, 2 other tools, genetics / omics, 2 references
[3] doi:10.1038/s41380-026-03497-4 [code]
Transcriptome-informed brain cartography of polygenic risk and association with brain structure in major psychiatric disorders.
Journal: Molecular psychiatry
In common: statsmodels, pandas, SciPy, 2 other tools, genetics / omics, cellular / molecular, 2 references
[4] doi:10.3389/fgene.2026.1807347 [code]
Epigenetic regulators are preferentially coordinated with protocadherin gene expression across the human brain: a genome-wide co-expression analysis.
Journal: Frontiers in genetics
In common: statsmodels, pandas, SciPy, 2 other tools, genetics / omics, cellular / molecular, 2 references
[5] doi:10.1002/advs.202521254 [code]
Persistently Increased Expression of PKMzeta and Unbiased Gene Expression Profiles Identify Hippocampal Molecular Traces of a Long-Term Active Place Avoidance Memory and "Shadow" Proteins.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: statsmodels, pandas, SciPy, 2 other tools, genetics / omics, cellular / molecular, 2 references
[6] doi:10.1038/s41467-026-73428-y [code]
Regional heterogeneity in phenotypic and genetic associations between bone and brain in humans.
Journal: Nature communications
In common: pandas, SciPy, NumPy, genetics / omics, cellular / molecular, 3 references
[7] doi:10.1016/j.xgen.2026.101280 [code]
BMI-genome interactions regulate global gene expression with emphasis in brain and gut.
Journal: Cell genomics
In common: statsmodels, pandas, SciPy, 1 other tool, genetics / omics, cellular / molecular, 2 references
[8] doi:10.1038/s41467-026-75193-4 [code]
Multi-ancestry gene expression models amplify transcriptome-wide association study discovery and validation.
Journal: Nature communications
In common: statsmodels, pandas, SciPy, 1 other tool, genetics / omics, cellular / molecular, 2 references
[9] doi:10.1038/s41588-026-02646-3 [code]
Co-expression-based models improve eQTL predictions for transcriptome-wide association studies and highlight new schizophrenia-associated genes.
Journal: Nature genetics
In common: statsmodels, pandas, SciPy, 1 other tool, genetics / omics, cellular / molecular, 2 references
[10] doi:10.1038/s41380-026-03571-x [code]
Convergent coexpression reveals shared biological mechanisms underlying common and rare variant risk in six neuropsychiatric disorders.
Journal: Molecular psychiatry
In common: statsmodels, pandas, SciPy, 1 other tool, 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.