OSCR

Decoding the role of transcriptomic clocks in the human prefrontal cortex.

Code ↔ Paper

24 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 24 matches
  1. [1] § STAR★Methods › Method details › Functional analyses › MENTOR network embedding of functional gene connections ↔ 10other-07302025/03plots-08272025/00plots-08282025.ipynb, lines 139–234 · score 0.96 · deMagalhaes, DESeq2, StocH, StocP, StocZ, GenAge
  2. [2] § STAR★Methods › Method details › Stochastic modeling of transcriptomic clocks ↔ 09stochastic-08172025/0stochastic-0817025.ipynb, lines 424–463 · score 0.92 · small jitter, standard deviations, reduce perfect correlations, multivariate normal, covariance matrices, gene correlation
  3. [3] § STAR★Methods › Method details › Transcriptomic clocks analyses ↔ 10other-07302025/04plots-0828205/00plots-08282025.ipynb, lines 79–100 · score 0.89 · de Magalhaes, DESeq2, GenAge, RNA seq, GTExAge, brain tissue
  4. [4] § STAR★Methods › Method details › Transcriptomic clocks analyses ↔ 11functional-08202025/00functional_analysis-08202025.ipynb, lines 78–97 · score 0.89 · de Magalhaes, DESeq2, GenAge, RNA seq, GTExAge, brain tissue
  5. [5] § STAR★Methods › Method details › Deep learning age prediction ↔ 08training-08152025/01deeplearning/00deeplearnin-08152025.ipynb, lines 271–312 · score 0.82 · dropout rate, Keras, Glorot, loss, Adam, MAE
  6. [6] § STAR★Methods › Method details › Deep learning age prediction ↔ 09stochastic-08172025/0stochastic-0817025.ipynb, lines 946–987 · score 0.82 · dropout rate, Keras, Glorot, loss, Adam, MAE
  7. [7] § STAR★Methods › Method details › Elastic net age prediction ↔ 08training-08152025/00elasticnet/00elasticnet-08152025.ipynb, lines 374–413 · score 0.79 · ElasticNetCV, l1_ratio, fold predictions, split, metrics, scores
  8. [8] § STAR★Methods › Method details › Elastic net age prediction ↔ 09stochastic-08172025/0stochastic-0817025.ipynb, lines 786–825 · score 0.78 · ElasticNetCV, l1_ratio, fold predictions, split, metrics, scores
  9. [9] § STAR★Methods › Method details › Functional analyses › MENTOR network embedding of functional gene connections ↔ src/mentor/cli.py, lines 195–340 · score 0.74 · pairwise distance, MENTOR, Spearman, linkage, singleton, subclustering
  10. [10] § Results › Comparison of transcriptomic clocks and epigenetic clocks ↔ 10other-07302025/03plots-08272025/00plots-08282025.ipynb, lines 139–234 · score 0.73 · DamAge, PhenoAge, RNAAgeCalc, retrained transcriptomic clock, deep learning, Epigenetic clocks
  11. [11] § STAR★Methods › Method details › Gene annotation of epigenetic clock ↔ 10other-07302025/00plots-07302025/01rnaagecalc-08032025.ipynb, lines 264–321 · score 0.71 · StocH, StocP, StocZ, PhenoAge, epigenetic clock, Zhang2019
  12. [12] § Results › Performance of the transcriptomic clocks in predicting age ↔ 01GEOcohorts/01rnagecalc-08052025.R, lines 50–103 · score 0.69 · DESeq2, GenAge, GTExAge, RNAAgeCalc, Dev, Peters
  13. [13] § Results › Performance of the transcriptomic clocks in predicting age ↔ 01ragecalc-07272025/00rnagecalc-07272025.R, lines 41–106 · score 0.69 · DESeq2, GenAge, GTExAge, RNAAgeCalc, Dev, Peters
  14. [14] § Results › Stochastic effects of the transcriptomic clocks ↔ 10other-07302025/02plots-08272025/00plots-08272025.ipynb, lines 256–358 · score 0.65 · R2 stochastic, R2 deterministic, variance explained, deep learning, RR2, Elastic
  15. [15] § STAR★Methods › Experimental model and study participant details › University of Texas health science center at Houston (UTHealth) brain collection ↔ 07assoc-08052025/00uthealth-08052025.ipynb, lines 180–214 · score 0.65 · consensus diagnoses, UTHealth Brain, autopsy, ethnicity, drug, alcohol
  16. [16] § STAR★Methods › Experimental model and study participant details › University of Texas health science center at Houston (UTHealth) brain collection ↔ 07assoc-08052025/02uthealth-08282025.ipynb, lines 195–229 · score 0.65 · consensus diagnoses, UTHealth Brain, autopsy, ethnicity, drug, alcohol
  17. [17] § Results › Associations of transcriptomic age with psychiatric-related phenotypes ↔ 07assoc-08052025/04vabb-09012025.ipynb, lines 388–467 · score 0.64 · anticonvulsants, antidepressants, hallucinogen, polysubstance, antipsychotics, lifetime
  18. [18] § Results › Associations of transcriptomic age with psychiatric-related phenotypes ↔ 07assoc-08052025/01vabb-08062025.ipynb, lines 198–290 · score 0.63 · anticonvulsants, antidepressants, hallucinogen, polysubstance, antipsychotics, lifetime
  19. [19] § STAR★Methods › Quantification and statistical analysis ↔ 10other-07302025/01plots-08242025/00plots-08242025.ipynb, lines 618–742 · score 0.60 · squared error, Pearson correlation, Correlation matrices, R2, RMSE, retrained
  20. [20] § STAR★Methods › Quantification and statistical analysis ↔ 10other-07302025/02plots-08272025/00plots-08272025.ipynb, lines 256–358 · score 0.56 · variance explained, Pearson correlation, RR2, squared, stochastic, Prediction
  21. [21] § Results › Performance of transcriptomic clocks in predicting delta of age ↔ 10other-07302025/01plots-08242025/01plots-08262025.ipynb, lines 508–585 · score 0.55 · concordant decrease, delta age, density, deep learning, retrained, elastic
  22. [22] § STAR★Methods › Method details › Retraining transcriptomic clocks ↔ 08training-08152025/00elasticnet/00elasticnet-08152025.ipynb, lines 252–306 · score 0.54 · Benjamini Hochberg, FDR, sex, fit, models, gene
  23. [23] § STAR★Methods › Method details › DNA methylation analysis ↔ 11functional-08202025/01CpG2genelistepigenetics-08202025.R, lines 1–71 · score 0.53 · Human Methylation, EPIC, San, Illumina
  24. [24] § Results › Comparison of transcriptomic clocks and epigenetic clocks ↔ 10other-07302025/00plots-07302025/01rnaagecalc-08032025.ipynb, lines 264–321 · score 0.51 · DamAge, PhenoAge, Epigenetic clocks, cortical, Zhang2019, transcriptomic clock

Paper

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

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,144 lines · 38 KB · no license · 3 matches

  1. # %%
  2. import os, random
  3. import re
  4. import pyreadr
  5. import numpy as np
  6. import pandas as pd
  7. import warnings
  8. import optuna
  9. from scipy import stats
  10. import tensorflow as tf
  11. from typing import List, Tuple
  12. from sklearn.decomposition import PCA
  13. from sklearn.preprocessing import StandardScaler
  14. from sklearn.exceptions import ConvergenceWarning
  15. from scipy.interpolate import interp1d
  16. from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
  17. from sklearn.model_selection import KFold, cross_val_score
  18. from sklearn.linear_model import ElasticNetCV
  19. import statsmodels.api as sm
  20. import seaborn as sns
  21. import matplotlib.pyplot as plt
  22. from sklearn.covariance import LedoitWolf
  23. from statsmodels.tools.sm_exceptions import ConvergenceWarning, PerfectSeparationWarning
  24. from tensorflow.keras.models import Sequential
  25. from tensorflow.keras.models import load_model
  26. from tensorflow.keras import losses # For mse
  27. from tensorflow.keras.losses import MeanSquaredError
  28. from tensorflow.keras.layers import Dense, Dropout
  29. from tensorflow.keras.optimizers import Adam
  30. from tensorflow.keras.callbacks import EarlyStopping
  31. from sklearn.model_selection import train_test_split
  32. from sklearn.feature_selection import SelectFromModel
  33. from tensorflow.keras.models import load_model
  34. import tensorflow.keras.metrics as metrics
  35. import matplotlib.pyplot as plt
  36. from tensorflow.keras.models import Sequential
  37. from tensorflow.keras.layers import Dense, Dropout
  38. from tensorflow.keras.optimizers import Adam
  39. from tensorflow.keras.initializers import GlorotUniform
  41. # Shutdown warnings
  42. warnings.simplefilter("ignore", ConvergenceWarning)
  43. warnings.simplefilter("ignore", PerfectSeparationWarning)
  44. # Suppress all sklearn convergence warnings globally
  45. warnings.filterwarnings("ignore", category=ConvergenceWarning)
  46. warnings.filterwarnings("ignore", category=RuntimeWarning)
  48. # %%
  49. # -----------------------------------------------------------------
  50. # Set working directory
  51. wkdir = r"C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/00databases/01transcriptome/00data/"
  52. os.chdir(wkdir)
  54. # -----------------------------------------------------------------
  55. # File paths
  56. #uthbCounts = r"00uthb_genecounts-08152025.csv.gz"
  57. vabbCounts = r"00vabb_genecounts-08152025.csv.gz"
  59. # Reading annotation
  60. #uthbAnnot = r"00uthb_annotation-08152025.csv.gz"
  61. vabbAnnot = r"00vabb_annotation-08152025.csv.gz"
  63. # -----------------------------------------------------------------
  64. # Read gzipped CSVs
  65. #expr_uthb = pd.read_csv(uthbCounts, index_col=0, compression='gzip')
  66. expr_vabb = pd.read_csv(vabbCounts, index_col=0, compression='gzip')
  67. #annot_uthb = pd.read_csv(uthbAnnot, index_col=0, compression='gzip')
  68. annot_vabb = pd.read_csv(vabbAnnot, index_col=0, compression='gzip')
  70. # -----------------------------------------------------------------
  72. # %%
  73. # -----------------------------------------------------------------
  74. # Loading metadata
  75. # loading phenotype data
  76. phenoPath = "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/00databases/00phenotype/00vabb/00PhenoMeta_trainclock-07272025.csv"
  77. phenoData = pd.read_csv(phenoPath)
  78. # Create suffixes
  79. suffixes = ["_9", "_24", "_25", "_11"]
  80. # Replicate and modify pheno_data
  81. phenoData = pd.concat([
  82. phenoData.assign(SampleID="Sample" + phenoData["SampleID"].astype(str) + suffix)
  83. for suffix in suffixes
  84. ], ignore_index=True)
  85. phenoData.shape
  86. # -----------------------------------------------------------------
  88. # %%
  89. # Load cell types and SVAs
  90. cells = pd.read_csv("C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/05celltype-08042025/00vabb_cellprop_svas-08042025.csv")
  91. cells = cells.rename(columns={'name': 'SampleID'})
  92. cells.shape
  94. # %%
  95. # Merge in sequence after conversion
  96. merged_df_com = pd.merge(phenoData, cells, on='SampleID')
  97. meta = merged_df_com
  98. meta.shape
  100. # %%
  101. # -----------------------------------------------------------------
  102. # Subset meta to only required columns
  103. meta = meta[['SampleID', 'Sex', 'PMI', 'RIN', 'AgeDeath', 'ast', 'end', 'mic', 'neu', 'oli', 'opc', 'W_1', 'W_2']]
  104. # Columns not to scale
  105. exclude_cols = ['SampleID', 'Sex', 'AgeDeath']
  107. # Identify columns to scale
  108. cols_to_scale = [col for col in meta.columns if col not in exclude_cols]
  110. # Initialize scaler
  111. scaler = StandardScaler(with_mean=True, with_std=True)
  113. # Scale selected columns
  114. meta_scaled = meta.copy()
  115. meta_scaled[cols_to_scale] = scaler.fit_transform(meta_scaled[cols_to_scale])
  116. # Set SampleID as index
  117. meta_scaled = meta_scaled.set_index('SampleID')
  118. meta_scaled.head()
  120. # %%
  121. # -----------------------------------------------------------------
  122. # Function to clean sample names
  123. def clean_sample_name(name):
  124. # Match "Sample<number><letters><number>" and replace with "Sample<number>_<number>"
  125. return re.sub(r"(Sample\d+)[A-Za-z]+(\d+)", r"\1_\2", name)
  127. # Apply to a DataFrame
  128. # Example: expr_uthb.columns = [clean_sample_name(c) for c in expr_uthb.columns]
  129. expr_vabb.columns = [clean_sample_name(c) for c in expr_vabb.columns]
  130. expr_vabb.head()
  131. # -----------------------------------------------------------------
  133. # %%
  134. # -----------------------------------------------------------------
  135. # Define a function for the estimation FPKM of transcriptome data
  136. def counts_to_fpkm(counts: pd.DataFrame, gene_lengths: pd.Series, lengths_in: str = "bp") -> pd.DataFrame:
  137. """
  138. Convert raw counts to FPKM.
  139. counts: genes x samples
  140. gene_lengths: length per gene (index aligned to counts)
  141. lengths_in: 'bp' or 'kb'
  142. """
  143. assert set(counts.index).issubset(set(gene_lengths.index)), "Gene lengths missing for some genes"
  144. # Convert lengths to kilobases
  145. L = gene_lengths.loc[counts.index].astype(float)
  146. if lengths_in == "bp":
  147. L = L / 1e3
  148. elif lengths_in != "kb":
  149. raise ValueError("lengths_in must be 'bp' or 'kb'")
  151. # RPK: counts / length_kb
  152. rpk = counts.divide(L, axis=0)
  154. # FPKM: RPK / (total_mapped_reads_in_millions)
  155. per_sample_scaler = counts.sum(axis=0) / 1e6
  156. fpkm = rpk.divide(per_sample_scaler, axis=1)
  158. return fpkm
  160. # -----------------------------------------------------------------
  161. # Define a function for the estimation of gene lengths
  162. def extract_gene_bounds(df: pd.DataFrame, start_col='Start', end_col='End'):
  163. """
  164. Extract first start and last end for genes where Start/End columns may contain multiple positions separated by ';'.
  166. Args:
  167. df: DataFrame with at least 'GeneID', start_col, end_col
  168. start_col: name of the start column
  169. end_col: name of the end column
  171. Returns:
  172. DataFrame with 'GeneID', 'Start', 'End', 'GeneLength'
  173. """
  174. starts = []
  175. ends = []
  176. for s, e in zip(df[start_col], df[end_col]):
  177. # Split by ',' and convert to int
  178. start_vals = [int(x) for x in str(s).split(';') if x.strip().isdigit()]
  179. end_vals = [int(x) for x in str(e).split(';') if x.strip().isdigit()]
  180. if start_vals and end_vals:
  181. starts.append(min(start_vals)) # first start = smallest
  182. ends.append(max(end_vals)) # last end = largest
  183. else:
  184. starts.append(None)
  185. ends.append(None)
  187. result = pd.DataFrame({
  188. 'GeneID': df['GeneID'],
  189. 'Start': starts,
  190. 'End': ends
  191. })
  193. # Compute gene length
  194. result['GeneLength'] = result['End'] - result['Start'] + 1
  195. return result
  197. # %%
  198. # -----------------------------------------------------------------
  199. # Transforming to FPKM
  200. gene_info_vabb = extract_gene_bounds(annot_vabb)
  201. # Make Series: index = GeneID, values = GeneLength
  202. gene_lengths = pd.Series(gene_info_vabb["GeneLength"].values, index=gene_info_vabb["GeneID"])
  203. fpkm_matrix = counts_to_fpkm(expr_vabb, gene_lengths, lengths_in="bp")
  204. fpkm_matrix.shape
  206. # %%
  207. # -----------------------------------------------------------------
  208. # Filtering to FPKM
  209. # To avoid the influence of low count genes on the analysis result, genes with more than 30% samples having count per million (CPM) less than one were filtered out
  210. mfilter = (fpkm_matrix < 1).sum(axis=1) / fpkm_matrix.shape[1] <= 0.3
  211. fpkm_matrix = fpkm_matrix.loc[mfilter] # assign back
  212. fpkm_matrix.shape
  214. # %%
  215. # -----------------------------------------------------------------
  216. # Standardize by gene (rows)
  217. scaler = StandardScaler(with_mean=True, with_std=True)
  218. expr_z = pd.DataFrame(
  219. scaler.fit_transform(fpkm_matrix.T).T, # transpose -> scale -> transpose back
  220. index=fpkm_matrix.index, # keep gene names
  221. columns=fpkm_matrix.columns # keep sample names
  222. )
  223. expr_z.head()
  225. # %%
  226. # -----------------------------------------------------------------
  227. # Subsetting to only unique samples
  228. # Example: columns in your expression matrix
  229. columns = expr_z.columns # expr is your genes x samples DataFrame
  230. # Extract the prefix before the underscore
  231. prefixes = [col.split('_')[0] for col in columns]
  232. # Create a DataFrame to associate columns with prefixes
  233. col_df = pd.DataFrame({'col': columns, 'prefix': prefixes})
  235. # Set random seed for reproducibility
  236. seed = 42
  237. np.random.seed(seed)
  238. # Sample one column per prefix
  239. unique_cols = col_df.groupby('prefix')['col'].apply(lambda x: np.random.choice(x)).values
  240. # Subset your expression matrix
  241. expr_unique = expr_z[unique_cols]
  242. # Check
  243. print(expr_unique.shape)
  244. print(expr_z.shape)
  246. # %%
  247. # -----------------------------------------------------------------
  248. # Columns that were selected
  249. unique_cols_set = set(unique_cols)
  250. # Columns that were not selected
  251. excluded_cols = [col for col in expr_z.columns if col not in unique_cols_set]
  252. # Check
  253. print("Number of selected columns:", len(unique_cols))
  254. print("Number of excluded columns:", len(excluded_cols))
  255. # Optional: print first 10 excluded columns
  256. print("Excluded columns:", excluded_cols[:10])
  258. # %%
  259. # ------------------------------
  260. # Align samples
  261. common_samples = expr_unique.columns.intersection(meta_scaled.index)
  262. expr_unique = expr_unique[common_samples]
  263. meta_scaled2 = meta_scaled.loc[common_samples]
  264. print(expr_unique.shape)
  265. print(meta_scaled2.shape)
  267. # Plotting Age
  268. plt.hist(meta_scaled2['AgeDeath'], bins=30, edgecolor='black')
  269. plt.xlabel("Age")
  270. plt.ylabel("Count")
  271. plt.title("Age distribution")
  272. plt.show()
  274. # %%
  275. # ------------------------------
  276. # Keep only samples <30 or >60
  277. meta_subset = meta_scaled2[(meta_scaled2['AgeDeath'] < 30) | (meta_scaled2['AgeDeath'] > 60)].copy()
  278. # Add new variable: 0 = young (<30), 1 = old (>60)
  279. meta_subset['age_group'] = (meta_subset['AgeDeath'] > 60).astype(int)
  280. # Subset expression matrix to the same samples
  281. expr_subset = expr_unique[meta_subset.index]
  282. print(meta_subset.shape, expr_subset.shape)
  283. meta_subset[['AgeDeath', 'age_group']].head()
  285. # %%
  286. # ------------------------------
  287. # Keep only genes in clock
  288. coefPath = "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/06tranningclock-08152025/00elasticnet-08152025/02elastic_net_nonzero_coefficients-08152025.csv"
  289. # Load coefficients
  290. df_coef = pd.read_csv(coefPath, index_col=0)
  291. # Convert the single-column DataFrame to a Series
  292. coef_series = pd.Series(df_coef['coef'].values, index=df_coef.index, name='coef')
  293. # Ensure it's a Series
  294. coef_series = pd.Series(coef_series, name='coef')
  295. coef = df_coef
  296. genes_to_keep = coef.index
  297. expr_subset = expr_subset.loc[genes_to_keep]
  298. expr_subset.shape
  300. # %%
  301. # ------------------------------
  302. # Prepare data
  303. # ------------------------------
  304. # Keep covariates of interest
  305. covariates = meta_subset[['Sex', 'PMI', 'RIN']].copy()
  307. # One-hot encode Sex if necessary
  308. covariates = pd.get_dummies(covariates, columns=['Sex'], drop_first=True)
  310. # Add intercept
  311. Xcov_only = sm.add_constant(covariates).astype(float)
  313. # Outcome variable: 0 = young, 1 = old
  314. y = meta_subset['age_group'].values
  316. # Align samples with expression matrix
  317. expr_use = expr_subset
  318. common_samples = expr_use.columns.intersection(meta_subset.index)
  319. expr_use = expr_use[common_samples]
  320. Xcov_only = Xcov_only.loc[common_samples]
  321. y = y[np.isin(meta_subset.index, common_samples)]
  323. # ------------------------------
  324. # Fit null logistic model
  325. # ------------------------------
  326. null_model = sm.Logit(y, Xcov_only).fit(disp=0)
  327. p_hat = null_model.predict()
  329. # W matrix: diagonal of p*(1-p)
  330. W = np.diag(p_hat * (1 - p_hat))
  332. # ------------------------------
  333. # Compute z-scores for each gene
  334. # ------------------------------
  335. results = []
  337. for gene_id in expr_use.index:
  338. x = expr_use.loc[gene_id].values.reshape(-1,1) # already standardized
  340. # If gene has zero variance, set z=0
  341. if np.var(x) == 0:
  342. results.append((gene_id, 0.0))
  343. continue
  345. # Wald z-score: z = (x^T (y - p_hat)) / sqrt(x^T W x)
  346. numerator = x.T @ (y - p_hat)
  347. denominator = np.sqrt(x.T @ W @ x)
  348. z = float(numerator / denominator)
  350. results.append((gene_id, z))
  352. # ------------------------------
  353. # Compile results
  354. # ------------------------------
  355. deg_z = pd.DataFrame(results, columns=['gene_id','z_expr']).set_index('gene_id')
  356. deg_z.head()
  358. # %%
  359. # Plot histogram of the first few genes (deg_z.head())
  360. plt.figure(figsize=(6,4))
  361. plt.hist(deg_z['z_expr'], bins=10, color='skyblue', edgecolor='black')
  362. plt.xlabel('Z-score')
  363. plt.ylabel('Frequency')
  364. plt.title('Histogram of Z-scores')
  365. plt.show()
  367. # %%
  368. # ------------------------------
  369. # Start stochastic modelling
  370. # ------------------------------
  371. # Define helper functions
  372. def cov_from_corr_sd(corr: pd.DataFrame, sd: pd.Series) -> pd.DataFrame:
  373. """Covariance = D * R * D, with D=diag(sd)."""
  374. D = np.diag(sd.values)
  375. cov = D @ corr.values @ D
  376. return pd.DataFrame(cov, index=corr.index, columns=corr.columns)
  379. def shrink_corr_from_matrix(X: pd.DataFrame) -> pd.DataFrame:
  380. """
  381. Estimate correlation via Ledoit-Wolf shrinkage
  382. from an expression matrix (samples x genes).
  383. X = samples x genes
  384. """
  385. Xc = X - X.mean(axis=0)
  386. lw = LedoitWolf().fit(Xc.values)
  387. cov = pd.DataFrame(lw.covariance_, index=X.columns, columns=X.columns)
  388. d = np.sqrt(np.diag(cov.values))
  389. corr = cov / d[:, None] / d[None, :]
  390. return corr
  392. # %%
  393. # ------------------------------
  394. # Convert z-scores → age-specific means
  395. # ------------------------------
  396. def z_to_means(genes, z_age80_vs_25, space="z",
  397. baseline_mean_age25=None, baseline_sd_age25=None, se_per_gene=None):
  398. """
  399. space = "z": assume standardized space
  400. - age 25 mean = 0, sd = 1
  401. - age 80 mean = z
  402. space = "raw": use delta = z * SE
  403. """
  404. z_age80_vs_25 = z_age80_vs_25.reindex(genes)
  406. if space == "z":
  407. mean25 = pd.Series(0.0, index=genes)
  408. sd25 = pd.Series(1.0, index=genes) if baseline_sd_age25 is None else baseline_sd_age25.reindex(genes)
  409. mean80 = z_age80_vs_25.copy()
  410. sd80 = sd25.copy()
  411. elif space == "raw":
  412. if se_per_gene is None or baseline_mean_age25 is None:
  413. raise ValueError("Need baseline_mean_age25 and se_per_gene for raw mode")
  414. mean25 = baseline_mean_age25.reindex(genes)
  415. sd25 = pd.Series(1.0, index=genes) if baseline_sd_age25 is None else baseline_sd_age25.reindex(genes)
  416. delta = z_age80_vs_25 * se_per_gene.reindex(genes)
  417. mean80 = mean25 + delta
  418. sd80 = sd25
  419. else:
  420. raise ValueError("space must be 'z' or 'raw'")
  422. return mean25, sd25, mean80, sd80
  424. # %%
  425. # ------------------------------
  426. # Adjusted MVN simulation to reduce excessive collinearity
  427. # ------------------------------
  428. def simulate_mvn(n, mean, corr=None, sd=None, jitter=1e-3):
  429. """
  430. Simulate multivariate normal samples for genes with slight jitter to reduce
  431. perfect correlations for better convergence in ElasticNet.
  433. Parameters
  434. ----------
  435. n : int
  436. Number of samples.
  437. mean : pd.Series
  438. Mean per gene.
  439. corr : pd.DataFrame, optional
  440. Gene-gene correlation matrix.
  441. sd : pd.Series, optional
  442. Standard deviation per gene.
  443. jitter : float, optional
  444. Small value to add to diagonal of covariance to reduce collinearity.
  445. """
  446. genes = mean.index
  447. p = len(genes)
  449. # Default identity correlation / unit SD
  450. if corr is None:
  451. corr = pd.DataFrame(np.eye(p), index=genes, columns=genes)
  452. if sd is None:
  453. sd = pd.Series(1.0, index=genes)
  455. # Covariance matrix from correlation and SD
  456. cov = cov_from_corr_sd(corr, sd)
  458. # Add small jitter to diagonal to reduce perfect correlation
  459. cov += np.eye(p) * jitter
  461. # Draw samples
  462. X = np.random.multivariate_normal(mean.values, cov.values, size=n)
  463. return pd.DataFrame(X, columns=genes)
  465. # %%
  466. # ------------------------------
  467. # Define a Multivariate Normal simulation
  468. # ------------------------------
  469. def negbin_params_from_mean_k(mu, k):
  470. n_param = k
  471. p_param = k / (k + mu)
  472. return n_param, p_param
  474. def gaussian_copula_transform(n, corr):
  475. R = corr.values
  476. Z = np.random.multivariate_normal(mean=np.zeros(len(R)), cov=R, size=n)
  477. U = stats.norm.cdf(Z)
  478. return U
  480. def simulate_gaussian_copula_negbin(n, corr, mean_counts, k_overdisp):
  481. genes = mean_counts.index
  482. U = gaussian_copula_transform(n, corr)
  483. out = np.empty_like(U, dtype=np.int64)
  485. for j, g in enumerate(genes):
  486. mu = float(mean_counts[g])
  487. k = float(k_overdisp[g])
  488. n_param, p_param = negbin_params_from_mean_k(mu, k)
  489. out[:, j] = stats.nbinom.ppf(U[:, j], n=n_param, p=p_param).astype(int)
  491. return pd.DataFrame(out, columns=genes)
  493. # %%
  494. # ------------------------------
  495. # Running simulations for ages 20 and 60
  496. # ------------------------------
  497. # This generates
  498. genes = deg_z.index
  499. z = deg_z
  501. # Convert z-scores → age-specific means in z-space
  502. mean20, sd20, mean60, sd60 = z_to_means(genes, z, space="z")
  504. # Correlation of genes
  505. expr_use_T = expr_use.T # shape = (samples, genes)
  506. corr = shrink_corr_from_matrix(expr_use_T)
  508. # Ensure 1D Series and alignment
  509. mean20 = mean20.squeeze()
  510. mean60 = mean60.squeeze()
  511. sd20 = sd20.squeeze()
  512. sd60 = sd60.squeeze()
  514. genes_common = mean20.index.intersection(corr.index)
  515. mean20 = mean20.loc[genes_common]
  516. mean60 = mean60.loc[genes_common]
  517. sd20 = sd20.loc[genes_common]
  518. sd60 = sd60.loc[genes_common]
  519. corr = corr.loc[genes_common, genes_common]
  521. # ---- MVN samples ----
  522. samples20 = simulate_mvn(2000, mean20, corr=corr, sd=sd20)
  523. samples60 = simulate_mvn(2000, mean60, corr=corr, sd=sd60)
  525. samples20["age"] = 20
  526. samples60["age"] = 60
  527. mvn_samples = pd.concat([samples20, samples60], ignore_index=True)
  529. mvn_samples.head()
  531. # %%
  532. # ------------------------------
  533. # Running Gaussian Copula–NB simulations for ages 20 and 60
  534. # ------------------------------
  535. # This generates COUNTS
  536. # If z is DataFrame with one column
  537. z = z.squeeze() # now a Series
  539. # Example: map z-scores to mean counts (toy example)
  540. # You can adjust this to your real counts
  541. baseline_mu = pd.Series(100.0, index=genes_common) # baseline counts
  542. fold_change = np.exp(z[genes_common] * 0.2) # simple mapping z -> fold change
  543. mu20 = baseline_mu # age 20
  544. mu60 = baseline_mu * fold_change # age 60
  546. # Overdispersion parameter k (can be gene-specific)
  547. k20 = pd.Series(10.0, index=genes_common)
  548. k60 = pd.Series(10.0, index=genes_common)
  550. # Simulate
  551. samples20_nb = simulate_gaussian_copula_negbin(2000, corr, mu20, k20)
  552. samples60_nb = simulate_gaussian_copula_negbin(2000, corr, mu60, k60)
  554. # Add age labels
  555. samples20_nb["age"] = 20
  556. samples60_nb["age"] = 60
  558. # Concatenate
  559. copula_nb_samples = pd.concat([samples20_nb, samples60_nb], ignore_index=True)
  560. copula_nb_samples.head()
  562. # %%
  563. # ------------------------------
  564. # PCA analysis of stochastic genes
  565. # Extract gene expression only
  566. X = mvn_samples.drop(columns=["age"]).values
  567. ages = mvn_samples["age"].values
  569. pca = PCA(n_components=2)
  570. pcs = pca.fit_transform(X)
  572. # Wrap into DataFrame
  573. pca_df = pd.DataFrame(pcs, columns=["PC1", "PC2"])
  574. pca_df["age"] = ages
  576. # Plotting
  577. plt.figure(figsize=(8,6))
  578. sns.scatterplot(data=pca_df, x="PC1", y="PC2", hue="age", palette="coolwarm", alpha=0.3)
  579. plt.title("PCA trajectory of ages 20 → 60")
  580. plt.xlabel(f"PC1 ({pca.explained_variance_ratio_[0]*100:.1f}% variance)")
  581. plt.ylabel(f"PC2 ({pca.explained_variance_ratio_[1]*100:.1f}% variance)")
  582. plt.show()
  584. # %%
  585. # ------------------------------
  586. # Compute mean expression per gene for each age group
  587. mean_expr = mvn_samples.groupby(['age']).mean().T # Transpose to have genes as rows
  588. mean_expr = mean_expr.loc[:, [20, 60]] # Select only ages 20 and 60
  590. # Scatter plots of stochastic genes
  591. plt.figure(figsize=(8, 6))
  593. # Draw a line for each gene
  594. for gene in mean_expr.index:
  595. plt.plot([20, 60], [mean_expr.loc[gene, 20], mean_expr.loc[gene, 60]],
  596. color='gray', alpha=0.5, lw=0.8)
  598. # Add points on top
  599. plt.scatter([20]*len(mean_expr), mean_expr[20], color='blue', label='Age 20')
  600. plt.scatter([60]*len(mean_expr), mean_expr[60], color='red', label='Age 60')
  602. plt.xlabel('Age')
  603. plt.ylabel('Mean Gene Expression')
  604. plt.title('Gene-level trajectories: Age 20 → 60')
  605. plt.legend()
  606. plt.grid(True)
  607. plt.show()
  609. # %%
  610. # ------------------------------
  611. # Scatter plots of expression per sample between ages 20 and 60
  612. # ------------------------------
  613. # ------------------------------
  614. # Scatter plots of expression per sample (no sample_id)
  615. # ------------------------------
  616. # Keep only age 20 and 60
  617. samples_20_60 = mvn_samples[mvn_samples['age'].isin([20, 60])]
  619. # Split by age
  620. samples_20 = samples_20_60[samples_20_60['age'] == 20].drop(columns=['age'])
  621. samples_60 = samples_20_60[samples_20_60['age'] == 60].drop(columns=['age'])
  623. plt.figure(figsize=(8, 6))
  625. # Loop over rows (each row = one sample’s expression vector)
  626. for i in range(min(len(samples_20), len(samples_60))):
  627. plt.plot([20, 60],
  628. [samples_20.iloc[i].mean(), samples_60.iloc[i].mean()],
  629. color='gray', alpha=0.5, lw=0.8)
  631. # Scatter points
  632. plt.scatter([20]*len(samples_20), samples_20.mean(axis=1),
  633. color='blue', label='Age 20')
  634. plt.scatter([60]*len(samples_60), samples_60.mean(axis=1),
  635. color='red', label='Age 60')
  637. plt.xlabel('Age')
  638. plt.ylabel('Sample-level Mean Gene Expression')
  639. plt.title('Sample-level trajectories: Age 20 → 60')
  640. plt.legend()
  641. plt.grid(True)
  642. plt.show()
  644. # %%
  645. # ------------------------------
  646. # Age range and interpolation
  647. # ------------------------------
  648. ages = [20, 30, 40, 50, 60]
  650. # Interpolate gene means for each age
  651. mean_df = pd.DataFrame(index=genes_common, columns=ages, dtype=float)
  652. for gene in genes_common:
  653. mean_df.loc[gene] = np.linspace(mean20[gene], mean60[gene], num=len(ages))
  655. # Optional: interpolate SDs across ages
  656. sd_df = pd.DataFrame(index=genes_common, columns=ages, dtype=float)
  657. for gene in genes_common:
  658. sd_df.loc[gene] = np.linspace(sd20[gene], sd60[gene], num=len(ages))
  660. # ------------------------------
  661. # Generate MVN samples for each age
  662. # ------------------------------
  663. n_samples_per_age = 2000
  664. all_samples = []
  666. for age in ages:
  667. mean_vec = mean_df[age].astype(float)
  668. sd_vec = sd_df[age].astype(float)
  670. samples_age = simulate_mvn(n_samples_per_age, mean_vec, corr=corr, sd=sd_vec)
  671. samples_age["age"] = age
  672. all_samples.append(samples_age)
  674. # Combine all ages
  675. mvn_samples_multi_age = pd.concat(all_samples, ignore_index=True)
  676. mvn_samples_multi_age.head()
  678. # %%
  679. # Compute mean expression per gene for each age group
  680. mean_expr = mvn_samples_multi_age.groupby('age').mean().T # genes as rows
  682. # Optional: select a subset of ages to plot
  683. ages_to_plot = [20, 30, 40, 50, 60]
  684. mean_expr = mean_expr.loc[:, ages_to_plot]
  686. # Plot
  687. plt.figure(figsize=(10, 6))
  689. # Draw a line for each gene
  690. for gene in mean_expr.index:
  691. plt.plot(ages_to_plot, mean_expr.loc[gene, ages_to_plot],
  692. color='gray', alpha=0.5, lw=0.8)
  694. # Add points on top
  695. for age in ages_to_plot:
  696. plt.scatter([age]*len(mean_expr), mean_expr[age],
  697. label=f'Age {age}' if age in [20, 60] else "", # label only 20 & 60
  698. s=20,
  699. color='blue' if age==20 else 'red' if age==60 else 'green')
  701. plt.xlabel('Age')
  702. plt.ylabel('Mean Gene Expression')
  703. plt.title('Gene-level trajectories: Age 20 → 60')
  704. plt.legend()
  705. plt.grid(True)
  706. plt.show()
  708. # %%
  709. # ------------------------------
  710. # Sample-level trajectories across ages
  711. # ------------------------------
  712. ages_to_plot = [20, 30, 40, 50, 60]
  713. samples_multi_age = mvn_samples_multi_age[mvn_samples_multi_age['age'].isin(ages_to_plot)]
  715. plt.figure(figsize=(10, 6))
  717. # Pivot: rows = sample index, cols = ages, values = mean expression
  718. sample_means = samples_multi_age.drop(columns=['age']).mean(axis=1) # per-sample mean
  719. sample_means = pd.DataFrame({
  720. "age": samples_multi_age['age'].values,
  721. "mean_expr": sample_means.values
  722. })
  724. # Group by index position (pseudo-IDs) across ages
  725. for idx in range(sample_means.shape[0] // len(ages_to_plot)):
  726. subset = sample_means.iloc[idx::(sample_means.shape[0] // len(ages_to_plot))]
  727. plt.plot(subset["age"], subset["mean_expr"], color="gray", alpha=0.5, lw=0.8)
  729. # Scatter points on top
  730. for age in ages_to_plot:
  731. subset = sample_means[sample_means['age'] == age]
  732. plt.scatter([age]*len(subset), subset['mean_expr'],
  733. color='blue' if age==20 else 'red' if age==60 else 'green',
  734. alpha=0.7, s=20, label=f'Age {age}' if age in [20, 60] else "")
  736. plt.xlabel('Age')
  737. plt.ylabel('Sample-level Mean Gene Expression')
  738. plt.title('Sample-level gene expression trajectories')
  739. plt.legend()
  740. plt.grid(True)
  741. plt.show()
  743. # %%
  744. # ------------------------------
  745. # Align samples
  746. # ------------------------------
  747. common_samples = expr_use.columns.intersection(meta_subset.index)
  748. expr_use_aligned = expr_use[common_samples]
  749. age_series = meta_subset.loc[common_samples, 'AgeDeath']
  751. # ------------------------------
  752. # Define age bins
  753. # ------------------------------
  754. bins = [20, 30, 40, 50, 60, 70, 80] # adjust as needed
  755. age_bins = pd.cut(age_series, bins=bins, right=False)
  757. # Compute mean expression per gene per age bin
  758. mean_expr = expr_use_aligned.T.groupby(age_bins).mean().T # genes x bins
  760. # Compute numeric bin centers
  761. bin_centers = [interval.left + (interval.right - interval.left)/2 for interval in mean_expr.columns]
  763. # Ensure numeric values and fill missing
  764. mean_expr_numeric = mean_expr.copy()
  765. mean_expr_numeric = mean_expr_numeric.astype(float).interpolate(axis=1) # fill NaNs if any
  767. # Plot
  768. plt.figure(figsize=(10, 6))
  770. for gene in mean_expr_numeric.index:
  771. y = mean_expr_numeric.loc[gene, :].values # numeric values across bins
  772. plt.plot(bin_centers, y, color='gray', alpha=0.5, lw=0.8)
  774. # Overlay points for each bin
  775. for i, center in enumerate(bin_centers):
  776. plt.scatter([center]*len(mean_expr_numeric),
  777. mean_expr_numeric.iloc[:, i].values,
  778. color='blue', s=15)
  780. plt.xlabel('Age')
  781. plt.ylabel('Mean Gene Expression')
  782. plt.title('Gene-level trajectories across AgeDeath')
  783. plt.grid(True)
  784. plt.show()
  786. # %%
  787. # ------------------------------
  788. # Model training
  789. # ------------------------------
  790. def train_elastic_net(
  791. z_expr: pd.DataFrame,
  792. y_age: pd.Series,
  793. selected_genes: List[str],
  794. n_splits: int = 5,
  795. random_state: int = 42,
  796. l1_ratios: Tuple[float, ...] = (0.1, 0.3, 0.5, 0.7, 0.9)
  797. ) -> Tuple[ElasticNetCV, dict]:
  798. """Fit ElasticNetCV on selected genes, return model and CV metrics."""
  800. # Align genes and transpose so samples are rows
  801. X = z_expr.loc[selected_genes].T
  802. y = y_age.loc[X.index] # ensure samples match
  804. # Initialize CV
  805. cv = KFold(n_splits=n_splits, shuffle=True, random_state=random_state)
  807. # Fit ElasticNetCV on all data
  808. model = ElasticNetCV(l1_ratio=l1_ratios, alphas=None, cv=cv, max_iter=1000)
  809. model.fit(X, y)
  811. # Out-of-fold predictions
  812. preds = np.zeros(len(y), dtype=float)
  813. for train_idx, test_idx in cv.split(X):
  814. Xtr, Xte = X.iloc[train_idx], X.iloc[test_idx]
  815. ytr = y.iloc[train_idx]
  816. m = ElasticNetCV(l1_ratio=l1_ratios, alphas=None, cv=cv, max_iter=1000)
  817. m.fit(Xtr, ytr)
  818. preds[test_idx] = m.predict(Xte)
  820. # CV metrics
  821. r2 = r2_score(y, preds)
  822. mae = mean_absolute_error(y, preds)
  823. metrics = {"cv_r2": float(r2), "cv_mae": float(mae)}
  825. return model, metrics
  827. # %%
  828. # ------------------------------
  829. # Model training
  830. # ------------------------------
  831. # Features: genes x samples
  832. # Drop the 'age' column from simulated data
  833. X_sim = mvn_samples_multi_age.drop(columns=['age'])
  834. y_sim = mvn_samples_multi_age['age']
  836. # Features: genes x samples
  837. X_sim_selected = X_sim[genes] # genes are currently columns
  838. X_sim_selected = X_sim_selected.T # now genes are rows, samples are columns
  840. # Train the model
  841. model, metrics = train_elastic_net(
  842. z_expr=X_sim_selected, # or X_sim if using all genes
  843. y_age=y_sim,
  844. selected_genes=genes, # can be None to use all genes
  845. n_splits=5,
  846. random_state=42,
  847. l1_ratios=(0.1, 0.3, 0.5, 0.7, 0.9)
  848. )
  850. # %%
  851. # ------------------------------
  852. # Plotting metrics
  853. # ------------------------------
  854. # Convert to lists for plotting
  855. names = list(metrics.keys())
  856. values = list(metrics.values())
  858. # Create bar plot
  859. plt.figure(figsize=(5, 4))
  860. plt.bar(names, values, color=['skyblue', 'salmon'])
  861. plt.ylabel("Value")
  862. plt.title("ElasticNet CV Metrics")
  863. plt.ylim(0, max(values)*1.2) # add some space on top
  864. plt.show()
  866. # %%
  867. # ------------------------------
  868. # Saving Coefficients
  869. # ------------------------------
  870. # Coefficients as Series (gene names as index)
  871. coef_series = pd.Series(model.coef_, index=genes, name='coef')
  872. coef_series.to_csv("C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/00elasticnet-08152025/01elastic_net_coefficients-08152025.csv", header=True)
  874. # Optional only non-zero coeficients
  875. nonzero_coef = coef_series[coef_series != 0]
  876. nonzero_coef.to_csv("C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/00elasticnet-08152025/02elastic_net_nonzero_coefficients-08152025.csv", header=True)
  878. # Print intercept
  879. print("Intercept:", model.intercept_)
  881. # Save intercept
  882. pd.Series({"intercept": model.intercept_}).to_csv(
  883. "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/00elasticnet-08152025/03elastic_net_intercept-08152025.csv"
  884. )
  886. # %%
  887. # ------------------------------
  888. # Predicting with Coefficients
  889. # ------------------------------
  890. # Expression matrix for excluded samples
  891. expr_excluded = expr_z
  892. # Make sure the same genes are used (selected_genes)
  893. X_excluded = expr_excluded.loc[genes]
  895. # Load coefficients
  896. df_coef = pd.read_csv("C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/00elasticnet-08152025/02elastic_net_nonzero_coefficients-08152025.csv", index_col=0)
  897. # Convert the single-column DataFrame to a Series
  898. coef_series = pd.Series(df_coef['coef'].values, index=df_coef.index, name='coef')
  899. # Ensure it's a Series
  900. coef_series = pd.Series(coef_series, name='coef')
  901. # Subset new expression matrix to selected genes
  902. X_new = X_excluded.loc[coef_series.index].T # samples x genes
  904. # Load the intercept CSV
  905. df_intercept = pd.read_csv(
  906. "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/00elasticnet-08152025/03elastic_net_intercept-08152025.csv",
  907. index_col=0
  908. )
  910. # Extract the value
  911. intercept = df_intercept.loc["intercept"].values[0]
  913. # Dot product to get predicted age
  914. predicted_age = X_new.dot(coef_series) + intercept
  916. # Convert to Series
  917. predicted_age = pd.Series(predicted_age, index=X_new.index, name='Predicted_Age')
  919. # %%
  920. # If predicted_age is a Pandas Series
  921. plt.figure(figsize=(6,4))
  922. plt.hist(predicted_age, bins=30, edgecolor="black")
  923. plt.xlabel("Predicted Age")
  924. plt.ylabel("Frequency")
  925. plt.title("Histogram of Predicted Age")
  926. plt.show()
  928. # %%
  929. # Now align with predicted_age
  930. common_ids = predicted_age.index.intersection(meta_scaled.index)
  931. ages = meta_scaled.loc[common_ids, "AgeDeath"]
  932. pred = predicted_age.loc[common_ids]
  934. # Correlation
  935. corr = pred.corr(ages)
  936. print(f"Correlation between predicted and AgeDeath: {corr:.3f}")
  938. # Scatter plot
  939. plt.figure(figsize=(5,5))
  940. plt.scatter(ages, pred, alpha=0.6, edgecolor="k")
  941. plt.xlabel("Age at Death")
  942. plt.ylabel("Predicted Age")
  943. plt.title(f"Predicted vs Actual Age at Death\nCorrelation = {corr:.3f}")
  944. plt.show()
  946. # %%
  947. # ------------------------------
  948. # Training Deep Learning Model
  949. # ------------------------------
  950. seed = 42
  951. os.environ['PYTHONHASHSEED'] = str(seed)
  952. os.environ['TF_DETERMINISTIC_OPS'] = '1'
  953. np.random.seed(seed)
  954. random.seed(seed)
  955. tf.random.set_seed(seed)
  957. # Define a flexible Keras model
  958. def create_model(trial, input_dim):
  959. n_layers = trial.suggest_int("n_layers", 1, 4)
  960. model = Sequential()
  962. for i in range(n_layers):
  963. units = trial.suggest_int(f"units_l{i}", 32, 512, step=32)
  964. dropout_rate = trial.suggest_float(f"dropout_l{i}", 0.0, 0.5)
  966. if i == 0:
  967. model.add(Dense(
  968. units, activation='relu', input_dim=input_dim,
  969. kernel_initializer=GlorotUniform(seed=42)
  970. ))
  971. else:
  972. model.add(Dense(
  973. units, activation='relu',
  974. kernel_initializer=GlorotUniform(seed=42)
  975. ))
  977. # Fix dropout randomness
  978. model.add(Dropout(dropout_rate, seed=42))
  980. # Output layer
  981. model.add(Dense(1, activation='linear', kernel_initializer=GlorotUniform(seed=42)))
  983. lr = trial.suggest_float("learning_rate", 1e-4, 1e-2, log=True)
  984. optimizer = Adam(learning_rate=lr)
  985. model.compile(optimizer=optimizer, loss='mse', metrics=['mae'])
  987. return model
  989. # %%
  990. # ------------------------------
  991. # Training deep models from the Elastic Net with cross-validation
  992. # ------------------------------
  993. # Setting seeds for replicability
  994. seed = 42
  995. np.random.seed(seed)
  997. # Adjust shape of samples
  998. X_sim_selected = X_sim[genes] # genes are currently columns
  1000. # Specify folder to save models
  1001. # change this to your desired path
  1002. MODEL_SAVE_PATH = "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/01deeplearning-08162025/"
  1003. os.makedirs(MODEL_SAVE_PATH, exist_ok=True)
  1005. # Define a function
  1006. def objective(trial):
  1007. # Create model with trial hyperparameters
  1008. model = create_model(trial, input_dim=X_sim_selected.shape[1])
  1010. # Early stopping
  1011. es = EarlyStopping(monitor='val_loss', patience=10, restore_best_weights=True)
  1013. # Train/validation split
  1014. X_train, X_val, y_train, y_val = train_test_split(X_sim_selected, y_sim, test_size=0.2, random_state=42)
  1016. # Train model
  1017. history = model.fit(
  1018. X_train, y_train,
  1019. validation_data=(X_val, y_val),
  1020. epochs=100,
  1021. batch_size=trial.suggest_categorical("batch_size", [16, 32, 64, 128]),
  1022. callbacks=[es],
  1023. verbose=0
  1024. )
  1026. # Save the trained model in the specified folder
  1027. model_filename = os.path.join(MODEL_SAVE_PATH, f"trial_{trial.number}_model.h5")
  1028. model.save(model_filename)
  1029. print(f"Saved model for trial {trial.number} as {model_filename}")
  1031. # Return validation metric for optimization
  1032. val_mae = min(history.history['val_mae'])
  1033. return val_mae
  1035. # Run Optuna optimization
  1036. study = optuna.create_study(direction="minimize", sampler=optuna.samplers.TPESampler(seed=42))
  1037. study.optimize(objective, n_trials=50)
  1039. # Print best trial
  1040. print("Best trial:")
  1041. trial = study.best_trial
  1042. print(trial.params)
  1044. # %%
  1045. # Build path to best model
  1046. best_model_filename = os.path.join(MODEL_SAVE_PATH, f"trial_{study.best_trial.number}_model.h5")
  1047. print(best_model_filename)
  1048. # Print best trial
  1049. print("Best trial:")
  1050. trial = study.best_trial
  1051. print(trial.params)
  1053. # %%
  1054. # ------------------------------
  1055. # Prepare data for prediction
  1056. # ------------------------------
  1057. # Filter expr_excluded to only these genes
  1058. expr_excluded_filtered = expr_z[excluded_cols].loc[genes.intersection(expr_z.index)]
  1060. # Transpose to samples x genes
  1061. X_excluded = expr_excluded_filtered.T
  1062. X_excluded.shape
  1064. # %%
  1065. # ------------------------------
  1066. # Load the best model (without compiling)
  1067. # ------------------------------
  1068. best_model_filename = "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/01deeplearning-08162025/trial_47_model.h5"
  1069. best_model = load_model(best_model_filename, compile=False)
  1070. print(f"Loaded best model from {best_model_filename}")
  1072. # ------------------------------
  1073. # Re-compile the model with standard loss/metrics
  1074. # ------------------------------
  1075. best_model.compile(optimizer='adam', loss='mse', metrics=['mae'])
  1077. # ------------------------------
  1078. # Save under a new name
  1079. # ------------------------------
  1080. new_model_filename = "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/00Stochasticdeepclock-08192025.h5"
  1081. best_model.save(new_model_filename)
  1082. print(f"Saved model as {new_model_filename}")
  1084. # ------------------------------
  1085. # Predict on new data
  1086. # ------------------------------
  1087. # X_excluded should have the same shape/features as the training data
  1088. predicted_values = best_model.predict(X_excluded)
  1089. print("Predictions shape:", predicted_values.shape)
  1091. # ------------------------------
  1092. # Predict
  1093. # ------------------------------
  1094. predicted_values = best_model.predict(X_excluded)
  1095. print(predicted_values.shape)
  1097. # %%
  1098. # If predicted_age is a Pandas Series
  1099. plt.figure(figsize=(6,4))
  1100. plt.hist(predicted_values, bins=30, edgecolor="black")
  1101. plt.xlabel("Predicted Age")
  1102. plt.ylabel("Frequency")
  1103. plt.title("Histogram of Predicted Age")
  1104. plt.show()
  1106. # %%
  1107. # Assume X_excluded has samples as rows and their index contains the sample IDs
  1108. pred = pd.Series(predicted_values.flatten(), index=X_excluded.index, name="PredictedAge")
  1110. # Align with actual ages
  1111. common_ids = pred.index.intersection(meta_scaled.index)
  1112. ages = meta_scaled.loc[common_ids, "AgeDeath"]
  1113. pred = pred.loc[common_ids]
  1116. # Correlation
  1117. corr = pred.corr(ages)
  1118. print(f"Correlation between predicted and AgeDeath: {corr:.3f}")
  1120. plt.figure(figsize=(5,5))
  1121. plt.scatter(ages, pred, alpha=0.6, edgecolor="k")
  1122. plt.xlabel("Age at Death")
  1123. plt.ylabel("Predicted Age")
  1124. plt.title(f"Predicted vs Actual Age at Death\nCorrelation = {corr:.3f}")
  1125. plt.show()
  1127. # %%
  1128. # Saving predictions
  1129. # Create a DataFrame with sample IDs and predicted age
  1130. common_index = predicted_age.index.intersection(pred.index)
  1132. df_pred = pd.DataFrame({
  1133. "SampleID": common_index,
  1134. "Predicted_Age_stochastic_elastic": predicted_age.loc[common_index].values,
  1135. "Predicted_Age_stochastic_deep": pred.loc[common_index].values
  1136. })
  1138. # Save to CSV
  1139. df_pred.to_csv(
  1140. "C:/Users/jjm262/OneDrive - Yale University/Documents/Documents/00yale/04fourthyear/01projects/03tclock/02results/07stochasticclock-08172025/04vabb_prediction-08152025.csv",
  1141. index=False
  1142. )
  1144. print("Predicted ages saved to 'predicted_age_samples.csv'.")

0stochastic-0817025.ipynb at commit a4dc4b1, no license · at the source

Overview

Authors: José J Martínez-Magaña1,2, Anna HC Vlot3, Kyle A Sullivan3, Daniel A Jacobson3, John H Krystal1,2,4, Matthew J Girgenti1,2, Diana L Núnez-Ríos1,2, Sheila T Nagamatsu1,2, Diego E Andrade-Brito1,2, Jean Merlet3,5, Alice Townsend3,5, Alana Wells6, Christiane Alvarez5, Matthew Lane5,6, Paul E Holtzheimer7, Traumatic Stress Brain Research Group, Consuelo Walss-Bass8,9, Janitza L Montalvo-Ortiz1,2,3
ORCID iDs: José J Martínez-Magaña, Janitza L Montalvo-Ortiz
  1. Division of Human Genetics, Department of Psychiatry, Yale University School of Medicine, New Haven, CT, USA
  2. National Center for PTSD, US Department of Veterans Affairs, West Haven, CT, USA
  3. Biosciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA
  4. Psychiatry Service, VA Connecticut Health Care System, West Haven, CT, USA
  5. The University of Tennessee, Knoxville, TN, USA
  6. Rhodes College, Memphis, TN, USA
  7. Department of Psychiatry, Geisel School of Medicine at Dartmouth, Lebanon, NH 03756, USA
  8. Louis A. Faillace, MD, Department of Psychiatry and Behavioral Sciences, McGovern Medical School, University of Texas Health Science Center at Houston, Houston, TX, USA
  9. MD Anderson Cancer Center, University of Texas Health Science Center at Houston Graduate School of Biomedical Sciences, Houston, TX, USA
Institutions: United States Department of Veterans Affairs (United States); Yale University (United States); National Center for PTSD (United States); National Center for Post Traumatic Stress Disorder (United States); Oak Ridge National Laboratory (United States); VA Connecticut Healthcare System (United States); University of Tennessee at Knoxville (United States); Rhodes College (United States); Dartmouth College (United States); The University of Texas MD Anderson Cancer Center (United States); The University of Texas Health Science Center at Houston (United States)
Journal: iScience, volume 29, issue 7, article 116439
Dates: received 1 May 2024; accepted 2 June 2026; published online 2 July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.isci.2026.116439 · PMID 42491807 · PMCID PMC13378135 · OpenAlex W4366998610
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2023.04.19.23288765
Status: code verified
Categories: genetics / omics (modality), human (organism)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Connectivity
Keywords: Biological sciences, Molecular biology, Neuroscience, Molecular neuroscience, Bioinformatics, Computational bioinformatics, Genomic analysis
Topic: Bioinformatics and Genomic Networks (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: CSRD VA (IK2 CX002095); National Institute of Mental Health (MH120170); National Institute on Drug Abuse (DP1DA058737, R21DA050160); National Center for PTSD, U.S. Department of Veterans Affairs; NIDA NIH HHS (R21 DA050160, DP1 DA058737); Office of Science; U.S. Department of Energy (DE-AC05-00OR22725); US Department of Veterans Affairs (1IK2CX002095-01A1); NIMH NIH HHS (R01 MH120170)
Citations: cited by 1 paper (Europe PMC); 125 references in the paper

Abstract

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

Repositories

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

martinezjaime/PFC-Transcriptomic

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: a4dc4b14a9abe71a0592da6aeb67085ad3e86026, 20 April 2026
Languages: Jupyter (26), R (11), Shell (6)
Size: 51 files, 43 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, documentation, 26 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration
Tools: pandas (23 files), NumPy (22 files), scikit-learn (20 files), Matplotlib (17 files), SciPy (17 files), seaborn (14 files), statsmodels (7 files), data.table (6 files), tidyverse (6 files), glmnet (4 files), Keras (4 files), TensorFlow (4 files), NetworkX (2 files), ggplot2 (1 file), patchwork (1 file), STAR (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
44 files
At the source: github.com/martinezjaime/PFC-Transcriptomic

Jacobson-CompSysBio/MENTOR-py

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 9bea8a4963bccf5d5842c72dcc6f6207b164a752, 17 July 2024
Languages: Python (14), R (8)
Size: 32 files, 22 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Holds: README, license file, CITATION.cff, environment (environment.yml, pyproject.toml, setup.cfg), tests
Not found: continuous integration, documentation
Tools: NumPy (12 files), pandas (9 files), scikit-learn (8 files), SciPy (8 files), tidyverse (7 files), data.table (2 files), igraph (2 files), circlize (1 file), ComplexHeatmap (1 file), cowplot (1 file), patchwork (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
24 files
At the source: github.com/Jacobson-CompSysBio/MENTOR-py

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:

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

Code and data availability statement

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

Read it in the paper: doi.org/10.1016/j.isci.2026.116439.

Versions

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

Version 2, 28 September 2026

  • Authors: added José J Martínez-Magaña (0000-0003-0390-8252); Janitza L Montalvo-Ortiz (0000-0002-7657-8365); removed José J Martínez-Magaña; Janitza L Montalvo-Ortiz

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 18 authors, 7 keywords, 9 funders, 125 references.

Cite

This paper

Martínez-Magaña, J. J., Vlot, A. H., Sullivan, K. A., Jacobson, D. A., Krystal, J. H., Girgenti, M. J., Núnez-Ríos, D. L., Nagamatsu, S. T., Andrade-Brito, D. E., Merlet, J., Townsend, A., Wells, A., Alvarez, C., Lane, M., Holtzheimer, P. E., Traumatic Stress Brain Research Group, Walss-Bass, C., & Montalvo-Ortiz, J. L. (2026). Decoding the role of transcriptomic clocks in the human prefrontal cortex. iScience, 29(7), 116439. https://doi.org/10.1016/j.isci.2026.116439

BibTeX

@article{martinezmagana2026decoding,
author = {Martínez-Magaña, José J and Vlot, Anna HC and Sullivan, Kyle A and Jacobson, Daniel A and Krystal, John H and Girgenti, Matthew J and Núnez-Ríos, Diana L and Nagamatsu, Sheila T and Andrade-Brito, Diego E and Merlet, Jean and Townsend, Alice and Wells, Alana and Alvarez, Christiane and Lane, Matthew and Holtzheimer, Paul E and {Traumatic Stress Brain Research Group} and Walss-Bass, Consuelo and Montalvo-Ortiz, Janitza L},
title = {{Decoding the role of transcriptomic clocks in the human prefrontal cortex}},
journal = {iScience},
year = {2026},
month = jul,
volume = {29},
number = {7},
pages = {116439},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.116439},
url = {https://doi.org/10.1016/j.isci.2026.116439},
pmid = {42491807},
pmcid = {PMC13378135}
}

RIS

TY - JOUR
AU - Martínez-Magaña, José J
AU - Vlot, Anna HC
AU - Sullivan, Kyle A
AU - Jacobson, Daniel A
AU - Krystal, John H
AU - Girgenti, Matthew J
AU - Núnez-Ríos, Diana L
AU - Nagamatsu, Sheila T
AU - Andrade-Brito, Diego E
AU - Merlet, Jean
AU - Townsend, Alice
AU - Wells, Alana
AU - Alvarez, Christiane
AU - Lane, Matthew
AU - Holtzheimer, Paul E
AU - Traumatic Stress Brain Research Group
AU - Walss-Bass, Consuelo
AU - Montalvo-Ortiz, Janitza L
TI - Decoding the role of transcriptomic clocks in the human prefrontal cortex
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/07/02
VL - 29
IS - 7
SP - 116439
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.116439
UR - https://doi.org/10.1016/j.isci.2026.116439
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.116439",
"type": "article-journal",
"title": "Decoding the role of transcriptomic clocks in the human prefrontal cortex",
"container-title": "iScience",
"author": [
{
"family": "Martínez-Magaña",
"given": "José J"
},
{
"family": "Vlot",
"given": "Anna HC"
},
{
"family": "Sullivan",
"given": "Kyle A"
},
{
"family": "Jacobson",
"given": "Daniel A"
},
{
"family": "Krystal",
"given": "John H"
},
{
"family": "Girgenti",
"given": "Matthew J"
},
{
"family": "Núnez-Ríos",
"given": "Diana L"
},
{
"family": "Nagamatsu",
"given": "Sheila T"
},
{
"family": "Andrade-Brito",
"given": "Diego E"
},
{
"family": "Merlet",
"given": "Jean"
},
{
"family": "Townsend",
"given": "Alice"
},
{
"family": "Wells",
"given": "Alana"
},
{
"family": "Alvarez",
"given": "Christiane"
},
{
"family": "Lane",
"given": "Matthew"
},
{
"family": "Holtzheimer",
"given": "Paul E"
},
{
"literal": "Traumatic Stress Brain Research Group"
},
{
"family": "Walss-Bass",
"given": "Consuelo"
},
{
"family": "Montalvo-Ortiz",
"given": "Janitza L"
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "7",
"page": "116439",
"DOI": "10.1016/j.isci.2026.116439",
"PMID": "42491807",
"PMCID": "PMC13378135",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.116439",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
2
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41514-026-00391-9 [code]
Region-specific transcriptional signatures of brain aging in the absence of neuropathology at the single-cell level.
Journal: npj aging
In common: circlize, ComplexHeatmap, cowplot, 12 other tools, genetics / omics, 5 references
[2] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: STAR, igraph, circlize, 14 other tools, genetics / omics, 3 references
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, Keras, igraph, 16 other tools
[4] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Keras, igraph, circlize, 15 other tools, genetics / omics, 1 reference
[5] doi:10.1038/s41598-026-48613-0 [code]
An snRNA-seq aging clock for the fruit fly head sheds light on sex-biased aging.
Journal: Scientific reports
In common: Keras, TensorFlow, seaborn, 5 other tools, genetics / omics, 8 references
[6] doi:10.1126/sciadv.aeg3223 [code]
The extreme diversity of retinal amacrine cells has deep evolutionary roots.
Journal: Science advances
In common: glmnet, igraph, circlize, 13 other tools, genetics / omics, 1 reference
[7] doi:10.1093/neuonc/noag128 [code]
Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response.
Journal: Neuro-oncology
In common: glmnet, Keras, circlize, 12 other tools, genetics / omics
[8] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: igraph, circlize, NetworkX, 12 other tools, genetics / omics, 2 references
[9] doi:10.1038/s41467-026-76676-0 [code]
Determinants of functional burden pleiotropy and gene dosage responses across human traits.
Journal: Nature communications
In common: igraph, circlize, ComplexHeatmap, 13 other tools, genetics / omics
[10] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: STAR, igraph, NetworkX, 12 other tools, genetics / omics, 1 reference

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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