OSCR

Network-based discovery of regulatory drivers of cognitive decline in alzheimer's disease.

Code ↔ Paper

7 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 7 matches
  1. [1] § Results › Differential expression analysis ↔ src/code_for_lioness_networks_analysis/Sankey_Plot/plots(box_plot_sankey_plots_network_plots).ipynb, lines 1261–1262 · score 0.90 · TNF alpha, sodium symporter activity, hormone activity, extracellular matrix, amino acid, microvillus
  2. [2] § Results › Network metrics differences in AD-specific and NCI-specific GRNs ↔ src/code_for_lioness_networks_analysis/calculating_graph_features/nci_netwoks_features_calculation.ipynb, lines 35–83 · score 0.71 · network properties, network metrics, connected component, largest, regulatory networks, modularity
  3. [3] § Results › Network metrics differences in AD-specific and NCI-specific GRNs ↔ src/code_for_lioness_networks_analysis/calculating_graph_features/ad_netwoks_features_calculation.ipynb, lines 58–107 · score 0.71 · network properties, network metrics, connected component, largest, regulatory networks, modularity
  4. [4] § Methods › Pathway and cell type-specific enrichment analysis ↔ src/code_for_lioness_networks_analysis/Sankey_Plot/plots(box_plot_sankey_plots_network_plots).ipynb, lines 389–436 · score 0.67 · GO Cellular Component, GO Molecular Function, KEGG, Hallmark, Human, genes
  5. [5] § Methods › Network analysis of AD-specific and NCI-specific GRNs ↔ src/code_for_lioness_networks_analysis/calculating_graph_features/nci_netwoks_features_calculation.ipynb, lines 35–83 · score 0.62 · network metrics, properties, largest, regulatory networks, modularity, nestedness
  6. [6] § Methods › Network analysis of AD-specific and NCI-specific GRNs ↔ src/code_for_lioness_networks_analysis/calculating_graph_features/ad_netwoks_features_calculation.ipynb, lines 58–107 · score 0.62 · network metrics, properties, largest, regulatory networks, modularity, nestedness
  7. [7] § Methods › Data preprocessing ↔ src/code_for_data_preprocessing/combat_seq_batch_removal.ipynb, lines 1–15 · score 0.53 · Combat Seq, batch, preprocessing

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,329 lines · 36 KB · no license · 2 matches

  1. # %%
  2. import pandas as pd
  3. import os
  4. import networkx as nx
  5. from tqdm import tqdm
  7. # %%
  8. preprocessed_ad_single_sample_network_dir_path = '../../../data_preprocessed/grn_output/ad_nci_seperate_network/lioness/ad_preprocessed/'
  9. preprocessed_nci_single_sample_network_dir_path = '../../../data_preprocessed/grn_output/ad_nci_seperate_network/lioness/nci_preprocessed/'
  12. # %%
  13. AD_data = pd.read_csv(preprocessed_ad_single_sample_network_dir_path +'lioness.427_120507.58.csv').rename(columns={'2':'lioness.427_120507.58.csv'.split('.')[1]})
  15. for i in tqdm(os.listdir(preprocessed_ad_single_sample_network_dir_path)):
  16. if i.split('.')[1] in AD_data.columns:
  17. continue
  18. if 'csv' in i:
  19. df = pd.read_csv(preprocessed_ad_single_sample_network_dir_path + i).rename(columns={'2':i.split('.')[1]})
  20. AD_data = pd.merge(AD_data,df, on=['0','1'])
  23. # %%
  24. NCI_data = pd.read_csv(preprocessed_nci_single_sample_network_dir_path +'lioness.713_120531.47.csv').rename(columns={'2':'lioness.713_120531.47.csv'.split('.')[1]})
  26. for i in tqdm(os.listdir(preprocessed_nci_single_sample_network_dir_path)):
  27. if i.split('.')[1] in NCI_data.columns:
  28. continue
  29. if 'csv' in i:
  30. df = pd.read_csv(preprocessed_nci_single_sample_network_dir_path + i).rename(columns={'2':i.split('.')[1]})
  31. NCI_data = pd.merge(NCI_data,df, on=['0','1'])
  33. # %%
  34. AD_data.index = AD_data['0'] + '_' + AD_data['1']
  35. AD_data.drop(columns=['0','1'],inplace = True)
  36. AD_data
  38. # %%
  39. NCI_data.index = NCI_data['0'] + '_' + NCI_data['1']
  40. NCI_data.drop(columns=['0','1'],inplace = True)
  41. NCI_data
  43. # %%
  44. common_edge = set(NCI_data.index).intersection(set(AD_data.index))
  45. AD_data = AD_data.loc[list(common_edge),:]
  46. NCI_data = NCI_data.loc[list(common_edge),:]
  48. # %%
  49. AD_data = AD_data.T
  50. AD_data
  52. # %%
  53. NCI_data = NCI_data.T
  54. NCI_data
  56. # %%
  57. combined = pd.concat([AD_data,NCI_data])
  59. # %%
  60. Label = ['AD' for i in range(len(AD_data))] + ['NCI' for i in range(len(NCI_data))]
  61. Label
  63. # %%
  64. import pandas as pd
  65. from sklearn.decomposition import PCA
  66. import matplotlib.pyplot as plt
  67. import seaborn as sns
  69. # Generate sample data (replace this with your own DataFrame)
  71. combined['Label'] = Label
  72. # Assuming your DataFrame contains only numeric values
  73. # If not, you might need to preprocess your data accordingly
  75. # Perform PCA
  76. pca = PCA(n_components=2)
  77. principal_components = pca.fit_transform(combined.drop(columns = 'Label'))
  79. # Create a new DataFrame with the principal components
  80. pc_df = pd.DataFrame(data=principal_components, columns=['PC1', 'PC2'])
  82. # Concatenate the principal components DataFrame with the original DataFrame
  83. final_df = pd.concat([pc_df], axis=1)
  84. final_df['Label'] = Label
  85. final_df.index = combined.index
  88. # Plot the PCA results
  89. plt.figure(figsize=(8, 8))
  90. sns.scatterplot(x='PC1', y='PC2', data=final_df, hue = Label)
  91. plt.title('PCA Plot')
  92. plt.xlabel('Principal Component 1 (PC1)')
  93. plt.ylabel('Principal Component 2 (PC2)')
  94. plt.show()
  97. # %%
  98. remove = final_df[final_df.PC1> 80].index.tolist()
  99. remove
  101. # %%
  102. outlier_corrected_combined = combined[~combined.index.isin(remove)]
  104. # %%
  105. import pandas as pd
  106. from sklearn.decomposition import PCA
  107. import matplotlib.pyplot as plt
  108. import seaborn as sns
  110. # Generate sample data (replace this with your own DataFrame)
  113. # Assuming your DataFrame contains only numeric values
  114. # If not, you might need to preprocess your data accordingly
  116. # Perform PCA
  117. pca = PCA(n_components=2)
  118. principal_components = pca.fit_transform(outlier_corrected_combined .drop(columns = 'Label'))
  120. # Create a new DataFrame with the principal components
  121. pc_df = pd.DataFrame(data=principal_components, columns=['PC1', 'PC2'])
  123. # Concatenate the principal components DataFrame with the original DataFrame
  124. final_df = pd.concat([pc_df], axis=1)
  125. final_df['Label'] = outlier_corrected_combined.Label.tolist()
  126. final_df.index = outlier_corrected_combined.index
  127. Label = outlier_corrected_combined.Label.tolist()
  129. # Plot the PCA results
  130. plt.figure(figsize=(8, 8))
  131. sns.scatterplot(x='PC1', y='PC2', data=final_df, hue = Label )
  132. plt.title('PCA Plot')
  133. plt.xlabel('Principal Component 1 (PC1)')
  134. plt.ylabel('Principal Component 2 (PC2)')
  135. plt.show()
  138. # %%
  139. l = [
  140. 'ENSG00000166823_ENSG00000178199',
  141. 'ENSG00000157557_ENSG00000121903',
  143. 'ENSG00000028277_ENSG00000039139',
  144. 'ENSG00000081059_ENSG00000255561',
  145. 'ENSG00000167981_ENSG00000143222',
  146. 'ENSG00000198482_ENSG00000156219',
  147. 'ENSG00000188786_ENSG00000083799',
  148. 'ENSG00000168874_ENSG00000139197',
  149. 'ENSG00000139800_ENSG00000150175',
  150. 'ENSG00000079432_ENSG00000214425',
  151. 'ENSG00000162992_ENSG00000139197',
  152. 'ENSG00000256294_ENSG00000272576',
  153. 'ENSG00000169635_ENSG00000234841',
  154. 'ENSG00000125850_ENSG00000198783',
  155. 'ENSG00000078900_ENSG00000062096',
  156. 'ENSG00000120837_ENSG00000123179',
  157. 'ENSG00000079432_ENSG00000188234',
  158. 'ENSG00000188786_ENSG00000162384',
  159. 'ENSG00000251369_ENSG00000265018',
  160. 'ENSG00000234284_ENSG00000149089',
  161. 'ENSG00000198429_ENSG00000107262',
  162. 'ENSG00000188785_ENSG00000122641']
  164. # %%
  165. import pickle
  167. # Path to your pickle file
  168. pickle_file_path = '../../../data_preprocessed/pickle_file/ensemble_to_gene_name.pickle'
  170. # Load the pickle file
  171. with open(pickle_file_path, 'rb') as file:
  172. ensemble_to_gene_name = pickle.load(file)
  175. # %%
  177. def return_gene_name(gene):
  179. try:
  180. if ensemble_to_gene_name[gene] == 'Gene name not found':
  181. return gene
  182. return ensemble_to_gene_name[gene]
  183. except:
  184. return gene
  187. # %%
  188. n = []
  189. for i in l:
  190. m = []
  191. for j in i.split('_'):
  193. m.append(return_gene_name(j))
  195. n.append(m[0] + '-' + m[1])
  199. # %%
  200. n
  202. # %%
  203. l.append('Label')
  204. l
  206. # %%
  207. important_gene_dataframe = combined[l]
  208. important_gene_dataframe.head()
  210. # %%
  211. import pandas as pd
  212. import seaborn as sns
  213. import matplotlib.pyplot as plt
  214. from sklearn.preprocessing import MinMaxScaler
  216. # Set the Seaborn theme for better aesthetics
  217. sns.set_theme(style="whitegrid", palette="muted")
  219. # Assuming 'important_gene_dataframe' is your DataFrame and it has a 'Label' column
  221. # Normalize the entire DataFrame (excluding the 'Label' column) to a 0-1 range
  222. scaler = MinMaxScaler()
  223. important_gene_dataframe_normalized = important_gene_dataframe.copy()
  225. # Apply normalization (only to the gene expression columns, not the 'Label' column)
  226. gene_columns = important_gene_dataframe_normalized.drop(columns=['Label']).columns
  227. important_gene_dataframe_normalized[gene_columns] = scaler.fit_transform(important_gene_dataframe_normalized[gene_columns])
  229. # Filter DataFrames for 'AD' and 'NCI'
  230. ad_data = important_gene_dataframe_normalized[important_gene_dataframe_normalized['Label'] == 'AD']
  231. nci_data = important_gene_dataframe_normalized[important_gene_dataframe_normalized['Label'] == 'NCI']
  233. # Drop the 'Label' column and only keep the gene expression values
  234. ad_data = ad_data.drop(columns=['Label'])
  235. nci_data = nci_data.drop(columns=['Label'])
  237. # List of genes (columns) to iterate over
  238. genes = ad_data.columns
  240. # Number of genes
  241. num_genes = len(genes)
  243. # Determine the number of rows and columns for subplots (depending on the number of genes)
  244. ncols = 12 # Set the number of columns (adjust based on your preference)
  245. nrows = (num_genes + ncols - 1) // ncols # Automatically adjust the number of rows
  247. # Create a figure and axes for subplots
  248. fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(50, 10 * nrows))
  249. axes = axes.flatten() # Flatten axes array for easy indexing
  251. # Loop through each gene and create a boxplot on a subplot
  252. for idx, gene in enumerate(genes):
  253. ax = axes[idx]
  255. # Prepare the data for the current gene (AD vs NCI)
  256. data = pd.DataFrame({
  257. 'Group': ['AD'] * len(ad_data[gene]) + ['NCI'] * len(nci_data[gene]),
  258. 'Expression': list(ad_data[gene]) + list(nci_data[gene])
  259. })
  261. # Create a boxplot for the current gene
  262. sns.boxplot(data=data, x='Group', y='Expression', palette='pastel', ax=ax)
  263. gene = n[idx]
  264. # Add a title and custom labels
  265. ax.set_title(f'{gene}', fontsize=14, weight='bold')
  266. ax.set_xlabel('Group', fontsize=12)
  267. ax.set_ylabel('Normalized Expression (0-1)', fontsize=12)
  269. # Add gridlines for better readability
  270. ax.grid(True, axis='y', linestyle='--', alpha=0.6)
  272. # Remove any empty subplots if the number of genes is not a perfect multiple of ncols
  273. for idx in range(num_genes, len(axes)):
  274. fig.delaxes(axes[idx])
  276. # Adjust layout for better spacing between subplots
  277. plt.tight_layout()
  279. # Save the figure as an SVG with high resolution
  280. fig.savefig('gene_expression_boxplots.svg', format='svg', dpi=600)
  282. # Show the plot
  283. plt.show()
  286. # %%
  287. l.remove('Label')
  288. plt.figure(figsize=(8, 6)) # Optional: adjust the figure size
  289. sns.heatmap(important_gene_dataframe[l], annot=False, cmap='YlGnBu') # 'annot=True' adds annotations inside the heatmap cells
  290. plt.title('Heatmap of DataFrame')
  291. plt.show()
  294. # %%
  295. AD_panda = pd.read_csv('../../../data_preprocessed/grn_output/ad_nci_seperate_network/panda/output_panda_ad_csv.csv')
  297. # %%
  298. AD_panda.sort_values(by = 'force', ascending=False)
  300. # %%
  302. import pandas as pd
  303. import seaborn as sns
  304. import matplotlib.pyplot as plt
  306. import numpy as np
  307. # Set the plot size
  308. plt.figure(figsize=(8, 6))
  310. # Plot percentile plot using seaborn
  311. sns.ecdfplot(AD_panda['force'], color='skyblue')
  314. percentile_95 = np.percentile(AD_panda['force'], 95)
  316. plt.axvline(x=percentile_95, color='red', linestyle='--', label='90th Percentile')
  317. plt.axhline(y=0.9, color='green', linestyle='--')
  319. # Annotate the point on the plot
  320. plt.annotate(f'({percentile_95:.2f}, 0.95)', xy=(percentile_95, 0.9), xytext=(percentile_95 + 0.1, 0.9),
  321. arrowprops=dict(facecolor='black', arrowstyle='->'))
  322. # Display the plot
  323. plt.title('Percentile Plot (CDF)')
  324. plt.xlabel('Data')
  325. plt.ylabel('Cumulative Probability')
  326. plt.show()
  328. # %%
  329. AD_panda_top_5 = AD_panda[AD_panda.force>=3.71]
  330. AD_panda_top_5
  332. # %%
  333. list_of_tf = list(set(AD_panda_top_5.tf))
  334. AD_panda_top_5_bipartite_complete = AD_panda_top_5[~AD_panda_top_5.gene.isin(list_of_tf)]
  336. # %%
  337. import pandas as pd
  338. import networkx as nx
  339. import matplotlib.pyplot as plt
  340. from tqdm import tqdm
  343. # Create a weighted bipartite graph
  344. AD_BG_bipartite = nx.Graph()
  345. AD_BG_bipartite.add_nodes_from(AD_panda_top_5_bipartite_complete['tf'], bipartite = 0)
  346. AD_BG_bipartite.add_nodes_from(AD_panda_top_5_bipartite_complete['gene'], bipartite = 1)
  348. # Add weighted edges
  349. edges = [(row['tf'], row['gene'], {'weight': row['force']}) for idx, row in tqdm(AD_panda_top_5_bipartite_complete.iterrows())]
  350. AD_BG_bipartite.add_edges_from(edges)
  353. # %%
  354. relevant_genes = l.copy()
  356. # %%
  357. relevant_genes
  359. # %%
  360. l = []
  361. for i in relevant_genes:
  362. l.append(i.split('_')[0])
  363. l.append(i.split('_')[1])
  365. # %%
  366. import gseapy as gp
  367. from gseapy import Msigdb
  368. names = gp.get_library_name(organism='Human')
  370. # %%
  371. l[0:5]
  373. # %%
  374. def return_gene_name(gene):
  376. try:
  377. if ensemble_to_gene_name[gene] == 'Gene name not found':
  378. return gene
  379. return ensemble_to_gene_name[gene]
  380. except:
  381. return gene
  383. # %%
  384. 0%2 != 0
  386. # %%
  387. 2%2 != 0
  389. # %%
  390. df = pd.DataFrame()
  391. data_dict = dict()
  392. for i,k in tqdm(enumerate(l)):
  393. if i%2 != 0:
  394. continue
  396. m = []
  397. try:
  398. #print(k)
  399. #print(l[i+1])
  400. m.append(k)
  401. m.extend(list(AD_BG_bipartite[k]))
  402. m.extend(list(AD_BG_bipartite[l[i+1]]))
  403. except:
  404. print(k)
  407. m = set(m)
  408. print(len(m))
  410. gene_set = [return_gene_name(gene_id) for gene_id in m]
  411. gene_set = list(set(gene_set))
  413. try:
  414. gene_set.remove('Gene name not found')
  415. except:
  416. print('skip')
  418. filtered_gene_set = []
  419. k = return_gene_name(k)
  420. h = return_gene_name(l[i+1])
  422. for p,j in enumerate(gene_set):
  424. if type(j) == type(gene_set[0]):
  425. filtered_gene_set.append(j)
  426. else:
  427. continue
  429. print(len(filtered_gene_set))
  430. enr = gp.enrichr(gene_list= filtered_gene_set, gene_sets=['GO_Molecular_Function_2023','MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Cellular_Component_2023'], organism='human')
  431. df2 = enr.results.sort_values(by = 'Adjusted P-value')[enr.results.sort_values(by = 'Adjusted P-value')['Adjusted P-value']<=0.05]
  432. df2 = df2[df2['Adjusted P-value']<0.05]
  434. data_dict[f'{k}-{h}'] = list(df2.Term)
  435. print(f'{k}-{h}')
  436. df = pd.concat([df,df2])
  438. # %%
  439. df
  441. # %%
  442. data_dict.keys()
  444. # %%
  445. import pandas as pd
  446. from tabulate import tabulate
  447. from openpyxl import Workbook, load_workbook
  448. from openpyxl.styles import Border, Side, Font, PatternFill
  450. # Input dictionary
  451. data = data_dict
  452. # Extract unique pathways
  453. pathways = sorted(set(p for v in data.values() for p in v))
  455. columns = list(data.keys())
  456. df = pd.DataFrame(index=pathways, columns=columns)
  457. column_widths = {
  458. "A": 80, # Increase column A width
  459. "B": 25, # Increase column B width
  460. "C": 25, # Increase column C width
  461. "D": 25,
  462. "E":25,
  464. "F": 25, # Increase column A width
  465. "G": 25, # Increase column B width
  466. "H": 25, # Increase column C width
  467. "I": 25,
  468. "K":25,
  470. "L": 25, # Increase column A width
  471. "M": 25, # Increase column B width
  472. "N": 25, # Increase column C width
  473. "O": 25,
  474. "P":25,
  476. "Q": 25, # Increase column A width
  477. "R": 25, # Increase column B width
  478. "S": 25, # Increase column C width
  479. "T": 25,
  480. "U":25,
  482. "V": 25, # Increase column A width
  483. "W": 25 #Increase column D width
  484. }
  486. # Function to insert tick mark
  487. def tick():
  488. return "✔"
  490. # Fill DataFrame with tick marks
  491. for key, values in data.items():
  492. for value in values:
  493. df.loc[value, key] = tick()
  495. df.fillna("", inplace=True)
  497. # Define styles
  498. thin_border = Border(left=Side(style='thin'),
  499. right=Side(style='thin'),
  500. top=Side(style='thin'),
  501. bottom=Side(style='thin'))
  502. bold_font = Font(bold=True)
  503. gray_fill = PatternFill(start_color="D3D3D3", end_color="D3D3D3", fill_type="solid")
  505. # Set chunk sizes
  506. columns_per_sheet = 25
  507. rows_per_sheet = 10000
  509. # Split into chunks of 8 columns and 28 rows
  510. total_col_chunks = (len(df.columns) // columns_per_sheet) + (1 if len(df.columns) % columns_per_sheet else 0)
  511. total_row_chunks = (len(df) // rows_per_sheet) + (1 if len(df) % rows_per_sheet else 0)
  513. file_count = 1
  515. for col_chunk in range(total_col_chunks):
  516. for row_chunk in range(total_row_chunks):
  517. start_col = col_chunk * columns_per_sheet
  518. end_col = start_col + columns_per_sheet
  519. start_row = row_chunk * rows_per_sheet
  520. end_row = start_row + rows_per_sheet
  522. df_chunk = df.iloc[start_row:end_row, start_col:end_col]
  524. # Ensure the chunk has the correct dimensions
  525. while len(df_chunk.columns) < columns_per_sheet:
  526. df_chunk[f"Extra_{len(df_chunk.columns)+1}"] = ""
  528. excel_path = f"/12tb_dsk2/danish/GRN_new_21_feb_2025/data_preprocessed/grn_analysis_output/gsea_excel_sheets/pathway_table_part_{file_count}.xlsx"
  529. df_chunk.to_excel(excel_path)
  531. # Load the workbook and apply formatting
  532. wb = load_workbook(excel_path)
  533. ws = wb.active
  535. # Freeze headers
  536. ws.freeze_panes = "B2" # Freeze first row and column
  538. # Apply styles
  539. for row in ws.iter_rows():
  540. for cell in row:
  541. cell.border = thin_border
  542. if cell.value == "✔":
  543. cell.font = bold_font
  544. cell.fill = gray_fill
  546. # Increase row height
  547. for row in ws.iter_rows():
  548. ws.row_dimensions[row[0].row].height = 20
  550. for col_letter, width in column_widths.items():
  551. ws.column_dimensions[col_letter].width = width
  553. wb.save(excel_path)
  555. file_count += 1
  557. print("Excel files created successfully! ✅")
  560. # %%
  561. df
  563. # %%
  564. # Convert DataFrame to HTML
  565. html_path = "temp.html"
  566. df.to_html(html_path, index=False)
  568. # Convert HTML to PDF
  569. pdf_path = "output.pdf"
  570. pdfkit.from_file(html_path, pdf_path)
  572. print(f"Excel file converted to PDF at {pdf_path}")
  574. # %%
  575. # Convert DataFrame to HTML
  576. html_path = "temp.html"
  577. df.to_html(html_path, index=False)
  579. # Convert HTML to PDF
  580. pdf_path = "output.pdf"
  581. pdfkit.from_file(html_path, pdf_path)
  583. print(f"Excel file converted to PDF at {pdf_path}")
  585. # %%
  586. hub_tf= ['ENSG00000186300','ENSG00000256294','ENSG00000234284','ENSG00000188785']
  588. # %%
  589. hub_tf = set(l).intersection(set(hub_tf))
  591. # %%
  592. hub_tf = list(hub_tf)
  593. hub_tf
  595. # %%
  596. relevant_genes
  598. # %%
  599. hub_tf = ['ENSG00000188785', 'ENSG00000234284', 'ENSG00000256294']
  601. # %%
  602. hub_attached_genes = ['ENSG00000122641','ENSG00000149089','ENSG00000272576'] # manually checked
  604. # %%
  605. m = []
  606. for k,i in enumerate(hub_tf):
  607. print(return_gene_name(hub_tf[k]))
  608. print(return_gene_name(hub_attached_genes[k]))
  610. m.append(f'{return_gene_name(hub_tf[k])}-{return_gene_name(hub_attached_genes[k])}')
  612. # %%
  613. m
  615. # %%
  616. m.append('Label')
  618. # %%
  619. n
  621. # %%
  622. n.append('Label')
  623. n
  625. # %%
  626. important_gene_dataframe.columns = n
  628. # %%
  629. important_gene_dataframe
  631. # %%
  633. df = important_gene_dataframe[m].copy()
  635. # Splitting data by labels
  636. df
  638. # %%
  639. important_gene_dataframe.columns
  641. # %%
  642. ad_data = df[df['Label'] == 'AD'].drop(columns='Label')
  643. nci_data = df[df['Label'] == 'NCI'].drop(columns='Label')
  644. from scipy import stats
  646. # Initialize dictionary to store p-values for each column
  647. p_values = {}
  649. # Perform median test for each column
  650. for column in ad_data.columns:
  651. # Get the data for the current column in each group
  652. group1 = ad_data[column]
  653. group2 = nci_data[column]
  655. # Perform the median test
  656. stat, p_value, med, table = stats.median_test(group1, group2)
  658. # Store the p-value for the column
  659. p_values[column] = p_value
  661. # Output the p-values
  662. print(p_values)
  664. # %%
  665. ad_data[['ZNF548-INHBA','ZNF879-APIP']]
  667. # %%
  670. # %%
  671. import pandas as pd
  672. import seaborn as sns
  673. import matplotlib.pyplot as plt
  674. from sklearn.preprocessing import MinMaxScaler
  676. # Set the Seaborn theme for better aesthetics
  677. sns.set_theme(style="whitegrid", palette="muted")
  679. # Assuming 'important_gene_dataframe' is your DataFrame and it has a 'Label' column
  681. # Normalize the entire DataFrame (excluding the 'Label' column) to a 0-1 range
  682. scaler = MinMaxScaler()
  683. important_gene_dataframe_normalized = important_gene_dataframe[m].copy()
  685. # Apply normalization (only to the gene expression columns, not the 'Label' column)
  686. gene_columns = important_gene_dataframe_normalized.drop(columns=['Label']).columns
  687. important_gene_dataframe_normalized[gene_columns] = scaler.fit_transform(important_gene_dataframe_normalized[gene_columns])
  689. # Filter DataFrames for 'AD' and 'NCI'
  690. ad_data = important_gene_dataframe_normalized[important_gene_dataframe_normalized['Label'] == 'AD']
  691. nci_data = important_gene_dataframe_normalized[important_gene_dataframe_normalized['Label'] == 'NCI']
  693. # Drop the 'Label' column and only keep the gene expression values
  694. ad_data = ad_data.drop(columns=['Label'])
  695. nci_data = nci_data.drop(columns=['Label'])
  697. # List of genes (columns) to iterate over
  698. genes = ad_data.columns
  700. # Number of genes
  701. num_genes = len(genes)
  702. ad_data = remove_outliers(ad_data)
  703. nci_data = remove_outliers(nci_data)
  705. # Determine the number of rows and columns for subplots (depending on the number of genes)
  706. ncols = 12 # Set the number of columns (adjust based on your preference)
  707. nrows = (num_genes + ncols - 1) // ncols # Automatically adjust the number of rows
  709. # Create a figure and axes for subplots
  710. fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(50, 10 * nrows))
  711. axes = axes.flatten() # Flatten axes array for easy indexing
  713. # Loop through each gene and create a boxplot on a subplot
  714. for idx, gene in enumerate(genes):
  715. ax = axes[idx]
  717. # Prepare the data for the current gene (AD vs NCI)
  718. data = pd.DataFrame({
  719. 'Group': ['AD'] * len(ad_data[gene]) + ['NCI'] * len(nci_data[gene]),
  720. 'Expression': list(ad_data[gene]) + list(nci_data[gene])
  721. })
  723. # Create a boxplot for the current gene
  724. sns.boxplot(data=data, x='Group', y='Expression', palette='pastel', ax=ax)
  725. gene = m[idx]
  726. # Add a title and custom labels
  727. ax.set_title(f'{gene}', fontsize=14, weight='bold')
  728. ax.set_xlabel('Group', fontsize=12)
  729. ax.set_ylabel('Normalized Expression (0-1)', fontsize=12)
  731. # Add gridlines for better readability
  732. max_y = max(data['Expression']) + 0.05
  733. line_y = max_y + 0.02 # Line position
  734. cap_height = 0.01 # Height of vertical caps
  736. # Add horizontal line for significance comparison
  737. ax.hlines(y=line_y, xmin=0, xmax=1, color='black', linewidth=1.5)
  739. # Add vertical caps at both ends of the horizontal line
  740. ax.vlines(x=0, ymin=line_y - cap_height, ymax=line_y, color='black', linewidth=1.5)
  741. ax.vlines(x=1, ymin=line_y - cap_height, ymax=line_y, color='black', linewidth=1.5)
  743. # Display p-value above the plot
  744. ax.text(0.5, line_y + 0.02, f"p = {p_value:.3f}", ha='center', va='bottom', fontsize=12, color='black')
  746. # Gridlines for better readability
  747. ax.grid(True, axis='y', linestyle='--', alpha=0.6)
  749. # Remove any empty subplots if the number of genes is not a perfect multiple of ncols
  750. for idx in range(num_genes, len(axes)):
  751. fig.delaxes(axes[idx])
  753. # Adjust layout for better spacing between subplots
  754. plt.tight_layout()
  756. # Save the figure as an SVG with high resolution
  757. fig.savefig('/12tb_dsk2/danish/GRN_new_21_feb_2025/data_preprocessed/grn_analysis_output/plots/gene_expression_boxplots_subset_after_removing_outliers.svg', format='svg', dpi=300)
  759. # Show the plot
  760. plt.show()
  763. # %%
  764. import pandas as pd
  765. import seaborn as sns
  766. import matplotlib.pyplot as plt
  767. from sklearn.preprocessing import MinMaxScaler
  769. # Set the Seaborn theme for better aesthetics
  770. sns.set_theme(style="whitegrid", palette="muted")
  772. # Assuming 'important_gene_dataframe' is your DataFrame and it has a 'Label' column
  773. scaler = MinMaxScaler()
  775. # Copy the DataFrame and exclude the 'Label' column from normalization
  776. important_gene_dataframe_normalized = important_gene_dataframe.copy()
  777. gene_columns = important_gene_dataframe.drop(columns=['Label']).columns
  778. important_gene_dataframe_normalized[gene_columns] = scaler.fit_transform(important_gene_dataframe[gene_columns])
  780. # Filter DataFrames for 'AD' and 'NCI'
  781. ad_data = important_gene_dataframe_normalized[important_gene_dataframe_normalized['Label'] == 'AD'].copy()
  782. nci_data = important_gene_dataframe_normalized[important_gene_dataframe_normalized['Label'] == 'NCI'].copy()
  784. # Remove the 'Label' column to work only with gene expression values
  785. ad_data = ad_data.drop(columns=['Label'])
  786. nci_data = nci_data.drop(columns=['Label'])
  788. # Function to remove outliers using IQR
  789. def remove_outliers(df):
  790. Q1 = df.quantile(0.25) # First quartile
  791. Q3 = df.quantile(0.75) # Third quartile
  792. IQR = Q3 - Q1 # Interquartile range
  793. lower_bound = Q1 - 1.5 * IQR
  794. upper_bound = Q3 + 1.5 * IQR
  795. return df[~((df < lower_bound) | (df > upper_bound)).any(axis=1)]
  797. # Apply outlier removal
  798. ad_data_clean = remove_outliers(ad_data)
  799. nci_data_clean = remove_outliers(nci_data)
  801. # List of genes to iterate over
  802. genes = ad_data_clean.columns
  803. num_genes = len(genes)
  805. # Determine subplot layout
  806. ncols = 12
  807. nrows = (num_genes + ncols - 1) // ncols # Auto-adjust rows
  809. # Create figure
  810. fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(50, 10 * nrows))
  811. axes = axes.flatten()
  813. # Loop through each gene to create boxplots
  814. for idx, gene in enumerate(genes):
  815. ax = axes[idx]
  817. # Prepare data for plotting
  818. data = pd.DataFrame({
  819. 'Group': ['AD'] * len(ad_data_clean[gene]) + ['NCI'] * len(nci_data_clean[gene]),
  820. 'Expression': list(ad_data_clean[gene]) + list(nci_data_clean[gene])
  821. })
  823. # Create boxplot
  824. sns.boxplot(data=data, x='Group', y='Expression', palette='pastel', ax=ax)
  826. # Set title and labels
  827. ax.set_title(f'{gene}', fontsize=14, weight='bold')
  828. ax.set_xlabel('Group', fontsize=12)
  829. ax.set_ylabel('Normalized Expression (0-1)', fontsize=12)
  831. # Grid for better readability
  832. ax.grid(True, axis='y', linestyle='--', alpha=0.6)
  834. # Remove empty subplots
  835. for idx in range(num_genes, len(axes)):
  836. fig.delaxes(axes[idx])
  838. # Adjust layout
  839. plt.tight_layout()
  841. # Show plot
  842. plt.show()
  845. # %%
  846. m
  848. # %%
  849. df = pd.DataFrame()
  850. data_dict = dict()
  851. for i,k in tqdm(enumerate(l)):
  852. if i%2 != 0:
  853. continue
  855. m = []
  856. try:
  857. #print(k)
  858. #print(l[i+1])
  859. m.append(k)
  860. m.extend(list(AD_BG_bipartite[k]))
  861. m.extend(list(AD_BG_bipartite[l[i+1]]))
  862. except:
  863. print(k)
  866. m = set(m)
  867. print(len(m))
  869. gene_set = [return_gene_name(gene_id) for gene_id in m]
  870. gene_set = list(set(gene_set))
  872. try:
  873. gene_set.remove('Gene name not found')
  874. except:
  875. print('skip')
  877. filtered_gene_set = []
  878. k = return_gene_name(k)
  879. h = return_gene_name(l[i+1])
  881. for p,j in enumerate(gene_set):
  883. if type(j) == type(gene_set[0]):
  884. filtered_gene_set.append(j)
  885. else:
  886. continue
  888. print(len(filtered_gene_set))
  889. enr = gp.enrichr(gene_list= filtered_gene_set, gene_sets=['GO_Molecular_Function_2023','MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Cellular_Component_2023'], organism='human')
  890. df2 = enr.results.sort_values(by = 'Adjusted P-value')[enr.results.sort_values(by = 'Adjusted P-value')['Adjusted P-value']<=0.05]
  891. df2 = df2[df2['Adjusted P-value']<0.005]
  893. data_dict[f'{k}-{h}'] = list(df2.Term)
  894. print(f'{k}-{h}')
  895. df = pd.concat([df,df2])
  897. # %%
  898. # 'GO_Biological_Process_2023'
  900. # %%
  901. df.Gene_set.value_counts()
  903. # %%
  904. df
  906. # %%
  907. import plotly.graph_objects as go
  909. # Data preparation as before
  910. sources = []
  911. targets = []
  912. for gene, pathways in data_dict.items():
  913. for pathway in pathways:
  914. sources.append(gene)
  915. targets.append(pathway)
  917. # Create a mapping of genes and pathways to unique indices
  918. unique_items = list(set(sources + targets))
  919. item_indices = {item: idx for idx, item in enumerate(unique_items)}
  921. # Create lists of indices for the Sankey diagram
  922. source_indices = [item_indices[gene] for gene in sources]
  923. target_indices = [item_indices[pathway] for pathway in targets]
  925. # Define the Sankey diagram figure
  926. fig = go.Figure(go.Sankey(
  927. node=dict(
  928. pad=15,
  929. thickness=20,
  930. line=dict(color="black", width=0.5),
  931. label=unique_items,
  932. align='right'
  933. ),
  934. link=dict(
  935. source=source_indices,
  936. target=target_indices,
  937. value=[1] * len(sources), # All links have equal value
  938. )
  939. ))
  941. # Update layout with increased figure size
  942. fig.update_layout(
  943. title_text="Gene to Pathway Sankey Diagram",
  944. font_size=12,
  945. width=900, # Increase width
  946. height=1550, # Increase height
  947. )
  949. # Save the figure as SVG
  950. fig.write_image("/12tb_dsk2/danish/GRN_new_21_feb_2025/data_preprocessed/grn_analysis_output/plots/sankey_enriched_plots.svg", format="svg", scale=1)# Show the figure
  951. fig.show()
  954. # %%
  955. set(df.Term)
  957. # %%
  958. df = pd.DataFrame()
  959. data_dict = dict()
  960. for i,k in tqdm(enumerate(l)):
  961. if i%2 != 0:
  962. continue
  964. m = []
  965. try:
  966. #print(k)
  967. #print(l[i+1])
  968. m.append(k)
  969. m.extend(list(AD_BG_bipartite[k]))
  970. m.extend(list(AD_BG_bipartite[l[i+1]]))
  971. except:
  972. print(k)
  975. m = set(m)
  976. print(len(m))
  978. gene_set = [return_gene_name(gene_id) for gene_id in m]
  979. gene_set = list(set(gene_set))
  981. try:
  982. gene_set.remove('Gene name not found')
  983. except:
  984. print('skip')
  986. filtered_gene_set = []
  987. k = return_gene_name(k)
  988. h = return_gene_name(l[i+1])
  990. for p,j in enumerate(gene_set):
  992. if type(j) == type(gene_set[0]):
  993. filtered_gene_set.append(j)
  994. else:
  995. continue
  997. print(len(filtered_gene_set))
  998. enr = gp.enrichr(gene_list= filtered_gene_set, gene_sets=['GO_Molecular_Function_2023','MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Cellular_Component_2023'], organism='human')
  999. df2 = enr.results.sort_values(by = 'Adjusted P-value')[enr.results.sort_values(by = 'Adjusted P-value')['Adjusted P-value']<=0.05]
  1000. df2 = df2[df2['Adjusted P-value']<0.05]
  1002. data_dict[f'{k}-{h}'] = list(df2.Term)
  1003. print(f'{k}-{h}')
  1004. df = pd.concat([df,df2])
  1006. # %%
  1007. data_dict.keys()
  1009. # %%
  1010. data_dict['MTF1-CYLD']
  1012. # %%
  1013. relevant_genes
  1015. # %%
  1016. AD_related_genes = pd.read_csv('../../../data_preprocessed/previously_researched_data/gene-list.csv')
  1017. AD_related_genes
  1019. # %%
  1020. gene_name_to_ensemble = {v: k for k, v in ensemble_to_gene_name.items()}
  1023. # %%
  1024. def return_ensemble_name(gene_name):
  1026. try:
  1028. return gene_name_to_ensemble[gene_name]
  1029. except:
  1030. return gene_name
  1032. # %%
  1033. AD_related_genes_list = [return_ensemble_name(gene_name) for gene_name in AD_related_genes['Gene Symbol'].tolist()]
  1035. # %%
  1036. AD_related_genes_list[0:10]
  1038. # %%
  1039. s = set()
  1041. for n in range(0,len(l),2):
  1042. k = l[n]
  1043. connected_nodes = list(AD_BG_bipartite.neighbors(k))
  1044. genes_l = set(connected_nodes).intersection(set(AD_related_genes_list))
  1045. x = len(set(connected_nodes).intersection(set(AD_related_genes_list)))
  1046. if x>0:
  1047. print(k)
  1048. print(len(connected_nodes))
  1049. print(len(genes_l))
  1050. print(len(genes_l)/len(connected_nodes))
  1051. print(genes_l)
  1052. s.update(genes_l)
  1054. # %%
  1055. len(s)
  1057. # %%
  1058. len(s)/948
  1060. # %%
  1061. import networkx as nx
  1063. # Assuming l is your list of nodes and AD_BG_bipartite is your bipartite graph
  1065. # Initialize a set to store the nodes to be included in the subgraph
  1066. nodes_to_include = set(l)
  1068. for n in range(0,len(l),2):
  1069. k = l[n]
  1070. connected_nodes = list(AD_BG_bipartite.neighbors(k))
  1071. nodes_to_include.update(connected_nodes)
  1073. # Now create the subgraph containing only the nodes in nodes_to_include
  1074. subgraph = AD_BG_bipartite.subgraph(nodes_to_include).copy()
  1076. # subgraph now contains the nodes in 'l' and their neighbors
  1078. # %%
  1079. highlighted_edges = []
  1080. for i in relevant_genes:
  1081. tf = i.split('_')[0]
  1082. gene = i.split('_')[1]
  1083. highlighted_edges.append((tf, gene))
  1085. # %%
  1086. import networkx as nx
  1087. import matplotlib.pyplot as plt
  1088. from adjustText import adjust_text
  1090. # Assuming l is your list of nodes and AD_BG_bipartite is your bipartite graph
  1092. # Initialize a set to store the nodes to be included in the subgraph
  1093. nodes_to_include = set()
  1094. s = []
  1095. # Add neighbors to nodes_to_include
  1096. for n in range(0, len(l), 2):
  1097. k = l[n]
  1098. nodes_to_include.update(set([k]))
  1099. s.append(l[n])
  1100. connected_nodes = list(AD_BG_bipartite.neighbors(k))
  1101. genes_l = set(connected_nodes).intersection(set(AD_related_genes_list))
  1102. genes_l.update(l)
  1103. nodes_to_include.update(genes_l)
  1106. # Now create the subgraph containing only the nodes in nodes_to_include
  1107. subgraph = AD_BG_bipartite.subgraph(nodes_to_include).copy()
  1109. # Position nodes using a layout (for example, spring layout)
  1110. pos = nx.spring_layout(subgraph)
  1112. # Draw nodes and edges
  1113. plt.figure(figsize=(20, 20))
  1114. node_sizes = []
  1115. # Draw the nodes
  1116. # Nodes from list 'l' will be squares, others will be circles
  1117. node_shapes = []
  1118. for node in subgraph.nodes():
  1120. if node in s:
  1121. node_shapes.append('s')
  1122. node_sizes.append(380) # Square for nodes in 'l'
  1123. else:
  1124. node_shapes.append('o')
  1125. node_sizes.append(250) # Circle for connected nodes
  1127. # Draw the graph with node shapes based on 'l' membership
  1128. nx.draw_networkx_edges(subgraph, pos, width=1.0, alpha=0.7, edge_color='blue')
  1130. # Draw the nodes with different shapes (squares for 'l' and circles for others)
  1131. for i, node in enumerate(subgraph.nodes()):
  1132. shape = node_shapes[i]
  1133. size = node_sizes[i]
  1134. nx.draw_networkx_nodes(subgraph, pos, nodelist=[node], node_size=size, node_color='blue', node_shape=shape)
  1135. nx.draw_networkx_nodes(subgraph, pos, nodelist=s, node_size=size, node_color='red',node_shape=shape)
  1137. gene_imp = set(l) - set(s)
  1138. gene_imp = list(gene_imp)
  1139. nx.draw_networkx_nodes(subgraph, pos, nodelist=gene_imp, node_size=450, node_color='red',node_shape='o')
  1140. nx.draw_networkx_edges(subgraph, pos, edgelist=highlighted_edges, width=3, alpha=0.8, edge_color='red')
  1142. # Labeling (only labels for nodes in 'l')
  1143. labels = {node: return_gene_name(node) for node in l} # Label nodes from 'l' only
  1144. texts = []
  1146. for node, label in labels.items():
  1147. text = plt.text(pos[node][0], pos[node][1], label, fontsize=6, fontweight='bold', color='black', ha='left')
  1148. texts.append(text)
  1150. # Adjust label positions to avoid overlap
  1151. adjust_text(texts, arrowprops=dict(arrowstyle="->", color='black', lw=1))
  1153. # Finalize plot
  1154. plt.title("Subgraph with Adjusted Labels for 'l' Nodes", fontsize=16)
  1155. plt.axis('off') # Remove axis
  1156. plt.savefig('../../../data_preprocessed/grn_analysis_output/plots/network.svg', format='svg', dpi=600) # Corrected .savefig -> plt.savefig
  1157. plt.show()
  1161. # %%
  1162. labels
  1164. # %% [markdown]
  1165. # # AD-targeted Genes-collected: https://agora.adknowledgeportal.org/genes/nominated-targets
  1167. # %%
  1170. # %%
  1173. # %%
  1176. # %%
  1179. # %%
  1182. # %%
  1185. # %%
  1188. # %%
  1191. # %%
  1194. # %%
  1197. # %%
  1198. TNF-alpha signaling via NF-kB, structural components like collagen-containing extracellular matrix (GO:0062023) and microvillus (GO:0005902), as well as molecular activities such as hormone activity (GO:0005179)
  1200. # %%
  1201. 'collagen-containing extracellular matrix (GO:0062023)' in df.Term.tolist()
  1203. # %%
  1204. df.Term.tolist()
  1206. # %%
  1207. df = pd.DataFrame()
  1208. data_dict = dict()
  1209. for i,k in tqdm(enumerate(l)):
  1210. if i%2 != 0:
  1211. continue
  1213. m = []
  1214. try:
  1215. #print(k)
  1216. #print(l[i+1])
  1217. m.append(k)
  1218. m.extend(list(AD_BG_bipartite[k]))
  1219. m.extend(list(AD_BG_bipartite[l[i+1]]))
  1220. except:
  1221. print(k)
  1224. m = set(m)
  1225. print(len(m))
  1227. gene_set = [return_gene_name(gene_id) for gene_id in m]
  1228. gene_set = list(set(gene_set))
  1230. try:
  1231. gene_set.remove('Gene name not found')
  1232. except:
  1233. print('skip')
  1235. filtered_gene_set = []
  1236. k = return_gene_name(k)
  1237. h = return_gene_name(l[i+1])
  1239. for p,j in enumerate(gene_set):
  1241. if type(j) == type(gene_set[0]):
  1242. filtered_gene_set.append(j)
  1243. else:
  1244. continue
  1246. print(len(filtered_gene_set))
  1247. enr = gp.enrichr(gene_list= filtered_gene_set, gene_sets=['GO_Molecular_Function_2023','MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Cellular_Component_2023'], organism='human')
  1248. df2 = enr.results.sort_values(by = 'Adjusted P-value')[enr.results.sort_values(by = 'Adjusted P-value')['Adjusted P-value']<=0.05]
  1249. df2 = df2[df2['Adjusted P-value']<0.05]
  1251. data_dict[f'{k}-{h}'] = list(df2.Term)
  1252. print(f'{k}-{h}')
  1253. df = pd.concat([df,df2])
  1255. # %%
  1256. df
  1258. # %%
  1259. 'collagen-containing extracellular matrix (GO:0062023)' in df.Term.tolist()
  1261. # %%
  1262. diff_exp_enriched = set(['TNF-alpha Signaling via NF-kB', 'Collagen-Containing Extracellular Matrix (GO:0062023)','Microvillus (GO:0005902)', 'Hormone Activity (GO:0005179)', 'Neutral L-amino Acid:Sodium Symporter Activity (GO:0005295)'])
  1264. # %%
  1265. diff_exp_enriched - set(df.Term.tolist())
  1267. # %%
  1268. diff_exp_enriched.
  1270. # %%
  1271. df = pd.DataFrame()
  1272. data_dict = dict()
  1273. for i,k in tqdm(enumerate(l)):
  1274. if i%2 != 0:
  1275. continue
  1277. m = []
  1278. try:
  1279. #print(k)
  1280. #print(l[i+1])
  1281. m.append(k)
  1282. m.extend(list(AD_BG_bipartite[k]))
  1283. m.extend(list(AD_BG_bipartite[l[i+1]]))
  1284. except:
  1285. print(k)
  1288. m = set(m)
  1289. print(len(m))
  1291. gene_set = [return_gene_name(gene_id) for gene_id in m]
  1292. gene_set = list(set(gene_set))
  1294. try:
  1295. gene_set.remove('Gene name not found')
  1296. except:
  1297. print('skip')
  1299. filtered_gene_set = []
  1300. k = return_gene_name(k)
  1301. h = return_gene_name(l[i+1])
  1303. for p,j in enumerate(gene_set):
  1305. if type(j) == type(gene_set[0]):
  1306. filtered_gene_set.append(j)
  1307. else:
  1308. continue
  1310. print(len(filtered_gene_set))
  1311. enr = gp.enrichr(gene_list= filtered_gene_set, gene_sets=['GO_Molecular_Function_2023','MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Cellular_Component_2023'], organism='human')
  1312. df2 = enr.results.sort_values(by = 'Adjusted P-value')[enr.results.sort_values(by = 'Adjusted P-value')['Adjusted P-value']<=0.05]
  1313. df2 = df2[df2['Adjusted P-value']<0.005]
  1315. data_dict[f'{k}-{h}'] = list(df2.Term)
  1316. print(f'{k}-{h}')
  1317. df = pd.concat([df,df2])
  1319. # %%
  1320. diff_exp_enriched - set(df.Term.tolist())
  1322. # %%
  1325. # %%
  1328. # %%

plots(box_plot_sankey_plots_network_plots).ipynb at commit 3f3f8d9, no license · at the source

Overview

Authors: Danish Anwer1, Arina A1,2, Agata Marchi1, Eduard Kerkhoven1, Annikka Polster1,2
  1. Division of Systems and Synthetic Biology, Department of Life Sciences, Chalmers University of Technology,Gothenburg, Sweden
  2. Department of Microbiology, Oslo University Hospital and University of Oslo,Oslo, Norway
Institutions: Chalmers University of Technology (Sweden); Oslo University Hospital (Norway)
Journal: npj aging, volume 12, issue 1, article 95
Dates: received 16 June 2025; accepted 7 July 2026; published online 16 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41514-026-00443-0 · PMID 42463670 · PMCID PMC13376175 · OpenAlex W7168694839
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: Alzheimer's / dementia (population)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Connectivity, Graphs
Keywords: Neuroscience, Diseases
Topic: Bioinformatics and Genomic Networks (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Area of Advance Health Engineering at Chalmers University of Technology
Citations: not cited yet (Europe PMC); 48 references in the paper

Abstract

Alzheimer’s disease (AD) is a multifactorial neurodegenerative disorder marked by progressive cognitive decline, yet its transcriptional regulatory architecture remains poorly understood. Here, we model sample-specific gene regulatory networks (GRNs) from dorsolateral prefrontal cortex transcriptomes of 87 individuals with AD and 67 non-cognitively impaired (NCI) controls and use a machine learning classifier to detect consistent disease-specific network features. This sample-specific network approach captures inter-individual variation in transcriptional regulation and revealed 22 key transcription factor-gene regulations that distinguish AD from NCI with 96% weighted accuracy. The key transcription factor-gene interactions were enriched in pathways central to AD pathology, including synaptic signalling, mitochondrial function, proteostasis, and neuroinflammation. Network analysis uncovered significant differences in regulatory connectivity between AD and controls, with ZNF225, ZNF849, and ZNF548 emerging as AD-specific regulatory hubs. Moreover, several key regulatory edges showed significant correlations with longitudinal cognitive decline, supporting their clinical relevance. Our findings highlight pervasive transcriptional dysregulation in AD, emphasizing sample-specific GRN modelling’s value in uncovering regulatory mechanisms.

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

Repository

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

Polster-lab/Analysis_GRN_ROSMAP

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 3f3f8d9e1bb54c67b99042add4a411fbae4b0127, 3 April 2025
Languages: Jupyter (55)
Size: 62 files, 55 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: 27 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: pandas (32 files), NumPy (24 files), Matplotlib (12 files), seaborn (12 files), SciPy (10 files), NetworkX (7 files), Plotly (6 files), scikit-learn (6 files), DESeq2 (1 file), ggplot2 (1 file), ggpubr (1 file), reshape2 (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
39 files
At the source: github.com/Polster-lab/Analysis_GRN_ROSMAP

Code availability

GitHub: https://github.com/Polster-lab/Analysis_GRN_ROSMAP,.

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

Tracing map

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

What the map holds:

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

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

Data

Datasets cited

Data availability

RNA-Seq data: https://www.synapse.org/Synapse:syn3219045 The datasets generated during this study are available in the Zenodo repository: 10.5281/zenodo.19543885.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 2 keywords, 1 funder, 44 references.

Cite

This paper

Anwer, D., A, A., Marchi, A., Kerkhoven, E., & Polster, A. (2026). Network-based discovery of regulatory drivers of cognitive decline in alzheimer's disease. npj aging, 12(1), 95. https://doi.org/10.1038/s41514-026-00443-0

BibTeX

@article{anwer2026network,
author = {Anwer, Danish and A, Arina and Marchi, Agata and Kerkhoven, Eduard and Polster, Annikka},
title = {{Network-based discovery of regulatory drivers of cognitive decline in alzheimer's disease}},
journal = {npj aging},
year = {2026},
month = jul,
volume = {12},
number = {1},
pages = {95},
publisher = {Nature Publishing Group},
issn = {2731-6068},
doi = {10.1038/s41514-026-00443-0},
url = {https://doi.org/10.1038/s41514-026-00443-0},
pmid = {42463670},
pmcid = {PMC13376175}
}

RIS

TY - JOUR
AU - Anwer, Danish
AU - A, Arina
AU - Marchi, Agata
AU - Kerkhoven, Eduard
AU - Polster, Annikka
TI - Network-based discovery of regulatory drivers of cognitive decline in alzheimer's disease
T2 - npj aging
J2 - NPJ Aging
PY - 2026
DA - 2026/07/16
VL - 12
IS - 1
SP - 95
SN - 2731-6068
PB - Nature Publishing Group
DO - 10.1038/s41514-026-00443-0
UR - https://doi.org/10.1038/s41514-026-00443-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41514-026-00443-0",
"type": "article-journal",
"title": "Network-based discovery of regulatory drivers of cognitive decline in alzheimer's disease",
"container-title": "npj aging",
"author": [
{
"family": "Anwer",
"given": "Danish"
},
{
"family": "A",
"given": "Arina"
},
{
"family": "Marchi",
"given": "Agata"
},
{
"family": "Kerkhoven",
"given": "Eduard"
},
{
"family": "Polster",
"given": "Annikka"
}
],
"container-title-short": "NPJ Aging",
"volume": "12",
"issue": "1",
"page": "95",
"DOI": "10.1038/s41514-026-00443-0",
"PMID": "42463670",
"PMCID": "PMC13376175",
"ISSN": "2731-6068",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41514-026-00443-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
16
]
]
}
}

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.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: DESeq2, NetworkX, Plotly, 10 other tools, 2 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: DESeq2, NetworkX, Plotly, 10 other tools, 1 reference
[3] doi:10.1038/s43587-026-01207-x [code]
A microprotein atlas of the human frontal cortex in Alzheimer's disease.
Journal: Nature aging
In common: DESeq2, Plotly, reshape2, 6 other tools, synapse.org/synapse:syn3219045, Alzheimer's / dementia
[4] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: DESeq2, NetworkX, Plotly, 10 other tools
[5] doi:10.1038/s44318-026-00818-9 [code]
FAM134B-mediated ER-phagy degrades APP and suppresses Alzheimer's disease pathology.
Journal: The EMBO journal
In common: DESeq2, Plotly, reshape2, 8 other tools, Alzheimer's / dementia, 1 reference
[6] doi:10.1016/j.xcrm.2026.102651 [code]
Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: DESeq2, NetworkX, Plotly, 9 other tools
[7] doi:10.1371/journal.pcbi.1014327 [code]
Supervised deep learning with gene functional annotation for cell classification.
Journal: PLoS computational biology
In common: DESeq2, reshape2, ggpubr, 8 other tools, Alzheimer's / dementia, 1 reference
[8] doi:10.1038/s41467-026-72598-z [code]
Functional impact of genetic background on variable expressivity in neurodevelopmental disorders.
Journal: Nature communications
In common: DESeq2, NetworkX, ggpubr, 7 other tools, 2 references
[9] 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: DESeq2, NetworkX, reshape2, 9 other tools
[10] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: DESeq2, NetworkX, reshape2, 9 other tools

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.