OSCR

MLMarker: a machine learning framework for tissue inference and biomarker discovery.

Code ↔ Paper

34 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 34 matches
  1. [1] § Reuse of pan-cancer dataset (MSV000095036) ↔ Cancer_datasets_Kuster/baseline_comparison_clean.ipynb, lines 79–138 · score 0.85 · parotid gland, small intestine, pituitary gland, colon, duodenum, rectum
  2. [2] § Baseline atlas comparisons and alternative classifiers › Interpretation analyses ↔ custom_functions.py, lines 248–305 · score 0.84 · cellular components, Human Protein Atlas, molecular functions, biological processes, KEGG, pathways
  3. [3] § Results › Reprocessed data from cerebral melanoma (PXD007592) ↔ pages/1_QC.py, lines 187–249 · score 0.84 · intensity variation, signal proportion, global sample, intensity distribution, protein intensity, feature space
  4. [4] § Reuse of pan-cancer dataset (MSV000095036) ↔ Home.py, lines 243–243 · score 0.77 · parotid gland, small intestine, pituitary gland, colon, duodenum, rectum
  5. [5] § Baseline atlas comparisons and alternative classifiers › Reference atlases ↔ Cancer_datasets_Kuster/baseline_comparison_clean.ipynb, lines 54–77 · score 0.76 · excluding cancer, normalised intensities, UniProt, tissue matrix, training atlas, HPA
  6. [6] § Results › Circularity and double-dipping ↔ PXD007592_responders_melanoma/03_shap_vs_dea_correlation.ipynb, lines 630–680 · score 0.74 · circular reasoning, reflects cohort, capture tissue, prone, complementary, leakage
  7. [7] § Results › Penalty factor › Performance across sample types: dense vs. sparse ↔ Missingnes_testing/02_penalty_factor_comparison.ipynb, lines 550–581 · score 0.72 · high dropout, low dropout, sparse biofluid, Penalty benefit, dense tissue, beneficial
  8. [8] § Penalty factor evaluation by simulated protein missingness › Penalty strategies ↔ Missingnes_testing/penalty_factor_analysis.py, lines 128–144 · score 0.71 · piecewise linear interpolation, Adaptive penalties, low coverage
  9. [9] § Reprocessed data from cerebrospinal fluid (PXD008029) ↔ Home.py, lines 243–243 · score 0.70 · bone marrow, adipose tissue, pituitary gland, testis, lung, monocytes
  10. [10] § Reprocessed data from cerebrospinal fluid (PXD008029) ↔ test.ipynb, lines 140–141 · score 0.70 · bone marrow, adipose tissue, pituitary gland, testis, lung, monocytes
  11. [11] § Results › Penalty factor › Performance across sample types: dense vs. sparse ↔ Missingnes_testing/02_penalty_factor_comparison.ipynb, lines 405–548 · score 0.69 · targeted missingness, sparse samples, dense tissue, loss, missing proteins, MNAR
  12. [12] § Materials and methods › The MLMarker package › Penalty factor for missing feature effects ↔ old_modules/mlmarker.py, lines 116–165 · score 0.64 · identifies absent proteins, remain unchanged, penalty factor, SHAP, zero, MLMarker
  13. [13] § Results › Reprocessed data from cerebral melanoma (PXD007592) ↔ PXD007592_responders_melanoma/classify_04_protein_analysis.ipynb, lines 41–127 · score 0.64 · Mann Whitney, good responders, higher brain, poor responders, brain prediction, scores
  14. [14] § Results › Reprocessed data from cerebral melanoma (PXD007592) ↔ PXD007592_responders_melanoma/05_clusteringgroups.ipynb, lines 224–297 · score 0.63 · response status, Mann Whitney, higher brain, Poor responders, UMAP, clustering
  15. [15] § Results › Penalty factor ↔ Missingnes_testing/02_penalty_factor_comparison.ipynb, lines 405–548 · score 0.62 · Targeted missingness, detection limits, loss, PXD009021, lost, MNAR
  16. [16] § Materials and methods › Reprocessing pipeline › Differential expression analysis ↔ PXD007592_responders_melanoma/msqrob_analysis.R, lines 106–138 · score 0.62 · msqrob2, Peptide intensities, aggregated, median, responders, models
  17. [17] § Materials and methods › The MLMarker package › Annotation utilities and visualisation ↔ custom_functions.py, lines 248–305 · score 0.62 · Human Protein Atlas, tissue expression, Profiler, enrichment, SHAP, score
  18. [18] § Results › Reprocessed data from cerebral melanoma (PXD007592) ↔ pages/1_QC.py, lines 20–28 · score 0.60 · predefined proteins, missing proteins, artifacts, quality, technical, intensity
  19. [19] § Baseline atlas comparisons and alternative classifiers › Alternative classification methods ↔ Cancer_datasets_Kuster/baseline_comparison_clean.ipynb, lines 368–415 · score 0.60 · nearest neighbours classification, sample profile, distance, correlation, training, tissue
  20. [20] § Reuse of pan-cancer dataset (MSV000095036) ↔ Cancer_datasets_Kuster/baseline_comparison_clean.ipynb, lines 79–138 · score 0.60 · healthy tissues, glioma, LN, OE, DLBCL, OSCC
  21. [21] § Results › Reprocessed data from cerebral melanoma (PXD007592) ↔ PXD007592_responders_melanoma/05_clusteringgroups.ipynb, lines 224–297 · score 0.59 · Mann Whitney, higher brain, poor responders, good responders, brain prediction, UMAP
  22. [22] § Baseline atlas comparisons and alternative classifiers › Interpretation analyses ↔ pages/2_Visualisations.py, lines 22–31 · score 0.58 · SHapley, exPlanations, Additive, predictions, tissue, proteins
  23. [23] § Results › Circularity and double-dipping ↔ PXD007592_responders_melanoma/03_shap_vs_dea_correlation.ipynb, lines 87–128 · score 0.57 · MLMarker feature space, proteins outside, log10 transformed, fold change, dipping, Circularity
  24. [24] § Reprocessed data from cerebrospinal fluid (PXD008029) ↔ PXD008029_CSF/classify_csf.ipynb, lines 211–256 · score 0.56 · adipose tissue, CSF samples, FDR, monocytes, plasma, PXD008029
  25. [25] § Materials and methods › The MLMarker package › Penalty factor for missing feature effects ↔ mlmarker/explainability.py, lines 146–208 · score 0.56 · optional penalty factor, Absent proteins, SHAP, zero, MLMarker, tissue
  26. [26] § Results › Circularity and double-dipping ↔ PXD007592_responders_melanoma/03_shap_vs_dea_correlation.ipynb, lines 630–680 · score 0.56 · indicates partial overlap, limited correlation, circular, SHAP, Spearman, MLMarker
  27. [27] § Results › Penalty factor › Performance across sample types: dense vs. sparse ↔ Missingnes_testing/rank_analysis_heatmap.py, lines 73–136 · score 0.56 · sparse biofluid, dense tissue, dropout rates, heatmap, adaptive, rank
  28. [28] § Penalty factor evaluation by simulated protein missingness ↔ Missingnes_testing/02_penalty_factor_comparison.ipynb, lines 14–46 · score 0.55 · pituitary gland, target tissue, PXD009021, Simulations, PXD008029, missingness
  29. [29] § Penalty factor evaluation by simulated protein missingness ↔ Missingnes_testing/penalty_factor_analysis.py, lines 193–205 · score 0.54 · pituitary gland, target tissue, NSAF, ionbot, PXD008029, missingness
  30. [30] § Results › Penalty factor › Performance across sample types: dense vs. sparse ↔ Missingnes_testing/rank_analysis_heatmap.py, lines 73–136 · score 0.54 · lower ranks, sparse biofluids, dense tissue, missingness, penalty
  31. [31] § Results › Circularity and double-dipping ↔ PXD007592_responders_melanoma/03_shap_vs_dea_correlation.ipynb, lines 87–128 · score 0.53 · MLMarker feature space, Double dipping, outside, circularity, external, protein
  32. [32] § Reprocessed data from cerebrospinal fluid (PXD008029) › SHAP analysis of CSF samples within the Streamlit application ↔ PXD008029_CSF/classify_csf.ipynb, lines 291–307 · score 0.52 · bone marrow, pituitary gland, CSF, brain, predicted
  33. [33] § Results › Differential expression analysis with MSqRob ↔ PXD007592_responders_melanoma/msqrob_analysis.R, lines 141–220 · score 0.51 · log fold change, log10 adjusted, volcano, MSqRob, responder, brain
  34. [34] § Results › Penalty factor › Performance across sample types: dense vs. sparse ↔ Missingnes_testing/penalty_factor_analysis.py, lines 322–359 · score 0.50 · adaptive piecewise, dropout rates, penalty factor, MNAR, simulating, missingness

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 · 2,059 lines · 86 KB · no license · 4 matches

  1. # %%
  2. # === 1. Setup and Imports ===
  3. import os
  4. import warnings
  5. warnings.filterwarnings('ignore')
  7. import pandas as pd
  8. import numpy as np
  9. import matplotlib.pyplot as plt
  10. import seaborn as sns
  11. from scipy.stats import spearmanr
  12. from tqdm import tqdm
  14. # Plot settings
  15. sns.set_style("whitegrid")
  16. plt.rcParams['figure.figsize'] = (12, 8)
  17. plt.rcParams['font.size'] = 12
  18. plt.rcParams['figure.dpi'] = 100
  20. print("Libraries loaded successfully.")
  22. # %% [markdown]
  23. # ## 2. Load Data
  24. #
  25. # Load protein expression data, MLMarker predictions, and reference atlases.
  27. # %%
  28. # Load protein expression data
  29. protein_data_path = 'search__Picked_Group_FDR_Outputs__Picked_Group_FDR_Outputs_-_FP_Booster__combined_protein.tsv'
  30. local_path = 'search__Picked_Group_FDR_Outputs__Picked_Group_FDR_Outputs_-_FP_Booster__combined_protein.tsv'
  32. protein_df = pd.read_csv(local_path if os.path.exists(local_path) else protein_data_path, sep='\t')
  34. # Load MLMarker predictions
  35. mlmarker_predictions = pd.read_csv('Entire_picked_FDR_MLMarker_predictions_dev_penalty1.csv', index_col=0)
  37. # Load metadata
  38. metadata = pd.read_csv('MassIVE_Mapping.txt', sep='\t')
  40. # Load reference atlases
  41. hpa = pd.read_csv('antibody_atlas.csv')
  42. pdb = pd.read_csv('proteomicsdb.csv')
  44. # Load MLMarker training data
  45. ml_train = pd.read_csv('training_atlas_7592%_10exp.csv.gz', index_col=0)
  47. print(f"Protein data: {protein_df.shape}")
  48. print(f"MLMarker predictions: {mlmarker_predictions.shape}")
  49. print(f"Metadata entries: {len(metadata)}")
  50. print(f"HPA: {hpa['Uniprot_id'].nunique()} proteins, {hpa['Tissue'].nunique()} tissues")
  51. print(f"PDB: {pdb['UniProt'].nunique()} proteins, {pdb['tissue'].nunique()} tissues")
  52. print(f"MLMarker training: {ml_train.shape[0]} samples, {ml_train['tissue_name'].nunique()} tissues")
  54. # %% [markdown]
  55. # ## 3. Create Atlas Matrices
  56. #
  57. # Convert all data sources to protein × tissue matrices.
  59. # %%
  60. # === HPA Atlas Matrix ===
  61. hpa_matrix = hpa.pivot_table(
  62. index='Uniprot_id', columns='Tissue', values='Level', aggfunc='mean'
  63. ).fillna(0)
  64. print(f"HPA matrix: {hpa_matrix.shape} (proteins × tissues)")
  66. # === PDB Atlas Matrix (exclude cancer/fluids) ===
  68. pdb_matrix = pdb.pivot_table(
  69. index='UniProt', columns='tissue', values='normalised intensity', aggfunc='mean'
  70. ).fillna(0)
  71. print(f"PDB matrix: {pdb_matrix.shape} (proteins × tissues)")
  73. # === MLMarker Training Atlas ===
  74. metadata_cols = ['tissue_name', 'disease_status', 'fluid', 'cell_type']
  75. protein_cols = [c for c in ml_train.columns if c not in metadata_cols]
  76. ml_atlas = ml_train.groupby('tissue_name')[protein_cols].mean().T
  77. print(f"MLMarker training atlas: {ml_atlas.shape} (proteins × tissues)")
  79. # %% [markdown]
  80. # ## 4. Define Tissue Mappings
  81. #
  82. # Map cancer cohorts to expected tissue of origin for each atlas.
  84. # %%
  85. # Cohorts to evaluate (excluding DLBCL+ as it's redundant with DLBCL)
  86. eval_cohorts = ['CRC', 'DLBCL', 'Glioma', 'OSCC', 'healthy LN', 'healthy OE']
  88. # === MLMarker Mappings ===
  89. cancer_to_mlmarker = {
  90. 'CRC': 'Colon', 'DLBCL': 'B-cells',
  91. 'Glioma': 'Brain', 'OSCC': 'Parotid gland',
  92. 'healthy LN': 'Lymph node', 'healthy OE': 'Parotid gland'
  93. }
  95. cancer_to_mlmarker_extended = {
  96. 'CRC': ['Colon', 'Small intestine', 'Rectum', 'Duodenum', 'Appendix'],
  97. 'DLBCL': ['B-cells', 'Lymph node', 'Spleen'],
  98. 'Glioma': ['Brain', 'Pituitary gland'],
  99. 'OSCC': ['Parotid gland', 'Tonsil'],
  100. 'healthy LN': ['Lymph node', 'B-cells', 'Spleen'],
  101. 'healthy OE': ['Parotid gland', 'Tonsil', 'Esophagus']
  102. }
  104. # === HPA Mappings ===
  105. cancer_to_hpa = {
  106. 'CRC': 'Colon', 'DLBCL': 'Lymph node',
  107. 'Glioma': 'Cerebral cortex', 'OSCC': 'Salivary gland',
  108. 'healthy LN': 'Lymph node', 'healthy OE': 'Oral mucosa'
  109. }
  111. cancer_to_hpa_extended = {
  112. 'CRC': ['Colon', 'Small intestine', 'Rectum', 'Duodenum', 'Appendix'],
  113. 'DLBCL': ['Lymph node', 'Spleen', 'Tonsil', 'Bone marrow'],
  114. 'Glioma': ['Cerebral cortex', 'Cerebellum', 'Hippocampus', 'Caudate', 'Pituitary gland'],
  115. 'OSCC': ['Salivary gland', 'Oral mucosa', 'Nasopharynx'],
  116. 'healthy LN': ['Lymph node', 'Spleen', 'Tonsil'],
  117. 'healthy OE': ['Oral mucosa', 'Salivary gland', 'Esophagus']
  118. }
  120. # === PDB Mappings ===
  121. cancer_to_pdb = {
  122. 'CRC': 'Colon', 'DLBCL': 'B-lymphocyte',
  123. 'Glioma': 'Brain', 'OSCC': 'Salivary gland',
  124. 'healthy LN': 'Lymph node', 'healthy OE': 'Mouth'
  125. }
  127. cancer_to_pdb_extended = {
  128. 'CRC': ['Colon', 'Colon muscle', 'Gut', 'Small intestine'],
  129. 'DLBCL': ['B-lymphocyte', 'Lymph node', 'Spleen'],
  130. 'Glioma': ['Brain', 'Cerebral cortex'],
  131. 'OSCC': ['Salivary gland', 'Mouth'],
  132. 'healthy LN': ['Lymph node', 'B-lymphocyte', 'Spleen'],
  133. 'healthy OE': ['Mouth', 'Salivary gland', 'Esophagus']
  134. }
  136. print(f"Cohorts to evaluate: {eval_cohorts}")
  137. print(f"\nExcluded: Melanoma, HELA, PDAC (no direct healthy tissue equivalent)")
  138. print(f"Excluded: DLBCL+ (redundant with DLBCL)")
  140. # %% [markdown]
  141. # ## 5. Prepare Sample Data
  143. # %%
  144. # Get expression columns
  145. expression_columns = [col for col in protein_df.columns if "MaxLFQ Intensity" in col]
  146. sample_names = [col.split(' MaxLFQ Intensity')[0] for col in expression_columns]
  148. # Create sample metadata
  149. sample_metadata = pd.DataFrame({'sample': sample_names, 'expression_col': expression_columns})
  150. sample_metadata = sample_metadata.merge(
  151. metadata[['Cohort', 'Folder name FP_DDA']],
  152. left_on='sample', right_on='Folder name FP_DDA', how='inner'
  153. ).drop_duplicates()
  154. sample_metadata = sample_metadata[sample_metadata['Cohort'].isin(eval_cohorts)]
  156. # Create expression matrix (log2 transformed)
  157. valid_cols = sample_metadata['expression_col'].tolist()
  158. sample_expr = protein_df[['Protein ID'] + valid_cols].set_index('Protein ID')
  159. sample_expr_log = np.log2(sample_expr + 1)
  161. # Create MLMarker column mapping
  162. sample_mapping_mlmarker = {}
  163. for sample in sample_metadata['sample']:
  164. matches = [c for c in mlmarker_predictions.columns if sample in c]
  165. if matches:
  166. sample_mapping_mlmarker[sample] = matches[0]
  168. print(f"Samples per cohort:")
  169. print(sample_metadata['Cohort'].value_counts())
  170. print(f"\nTotal samples: {len(sample_metadata)}")
  171. print(f"MLMarker mapping found: {len(sample_mapping_mlmarker)} samples")
  173. # %% [markdown]
  174. # ## 6. Define Classification Functions
  176. # %%
  177. def atlas_correlation_classify(sample_profile, atlas_matrix):
  178. """
  179. Classify a sample using Spearman correlation with atlas tissues.
  180. Returns: (predicted_tissue, score, all_scores_dict)
  181. """
  182. common_proteins = sample_profile.index.intersection(atlas_matrix.index)
  183. if len(common_proteins) < 10:
  184. return None, 0, {}
  186. sample_aligned = sample_profile.loc[common_proteins].values
  187. atlas_aligned = atlas_matrix.loc[common_proteins]
  189. correlations = {}
  190. for tissue in atlas_aligned.columns:
  191. tissue_profile = atlas_aligned[tissue].values
  192. if tissue_profile.std() == 0 or sample_aligned.std() == 0:
  193. correlations[tissue] = 0
  194. else:
  195. corr, _ = spearmanr(sample_aligned, tissue_profile)
  196. correlations[tissue] = corr if not np.isnan(corr) else 0
  198. if correlations:
  199. predicted = max(correlations, key=correlations.get)
  200. return predicted, correlations[predicted], correlations
  201. return None, 0, {}
  204. def get_topk_predictions(scores_dict, k=5):
  205. """Get top-k tissue predictions from scores dictionary."""
  206. if not scores_dict:
  207. return []
  208. return sorted(scores_dict.keys(), key=lambda x: scores_dict[x], reverse=True)[:k]
  211. print("Classification functions defined.")
  213. # %% [markdown]
  214. # ## 7. Run Classification (All 4 Methods)
  216. # %%
  217. # Initialize results storage
  218. results = []
  220. for _, row in tqdm(sample_metadata.iterrows(), total=len(sample_metadata), desc="Classifying samples"):
  221. sample_name = row['sample']
  222. cohort = row['Cohort']
  223. expr_col = row['expression_col']
  225. # Get sample profile (non-zero proteins only)
  226. sample_profile = sample_expr_log[expr_col]
  227. sample_profile = sample_profile[sample_profile > 0]
  229. result = {'sample': sample_name, 'cohort': cohort}
  231. # === 1. MLMarker Model ===
  232. if sample_name in sample_mapping_mlmarker:
  233. ml_col = sample_mapping_mlmarker[sample_name]
  234. ml_scores = mlmarker_predictions[ml_col].to_dict()
  235. result['mlmarker_pred'] = mlmarker_predictions[ml_col].idxmax()
  236. else:
  237. ml_scores = {}
  238. result['mlmarker_pred'] = None
  240. true_strict = cancer_to_mlmarker.get(cohort)
  241. true_extended = cancer_to_mlmarker_extended.get(cohort, [])
  243. for k in [1, 2, 3, 4, 5]:
  244. top_k = get_topk_predictions(ml_scores, k)
  245. result[f'mlmarker_top{k}_strict'] = true_strict in top_k if top_k else False
  246. result[f'mlmarker_top{k}_extended'] = any(t in top_k for t in true_extended) if top_k else False
  248. # === 2. MLMarker Training Atlas ===
  249. pred, _, scores = atlas_correlation_classify(sample_profile, ml_atlas)
  250. result['ml_training_pred'] = pred
  252. for k in [1, 2, 3, 4, 5]:
  253. top_k = get_topk_predictions(scores, k)
  254. result[f'ml_training_top{k}_strict'] = true_strict in top_k if top_k else False
  255. result[f'ml_training_top{k}_extended'] = any(t in top_k for t in true_extended) if top_k else False
  257. # === 3. HPA Atlas ===
  258. true_hpa_strict = cancer_to_hpa.get(cohort)
  259. true_hpa_extended = cancer_to_hpa_extended.get(cohort, [])
  261. pred, _, scores = atlas_correlation_classify(sample_profile, hpa_matrix)
  262. result['hpa_pred'] = pred
  264. for k in [1, 2, 3, 4, 5]:
  265. top_k = get_topk_predictions(scores, k)
  266. result[f'hpa_top{k}_strict'] = true_hpa_strict in top_k if top_k else False
  267. result[f'hpa_top{k}_extended'] = any(t in top_k for t in true_hpa_extended) if top_k else False
  269. # === 4. PDB Atlas ===
  270. true_pdb_strict = cancer_to_pdb.get(cohort)
  271. true_pdb_extended = cancer_to_pdb_extended.get(cohort, [])
  273. pred, _, scores = atlas_correlation_classify(sample_profile, pdb_matrix)
  274. result['pdb_pred'] = pred
  276. for k in [1, 2, 3, 4, 5]:
  277. top_k = get_topk_predictions(scores, k)
  278. result[f'pdb_top{k}_strict'] = true_pdb_strict in top_k if top_k else False
  279. result[f'pdb_top{k}_extended'] = any(t in top_k for t in true_pdb_extended) if top_k else False
  281. results.append(result)
  283. results_df = pd.DataFrame(results)
  284. print(f"\nClassification complete! Results: {results_df.shape}")
  286. # %% [markdown]
  287. # ## 8. Calculate Top-k Accuracy Summary
  289. # %%
  290. # Calculate overall accuracy for each method
  291. methods = ['mlmarker', 'ml_training', 'hpa', 'pdb']
  292. method_labels = ['MLMarker', 'MLMarker Training', 'HPA Atlas', 'PDB Atlas']
  294. summary_data = []
  295. for method, label in zip(methods, method_labels):
  296. for k in [1, 2, 3, 4, 5]:
  297. strict_acc = results_df[f'{method}_top{k}_strict'].mean() * 100
  298. extended_acc = results_df[f'{method}_top{k}_extended'].mean() * 100
  299. summary_data.append({
  300. 'Method': label, 'k': k,
  301. 'Strict (%)': strict_acc, 'Extended (%)': extended_acc
  302. })
  304. summary_df = pd.DataFrame(summary_data)
  306. # Display Top-1 and Top-5 summary
  307. print("=" * 70)
  308. print("OVERALL ACCURACY SUMMARY")
  309. print("=" * 70)
  310. print(f"\n{'Method':<20} {'Top-1 Strict':>12} {'Top-1 Extended':>15} {'Top-5 Extended':>15}")
  311. print("-" * 70)
  312. for method, label in zip(methods, method_labels):
  313. t1_strict = results_df[f'{method}_top1_strict'].mean() * 100
  314. t1_ext = results_df[f'{method}_top1_extended'].mean() * 100
  315. t5_ext = results_df[f'{method}_top5_extended'].mean() * 100
  316. print(f"{label:<20} {t1_strict:>11.1f}% {t1_ext:>14.1f}% {t5_ext:>14.1f}%")
  317. print(f"\nTotal samples: {len(results_df)}")
  319. # %% [markdown]
  320. # ## 9. Per-Cohort Accuracy Analysis
  322. # %%
  323. # Calculate per-cohort accuracy
  324. cohort_data = []
  325. for cohort in eval_cohorts:
  326. cohort_mask = results_df['cohort'] == cohort
  327. n_samples = cohort_mask.sum()
  329. for method, label in zip(methods, method_labels):
  330. for k in [1, 3, 5]:
  331. strict_acc = results_df.loc[cohort_mask, f'{method}_top{k}_strict'].mean() * 100
  332. extended_acc = results_df.loc[cohort_mask, f'{method}_top{k}_extended'].mean() * 100
  333. cohort_data.append({
  334. 'Cohort': cohort, 'Method': label, 'k': k, 'n': n_samples,
  335. 'Strict (%)': strict_acc, 'Extended (%)': extended_acc
  336. })
  338. cohort_df = pd.DataFrame(cohort_data)
  340. # Display per-cohort Top-1 Extended accuracy
  341. print("\nPer-Cohort Top-1 Extended Accuracy:")
  342. print("=" * 80)
  343. pivot = cohort_df[cohort_df['k'] == 1].pivot(index='Cohort', columns='Method', values='Extended (%)')
  344. pivot = pivot[method_labels] # Reorder columns
  345. print(pivot.round(1).to_string())
  347. # %% [markdown]
  348. # ## 10. K-Nearest Neighbors (KNN) Classification
  349. #
  350. # Instead of correlating with tissue centroids (means), KNN classifies based on the k nearest training samples. This preserves sample-level variability and may better capture tissue-specific patterns.
  352. # %%
  353. # Prepare MLMarker training data for KNN
  354. # ml_train has samples as rows, proteins as columns, with tissue_name column
  356. # Get protein columns (exclude metadata)
  357. ml_protein_cols = [c for c in ml_train.columns if c not in ['tissue_name', 'disease_status', 'fluid', 'cell_type']]
  358. ml_train_X = ml_train[ml_protein_cols]
  359. ml_train_y = ml_train['tissue_name']
  361. print(f"MLMarker Training data for KNN:")
  362. print(f" Samples: {len(ml_train_X)}")
  363. print(f" Proteins: {len(ml_protein_cols)}")
  364. print(f" Tissues: {ml_train_y.nunique()}")
  365. print(f" Samples per tissue (top 10):")
  366. print(ml_train_y.value_counts().head(10))
  368. # %%
  369. # KNN Classification Function
  370. from sklearn.neighbors import KNeighborsClassifier
  371. from collections import Counter
  373. def knn_classify_with_topk(sample_profile, train_X, train_y, k=5):
  374. """
  375. Classify sample using KNN and return top-k tissue predictions.
  377. Returns: (predicted_tissue, top_k_tissues, vote_counts)
  378. """
  379. # Find common proteins
  380. common_proteins = sample_profile.index.intersection(train_X.columns)
  381. if len(common_proteins) < 10:
  382. return None, [], {}
  384. # Align data
  385. sample_aligned = sample_profile.loc[common_proteins].values.reshape(1, -1)
  386. train_aligned = train_X[common_proteins].values
  387. train_labels = train_y.values
  389. # Handle NaN
  390. sample_aligned = np.nan_to_num(sample_aligned, 0)
  391. train_aligned = np.nan_to_num(train_aligned, 0)
  393. # Fit KNN with k neighbors
  394. knn = KNeighborsClassifier(n_neighbors=min(k, len(train_aligned)), metric='correlation')
  395. knn.fit(train_aligned, train_labels)
  397. # Get k nearest neighbors
  398. distances, indices = knn.kneighbors(sample_aligned)
  399. neighbor_labels = train_labels[indices[0]]
  401. # Count votes
  402. vote_counts = Counter(neighbor_labels)
  404. # Get top-k unique tissues by vote count
  405. sorted_tissues = [t for t, _ in vote_counts.most_common()]
  407. # Get probabilities for all classes
  408. proba = knn.predict_proba(sample_aligned)[0]
  409. class_proba = dict(zip(knn.classes_, proba))
  411. predicted = knn.predict(sample_aligned)[0]
  413. return predicted, sorted_tissues, class_proba
  415. print("KNN classification function defined.")
  417. # %%
  418. # Run KNN on MLMarker Training Data (sample-level)
  419. # Test different k values
  420. k_values = [1, 3, 5, 7, 11]
  422. knn_results = {f'knn_ml_k{k}_pred': [] for k in k_values}
  423. knn_results['sample'] = []
  424. knn_results['cohort'] = []
  426. # Add top-k accuracy tracking for BOTH strict and extended
  427. for k in k_values:
  428. for topk in [1, 2, 3, 4, 5]:
  429. knn_results[f'knn_ml_k{k}_top{topk}_strict'] = []
  430. knn_results[f'knn_ml_k{k}_top{topk}_extended'] = []
  432. print(f"Running KNN classification with k={k_values}...")
  434. for _, row in tqdm(sample_metadata.iterrows(), total=len(sample_metadata), desc="KNN MLMarker Training"):
  435. sample_name = row['sample']
  436. cohort = row['Cohort']
  437. expr_col = row['expression_col']
  439. sample_profile = sample_expr_log[expr_col]
  440. sample_profile = sample_profile[sample_profile > 0]
  442. knn_results['sample'].append(sample_name)
  443. knn_results['cohort'].append(cohort)
  445. # Get both strict and extended ground truth
  446. true_strict = cancer_to_mlmarker.get(cohort)
  447. true_extended = cancer_to_mlmarker_extended.get(cohort, [])
  449. for k in k_values:
  450. pred, top_tissues, class_proba = knn_classify_with_topk(sample_profile, ml_train_X, ml_train_y, k=k)
  451. knn_results[f'knn_ml_k{k}_pred'].append(pred)
  453. # Calculate top-k accuracy using class probabilities
  454. if class_proba:
  455. sorted_by_proba = sorted(class_proba.keys(), key=lambda x: class_proba[x], reverse=True)
  456. for topk in [1, 2, 3, 4, 5]:
  457. top_k_tissues = sorted_by_proba[:topk]
  458. # Strict: exact 1-to-1 mapping
  459. is_correct_strict = true_strict in top_k_tissues if true_strict else False
  460. knn_results[f'knn_ml_k{k}_top{topk}_strict'].append(is_correct_strict)
  461. # Extended: 1-to-many mapping
  462. is_correct_extended = any(t in top_k_tissues for t in true_extended)
  463. knn_results[f'knn_ml_k{k}_top{topk}_extended'].append(is_correct_extended)
  464. else:
  465. for topk in [1, 2, 3, 4, 5]:
  466. knn_results[f'knn_ml_k{k}_top{topk}_strict'].append(False)
  467. knn_results[f'knn_ml_k{k}_top{topk}_extended'].append(False)
  469. knn_ml_df = pd.DataFrame(knn_results)
  470. print(f"\nKNN results shape: {knn_ml_df.shape}")
  472. # %%
  473. # For HPA and PDB: Use tissue profiles as "samples" for KNN
  474. # This treats each tissue as a single sample
  476. # HPA KNN (using tissue profiles)
  477. # hpa_matrix is proteins × tissues, so transpose to get tissues × proteins (samples × features)
  478. hpa_train_X = hpa_matrix.T # Now: tissues × proteins (each tissue is a "sample")
  479. hpa_train_y = pd.Series(hpa_train_X.index, index=hpa_train_X.index) # Labels = tissue names
  481. # PDB KNN
  482. pdb_train_X = pdb_matrix.T # tissues × proteins
  483. pdb_train_y = pd.Series(pdb_train_X.index, index=pdb_train_X.index)
  485. print(f"HPA for KNN: {hpa_train_X.shape[0]} tissues (samples), {hpa_train_X.shape[1]} proteins (features)")
  486. print(f"PDB for KNN: {pdb_train_X.shape[0]} tissues (samples), {pdb_train_X.shape[1]} proteins (features)")
  488. # Run KNN on HPA and PDB
  489. knn_atlas_results = {'sample': [], 'cohort': [],
  490. 'knn_hpa_pred': [], 'knn_pdb_pred': []}
  492. # Initialize columns for BOTH strict and extended
  493. for topk in [1, 2, 3, 4, 5]:
  494. knn_atlas_results[f'knn_hpa_top{topk}_strict'] = []
  495. knn_atlas_results[f'knn_hpa_top{topk}_extended'] = []
  496. knn_atlas_results[f'knn_pdb_top{topk}_strict'] = []
  497. knn_atlas_results[f'knn_pdb_top{topk}_extended'] = []
  499. for _, row in tqdm(sample_metadata.iterrows(), total=len(sample_metadata), desc="KNN HPA/PDB"):
  500. sample_name = row['sample']
  501. cohort = row['Cohort']
  502. expr_col = row['expression_col']
  504. sample_profile = sample_expr_log[expr_col]
  505. sample_profile = sample_profile[sample_profile > 0]
  507. knn_atlas_results['sample'].append(sample_name)
  508. knn_atlas_results['cohort'].append(cohort)
  510. # HPA KNN - get both strict and extended ground truth
  511. true_hpa_strict = cancer_to_hpa.get(cohort)
  512. true_hpa_extended = cancer_to_hpa_extended.get(cohort, [])
  513. pred, _, class_proba = knn_classify_with_topk(sample_profile, hpa_train_X, hpa_train_y, k=1)
  514. knn_atlas_results['knn_hpa_pred'].append(pred)
  516. if class_proba:
  517. sorted_by_proba = sorted(class_proba.keys(), key=lambda x: class_proba[x], reverse=True)
  518. for topk in [1, 2, 3, 4, 5]:
  519. top_k_tissues = sorted_by_proba[:topk]
  520. is_correct_strict = true_hpa_strict in top_k_tissues if true_hpa_strict else False
  521. knn_atlas_results[f'knn_hpa_top{topk}_strict'].append(is_correct_strict)
  522. is_correct_extended = any(t in top_k_tissues for t in true_hpa_extended)
  523. knn_atlas_results[f'knn_hpa_top{topk}_extended'].append(is_correct_extended)
  524. else:
  525. for topk in [1, 2, 3, 4, 5]:
  526. knn_atlas_results[f'knn_hpa_top{topk}_strict'].append(False)
  527. knn_atlas_results[f'knn_hpa_top{topk}_extended'].append(False)
  529. # PDB KNN - get both strict and extended ground truth
  530. true_pdb_strict = cancer_to_pdb.get(cohort)
  531. true_pdb_extended = cancer_to_pdb_extended.get(cohort, [])
  532. pred, _, class_proba = knn_classify_with_topk(sample_profile, pdb_train_X, pdb_train_y, k=1)
  533. knn_atlas_results['knn_pdb_pred'].append(pred)
  535. if class_proba:
  536. sorted_by_proba = sorted(class_proba.keys(), key=lambda x: class_proba[x], reverse=True)
  537. for topk in [1, 2, 3, 4, 5]:
  538. top_k_tissues = sorted_by_proba[:topk]
  539. is_correct_strict = true_pdb_strict in top_k_tissues if true_pdb_strict else False
  540. knn_atlas_results[f'knn_pdb_top{topk}_strict'].append(is_correct_strict)
  541. is_correct_extended = any(t in top_k_tissues for t in true_pdb_extended)
  542. knn_atlas_results[f'knn_pdb_top{topk}_extended'].append(is_correct_extended)
  543. else:
  544. for topk in [1, 2, 3, 4, 5]:
  545. knn_atlas_results[f'knn_pdb_top{topk}_strict'].append(False)
  546. knn_atlas_results[f'knn_pdb_top{topk}_extended'].append(False)
  548. knn_atlas_df = pd.DataFrame(knn_atlas_results)
  549. print(f"\nKNN Atlas results shape: {knn_atlas_df.shape}")
  551. # %%
  552. # For HPA and PDB: Use tissue profiles as "samples" for KNN
  553. # This treats each tissue as a single sample
  555. # HPA KNN (using tissue profiles)
  556. hpa_train_X = hpa_matrix.T # tissues × proteins
  557. hpa_train_y = pd.Series(hpa_train_X.index, index=hpa_train_X.index)
  559. # PDB KNN
  560. pdb_train_X = pdb_matrix.T # tissues × proteins
  561. pdb_train_y = pd.Series(pdb_train_X.index, index=pdb_train_X.index)
  563. print(f"HPA for KNN: {hpa_train_X.shape[0]} tissues (samples), {hpa_train_X.shape[1]} proteins (features)")
  564. print(f"PDB for KNN: {pdb_train_X.shape[0]} tissues (samples), {pdb_train_X.shape[1]} proteins (features)")
  566. # Run KNN on HPA and PDB with multiple k values
  567. k_values_atlas = [1, 3, 5, 7, 11, 15, 20]
  568. knn_atlas_results = {'sample': [], 'cohort': []}
  570. # Initialize columns for all k values and both strict/extended
  571. for k in k_values_atlas:
  572. knn_atlas_results[f'knn_hpa_k{k}_pred'] = []
  573. knn_atlas_results[f'knn_pdb_k{k}_pred'] = []
  574. for topk in [1, 2, 3, 4, 5]:
  575. knn_atlas_results[f'knn_hpa_k{k}_top{topk}_strict'] = []
  576. knn_atlas_results[f'knn_hpa_k{k}_top{topk}_extended'] = []
  577. knn_atlas_results[f'knn_pdb_k{k}_top{topk}_strict'] = []
  578. knn_atlas_results[f'knn_pdb_k{k}_top{topk}_extended'] = []
  580. for _, row in tqdm(sample_metadata.iterrows(), total=len(sample_metadata), desc="KNN HPA/PDB"):
  581. sample_name = row['sample']
  582. cohort = row['Cohort']
  583. expr_col = row['expression_col']
  585. sample_profile = sample_expr_log[expr_col]
  586. sample_profile = sample_profile[sample_profile > 0]
  588. knn_atlas_results['sample'].append(sample_name)
  589. knn_atlas_results['cohort'].append(cohort)
  591. # Ground truth for HPA and PDB
  592. true_hpa_strict = cancer_to_hpa.get(cohort)
  593. true_hpa_extended = cancer_to_hpa_extended.get(cohort, [])
  594. true_pdb_strict = cancer_to_pdb.get(cohort)
  595. true_pdb_extended = cancer_to_pdb_extended.get(cohort, [])
  598. # Test different k values for HPA
  599. for k in k_values_atlas:
  600. pred_hpa, _, class_proba_hpa = knn_classify_with_topk(sample_profile, hpa_train_X, hpa_train_y, k=k)
  601. knn_atlas_results[f'knn_hpa_k{k}_pred'].append(pred_hpa)
  603. if class_proba_hpa:
  604. sorted_by_proba = sorted(class_proba_hpa.keys(), key=lambda x: class_proba_hpa[x], reverse=True)
  605. for topk in [1, 2, 3, 4, 5]:
  606. top_k_tissues = sorted_by_proba[:topk]
  607. is_correct_strict = true_hpa_strict in top_k_tissues if true_hpa_strict else False
  608. is_correct_extended = any(t in top_k_tissues for t in true_hpa_extended)
  609. knn_atlas_results[f'knn_hpa_k{k}_top{topk}_strict'].append(is_correct_strict)
  610. knn_atlas_results[f'knn_hpa_k{k}_top{topk}_extended'].append(is_correct_extended)
  611. else:
  612. for topk in [1, 2, 3, 4, 5]:
  613. knn_atlas_results[f'knn_hpa_k{k}_top{topk}_strict'].append(False)
  614. knn_atlas_results[f'knn_hpa_k{k}_top{topk}_extended'].append(False)
  616. # Test different k values for PDB
  617. for k in k_values_atlas:
  618. pred_pdb, _, class_proba_pdb = knn_classify_with_topk(sample_profile, pdb_train_X, pdb_train_y, k=k)
  619. knn_atlas_results[f'knn_pdb_k{k}_pred'].append(pred_pdb)
  621. if class_proba_pdb:
  622. sorted_by_proba = sorted(class_proba_pdb.keys(), key=lambda x: class_proba_pdb[x], reverse=True)
  623. for topk in [1, 2, 3, 4, 5]:
  624. top_k_tissues = sorted_by_proba[:topk]
  625. is_correct_strict = true_pdb_strict in top_k_tissues if true_pdb_strict else False
  626. is_correct_extended = any(t in top_k_tissues for t in true_pdb_extended)
  627. knn_atlas_results[f'knn_pdb_k{k}_top{topk}_strict'].append(is_correct_strict)
  628. knn_atlas_results[f'knn_pdb_k{k}_top{topk}_extended'].append(is_correct_extended)
  629. else:
  630. for topk in [1, 2, 3, 4, 5]:
  631. knn_atlas_results[f'knn_pdb_k{k}_top{topk}_strict'].append(False)
  632. knn_atlas_results[f'knn_pdb_k{k}_top{topk}_extended'].append(False)
  634. knn_atlas_df_multik = pd.DataFrame(knn_atlas_results)
  635. print(f"\nKNN Atlas results shape: {knn_atlas_df_multik.shape}")
  637. # %%
  638. # KNN Atlas: Top-1 Strict Accuracy vs k value
  639. fig, ax = plt.subplots(figsize=(10, 6))
  641. # HPA KNN - Top-1 Strict
  642. hpa_knn_strict_accs = [knn_atlas_df_multik[f'knn_hpa_k{k}_top5_extended'].mean() * 100 for k in k_values_atlas]
  643. ax.plot(k_values_atlas, hpa_knn_strict_accs, 'o-', color='#2A9D8F', linewidth=2.5,
  644. markersize=8, label='HPA (KNN)')
  646. # PDB KNN - Top-1 Strict
  647. pdb_knn_strict_accs = [knn_atlas_df_multik[f'knn_pdb_k{k}_top5_extended'].mean() * 100 for k in k_values_atlas]
  648. ax.plot(k_values_atlas, pdb_knn_strict_accs, 'D-', color='#F4A261', linewidth=2.5,
  649. markersize=8, label='PDB (KNN)')
  653. ax.set_xlabel('k (Number of Neighbors)', fontsize=12, fontweight='bold')
  654. ax.set_ylabel('Top-1 Strict Accuracy (%)', fontsize=12, fontweight='bold')
  655. ax.set_title('KNN Atlas Methods: Top-1 Strict Accuracy vs k', fontsize=14, fontweight='bold')
  656. ax.set_xticks(k_values_atlas)
  657. ax.set_ylim(0, 100)
  658. ax.grid(True, alpha=0.3)
  659. ax.legend(loc='best', fontsize=11)
  661. plt.tight_layout()
  662. # plt.savefig('knn_atlas_top5_extended_vs_k.png', dpi=300, bbox_inches='tight')
  663. # plt.savefig('knn_atlas_top5_extended_vs_k.svg', bbox_inches='tight')
  664. plt.show()
  666. print("KNN Atlas - Top-1 Strict Accuracy by k:")
  667. print("=" * 50)
  668. print(f"{'k':<5} {'HPA (%)':<12} {'PDB (%)':<12}")
  669. print("-" * 50)
  670. for k, hpa_acc, pdb_acc in zip(k_values_atlas, hpa_knn_strict_accs, pdb_knn_strict_accs):
  671. print(f"{k:<5} {hpa_acc:<11.1f} {pdb_acc:<11.1f}")
  672. print("-" * 50)
  673. # Spearman baselines computed in later cell (see cells 26-27)
  675. # %% [markdown]
  676. # ## 10. Figure 1: Top-k Accuracy Line Plot
  678. # %%
  679. # Color palette for methods - manuscript style
  680. # MLMarker: reddish, Training atlas: bluish, HPA: greenish, PDB: yellowish
  681. colors_spearman = {'MLMarker': '#E63946', 'Training Atlas': '#457B9D',
  682. 'HPA Atlas': '#2A9D8F', 'PDB Atlas': '#E9C46A'}
  684. # Extended colors for Spearman vs KNN comparison
  685. colors_all = {
  686. 'MLMarker': '#E63946',
  687. 'Training Atlas (Spearman)': '#457B9D',
  688. 'Training Atlas (KNN)': '#1D3557',
  689. 'HPA (Spearman)': '#2A9D8F',
  690. 'HPA (KNN)': '#264653',
  691. 'PDB (Spearman)': '#E9C46A',
  692. 'PDB (KNN)': '#F4A261'
  693. }
  695. # ============================================
  696. # Figure 1: Overall Top-k Accuracy - Spearman vs KNN Comparison
  697. # ============================================
  698. fig, axes = plt.subplots(1, 2, figsize=(14, 5))
  700. # Prepare data - need to recalculate with filtered cohorts
  701. # Filter results to only include eval_cohorts
  702. results_filtered = results_df[results_df['cohort'].isin(eval_cohorts)]
  703. knn_ml_filtered = knn_ml_df[knn_ml_df['cohort'].isin(eval_cohorts)]
  704. knn_atlas_filtered = knn_atlas_df[knn_atlas_df['cohort'].isin(eval_cohorts)]
  706. # Plot 1: Extended accuracy (Extended mapping - more meaningful)
  707. ax = axes[0]
  709. # MLMarker Model
  710. ml_topk_ext = [results_filtered[f'mlmarker_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  711. ax.plot([1,2,3,4,5], ml_topk_ext, 'o-', color=colors_all['MLMarker'],
  712. linewidth=2.5, markersize=8, label='MLMarker Model')
  714. # Training Atlas - Spearman
  715. ml_train_spear = [results_filtered[f'ml_training_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  716. ax.plot([1,2,3,4,5], ml_train_spear, 's-', color=colors_all['Training Atlas (Spearman)'],
  717. linewidth=2, markersize=7, label='Training Atlas (Spearman)')
  719. # Training Atlas - KNN k=5
  720. ml_train_knn = [knn_ml_filtered[f'knn_ml_k5_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  721. ax.plot([1,2,3,4,5], ml_train_knn, 's-', color=colors_all['Training Atlas (KNN)'],
  722. linewidth=2, markersize=7, label='Training Atlas (KNN)')
  724. # HPA - Spearman
  725. hpa_spear = [results_filtered[f'hpa_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  726. ax.plot([1,2,3,4,5], hpa_spear, '^-', color=colors_all['HPA (Spearman)'],
  727. linewidth=2, markersize=7, label='HPA (Spearman)')
  729. # HPA - KNN
  730. hpa_knn = [knn_atlas_filtered[f'knn_hpa_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  731. ax.plot([1,2,3,4,5], hpa_knn, '^-', color=colors_all['HPA (KNN)'],
  732. linewidth=2, markersize=7, label='HPA (KNN)')
  734. # PDB - Spearman
  735. pdb_spear = [results_filtered[f'pdb_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  736. ax.plot([1,2,3,4,5], pdb_spear, 'D-', color=colors_all['PDB (Spearman)'],
  737. linewidth=2, markersize=7, label='PDB (Spearman)')
  739. # PDB - KNN
  740. pdb_knn = [knn_atlas_filtered[f'knn_pdb_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  741. ax.plot([1,2,3,4,5], pdb_knn, 'D-', color=colors_all['PDB (KNN)'],
  742. linewidth=2, markersize=7, label='PDB (KNN)')
  744. ax.set_xlabel('k (Top-k)', fontsize=12)
  745. ax.set_ylabel('Extended Accuracy (%)', fontsize=12)
  746. ax.set_title('A) Overall Top-k Accuracy', fontsize=12, fontweight='bold')
  747. ax.set_xticks([1, 2, 3, 4, 5])
  748. ax.set_ylim(0, 100)
  749. ax.grid(True, alpha=0.3)
  751. # Plot 2: Top-1 accuracy bar comparison
  752. ax = axes[1]
  753. methods_bar = ['MLMarker', 'Training Atlas\n(Spearman)', 'Training Atlas\n(KNN)',
  754. 'HPA\n(Spearman)', 'HPA\n(KNN)', 'PDB\n(Spearman)', 'PDB\n(KNN)']
  755. accs_bar = [ml_topk_ext[0], ml_train_spear[0], ml_train_knn[0],
  756. hpa_spear[0], hpa_knn[0], pdb_spear[0], pdb_knn[0]]
  757. bar_colors = [colors_all['MLMarker'], colors_all['Training Atlas (Spearman)'], colors_all['Training Atlas (KNN)'],
  758. colors_all['HPA (Spearman)'], colors_all['HPA (KNN)'],
  759. colors_all['PDB (Spearman)'], colors_all['PDB (KNN)']]
  761. bars = ax.bar(methods_bar, accs_bar, color=bar_colors, edgecolor='black', linewidth=0.5)
  762. for bar in bars:
  763. height = bar.get_height()
  764. ax.annotate(f'{height:.1f}%', xy=(bar.get_x() + bar.get_width()/2, height),
  765. xytext=(0, 3), textcoords='offset points', ha='center', va='bottom', fontsize=9)
  767. ax.set_ylabel('Top-1 Extended Accuracy (%)', fontsize=12)
  768. ax.set_title('B) Top-1 Accuracy Comparison', fontsize=12, fontweight='bold')
  769. ax.set_ylim(-5, 100) # Start below 0 to see all bars clearly
  770. ax.tick_params(axis='x', rotation=45)
  771. ax.grid(True, alpha=0.3, axis='y')
  773. # Shared legend
  774. handles, labels = axes[0].get_legend_handles_labels()
  775. fig.legend(handles, labels, loc='center right', bbox_to_anchor=(1.18, 0.5), fontsize=9)
  777. plt.tight_layout()
  778. plt.savefig('baseline_spearman_vs_knn_overall.png', dpi=300, bbox_inches='tight')
  779. plt.savefig('baseline_spearman_vs_knn_overall.svg', bbox_inches='tight')
  780. plt.show()
  781. print("Saved: baseline_spearman_vs_knn_overall.png/svg")
  783. # ============================================
  784. # Figure 2: Per-Cohort Top-k Accuracy - Spearman vs KNN
  785. # ============================================
  786. fig, axes = plt.subplots(2, len(eval_cohorts), figsize=(18, 8), sharey=True)
  788. for col_idx, cohort in enumerate(eval_cohorts):
  789. # Filter data for this cohort
  790. cohort_results = results_filtered[results_filtered['cohort'] == cohort]
  791. cohort_knn_ml = knn_ml_filtered[knn_ml_filtered['cohort'] == cohort]
  792. cohort_knn_atlas = knn_atlas_filtered[knn_atlas_filtered['cohort'] == cohort]
  793. n_samples = len(cohort_results)
  795. # Row 1: Spearman methods
  796. ax = axes[0, col_idx]
  798. # MLMarker Model
  799. ml_topk = [cohort_results[f'mlmarker_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  800. ax.plot([1,2,3,4,5], ml_topk, 'o-', color=colors_all['MLMarker'], linewidth=2, markersize=6, label='MLMarker')
  802. # Training Atlas Spearman
  803. ml_train = [cohort_results[f'ml_training_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  804. ax.plot([1,2,3,4,5], ml_train, 's-', color=colors_all['Training Atlas (Spearman)'], linewidth=1.5, markersize=5)
  806. # HPA Spearman
  807. hpa = [cohort_results[f'hpa_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  808. ax.plot([1,2,3,4,5], hpa, '^-', color=colors_all['HPA (Spearman)'], linewidth=1.5, markersize=5)
  810. # PDB Spearman
  811. pdb = [cohort_results[f'pdb_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  812. ax.plot([1,2,3,4,5], pdb, 'D-', color=colors_all['PDB (Spearman)'], linewidth=1.5, markersize=5)
  814. ax.set_title(f'{cohort}\n(n={n_samples})', fontsize=10, fontweight='bold')
  815. ax.set_xticks([1, 2, 3, 4, 5])
  816. ax.set_ylim(-5, 105) # Start below 0
  817. ax.grid(True, alpha=0.3)
  818. if col_idx == 0:
  819. ax.set_ylabel('Spearman (%)', fontsize=10)
  821. # Row 2: KNN methods
  822. ax = axes[1, col_idx]
  824. # MLMarker Model (same as above for reference)
  825. ax.plot([1,2,3,4,5], ml_topk, 'o-', color=colors_all['MLMarker'], linewidth=2, markersize=6)
  827. # Training Atlas KNN
  828. ml_knn = [cohort_knn_ml[f'knn_ml_k5_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  829. ax.plot([1,2,3,4,5], ml_knn, 's-', color=colors_all['Training Atlas (KNN)'], linewidth=1.5, markersize=5)
  831. # HPA KNN
  832. hpa_k = [cohort_knn_atlas[f'knn_hpa_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  833. ax.plot([1,2,3,4,5], hpa_k, '^-', color=colors_all['HPA (KNN)'], linewidth=1.5, markersize=5)
  835. # PDB KNN
  836. pdb_k = [cohort_knn_atlas[f'knn_pdb_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  837. ax.plot([1,2,3,4,5], pdb_k, 'D-', color=colors_all['PDB (KNN)'], linewidth=1.5, markersize=5)
  839. ax.set_xlabel('k', fontsize=10)
  840. ax.set_xticks([1, 2, 3, 4, 5])
  841. ax.set_ylim(-5, 105) # Start below 0
  842. ax.grid(True, alpha=0.3)
  843. if col_idx == 0:
  844. ax.set_ylabel('KNN (%)', fontsize=10)
  846. # Row labels
  847. axes[0, 0].annotate('Spearman\nCorrelation', xy=(-0.45, 0.5), xycoords='axes fraction',
  848. fontsize=11, fontweight='bold', ha='center', va='center', rotation=90)
  849. axes[1, 0].annotate('KNN\nClassification', xy=(-0.45, 0.5), xycoords='axes fraction',
  850. fontsize=11, fontweight='bold', ha='center', va='center', rotation=90)
  852. # Legend
  853. legend_labels = ['MLMarker', 'Training Atlas', 'HPA Atlas', 'PDB Atlas']
  854. legend_colors = [colors_all['MLMarker'], colors_all['Training Atlas (Spearman)'],
  855. colors_all['HPA (Spearman)'], colors_all['PDB (Spearman)']]
  856. handles = [plt.Line2D([0], [0], color=c, linewidth=2, marker='o', markersize=6) for c in legend_colors]
  857. fig.legend(handles, legend_labels, loc='lower center', bbox_to_anchor=(0.5, -0.02), ncol=4, fontsize=11)
  859. plt.suptitle('Top-k Accuracy by Cohort: Spearman vs KNN', fontsize=14, fontweight='bold', y=1.02)
  860. plt.tight_layout()
  861. # plt.savefig('baseline_spearman_vs_knn_per_cohort.png', dpi=300, bbox_inches='tight')
  862. # plt.savefig('baseline_spearman_vs_knn_per_cohort.svg', bbox_inches='tight')
  863. plt.show()
  864. # print("Saved: baseline_spearman_vs_knn_per_cohort.png/svg")
  866. # %%
  867. # ============================================
  868. # STRICT (1-to-1) MAPPING VERSION
  869. # Same figures as above but using strict mapping instead of extended
  870. # ============================================
  872. # Color palette for methods - manuscript style
  873. colors_all = {
  874. 'MLMarker': '#E63946',
  875. 'Training Atlas (Spearman)': '#457B9D',
  876. 'Training Atlas (KNN)': '#1D3557',
  877. 'HPA (Spearman)': '#2A9D8F',
  878. 'HPA (KNN)': '#264653',
  879. 'PDB (Spearman)': '#E9C46A',
  880. 'PDB (KNN)': '#F4A261'
  881. }
  883. # ============================================
  884. # Figure 1: Overall Top-k Accuracy - Spearman vs KNN (STRICT 1-to-1)
  885. # ============================================
  886. fig, axes = plt.subplots(1, 2, figsize=(14, 5))
  888. # Prepare data - need to recalculate with filtered cohorts
  889. results_filtered = results_df[results_df['cohort'].isin(eval_cohorts)]
  890. knn_ml_filtered = knn_ml_df[knn_ml_df['cohort'].isin(eval_cohorts)]
  891. knn_atlas_filtered = knn_atlas_df[knn_atlas_df['cohort'].isin(eval_cohorts)]
  893. # Plot 1: Strict accuracy (1-to-1 mapping)
  894. ax = axes[0]
  896. # MLMarker Model - STRICT
  897. ml_topk_strict = [results_filtered[f'mlmarker_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  898. ax.plot([1,2,3,4,5], ml_topk_strict, 'o-', color=colors_all['MLMarker'],
  899. linewidth=2.5, markersize=8, label='MLMarker Model')
  901. # Training Atlas - Spearman STRICT
  902. ml_train_spear_strict = [results_filtered[f'ml_training_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  903. ax.plot([1,2,3,4,5], ml_train_spear_strict, 's-', color=colors_all['Training Atlas (Spearman)'],
  904. linewidth=2, markersize=7, label='Training Atlas (Spearman)')
  906. # Training Atlas - KNN k=5 STRICT
  907. ml_train_knn_strict = [knn_ml_filtered[f'knn_ml_k5_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  908. ax.plot([1,2,3,4,5], ml_train_knn_strict, 's-', color=colors_all['Training Atlas (KNN)'],
  909. linewidth=2, markersize=7, label='Training Atlas (KNN)')
  911. # HPA - Spearman STRICT
  912. hpa_spear_strict = [results_filtered[f'hpa_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  913. ax.plot([1,2,3,4,5], hpa_spear_strict, '^-', color=colors_all['HPA (Spearman)'],
  914. linewidth=2, markersize=7, label='HPA (Spearman)')
  916. # HPA - KNN STRICT
  917. hpa_knn_strict = [knn_atlas_filtered[f'knn_hpa_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  918. ax.plot([1,2,3,4,5], hpa_knn_strict, '^-', color=colors_all['HPA (KNN)'],
  919. linewidth=2, markersize=7, label='HPA (KNN)')
  921. # PDB - Spearman STRICT
  922. pdb_spear_strict = [results_filtered[f'pdb_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  923. ax.plot([1,2,3,4,5], pdb_spear_strict, 'D-', color=colors_all['PDB (Spearman)'],
  924. linewidth=2, markersize=7, label='PDB (Spearman)')
  926. # PDB - KNN STRICT
  927. pdb_knn_strict = [knn_atlas_filtered[f'knn_pdb_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  928. ax.plot([1,2,3,4,5], pdb_knn_strict, 'D-', color=colors_all['PDB (KNN)'],
  929. linewidth=2, markersize=7, label='PDB (KNN)')
  931. ax.set_xlabel('k (Top-k)', fontsize=12)
  932. ax.set_ylabel('Strict Accuracy (%)', fontsize=12)
  933. ax.set_title('A) Overall Top-k Accuracy (1-to-1 Mapping)', fontsize=12, fontweight='bold')
  934. ax.set_xticks([1, 2, 3, 4, 5])
  935. ax.set_ylim(0, 100)
  936. ax.grid(True, alpha=0.3)
  938. # Plot 2: Top-1 accuracy bar comparison (STRICT)
  939. ax = axes[1]
  940. methods_bar = ['MLMarker', 'Training Atlas\n(Spearman)', 'Training Atlas\n(KNN)',
  941. 'HPA\n(Spearman)', 'HPA\n(KNN)', 'PDB\n(Spearman)', 'PDB\n(KNN)']
  942. accs_bar_strict = [ml_topk_strict[0], ml_train_spear_strict[0], ml_train_knn_strict[0],
  943. hpa_spear_strict[0], hpa_knn_strict[0], pdb_spear_strict[0], pdb_knn_strict[0]]
  944. bar_colors = [colors_all['MLMarker'], colors_all['Training Atlas (Spearman)'], colors_all['Training Atlas (KNN)'],
  945. colors_all['HPA (Spearman)'], colors_all['HPA (KNN)'],
  946. colors_all['PDB (Spearman)'], colors_all['PDB (KNN)']]
  948. bars = ax.bar(methods_bar, accs_bar_strict, color=bar_colors, edgecolor='black', linewidth=0.5)
  949. for bar in bars:
  950. height = bar.get_height()
  951. ax.annotate(f'{height:.1f}%', xy=(bar.get_x() + bar.get_width()/2, height),
  952. xytext=(0, 3), textcoords='offset points', ha='center', va='bottom', fontsize=9)
  954. ax.set_ylabel('Top-1 Strict Accuracy (%)', fontsize=12)
  955. ax.set_title('B) Top-1 Accuracy Comparison (1-to-1 Mapping)', fontsize=12, fontweight='bold')
  956. ax.set_ylim(-5, 100)
  957. ax.tick_params(axis='x', rotation=45)
  958. ax.grid(True, alpha=0.3, axis='y')
  960. # Shared legend
  961. handles, labels = axes[0].get_legend_handles_labels()
  962. fig.legend(handles, labels, loc='center right', bbox_to_anchor=(1.18, 0.5), fontsize=9)
  964. plt.tight_layout()
  965. plt.savefig('baseline_spearman_vs_knn_overall_strict.png', dpi=300, bbox_inches='tight')
  966. plt.savefig('baseline_spearman_vs_knn_overall_strict.svg', bbox_inches='tight')
  967. plt.show()
  968. print("Saved: baseline_spearman_vs_knn_overall_strict.png/svg")
  970. # ============================================
  971. # Figure 2: Per-Cohort Top-k Accuracy - Spearman vs KNN (STRICT 1-to-1)
  972. # ============================================
  973. fig, axes = plt.subplots(2, len(eval_cohorts), figsize=(18, 8), sharey=True)
  975. for col_idx, cohort in enumerate(eval_cohorts):
  976. # Filter data for this cohort
  977. cohort_results = results_filtered[results_filtered['cohort'] == cohort]
  978. cohort_knn_ml = knn_ml_filtered[knn_ml_filtered['cohort'] == cohort]
  979. cohort_knn_atlas = knn_atlas_filtered[knn_atlas_filtered['cohort'] == cohort]
  980. n_samples = len(cohort_results)
  982. # Row 1: Spearman methods (STRICT)
  983. ax = axes[0, col_idx]
  985. # MLMarker Model
  986. ml_topk = [cohort_results[f'mlmarker_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  987. ax.plot([1,2,3,4,5], ml_topk, 'o-', color=colors_all['MLMarker'], linewidth=2, markersize=6, label='MLMarker')
  989. # Training Atlas Spearman
  990. ml_train = [cohort_results[f'ml_training_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  991. ax.plot([1,2,3,4,5], ml_train, 's-', color=colors_all['Training Atlas (Spearman)'], linewidth=1.5, markersize=5)
  993. # HPA Spearman
  994. hpa = [cohort_results[f'hpa_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  995. ax.plot([1,2,3,4,5], hpa, '^-', color=colors_all['HPA (Spearman)'], linewidth=1.5, markersize=5)
  997. # PDB Spearman
  998. pdb = [cohort_results[f'pdb_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  999. ax.plot([1,2,3,4,5], pdb, 'D-', color=colors_all['PDB (Spearman)'], linewidth=1.5, markersize=5)
  1001. ax.set_title(f'{cohort}\n(n={n_samples})', fontsize=10, fontweight='bold')
  1002. ax.set_xticks([1, 2, 3, 4, 5])
  1003. ax.set_ylim(-5, 105)
  1004. ax.grid(True, alpha=0.3)
  1005. if col_idx == 0:
  1006. ax.set_ylabel('Spearman (%)', fontsize=10)
  1008. # Row 2: KNN methods (STRICT)
  1009. ax = axes[1, col_idx]
  1011. # MLMarker Model (same as above for reference)
  1012. ax.plot([1,2,3,4,5], ml_topk, 'o-', color=colors_all['MLMarker'], linewidth=2, markersize=6)
  1014. # Training Atlas KNN
  1015. ml_knn = [cohort_knn_ml[f'knn_ml_k5_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  1016. ax.plot([1,2,3,4,5], ml_knn, 's-', color=colors_all['Training Atlas (KNN)'], linewidth=1.5, markersize=5)
  1018. # HPA KNN
  1019. hpa_k = [cohort_knn_atlas[f'knn_hpa_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  1020. ax.plot([1,2,3,4,5], hpa_k, '^-', color=colors_all['HPA (KNN)'], linewidth=1.5, markersize=5)
  1022. # PDB KNN
  1023. pdb_k = [cohort_knn_atlas[f'knn_pdb_top{k}_strict'].mean()*100 for k in [1,2,3,4,5]]
  1024. ax.plot([1,2,3,4,5], pdb_k, 'D-', color=colors_all['PDB (KNN)'], linewidth=1.5, markersize=5)
  1026. ax.set_xlabel('k', fontsize=10)
  1027. ax.set_xticks([1, 2, 3, 4, 5])
  1028. ax.set_ylim(-5, 105)
  1029. ax.grid(True, alpha=0.3)
  1030. if col_idx == 0:
  1031. ax.set_ylabel('KNN (%)', fontsize=10)
  1033. # Row labels
  1034. axes[0, 0].annotate('Spearman\nCorrelation', xy=(-0.45, 0.5), xycoords='axes fraction',
  1035. fontsize=11, fontweight='bold', ha='center', va='center', rotation=90)
  1036. axes[1, 0].annotate('KNN\nClassification', xy=(-0.45, 0.5), xycoords='axes fraction',
  1037. fontsize=11, fontweight='bold', ha='center', va='center', rotation=90)
  1039. # Legend
  1040. legend_labels = ['MLMarker', 'Training Atlas', 'HPA Atlas', 'PDB Atlas']
  1041. legend_colors = [colors_all['MLMarker'], colors_all['Training Atlas (Spearman)'],
  1042. colors_all['HPA (Spearman)'], colors_all['PDB (Spearman)']]
  1043. handles = [plt.Line2D([0], [0], color=c, linewidth=2, marker='o', markersize=6) for c in legend_colors]
  1044. fig.legend(handles, legend_labels, loc='lower center', bbox_to_anchor=(0.5, -0.02), ncol=4, fontsize=11)
  1046. plt.suptitle('Top-k Accuracy by Cohort: Spearman vs KNN (1-to-1 Mapping)', fontsize=14, fontweight='bold', y=1.02)
  1047. plt.tight_layout()
  1048. # plt.savefig('baseline_spearman_vs_knn_per_cohort_strict.png', dpi=300, bbox_inches='tight')
  1049. # plt.savefig('baseline_spearman_vs_knn_per_cohort_strict.svg', bbox_inches='tight')
  1050. plt.show()
  1051. # print("Saved: baseline_spearman_vs_knn_per_cohort_strict.png/svg")
  1053. # %%
  1054. # ============================================
  1055. # FIGURE 1: Grouped Bar Chart - Top-1 Accuracy (Extended vs Strict)
  1056. # ============================================
  1058. import matplotlib.gridspec as gridspec
  1059. import seaborn as sns
  1060. from matplotlib.lines import Line2D
  1061. from matplotlib.patches import Patch
  1063. # Color palette - manuscript colors
  1064. colors_atlas = {
  1065. 'MLMarker': '#E63946',
  1066. 'Training Atlas': '#457B9D',
  1067. 'HPA': '#2A9D8F',
  1068. 'PDB': '#E9C46A'
  1069. }
  1071. # Filter data
  1072. results_filtered = results_df[results_df['cohort'].isin(eval_cohorts)]
  1073. knn_ml_filtered = knn_ml_df[knn_ml_df['cohort'].isin(eval_cohorts)]
  1074. knn_atlas_filtered = knn_atlas_df[knn_atlas_df['cohort'].isin(eval_cohorts)]
  1076. # Calculate Top-1 accuracies for all methods (overall averages)
  1077. # Spearman methods
  1078. ml_ext = results_filtered['mlmarker_top1_extended'].mean() * 100
  1079. ml_strict = results_filtered['mlmarker_top1_strict'].mean() * 100
  1080. train_spear_ext = results_filtered['ml_training_top1_extended'].mean() * 100
  1081. train_spear_strict = results_filtered['ml_training_top1_strict'].mean() * 100
  1082. hpa_spear_ext = results_filtered['hpa_top1_extended'].mean() * 100
  1083. hpa_spear_strict = results_filtered['hpa_top1_strict'].mean() * 100
  1084. pdb_spear_ext = results_filtered['pdb_top1_extended'].mean() * 100
  1085. pdb_spear_strict = results_filtered['pdb_top1_strict'].mean() * 100
  1087. # KNN methods
  1088. train_knn_ext = knn_ml_filtered['knn_ml_k5_top1_extended'].mean() * 100
  1089. train_knn_strict = knn_ml_filtered['knn_ml_k5_top1_strict'].mean() * 100
  1090. hpa_knn_ext = knn_atlas_filtered['knn_hpa_top1_extended'].mean() * 100
  1091. hpa_knn_strict = knn_atlas_filtered['knn_hpa_top1_strict'].mean() * 100
  1092. pdb_knn_ext = knn_atlas_filtered['knn_pdb_top1_extended'].mean() * 100
  1093. pdb_knn_strict = knn_atlas_filtered['knn_pdb_top1_strict'].mean() * 100
  1095. # Method labels and data
  1096. method_labels_bar = ['MLMarker', 'Training\n(Spearman)', 'Training\n(KNN)', 'HPA\n(Spearman)', 'HPA\n(KNN)', 'PDB\n(Spearman)', 'PDB\n(KNN)']
  1097. ext_values = [ml_ext, train_spear_ext, train_knn_ext, hpa_spear_ext, hpa_knn_ext, pdb_spear_ext, pdb_knn_ext]
  1098. strict_values = [ml_strict, train_spear_strict, train_knn_strict, hpa_spear_strict, hpa_knn_strict, pdb_spear_strict, pdb_knn_strict]
  1099. bar_colors = [colors_atlas['MLMarker'], colors_atlas['Training Atlas'], colors_atlas['Training Atlas'],
  1100. colors_atlas['HPA'], colors_atlas['HPA'], colors_atlas['PDB'], colors_atlas['PDB']]
  1102. fig, ax = plt.subplots(figsize=(10, 5))
  1104. x = np.arange(len(method_labels_bar))
  1105. width = 0.35
  1107. bars_ext = ax.bar(x - width/2, ext_values, width, label='Extended (1-to-many)', color=bar_colors, alpha=1.0, edgecolor='black', linewidth=1)
  1108. bars_strict = ax.bar(x + width/2, strict_values, width, label='Strict (1-to-1)', color=bar_colors, alpha=0.5, edgecolor='black', linewidth=1, hatch='//')
  1110. ax.set_ylabel('Top-1 Accuracy (%)', fontsize=12, fontweight='bold')
  1111. ax.set_title('Top-1 Accuracy Comparison: Extended vs Strict Mapping', fontsize=13, fontweight='bold')
  1112. ax.set_xticks(x)
  1113. ax.set_xticklabels(method_labels_bar, fontsize=10)
  1114. ax.set_ylim(0, 105)
  1115. ax.grid(True, alpha=0.3, axis='y')
  1117. # Add value labels on bars
  1118. for bar in bars_ext:
  1119. height = bar.get_height()
  1120. ax.annotate(f'{height:.1f}', xy=(bar.get_x() + bar.get_width()/2, height),
  1121. xytext=(0, 3), textcoords="offset points", ha='center', va='bottom', fontsize=9)
  1122. for bar in bars_strict:
  1123. height = bar.get_height()
  1124. ax.annotate(f'{height:.1f}', xy=(bar.get_x() + bar.get_width()/2, height),
  1125. xytext=(0, 3), textcoords="offset points", ha='center', va='bottom', fontsize=9)
  1127. # Custom legend for bar chart
  1128. legend_elements = [Patch(facecolor='gray', edgecolor='black', alpha=1.0, label='Extended (1-to-many)'),
  1129. Patch(facecolor='gray', edgecolor='black', alpha=0.5, hatch='//', label='Strict (1-to-1)')]
  1130. ax.legend(handles=legend_elements, loc='upper right', fontsize=10)
  1132. plt.tight_layout()
  1133. # plt.savefig('../Figures_paper/fig11a.png', dpi=300, bbox_inches='tight')
  1134. # plt.savefig('baseline_barplot_extended_vs_strict.svg', bbox_inches='tight')
  1135. plt.show()
  1136. # print("Saved: ../Figures_paper/fig11a.png/svg")
  1138. # %%
  1139. # ============================================
  1140. # FIGURE 2: Per-Cohort Line Plots (Extended and Strict)
  1141. # ============================================
  1143. fig, axes = plt.subplots(2, len(eval_cohorts), figsize=(16, 7), sharey=True)
  1145. row_labels = ['Strict (1-to-1)', 'Extended (1-to-many)']
  1146. mapping_types = ['strict', 'extended']
  1148. for row_idx, mapping_type in enumerate(mapping_types):
  1149. for col_idx, cohort in enumerate(eval_cohorts):
  1150. ax = axes[row_idx, col_idx]
  1152. cohort_results = results_filtered[results_filtered['cohort'] == cohort]
  1153. cohort_knn_ml = knn_ml_filtered[knn_ml_filtered['cohort'] == cohort]
  1154. cohort_knn_atlas = knn_atlas_filtered[knn_atlas_filtered['cohort'] == cohort]
  1155. n_samples = len(cohort_results)
  1157. # MLMarker (only model - solid line)
  1158. ml_topk = [cohort_results[f'mlmarker_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1159. ax.plot([1,2,3,4,5], ml_topk, 'o-', color=colors_atlas['MLMarker'], linewidth=2, markersize=5, label='MLMarker')
  1161. # Training Atlas - Spearman (solid) and KNN (dashed)
  1162. ml_train_spear = [cohort_results[f'ml_training_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1163. ml_train_knn = [cohort_knn_ml[f'knn_ml_k5_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1164. ax.plot([1,2,3,4,5], ml_train_spear, 's-', color=colors_atlas['Training Atlas'], linewidth=1.5, markersize=4, label='Training (Spearman)')
  1165. ax.plot([1,2,3,4,5], ml_train_knn, 's--', color=colors_atlas['Training Atlas'], linewidth=1.5, markersize=4, label='Training (KNN)')
  1167. # HPA - Spearman (solid) and KNN (dashed)
  1168. hpa_spear = [cohort_results[f'hpa_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1169. hpa_knn = [cohort_knn_atlas[f'knn_hpa_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1170. ax.plot([1,2,3,4,5], hpa_spear, '^-', color=colors_atlas['HPA'], linewidth=1.5, markersize=4, label='HPA (Spearman)')
  1171. ax.plot([1,2,3,4,5], hpa_knn, '^--', color=colors_atlas['HPA'], linewidth=1.5, markersize=4, label='HPA (KNN)')
  1173. # PDB - Spearman (solid) and KNN (dashed)
  1174. pdb_spear = [cohort_results[f'pdb_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1175. pdb_knn = [cohort_knn_atlas[f'knn_pdb_top{k}_{mapping_type}'].mean()*100 for k in [1,2,3,4,5]]
  1176. ax.plot([1,2,3,4,5], pdb_spear, 'D-', color=colors_atlas['PDB'], linewidth=1.5, markersize=4, label='PDB (Spearman)')
  1177. ax.plot([1,2,3,4,5], pdb_knn, 'D--', color=colors_atlas['PDB'], linewidth=1.5, markersize=4, label='PDB (KNN)')
  1179. ax.set_xticks([1, 2, 3, 4, 5])
  1180. ax.set_ylim(-5, 105)
  1181. ax.grid(True, alpha=0.3)
  1183. # Column titles (cohort names) only on first row
  1184. if row_idx == 0:
  1185. ax.set_title(f'{cohort}\n(n={n_samples})', fontsize=11, fontweight='bold')
  1187. # Y-axis label only on first column
  1188. if col_idx == 0:
  1189. ax.set_ylabel(f'{row_labels[row_idx]}\nAccuracy (%)', fontsize=10, fontweight='bold')
  1191. # X-axis label only on last row
  1192. if row_idx == 1:
  1193. ax.set_xlabel('k (Top-k)', fontsize=10)
  1195. # Create legend outside the plots
  1196. style_handles = [
  1197. Line2D([0], [0], color='gray', linestyle='-', linewidth=2, label='Spearman'),
  1198. Line2D([0], [0], color='gray', linestyle='--', linewidth=2, label='KNN'),
  1199. Line2D([0], [0], marker='o', color='w', markerfacecolor=colors_atlas['MLMarker'], markersize=10, label='MLMarker'),
  1200. Line2D([0], [0], marker='s', color='w', markerfacecolor=colors_atlas['Training Atlas'], markersize=10, label='Training Atlas'),
  1201. Line2D([0], [0], marker='^', color='w', markerfacecolor=colors_atlas['HPA'], markersize=10, label='HPA'),
  1202. Line2D([0], [0], marker='D', color='w', markerfacecolor=colors_atlas['PDB'], markersize=10, label='PDB'),
  1203. ]
  1204. fig.legend(handles=style_handles, loc='center right', bbox_to_anchor=(1.12, 0.5), fontsize=10, frameon=True, title='Method')
  1206. plt.suptitle('Top-k Accuracy by Cohort', fontsize=14, fontweight='bold', y=1.02)
  1207. plt.tight_layout()
  1208. # plt.savefig('../Figures_paper/fig11c.png', dpi=300, bbox_inches='tight')
  1209. # plt.savefig('baseline_topk_per_cohort.svg', bbox_inches='tight')
  1210. plt.show()
  1211. print("Saved: ../Figures_paper/fig11c.png/svg")
  1213. # %%
  1214. # ============================================
  1215. # FIGURE 3: Heatmap - Top-1 Extended Accuracy by Cohort and Method
  1216. # ============================================
  1218. # Prepare heatmap data
  1219. heatmap_methods = ['MLMarker', 'Training (Spearman)', 'Training (KNN)', 'HPA (Spearman)', 'HPA (KNN)', 'PDB (Spearman)', 'PDB (KNN)']
  1220. heatmap_data = []
  1222. for cohort in eval_cohorts:
  1223. cohort_results = results_filtered[results_filtered['cohort'] == cohort]
  1224. cohort_knn_ml = knn_ml_filtered[knn_ml_filtered['cohort'] == cohort]
  1225. cohort_knn_atlas = knn_atlas_filtered[knn_atlas_filtered['cohort'] == cohort]
  1227. row = [
  1228. cohort_results['mlmarker_top1_extended'].mean() * 100,
  1229. cohort_results['ml_training_top1_extended'].mean() * 100,
  1230. cohort_knn_ml['knn_ml_k5_top1_extended'].mean() * 100,
  1231. cohort_results['hpa_top1_extended'].mean() * 100,
  1232. cohort_knn_atlas['knn_hpa_top1_extended'].mean() * 100,
  1233. cohort_results['pdb_top1_extended'].mean() * 100,
  1234. cohort_knn_atlas['knn_pdb_top1_extended'].mean() * 100,
  1235. ]
  1236. heatmap_data.append(row)
  1238. heatmap_df = pd.DataFrame(heatmap_data, index=eval_cohorts, columns=heatmap_methods)
  1240. # Colorblind-friendly colormap
  1241. cmap = sns.color_palette('YlGn', as_cmap=True)
  1243. fig, ax = plt.subplots(figsize=(10, 5))
  1245. sns.heatmap(heatmap_df, annot=True, fmt='.1f', cmap=cmap,
  1246. vmin=0, vmax=100, ax=ax, cbar_kws={'label': 'Accuracy (%)', 'shrink': 0.8},
  1247. linewidths=0.5, linecolor='white', annot_kws={'fontsize': 11})
  1248. ax.set_title('Top-1 Extended Accuracy by Cohort and Method', fontsize=13, fontweight='bold')
  1249. ax.set_xlabel('Method', fontsize=12)
  1250. ax.set_ylabel('Cohort', fontsize=12)
  1251. ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha='right', fontsize=10)
  1252. ax.set_yticklabels(ax.get_yticklabels(), rotation=0, fontsize=10)
  1254. plt.tight_layout()
  1255. # plt.savefig('../Figures_paper/fig11b.png', dpi=300, bbox_inches='tight')
  1256. # plt.savefig('baseline_heatmap_extended.svg', bbox_inches='tight')
  1257. plt.show()
  1258. print("Saved: ../Figures_paper/fig11b.png/svg")
  1260. # %%
  1263. # %% [markdown]
  1264. # ## 11. Figure 2: Per-Cohort Bar Plot (Top-1 Extended)
  1266. # %%
  1267. # Filter to Top-1 Extended
  1268. bar_data = cohort_df[(cohort_df['k'] == 1)].copy()
  1270. fig, ax = plt.subplots(figsize=(14, 6))
  1272. x = np.arange(len(eval_cohorts))
  1273. width = 0.2
  1275. for i, method in enumerate(method_labels):
  1276. method_data = bar_data[bar_data['Method'] == method].set_index('Cohort').loc[eval_cohorts]
  1277. bars = ax.bar(x + i*width, method_data['Extended (%)'], width,
  1278. label=method, color=colors[method], edgecolor='black', linewidth=0.5)
  1279. # Add value labels
  1280. for bar in bars:
  1281. height = bar.get_height()
  1282. ax.annotate(f'{height:.0f}', xy=(bar.get_x() + bar.get_width()/2, height),
  1283. xytext=(0, 3), textcoords='offset points', ha='center', va='bottom', fontsize=8)
  1285. ax.set_xlabel('Cancer Cohort', fontsize=12)
  1286. ax.set_ylabel('Top-1 Extended Accuracy (%)', fontsize=12)
  1287. ax.set_title('Classification Accuracy by Cohort (Extended Mapping)', fontsize=14, fontweight='bold')
  1288. ax.set_xticks(x + width * 1.5)
  1289. ax.set_xticklabels(eval_cohorts, rotation=45, ha='right')
  1290. ax.set_ylim(0, 115)
  1291. ax.legend(loc='upper right')
  1292. ax.grid(True, alpha=0.3, axis='y')
  1294. plt.tight_layout()
  1295. plt.savefig('baseline_cohort_barplot.png', dpi=300, bbox_inches='tight')
  1296. plt.savefig('baseline_cohort_barplot.svg', bbox_inches='tight')
  1297. plt.show()
  1298. print("Saved: baseline_cohort_barplot.png/svg")
  1300. # %% [markdown]
  1301. # ## 12. Figure 3: Per-Cohort Top-k Line Plots
  1303. # %%
  1304. # Calculate full k range per cohort
  1305. cohort_topk_data = []
  1306. for cohort in eval_cohorts:
  1307. cohort_mask = results_df['cohort'] == cohort
  1308. for method, label in zip(methods, method_labels):
  1309. for k in [1, 2, 3, 4, 5]:
  1310. extended_acc = results_df.loc[cohort_mask, f'{method}_top{k}_extended'].mean() * 100
  1311. cohort_topk_data.append({'Cohort': cohort, 'Method': label, 'k': k, 'Extended (%)': extended_acc})
  1313. cohort_topk_df = pd.DataFrame(cohort_topk_data)
  1315. # Create subplot grid
  1316. fig, axes = plt.subplots(2, 4, figsize=(16, 8))
  1317. axes = axes.flatten()
  1319. for idx, cohort in enumerate(eval_cohorts):
  1320. ax = axes[idx]
  1321. cohort_data = cohort_topk_df[cohort_topk_df['Cohort'] == cohort]
  1323. for method in method_labels:
  1324. method_data = cohort_data[cohort_data['Method'] == method]
  1325. ax.plot(method_data['k'], method_data['Extended (%)'], 'o-',
  1326. label=method, color=colors[method], linewidth=2, markersize=6)
  1328. n_samples = (results_df['cohort'] == cohort).sum()
  1329. ax.set_title(f'{cohort} (n={n_samples})', fontsize=11, fontweight='bold')
  1330. ax.set_xlabel('k')
  1331. ax.set_ylabel('Accuracy (%)')
  1332. ax.set_xticks([1, 2, 3, 4, 5])
  1333. ax.set_ylim(0, 105)
  1334. ax.grid(True, alpha=0.3)
  1336. # Hide empty subplot
  1337. axes[-1].axis('off')
  1339. # Add shared legend
  1340. handles, labels = axes[0].get_legend_handles_labels()
  1341. fig.legend(handles, labels, loc='lower right', bbox_to_anchor=(0.98, 0.12), fontsize=10)
  1343. plt.suptitle('Top-k Accuracy by Cohort (Extended Mapping)', fontsize=14, fontweight='bold', y=1.02)
  1344. plt.tight_layout()
  1345. plt.savefig('baseline_cohort_topk_lines.png', dpi=300, bbox_inches='tight')
  1346. plt.savefig('baseline_cohort_topk_lines.svg', bbox_inches='tight')
  1347. plt.show()
  1348. print("Saved: baseline_cohort_topk_lines.png/svg")
  1350. # %% [markdown]
  1351. # ## 13. Figure 4: Accuracy Heatmap
  1353. # %%
  1354. # Create heatmap data (Top-1 Extended)
  1355. heatmap_data = bar_data.pivot(index='Cohort', columns='Method', values='Extended (%)')
  1356. heatmap_data = heatmap_data.loc[eval_cohorts, method_labels] # Reorder
  1358. fig, ax = plt.subplots(figsize=(10, 6))
  1359. sns.heatmap(heatmap_data, annot=True, fmt='.1f', cmap='RdYlGn',
  1360. vmin=0, vmax=100, ax=ax, cbar_kws={'label': 'Accuracy (%)'},
  1361. linewidths=0.5, linecolor='white')
  1362. ax.set_title('Top-1 Extended Accuracy by Cohort and Method', fontsize=14, fontweight='bold')
  1363. ax.set_xlabel('Method', fontsize=12)
  1364. ax.set_ylabel('Cohort', fontsize=12)
  1365. plt.xticks(rotation=45, ha='right')
  1367. plt.tight_layout()
  1368. plt.savefig('baseline_accuracy_heatmap.png', dpi=300, bbox_inches='tight')
  1369. plt.savefig('baseline_accuracy_heatmap.svg', bbox_inches='tight')
  1370. plt.show()
  1371. print("Saved: baseline_accuracy_heatmap.png/svg")
  1373. # %% [markdown]
  1374. # ## 14. Confusion Matrices
  1376. # %%
  1377. from sklearn.metrics import confusion_matrix
  1379. def plot_confusion_matrix(y_true, y_pred, title, ax, normalize=True):
  1380. """Plot confusion matrix with cohort labels."""
  1381. labels = sorted(set(y_true) | set(y_pred))
  1382. cm = confusion_matrix(y_true, y_pred, labels=labels)
  1384. if normalize:
  1385. cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
  1386. cm = np.nan_to_num(cm) # Handle division by zero
  1388. sns.heatmap(cm, annot=False, fmt='.2f' if normalize else 'd', cmap='Blues',
  1389. xticklabels=labels, yticklabels=labels, ax=ax,
  1390. vmin=0, vmax=1 if normalize else None)
  1391. ax.set_title(title, fontsize=11, fontweight='bold')
  1392. ax.set_xlabel('Predicted')
  1393. ax.set_ylabel('True Cohort')
  1395. # Create confusion matrices for each method
  1396. fig, axes = plt.subplots(2, 2, figsize=(16, 14))
  1398. # MLMarker: Use predicted tissue grouped by cohort
  1399. ax = axes[0, 0]
  1400. y_true = results_df['cohort'].values
  1401. y_pred = results_df['mlmarker_pred'].fillna('Unknown').values
  1402. plot_confusion_matrix(y_true, y_pred, 'MLMarker Model', ax)
  1404. ax = axes[0, 1]
  1405. y_pred = results_df['ml_training_pred'].fillna('Unknown').values
  1406. plot_confusion_matrix(y_true, y_pred, 'MLMarker Training Atlas', ax)
  1408. ax = axes[1, 0]
  1409. y_pred = results_df['hpa_pred'].fillna('Unknown').values
  1410. plot_confusion_matrix(y_true, y_pred, 'HPA Atlas', ax)
  1412. ax = axes[1, 1]
  1413. y_pred = results_df['pdb_pred'].fillna('Unknown').values
  1414. plot_confusion_matrix(y_true, y_pred, 'PDB Atlas', ax)
  1416. plt.suptitle('Confusion Matrices: True Cohort vs Predicted Tissue', fontsize=14, fontweight='bold')
  1417. plt.tight_layout()
  1418. plt.savefig('baseline_confusion_matrices.png', dpi=300, bbox_inches='tight')
  1419. plt.savefig('baseline_confusion_matrices.svg', bbox_inches='tight')
  1420. plt.show()
  1421. print("Saved: baseline_confusion_matrices.png/svg")
  1423. # %% [markdown]
  1424. # ## 15. Summary Statistics and Save Results
  1426. # %%
  1427. # Final summary
  1428. print("=" * 80)
  1429. print("FINAL SUMMARY: MLMarker vs Atlas-Correlation Baselines")
  1430. print("=" * 80)
  1432. mlmarker_t1 = results_df['mlmarker_top1_extended'].mean() * 100
  1433. mltraining_t1 = results_df['ml_training_top1_extended'].mean() * 100
  1434. hpa_t1 = results_df['hpa_top1_extended'].mean() * 100
  1435. pdb_t1 = results_df['pdb_top1_extended'].mean() * 100
  1437. print(f"\nTop-1 Extended Accuracy:")
  1438. print(f" MLMarker Model: {mlmarker_t1:.1f}%")
  1439. print(f" MLMarker Training: {mltraining_t1:.1f}% (Δ = {mlmarker_t1 - mltraining_t1:+.1f}%)")
  1440. print(f" HPA Atlas: {hpa_t1:.1f}% (Δ = {mlmarker_t1 - hpa_t1:+.1f}%)")
  1441. print(f" PDB Atlas: {pdb_t1:.1f}% (Δ = {mlmarker_t1 - pdb_t1:+.1f}%)")
  1443. print(f"\nEvaluated on {len(results_df)} samples across {len(eval_cohorts)} cohorts")
  1444. print(f"Cohorts: {', '.join(eval_cohorts)}")
  1446. # Save results
  1447. summary_df.to_csv('baseline_comparison_summary.csv', index=False)
  1448. cohort_df.to_csv('baseline_comparison_cohort.csv', index=False)
  1449. results_df.to_csv('baseline_comparison_detailed.csv', index=False)
  1451. print("\nResults saved:")
  1452. print(" - baseline_comparison_summary.csv")
  1453. print(" - baseline_comparison_cohort.csv")
  1454. print(" - baseline_comparison_detailed.csv")
  1456. # %% [markdown]
  1457. # ## 16. Additional Baseline Methods
  1458. #
  1459. # ### 16.1 Nearest Centroid Classification
  1460. # Classify samples by Euclidean distance to tissue centroids.
  1462. # %%
  1463. from scipy.spatial.distance import euclidean, cdist
  1464. from sklearn.neighbors import KNeighborsClassifier
  1465. from sklearn.preprocessing import StandardScaler
  1467. def nearest_centroid_classify(sample_profile, atlas_matrix):
  1468. """Classify by Euclidean distance to tissue centroids."""
  1469. common_proteins = sample_profile.index.intersection(atlas_matrix.index)
  1470. if len(common_proteins) < 10:
  1471. return None, {}, {}
  1473. sample_aligned = sample_profile.loc[common_proteins].values
  1474. atlas_aligned = atlas_matrix.loc[common_proteins]
  1476. distances = {}
  1477. for tissue in atlas_aligned.columns:
  1478. tissue_profile = atlas_aligned[tissue].values
  1479. dist = euclidean(sample_aligned, tissue_profile)
  1480. distances[tissue] = -dist # Negative so higher = closer
  1482. if distances:
  1483. predicted = max(distances, key=distances.get)
  1484. return predicted, distances[predicted], distances
  1485. return None, 0, {}
  1487. # Run Nearest Centroid classification on all atlases
  1488. nc_results = {'sample': [], 'cohort': [],
  1489. 'nc_ml_pred': [], 'nc_hpa_pred': [], 'nc_pdb_pred': []}
  1491. for _, row in tqdm(sample_metadata.iterrows(), total=len(sample_metadata), desc="Nearest Centroid"):
  1492. sample_name = row['sample']
  1493. cohort = row['Cohort']
  1494. expr_col = row['expression_col']
  1495. sample_profile = sample_expr_log[expr_col]
  1496. sample_profile = sample_profile[sample_profile > 0]
  1498. nc_results['sample'].append(sample_name)
  1499. nc_results['cohort'].append(cohort)
  1501. # MLMarker Training Atlas
  1502. pred, _, scores = nearest_centroid_classify(sample_profile, ml_atlas)
  1503. nc_results['nc_ml_pred'].append(pred)
  1505. # HPA Atlas
  1506. pred, _, scores = nearest_centroid_classify(sample_profile, hpa_matrix)
  1507. nc_results['nc_hpa_pred'].append(pred)
  1509. # PDB Atlas
  1510. pred, _, scores = nearest_centroid_classify(sample_profile, pdb_matrix)
  1511. nc_results['nc_pdb_pred'].append(pred)
  1513. nc_df = pd.DataFrame(nc_results)
  1515. # Calculate accuracy
  1516. print("Nearest Centroid Classification Results (Top-1 Extended):")
  1517. print("=" * 60)
  1519. # MLMarker Training
  1520. correct = sum(nc_df.apply(lambda r: r['nc_ml_pred'] in cancer_to_mlmarker_extended.get(r['cohort'], []), axis=1))
  1521. print(f"MLMarker Training Atlas (NC): {correct/len(nc_df)*100:.1f}%")
  1523. # HPA
  1524. correct = sum(nc_df.apply(lambda r: r['nc_hpa_pred'] in cancer_to_hpa_extended.get(r['cohort'], []), axis=1))
  1525. print(f"HPA Atlas (NC): {correct/len(nc_df)*100:.1f}%")
  1527. # PDB
  1528. correct = sum(nc_df.apply(lambda r: r['nc_pdb_pred'] in cancer_to_pdb_extended.get(r['cohort'], []), axis=1))
  1529. print(f"PDB Atlas (NC): {correct/len(nc_df)*100:.1f}%")
  1531. # %%
  1532. # Figure: KNN vs Correlation Comparison
  1533. fig, axes = plt.subplots(1, 2, figsize=(14, 5))
  1535. # Plot 1: MLMarker Training - KNN k sensitivity
  1536. ax = axes[0]
  1537. knn_accs = [knn_ml_df[f'knn_ml_k{k}_top1_extended'].mean()*100 for k in k_values]
  1538. spearman_acc = results_df['ml_training_top1_extended'].mean()*100
  1539. mlmarker_acc = results_df['mlmarker_top1_extended'].mean()*100
  1541. ax.plot(k_values, knn_accs, 'o-', color='#457B9D', linewidth=2, markersize=8, label='KNN')
  1542. ax.axhline(y=spearman_acc, color='#2A9D8F', linestyle='--', linewidth=2, label='Spearman Correlation')
  1543. ax.axhline(y=mlmarker_acc, color='#E63946', linestyle='-', linewidth=2, label='MLMarker Model')
  1545. ax.set_xlabel('k (Number of Neighbors)', fontsize=12)
  1546. ax.set_ylabel('Top-1 Extended Accuracy (%)', fontsize=12)
  1547. ax.set_title('MLMarker Training Data: KNN vs Correlation', fontsize=14, fontweight='bold')
  1548. ax.set_xticks(k_values)
  1549. ax.set_ylim(0, 100)
  1550. ax.legend(loc='lower right')
  1551. ax.grid(True, alpha=0.3)
  1553. # Plot 2: Comparison across all methods
  1554. ax = axes[1]
  1555. methods_compare = ['MLMarker\nModel', 'ML Train\n(Spearman)', 'ML Train\n(KNN k=5)',
  1556. 'HPA\n(Spearman)', 'HPA\n(KNN)', 'PDB\n(Spearman)', 'PDB\n(KNN)']
  1557. accs_compare = [
  1558. mlmarker_acc,
  1559. spearman_acc,
  1560. knn_ml_df['knn_ml_k5_top1_extended'].mean()*100,
  1561. results_df['hpa_top1_extended'].mean()*100,
  1562. knn_atlas_df['knn_hpa_top1_extended'].mean()*100,
  1563. results_df['pdb_top1_extended'].mean()*100,
  1564. knn_atlas_df['knn_pdb_top1_extended'].mean()*100
  1565. ]
  1566. colors_compare = ['#E63946', '#457B9D', '#1D3557', '#2A9D8F', '#264653', '#E9C46A', '#F4A261']
  1568. bars = ax.bar(methods_compare, accs_compare, color=colors_compare, edgecolor='black', linewidth=0.5)
  1569. for bar in bars:
  1570. height = bar.get_height()
  1571. ax.annotate(f'{height:.1f}%', xy=(bar.get_x() + bar.get_width()/2, height),
  1572. xytext=(0, 3), textcoords='offset points', ha='center', va='bottom', fontsize=9)
  1574. ax.set_ylabel('Top-1 Extended Accuracy (%)', fontsize=12)
  1575. ax.set_title('KNN vs Spearman Correlation', fontsize=14, fontweight='bold')
  1576. ax.set_ylim(0, 100)
  1577. ax.tick_params(axis='x', rotation=45)
  1578. ax.grid(True, alpha=0.3, axis='y')
  1580. plt.tight_layout()
  1581. plt.savefig('baseline_knn_comparison.png', dpi=300, bbox_inches='tight')
  1582. plt.savefig('baseline_knn_comparison.svg', bbox_inches='tight')
  1583. plt.show()
  1584. print("Saved: baseline_knn_comparison.png/svg")
  1586. # %%
  1587. # Per-cohort KNN accuracy breakdown
  1588. print("\nPer-Cohort KNN Accuracy (Top-1 Extended):")
  1589. print("=" * 100)
  1590. print(f"{'Cohort':<12} {'MLMarker':>10} {'ML Spear':>10} {'ML KNN-5':>10} {'HPA Spear':>10} {'HPA KNN':>10} {'PDB Spear':>10} {'PDB KNN':>10}")
  1591. print("-" * 100)
  1593. for cohort in eval_cohorts:
  1594. mask = results_df['cohort'] == cohort
  1595. knn_mask = knn_ml_df['cohort'] == cohort
  1596. atlas_mask = knn_atlas_df['cohort'] == cohort
  1598. ml_model = results_df.loc[mask, 'mlmarker_top1_extended'].mean() * 100
  1599. ml_spear = results_df.loc[mask, 'ml_training_top1_extended'].mean() * 100
  1600. ml_knn = knn_ml_df.loc[knn_mask, 'knn_ml_k5_top1_extended'].mean() * 100
  1601. hpa_spear = results_df.loc[mask, 'hpa_top1_extended'].mean() * 100
  1602. hpa_knn = knn_atlas_df.loc[atlas_mask, 'knn_hpa_top1_extended'].mean() * 100
  1603. pdb_spear = results_df.loc[mask, 'pdb_top1_extended'].mean() * 100
  1604. pdb_knn = knn_atlas_df.loc[atlas_mask, 'knn_pdb_top1_extended'].mean() * 100
  1606. print(f"{cohort:<12} {ml_model:>9.1f}% {ml_spear:>9.1f}% {ml_knn:>9.1f}% {hpa_spear:>9.1f}% {hpa_knn:>9.1f}% {pdb_spear:>9.1f}% {pdb_knn:>9.1f}%")
  1608. # %%
  1609. # Top-k accuracy curves for KNN methods
  1610. fig, ax = plt.subplots(figsize=(10, 6))
  1612. # MLMarker Model
  1613. ml_topk = [results_df[f'mlmarker_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1614. ax.plot([1,2,3,4,5], ml_topk, 'o-', color='#E63946', linewidth=2, markersize=8, label='MLMarker Model')
  1616. # ML Training Spearman
  1617. ml_train_topk = [results_df[f'ml_training_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1618. ax.plot([1,2,3,4,5], ml_train_topk, 's--', color='#457B9D', linewidth=2, markersize=8, label='ML Training (Spearman)')
  1620. # ML Training KNN k=5
  1621. ml_knn_topk = [knn_ml_df[f'knn_ml_k5_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1622. ax.plot([1,2,3,4,5], ml_knn_topk, '^--', color='#1D3557', linewidth=2, markersize=8, label='ML Training (KNN k=5)')
  1624. # HPA Spearman
  1625. hpa_topk = [results_df[f'hpa_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1626. ax.plot([1,2,3,4,5], hpa_topk, 's--', color='#2A9D8F', linewidth=2, markersize=8, label='HPA (Spearman)')
  1628. # HPA KNN
  1629. hpa_knn_topk = [knn_atlas_df[f'knn_hpa_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1630. ax.plot([1,2,3,4,5], hpa_knn_topk, '^--', color='#264653', linewidth=2, markersize=8, label='HPA (KNN)')
  1632. # PDB Spearman
  1633. pdb_topk = [results_df[f'pdb_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1634. ax.plot([1,2,3,4,5], pdb_topk, 's--', color='#E9C46A', linewidth=2, markersize=8, label='PDB (Spearman)')
  1636. # PDB KNN
  1637. pdb_knn_topk = [knn_atlas_df[f'knn_pdb_top{k}_extended'].mean()*100 for k in [1,2,3,4,5]]
  1638. ax.plot([1,2,3,4,5], pdb_knn_topk, '^--', color='#F4A261', linewidth=2, markersize=8, label='PDB (KNN)')
  1640. ax.set_xlabel('k (Top-k)', fontsize=12)
  1641. ax.set_ylabel('Extended Accuracy (%)', fontsize=12)
  1642. ax.set_title('Top-k Accuracy: All Methods Including KNN', fontsize=14, fontweight='bold')
  1643. ax.set_xticks([1, 2, 3, 4, 5])
  1644. ax.set_ylim(0, 100)
  1645. ax.legend(loc='lower right', fontsize=9)
  1646. ax.grid(True, alpha=0.3)
  1648. plt.tight_layout()
  1649. # plt.savefig('baseline_topk_with_knn.png', dpi=300, bbox_inches='tight')
  1650. # plt.savefig('baseline_topk_with_knn.svg', bbox_inches='tight')
  1651. plt.show()
  1652. # print("Saved: baseline_topk_with_knn.png/svg")
  1654. # %%
  1655. # Save KNN results
  1656. knn_ml_df.to_csv('baseline_knn_mltraining.csv', index=False)
  1657. knn_atlas_df.to_csv('baseline_knn_atlas.csv', index=False)
  1659. print("KNN results saved:")
  1660. print(" - baseline_knn_mltraining.csv")
  1661. print(" - baseline_knn_atlas.csv")
  1663. # %% [markdown]
  1664. # ### 16.2 UMAP Visualization with Atlas Centroids
  1665. #
  1666. # Project cancer samples and atlas tissue centroids into a shared 2D space to visualize clustering.
  1668. # %%
  1669. import umap
  1671. # Use MLMarker Training Atlas for UMAP (most comparable to MLMarker model)
  1672. # Get common proteins between samples and atlas
  1673. common_proteins = sample_expr_log.index.intersection(ml_atlas.index)
  1674. print(f"Common proteins for UMAP: {len(common_proteins)}")
  1676. # Prepare sample data
  1677. sample_data = sample_expr_log.loc[common_proteins, valid_cols].T
  1678. sample_data.index = sample_metadata['sample'].values
  1680. # Prepare atlas centroids
  1681. atlas_centroids = ml_atlas.loc[common_proteins].T
  1683. # Combine for joint embedding
  1684. combined_data = pd.concat([sample_data, atlas_centroids])
  1685. combined_labels = list(sample_metadata['Cohort'].values) + list(atlas_centroids.index)
  1686. is_centroid = [False] * len(sample_data) + [True] * len(atlas_centroids)
  1688. # Fill NaN with 0 and scale
  1689. combined_filled = combined_data.fillna(0)
  1690. scaler = StandardScaler()
  1691. combined_scaled = scaler.fit_transform(combined_filled)
  1693. # Run UMAP
  1694. print("Running UMAP...")
  1695. reducer = umap.UMAP(n_neighbors=15, min_dist=0.1, random_state=42, n_components=2)
  1696. embedding = reducer.fit_transform(combined_scaled)
  1698. # Create embedding DataFrame
  1699. umap_df = pd.DataFrame({
  1700. 'UMAP1': embedding[:, 0],
  1701. 'UMAP2': embedding[:, 1],
  1702. 'Label': combined_labels,
  1703. 'Is_Centroid': is_centroid
  1704. })
  1706. print(f"UMAP embedding complete: {umap_df.shape}")
  1708. # %%
  1709. # Figure 5: UMAP with samples and expected tissue centroids highlighted
  1710. fig, ax = plt.subplots(figsize=(14, 10))
  1712. # Define colors for cohorts
  1713. cohort_colors = {'CRC': '#E63946', 'DLBCL': '#457B9D', 'DLBCL+': '#1D3557',
  1714. 'Glioma': '#2A9D8F', 'OSCC': '#E9C46A',
  1715. 'healthy LN': '#F4A261', 'healthy OE': '#264653'}
  1717. # Plot samples (smaller, transparent)
  1718. samples_df = umap_df[~umap_df['Is_Centroid']]
  1719. for cohort in eval_cohorts:
  1720. mask = samples_df['Label'] == cohort
  1721. ax.scatter(samples_df.loc[mask, 'UMAP1'], samples_df.loc[mask, 'UMAP2'],
  1722. c=cohort_colors[cohort], alpha=0.6, s=30, label=f'{cohort} samples')
  1724. # Plot expected tissue centroids with stars
  1725. expected_tissues = set()
  1726. for cohort in eval_cohorts:
  1727. expected_tissues.update(cancer_to_mlmarker_extended.get(cohort, []))
  1729. centroids_df = umap_df[umap_df['Is_Centroid']]
  1730. for tissue in expected_tissues:
  1731. if tissue in centroids_df['Label'].values:
  1732. row = centroids_df[centroids_df['Label'] == tissue]
  1733. ax.scatter(row['UMAP1'], row['UMAP2'], marker='*', s=500,
  1734. c='black', edgecolors='white', linewidth=2, zorder=10)
  1735. ax.annotate(tissue, (row['UMAP1'].values[0], row['UMAP2'].values[0]),
  1736. fontsize=9, fontweight='bold',
  1737. xytext=(5, 5), textcoords='offset points')
  1739. # Plot other centroids (grey, smaller)
  1740. other_tissues = set(centroids_df['Label']) - expected_tissues
  1741. for tissue in other_tissues:
  1742. row = centroids_df[centroids_df['Label'] == tissue]
  1743. ax.scatter(row['UMAP1'], row['UMAP2'], marker='*', s=150,
  1744. c='lightgray', edgecolors='gray', linewidth=1, zorder=5)
  1746. ax.set_xlabel('UMAP 1', fontsize=12)
  1747. ax.set_ylabel('UMAP 2', fontsize=12)
  1748. ax.set_title('UMAP: Cancer Samples with MLMarker Training Atlas Tissue Centroids',
  1749. fontsize=14, fontweight='bold')
  1750. ax.legend(loc='upper left', bbox_to_anchor=(1.02, 1), fontsize=9)
  1752. plt.tight_layout()
  1753. plt.savefig('baseline_umap_centroids.png', dpi=300, bbox_inches='tight')
  1754. plt.savefig('baseline_umap_centroids.svg', bbox_inches='tight')
  1755. plt.show()
  1756. print("Saved: baseline_umap_centroids.png/svg")
  1758. # %% [markdown]
  1759. # ### 16.3 Prediction Confidence Analysis
  1760. #
  1761. # Compare the "confidence margin" of correct predictions across methods. A higher margin indicates the model is more certain about its prediction.
  1763. # %%
  1764. # Calculate confidence margins (top1 - top2 score) for each method
  1765. confidence_data = []
  1767. for _, row in tqdm(sample_metadata.iterrows(), total=len(sample_metadata), desc="Confidence analysis"):
  1768. sample_name = row['sample']
  1769. cohort = row['Cohort']
  1770. expr_col = row['expression_col']
  1771. sample_profile = sample_expr_log[expr_col]
  1772. sample_profile = sample_profile[sample_profile > 0]
  1774. # MLMarker Model
  1775. if sample_name in sample_mapping_mlmarker:
  1776. ml_col = sample_mapping_mlmarker[sample_name]
  1777. scores = mlmarker_predictions[ml_col].sort_values(ascending=False)
  1778. margin = scores.iloc[0] - scores.iloc[1] if len(scores) > 1 else 0
  1779. is_correct = scores.index[0] in cancer_to_mlmarker_extended.get(cohort, [])
  1780. confidence_data.append({'Method': 'MLMarker', 'Margin': margin,
  1781. 'Correct': is_correct, 'Cohort': cohort})
  1783. # MLMarker Training Atlas (Spearman)
  1784. _, _, scores = atlas_correlation_classify(sample_profile, ml_atlas)
  1785. if scores:
  1786. sorted_scores = sorted(scores.values(), reverse=True)
  1787. margin = sorted_scores[0] - sorted_scores[1] if len(sorted_scores) > 1 else 0
  1788. pred = max(scores, key=scores.get)
  1789. is_correct = pred in cancer_to_mlmarker_extended.get(cohort, [])
  1790. confidence_data.append({'Method': 'MLMarker Training', 'Margin': margin,
  1791. 'Correct': is_correct, 'Cohort': cohort})
  1793. # HPA Atlas
  1794. _, _, scores = atlas_correlation_classify(sample_profile, hpa_matrix)
  1795. if scores:
  1796. sorted_scores = sorted(scores.values(), reverse=True)
  1797. margin = sorted_scores[0] - sorted_scores[1] if len(sorted_scores) > 1 else 0
  1798. pred = max(scores, key=scores.get)
  1799. is_correct = pred in cancer_to_hpa_extended.get(cohort, [])
  1800. confidence_data.append({'Method': 'HPA Atlas', 'Margin': margin,
  1801. 'Correct': is_correct, 'Cohort': cohort})
  1803. # PDB Atlas
  1804. _, _, scores = atlas_correlation_classify(sample_profile, pdb_matrix)
  1805. if scores:
  1806. sorted_scores = sorted(scores.values(), reverse=True)
  1807. margin = sorted_scores[0] - sorted_scores[1] if len(sorted_scores) > 1 else 0
  1808. pred = max(scores, key=scores.get)
  1809. is_correct = pred in cancer_to_pdb_extended.get(cohort, [])
  1810. confidence_data.append({'Method': 'PDB Atlas', 'Margin': margin,
  1811. 'Correct': is_correct, 'Cohort': cohort})
  1813. confidence_df = pd.DataFrame(confidence_data)
  1814. print(f"Confidence data collected: {len(confidence_df)} predictions")
  1816. # %%
  1817. # Figure 6: Confidence Margin Distribution
  1818. fig, axes = plt.subplots(1, 2, figsize=(14, 5))
  1820. # Plot 1: Box plot of margins by method
  1821. ax = axes[0]
  1822. order = ['MLMarker', 'MLMarker Training', 'HPA Atlas', 'PDB Atlas']
  1823. palette = {'MLMarker': '#E63946', 'MLMarker Training': '#457B9D',
  1824. 'HPA Atlas': '#2A9D8F', 'PDB Atlas': '#E9C46A'}
  1825. sns.boxplot(data=confidence_df, x='Method', y='Margin', order=order,
  1826. palette=palette, ax=ax)
  1827. ax.set_xlabel('Method', fontsize=12)
  1828. ax.set_ylabel('Confidence Margin (Top1 - Top2 Score)', fontsize=12)
  1829. ax.set_title('Prediction Confidence by Method', fontsize=14, fontweight='bold')
  1830. ax.tick_params(axis='x', rotation=45)
  1832. # Plot 2: Margin distribution for correct vs incorrect predictions
  1833. ax = axes[1]
  1834. correct_df = confidence_df[confidence_df['Correct']]
  1835. incorrect_df = confidence_df[~confidence_df['Correct']]
  1837. for i, method in enumerate(order):
  1838. correct_margins = correct_df[correct_df['Method'] == method]['Margin']
  1839. incorrect_margins = incorrect_df[incorrect_df['Method'] == method]['Margin']
  1841. positions = [i - 0.2, i + 0.2]
  1842. bp = ax.boxplot([correct_margins, incorrect_margins], positions=positions,
  1843. widths=0.35, patch_artist=True)
  1844. bp['boxes'][0].set_facecolor(palette[method])
  1845. bp['boxes'][0].set_alpha(0.8)
  1846. bp['boxes'][1].set_facecolor(palette[method])
  1847. bp['boxes'][1].set_alpha(0.3)
  1849. ax.set_xticks(range(len(order)))
  1850. ax.set_xticklabels(order, rotation=45, ha='right')
  1851. ax.set_xlabel('Method', fontsize=12)
  1852. ax.set_ylabel('Confidence Margin', fontsize=12)
  1853. ax.set_title('Confidence: Correct (dark) vs Incorrect (light)', fontsize=14, fontweight='bold')
  1855. # Add legend
  1856. from matplotlib.patches import Patch
  1857. legend_elements = [Patch(facecolor='gray', alpha=0.8, label='Correct'),
  1858. Patch(facecolor='gray', alpha=0.3, label='Incorrect')]
  1859. ax.legend(handles=legend_elements, loc='upper right')
  1861. plt.tight_layout()
  1862. plt.savefig('baseline_confidence_margins.png', dpi=300, bbox_inches='tight')
  1863. plt.savefig('baseline_confidence_margins.svg', bbox_inches='tight')
  1864. plt.show()
  1865. print("Saved: baseline_confidence_margins.png/svg")
  1867. # Print mean margins
  1868. print("\nMean Confidence Margins:")
  1869. print("=" * 50)
  1870. for method in order:
  1871. method_df = confidence_df[confidence_df['Method'] == method]
  1872. correct_margin = method_df[method_df['Correct']]['Margin'].mean()
  1873. incorrect_margin = method_df[~method_df['Correct']]['Margin'].mean()
  1874. print(f"{method:<20}: Correct={correct_margin:.3f}, Incorrect={incorrect_margin:.3f}")
  1876. # %% [markdown]
  1877. # ### 16.4 Statistical Significance: McNemar's Test
  1878. #
  1879. # Compare MLMarker vs each baseline using McNemar's test for paired nominal data.
  1881. # %%
  1882. from statsmodels.stats.contingency_tables import mcnemar
  1884. def run_mcnemar(correct_a, correct_b, name_a, name_b):
  1885. """Run McNemar's test comparing two classifiers."""
  1886. # Create contingency table
  1887. both_correct = sum(a and b for a, b in zip(correct_a, correct_b))
  1888. a_only = sum(a and not b for a, b in zip(correct_a, correct_b))
  1889. b_only = sum(not a and b for a, b in zip(correct_a, correct_b))
  1890. both_wrong = sum(not a and not b for a, b in zip(correct_a, correct_b))
  1892. table = [[both_correct, a_only], [b_only, both_wrong]]
  1893. result = mcnemar(table, exact=True)
  1895. return {
  1896. 'Comparison': f'{name_a} vs {name_b}',
  1897. 'Both Correct': both_correct,
  1898. f'{name_a} Only': a_only,
  1899. f'{name_b} Only': b_only,
  1900. 'Both Wrong': both_wrong,
  1901. 'p-value': result.pvalue,
  1902. 'Significant': 'Yes' if result.pvalue < 0.05 else 'No'
  1903. }
  1905. # Get correct/incorrect arrays
  1906. mlmarker_correct = results_df['mlmarker_top1_extended'].values
  1907. ml_training_correct = results_df['ml_training_top1_extended'].values
  1908. hpa_correct = results_df['hpa_top1_extended'].values
  1909. pdb_correct = results_df['pdb_top1_extended'].values
  1911. # Run McNemar's tests
  1912. mcnemar_results = []
  1913. mcnemar_results.append(run_mcnemar(mlmarker_correct, ml_training_correct, 'MLMarker', 'ML Training'))
  1914. mcnemar_results.append(run_mcnemar(mlmarker_correct, hpa_correct, 'MLMarker', 'HPA'))
  1915. mcnemar_results.append(run_mcnemar(mlmarker_correct, pdb_correct, 'MLMarker', 'PDB'))
  1916. mcnemar_results.append(run_mcnemar(ml_training_correct, hpa_correct, 'ML Training', 'HPA'))
  1917. mcnemar_results.append(run_mcnemar(ml_training_correct, pdb_correct, 'ML Training', 'PDB'))
  1918. mcnemar_results.append(run_mcnemar(hpa_correct, pdb_correct, 'HPA', 'PDB'))
  1920. mcnemar_df = pd.DataFrame(mcnemar_results)
  1921. print("McNemar's Test Results (Top-1 Extended Accuracy):")
  1922. print("=" * 80)
  1923. print(mcnemar_df.to_string(index=False))
  1925. # %% [markdown]
  1926. # ### 16.5 Random Baseline Comparison
  1927. #
  1928. # Add a random classifier to show all methods perform significantly above chance.
  1930. # %%
  1931. # Calculate random baseline (chance level)
  1932. # For MLMarker: 34 tissue classes
  1933. n_mlmarker_classes = len(mlmarker_predictions.index)
  1934. random_mlmarker_strict = 1 / n_mlmarker_classes * 100
  1935. # Extended: average number of acceptable tissues per cohort
  1936. avg_extended = np.mean([len(v) for v in cancer_to_mlmarker_extended.values()])
  1937. random_mlmarker_extended = avg_extended / n_mlmarker_classes * 100
  1939. # For HPA: 61 tissues
  1940. n_hpa_classes = len(hpa_matrix.columns)
  1941. avg_hpa_extended = np.mean([len(v) for v in cancer_to_hpa_extended.values()])
  1942. random_hpa_extended = avg_hpa_extended / n_hpa_classes * 100
  1944. # For PDB: 57 tissues
  1945. n_pdb_classes = len(pdb_matrix.columns)
  1946. avg_pdb_extended = np.mean([len(v) for v in cancer_to_pdb_extended.values()])
  1947. random_pdb_extended = avg_pdb_extended / n_pdb_classes * 100
  1949. print("Random Baseline (Chance Level):")
  1950. print("=" * 60)
  1951. print(f"MLMarker ({n_mlmarker_classes} classes):")
  1952. print(f" Strict: {random_mlmarker_strict:.2f}%")
  1953. print(f" Extended (avg {avg_extended:.1f} acceptable): {random_mlmarker_extended:.2f}%")
  1954. print(f"\nHPA ({n_hpa_classes} tissues):")
  1955. print(f" Extended (avg {avg_hpa_extended:.1f} acceptable): {random_hpa_extended:.2f}%")
  1956. print(f"\nPDB ({n_pdb_classes} tissues):")
  1957. print(f" Extended (avg {avg_pdb_extended:.1f} acceptable): {random_pdb_extended:.2f}%")
  1959. # %% [markdown]
  1960. # ## 17. Final Summary Figure: Combined Comparison
  1962. # %%
  1963. # Figure 7: Comprehensive comparison bar chart with random baseline
  1964. fig, ax = plt.subplots(figsize=(12, 7))
  1966. # Collect all method accuracies (Top-1 Extended)
  1967. method_accs = {
  1968. 'MLMarker\nModel': results_df['mlmarker_top1_extended'].mean() * 100,
  1969. 'MLMarker\nTraining Atlas': results_df['ml_training_top1_extended'].mean() * 100,
  1970. 'HPA Atlas\n(Spearman)': results_df['hpa_top1_extended'].mean() * 100,
  1971. 'PDB Atlas\n(Spearman)': results_df['pdb_top1_extended'].mean() * 100,
  1972. }
  1974. # Add random baselines
  1975. method_accs['Random\n(MLMarker)'] = random_mlmarker_extended
  1976. method_accs['Random\n(HPA)'] = random_hpa_extended
  1977. method_accs['Random\n(PDB)'] = random_pdb_extended
  1979. # Colors
  1980. bar_colors = ['#E63946', '#457B9D', '#2A9D8F', '#E9C46A', '#CCCCCC', '#CCCCCC', '#CCCCCC']
  1982. x = np.arange(len(method_accs))
  1983. bars = ax.bar(x, list(method_accs.values()), color=bar_colors, edgecolor='black', linewidth=0.5)
  1985. # Add value labels
  1986. for bar in bars:
  1987. height = bar.get_height()
  1988. ax.annotate(f'{height:.1f}%', xy=(bar.get_x() + bar.get_width()/2, height),
  1989. xytext=(0, 3), textcoords='offset points', ha='center', va='bottom',
  1990. fontsize=11, fontweight='bold')
  1992. ax.set_xticks(x)
  1993. ax.set_xticklabels(list(method_accs.keys()), fontsize=10)
  1994. ax.set_ylabel('Top-1 Extended Accuracy (%)', fontsize=12)
  1995. ax.set_title('Tissue Classification Accuracy: MLMarker vs Atlas Baselines',
  1996. fontsize=14, fontweight='bold')
  1997. ax.set_ylim(0, 100)
  1998. ax.axhline(y=50, color='gray', linestyle='--', alpha=0.5, label='50% threshold')
  1999. ax.grid(True, alpha=0.3, axis='y')
  2001. # Add annotation for improvement
  2002. best_baseline = max(method_accs['MLMarker\nTraining Atlas'],
  2003. method_accs['HPA Atlas\n(Spearman)'],
  2004. method_accs['PDB Atlas\n(Spearman)'])
  2005. improvement = method_accs['MLMarker\nModel'] - best_baseline
  2006. ax.annotate(f'+{improvement:.1f}% vs best baseline',
  2007. xy=(0, method_accs['MLMarker\nModel']),
  2008. xytext=(1.5, method_accs['MLMarker\nModel'] + 5),
  2009. fontsize=11, fontweight='bold', color='#E63946',
  2010. arrowprops=dict(arrowstyle='->', color='#E63946'))
  2012. plt.tight_layout()
  2013. plt.savefig('baseline_final_comparison.png', dpi=300, bbox_inches='tight')
  2014. plt.savefig('baseline_final_comparison.svg', bbox_inches='tight')
  2015. plt.show()
  2016. print("Saved: baseline_final_comparison.png/svg")
  2018. # %%
  2019. # Save all results
  2020. mcnemar_df.to_csv('baseline_mcnemar_tests.csv', index=False)
  2021. confidence_df.to_csv('baseline_confidence_margins.csv', index=False)
  2023. print("\n" + "=" * 80)
  2024. print("COMPLETE ANALYSIS SUMMARY")
  2025. print("=" * 80)
  2026. print(f"\nDataset: {len(results_df)} cancer samples across {len(eval_cohorts)} cohorts")
  2027. print(f"Cohorts: {', '.join(eval_cohorts)}")
  2028. print(f"\nMethods compared:")
  2029. print(f" 1. MLMarker Model (Random Forest, 34 tissues)")
  2030. print(f" 2. MLMarker Training Atlas (Spearman correlation)")
  2031. print(f" 3. HPA Atlas (Spearman correlation, {len(hpa_matrix.columns)} tissues)")
  2032. print(f" 4. PDB Atlas (Spearman correlation, {len(pdb_matrix.columns)} tissues)")
  2033. print(f" 5. Nearest Centroid (Euclidean distance)")
  2034. print(f" 6. Random baseline (chance level)")
  2036. print(f"\n{'Method':<25} {'Top-1 Extended':>15} {'vs Random':>12}")
  2037. print("-" * 55)
  2038. print(f"{'MLMarker Model':<25} {results_df['mlmarker_top1_extended'].mean()*100:>14.1f}% {'+' + str(round(results_df['mlmarker_top1_extended'].mean()*100 - random_mlmarker_extended, 1)):>11}%")
  2039. print(f"{'MLMarker Training':<25} {results_df['ml_training_top1_extended'].mean()*100:>14.1f}% {'+' + str(round(results_df['ml_training_top1_extended'].mean()*100 - random_mlmarker_extended, 1)):>11}%")
  2040. print(f"{'HPA Atlas':<25} {results_df['hpa_top1_extended'].mean()*100:>14.1f}% {'+' + str(round(results_df['hpa_top1_extended'].mean()*100 - random_hpa_extended, 1)):>11}%")
  2041. print(f"{'PDB Atlas':<25} {results_df['pdb_top1_extended'].mean()*100:>14.1f}% {'+' + str(round(results_df['pdb_top1_extended'].mean()*100 - random_pdb_extended, 1)):>11}%")
  2043. print(f"\nAll saved files:")
  2044. print(" - baseline_comparison_summary.csv")
  2045. print(" - baseline_comparison_cohort.csv")
  2046. print(" - baseline_comparison_detailed.csv")
  2047. print(" - baseline_mcnemar_tests.csv")
  2048. print(" - baseline_confidence_margins.csv")
  2049. print(" - baseline_topk_accuracy_lines.png/svg")
  2050. print(" - baseline_cohort_barplot.png/svg")
  2051. print(" - baseline_cohort_topk_lines.png/svg")
  2052. print(" - baseline_accuracy_heatmap.png/svg")
  2053. print(" - baseline_confusion_matrices.png/svg")
  2054. print(" - baseline_umap_centroids.png/svg")
  2055. print(" - baseline_confidence_margins.png/svg")
  2056. print(" - baseline_final_comparison.png/svg")
  2058. # %%

baseline_comparison_clean.ipynb at commit c164143, no license · at the source

Overview

Authors: Tine Claeys1,2, Sam van Puyenbroeck1,2, Kris Gevaert1,2, Lennart Martens1,2,3,4
ORCID iDs: Tine Claeys, Sam van Puyenbroeck, Kris Gevaert, Lennart Martens
  1. VIB-UGent Center for Medical Biotechnology, VIB,Ghent, Belgium
  2. Department of Biomolecular Medicine, Ghent University,Ghent, Belgium
  3. BioOrganic Mass Spectrometry Laboratory (LSMBO), UMR 7178, IPHC, University of Strasbourg, CNRS,Strasbourg, 67000 France
  4. Infrastructure Nationale de Protéomique, ProFI-UAR 2048, Strasbourg, France
Institutions: VIB-UGent Center for Medical Biotechnology (Belgium); Ghent University (Belgium); Université de Strasbourg (France)
Journal: Genome biology, volume 27, issue 1, article 207
Dates: received 9 July 2025; accepted 21 May 2026; published online 24 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1186/s13059-026-04125-8 · PMID 42343371 · PMCID PMC13292323 · OpenAlex W4411313251
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.06.10.658642
Status: code verified
Categories: genetics / omics (modality), human (organism), methods / tools (subfield)
Methods: Statistics, Machine learning
Keywords: Tissue prediction, Mass spectrometry-based proteomics, Machine learning, Public data reuse, AI
MeSH: Machine Learning*, Proteomics*, Software*, Biomarkers, Humans, Prediction Algorithms, Predictive Learning Models, Random Forest (* major topic)
Journal subjects: Methodology
Topic: AI in cancer detection (Artificial Intelligence, Computer Science), according to OpenAlex
Funding: Agence Nationale de la Recherche (G0GDV23N)
Citations: cited by 1 paper (Europe PMC); 47 references in the paper

Abstract

MLMarker is a machine learning tool that computes continuous tissue similarity scores for proteomics data, addressing the challenge of interpreting complex or sparse datasets. Trained on 34 healthy tissues, its Random Forest model generates probabilistic predictions with SHAP-based protein-level explanations. A penalty factor corrects for missing proteins, improving robustness for low-coverage samples. Across three public datasets, MLMarker revealed brain-like signatures in cerebral melanoma metastases, achieved high accuracy in a pan-cancer cohort, and identified brain and pituitary origins in biofluids. MLMarker provides an interpretable framework for tissue inference and hypothesis generation, available as a Python package and Streamlit app.

Supplementary Information: The online version contains supplementary material available at 10.1186/s13059-026-04125-8.

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

Repositories

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

TineClaeys/MLMarker-manuscript

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c164143cc2859301d5dff8ab50767401fcf72ba8, 13 July 2026
Languages: Jupyter (16), Python (6), R (2)
Size: 152 files, 24 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, 16 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (22 files), Matplotlib (20 files), pandas (20 files), seaborn (19 files), SciPy (10 files), Plotly (8 files), UMAP (7 files), scikit-learn (5 files), limma (2 files), tidyverse (2 files), statsmodels (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/TineClaeys/MLMarker-manuscript

TineClaeys/MLMarker

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 6444accee481886a09fd04c9b873ae1185c23420, 30 January 2026
Languages: Python (14), Jupyter (2)
Size: 72 files, 16 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, environment (pyproject.toml, requirements.txt, setup.py), tests, continuous integration, 2 notebooks
Not found: CITATION.cff, documentation
Tools: pandas (14 files), NumPy (10 files), SHAP (6 files), Matplotlib (4 files), Plotly (4 files), scikit-learn (4 files), XGBoost (4 files), imbalanced-learn (2 files), seaborn (2 files), PyTorch (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
18 files
At the source: github.com/TineClaeys/MLMarker

TineClaeys/MLMarker-streamlit

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 2bf8a28f94efe81f9d89ddb6223d255b1e3af9e0, 21 August 2026
Languages: Python (8), Jupyter (2)
Size: 50 files, 10 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, environment (requirements.txt), 2 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (8 files), pandas (8 files), Plotly (8 files), scikit-learn (4 files), SciPy (3 files), Matplotlib (2 files), seaborn (2 files), SHAP (2 files), Pillow (1 file), UMAP (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
12 files
At the source: github.com/TineClaeys/MLMarker-streamlit

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

Tracing map

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

What the map holds:

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

All datasets used in this study are publicly available via ProteomeXchange (PXD007592, PXD009021, PXD008029) and MassIVE (MSV000095036). MLMarker predictions, differential expression results, and SHAP-derived protein lists are available in the Supplementary Data & Github repository for the manuscript analyses, pip package and streamlit application respectively github.com/TineClaeys/MLMarker-manuscript; github.com/TineClaeys/MLMarker; github.com/TineClaeys/MLMarker-streamlit. The streamlit application is available through mlmarker.streamlit.app all falling under Apache-2.0 license.

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 2, 28 September 2026

  • Authors: added Tine Claeys (0000-0001-9408-488X); Sam van Puyenbroeck (0009-0001-5281-6891); Kris Gevaert (0000-0002-4237-0283); Lennart Martens (0000-0003-4277-658X); removed Tine Claeys; Sam van Puyenbroeck; Kris Gevaert; Lennart Martens
  • Funding: added CHIST-ERA: G0GDV23N

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 5 keywords, 8 MeSH terms, 41 references.

Cite

This paper

Claeys, T., van Puyenbroeck, S., Gevaert, K., & Martens, L. (2026). MLMarker: a machine learning framework for tissue inference and biomarker discovery. Genome biology, 27(1), 207. https://doi.org/10.1186/s13059-026-04125-8

BibTeX

@article{claeys2026mlmarker,
author = {Claeys, Tine and van Puyenbroeck, Sam and Gevaert, Kris and Martens, Lennart},
title = {{MLMarker: a machine learning framework for tissue inference and biomarker discovery}},
journal = {Genome biology},
year = {2026},
month = jun,
volume = {27},
number = {1},
pages = {207},
publisher = {BMC},
issn = {1474-7596},
doi = {10.1186/s13059-026-04125-8},
url = {https://doi.org/10.1186/s13059-026-04125-8},
pmid = {42343371},
pmcid = {PMC13292323}
}

RIS

TY - JOUR
AU - Claeys, Tine
AU - van Puyenbroeck, Sam
AU - Gevaert, Kris
AU - Martens, Lennart
TI - MLMarker: a machine learning framework for tissue inference and biomarker discovery
T2 - Genome biology
J2 - Genome Biol
PY - 2026
DA - 2026/06/24
VL - 27
IS - 1
SP - 207
SN - 1474-7596
PB - BMC
DO - 10.1186/s13059-026-04125-8
UR - https://doi.org/10.1186/s13059-026-04125-8
LA - en
ER -

CSL-JSON

{
"id": "10.1186/s13059-026-04125-8",
"type": "article-journal",
"title": "MLMarker: a machine learning framework for tissue inference and biomarker discovery",
"container-title": "Genome biology",
"author": [
{
"family": "Claeys",
"given": "Tine"
},
{
"family": "van Puyenbroeck",
"given": "Sam"
},
{
"family": "Gevaert",
"given": "Kris"
},
{
"family": "Martens",
"given": "Lennart"
}
],
"container-title-short": "Genome Biol",
"volume": "27",
"issue": "1",
"page": "207",
"DOI": "10.1186/s13059-026-04125-8",
"PMID": "42343371",
"PMCID": "PMC13292323",
"ISSN": "1474-7596",
"publisher": "BMC",
"URL": "https://doi.org/10.1186/s13059-026-04125-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
24
]
]
}
}

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/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: imbalanced-learn, SHAP, XGBoost, 13 other tools, 1 reference
[2] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: imbalanced-learn, SHAP, XGBoost, 11 other tools, genetics / omics
[3] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: imbalanced-learn, XGBoost, UMAP, 8 other tools, methods / tools
[4] doi:10.1038/s41467-026-76939-w [code]
HIPPIE: a generative model for electrophysiological analysis across species, technologies, and modalities.
Journal: Nature communications
In common: imbalanced-learn, UMAP, Pillow, 9 other tools, methods / tools
[5] doi:10.1016/j.isci.2026.115329 [code]
Brain metastases converge on shared geometric architecture and transcriptomic landscape yet remain distinct from gliomas.
Journal: iScience
In common: imbalanced-learn, SHAP, XGBoost, 8 other tools, genetics / omics
[6] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: imbalanced-learn, SHAP, XGBoost, 8 other tools
[7] doi:10.1038/s43856-026-01606-6 [code]
Validation of remote multimodal AI screening for Parkinson disease across diverse settings.
Journal: Communications medicine
In common: imbalanced-learn, SHAP, Plotly, 8 other tools
[8] doi:10.1126/sciadv.aed3650 [code]
Truthful visualizations for mass spectrometry imaging enable high-spatial-resolution interactive &lt;i&gt;m/z&lt;/i&gt; mapping and exploration.
Journal: Science advances
In common: SHAP, XGBoost, UMAP, 7 other tools, methods / tools, genetics / omics
[9] 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: limma, UMAP, Plotly, 8 other tools, genetics / omics
[10] doi:10.3390/s26175327 [code]
Subject Identity Confounds qEEG Emotion Recognition on DEAP and DREAMER.
Journal: Sensors (Basel, Switzerland)
In common: imbalanced-learn, SHAP, XGBoost, 7 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.