OSCR

Wnt activation and dual SMAD inhibition for induction and maintenance of hindbrain-like neural stem cell from hiPSCs.

Code ↔ Paper

3 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 3 matches
  1. [1] § Results › scRNA-seq of early- and late-passage number Hb-LiNSCs reveals distinct cell populations ↔ scRNAseq_HbL-iNSCs.ipynb, lines 702 to 730 · score 0.73 · stem cells, HOXB8, HES5, HOXB2, NES, HOXB7 · Report a mismatch
  2. [2] § Results › Induction of Hb-LiNSCs ↔ scRNAseq_HbL-iNSCs.ipynb, lines 427 to 470 · score 0.69 · cellular component, biological process, upregulated genes, Augmented, GO, PanglaoDB · Report a mismatch
  3. [3] § Results › Transcriptional differences between early- and late-passage Hb-LiNSCs ↔ scRNAseq_HbL-iNSCs.ipynb, lines 427 to 470 · score 0.60 · GO biological process, upregulated genes, Augmented, PanglaoDB, human, Hb · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,188 lines · 37 KB · CC-BY-4.0 · 3 matches

  1. # %%
  2. import scanpy as sc #important note: the source code for scanpy heatmap was modified locally to suite the visual preferences (Ziadoon Al-Akashi)
  3. import anndata as ad
  4. import numpy as np
  5. import pandas as pd
  6. import matplotlib.pyplot as plt
  7. import seaborn as sns
  8. from sklearn.preprocessing import StandardScaler
  9. import os
  10. from matplotlib.colors import Normalize
  11. import gseapy as gp
  12. from typing import Union, Optional, Tuple, Mapping
  13. from matplotlib.colors import Normalize
  14. import matplotlib.gridspec as gridspec
  15. from matplotlib.patches import Patch
  16. from matplotlib.lines import Line2D
  19. # %% [markdown]
  20. # # Reading and processing
  22. # %%
  23. data_folder = 'E:/scRNA seq Data/'
  25. # %%
  26. ribo_url = "http://software.broadinstitute.org/gsea/msigdb/download_geneset.jsp?geneSetName=KEGG_RIBOSOME&fileType=txt"
  27. ribo_genes = pd.read_table(ribo_url, skiprows=2, header = None, names=['gene'])
  28. ribo_genes
  30. # %%
  31. all_adata = []
  33. h5_paths = [
  34. data_folder + 'ACL_iNSCs_PN5',
  35. data_folder + 'ACL_iNSCs_PN21',
  36. data_folder + 'ACL_iNSCs_PN41',
  37. data_folder + 'ACL_iNSCs_PN59',
  38. ]
  40. # %%
  41. def read_and_append_data(path, all_adata):
  42. adata = sc.read_10x_mtx(path, var_names="gene_symbols", cache=True)
  43. sample_name = path.split('_')[2]
  44. adata.obs['Sample'] = sample_name
  45. all_adata.append(adata)
  47. def process_data(all_adata, ribo_genes):
  48. for i, adata in enumerate(all_adata):
  49. print(f"Before filtering - Number of cells: {adata.n_obs}")
  51. sc.pp.filter_cells(adata, min_genes=400)
  52. adata.var['mt'] = adata.var_names.str.startswith('MT-')
  53. adata.var['ribo'] = adata.var_names.isin(ribo_genes['gene'])
  54. sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo'], percent_top=None, log1p=False, inplace=True)
  56. upper_lim = np.quantile(adata.obs.n_genes_by_counts.values, 0.98)
  57. adata = adata[adata.obs.n_genes_by_counts < upper_lim]
  59. adata = adata[adata.obs.pct_counts_mt < 10, :].copy()
  60. all_adata[i] = adata
  62. print(f"After filtering - Number of cells: {adata.n_obs}")
  64. return all_adata
  66. # %%
  67. for path in h5_paths:
  68. read_and_append_data(path, all_adata)
  70. # %%
  71. processed_data = process_data(all_adata, ribo_genes)
  73. # %%
  74. all_adata[1]
  76. # %%
  77. concatenated_processed_data = ad.concat(all_adata)
  79. # %%
  80. concatenated_processed_data.obs_names_make_unique()
  82. # %%
  83. concatenated_processed_data.obs.sort_values('n_genes_by_counts')
  85. # %%
  86. sc.pl.violin(concatenated_processed_data,
  87. ['n_genes_by_counts', 'total_counts', 'pct_counts_mt', 'pct_counts_ribo'],
  88. jitter=0.4,
  89. multi_panel=True)
  91. # %%
  92. import scrublet as scr
  94. # %%
  95. adata_for_scrublet = concatenated_processed_data.copy()
  97. # %%
  98. scrub = scr.Scrublet(adata_for_scrublet.X)
  100. # %%
  101. doublet_scores, predicted_doublets = scrub.scrub_doublets()
  102. plt.hist(doublet_scores, bins=50, color='blue', edgecolor='black')
  103. plt.xlabel('Doublet Scores')
  104. plt.ylabel('Frequency')
  105. plt.title('Distribution of Doublet Scores')
  106. plt.show()
  108. # %%
  109. concatenated_processed_data.obs['doublet_scores'] = doublet_scores
  111. # %%
  112. concatenated_processed_data = concatenated_processed_data[~predicted_doublets, :]
  114. # %%
  115. adata = concatenated_processed_data.copy()
  117. # %% [markdown]
  118. # # Normalize and highy variable genes
  120. # %%
  121. sc.pp.normalize_total(adata, target_sum=1e4)
  122. sc.pp.log1p(adata)
  123. sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=2, min_disp=0.5)
  124. adata.raw = adata
  125. adata_raw = sc.AnnData(X=adata.raw.X, obs=adata.obs, var=adata.var, uns=adata.uns)
  126. adata = adata[:, adata.var.highly_variable]
  128. # %%
  129. sc.set_figure_params(transparent=True, fontsize=12)
  130. custom_palette = sns.color_palette("hls", 4)
  131. sns.set_palette('colorblind')
  133. # %%
  134. sc.pp.scale(adata, max_value=10)
  135. sc.tl.pca(adata, svd_solver='arpack')
  136. sc.pl.pca_variance_ratio(adata, log=False, n_pcs=50, save=True)
  138. # %%
  139. sc.pl.pca(adata, color='Sample', save='final_PCA_.png')
  141. # %%
  142. adata.obs
  144. # %%
  145. sc.pp.neighbors(adata, n_neighbors=16, n_pcs=20)
  147. # Compute UMAP
  148. sc.tl.umap(adata)
  149. sc.tl.leiden(adata, resolution=0.8)
  151. # %%
  152. file_path = "final_processed_adata.h5ad"
  154. # %%
  155. #adata.write(file_path)
  157. # %%
  158. adata = sc.read(file_path)
  160. # %%
  161. color_s = ['#800000', '#9A6324', '#808000', '#469990', '#000075', '#e6194B', '#f58231', '#ffe119', '#bfef45', '#3cb44b', '#42d4f4', '#4363d8', '#911eb4', '#f032e6', '#a9a9a9', '#fabed4', '#ffd8b1']
  162. legend_handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=10, label=f'{i}') for i, color in enumerate(color_s)]
  164. fig, ax = plt.subplots(figsize=(1, 1))
  165. ax.axis('off')
  166. ax.legend(handles=legend_handles, ncol=3, loc='center')
  167. plt.savefig('legend_clusters.png', bbox_inches='tight', transparent=True, dpi=300)
  168. plt.show()
  170. # %%
  171. # UMAP Leiden
  172. color_s = ['#800000', '#9A6324', '#808000', '#469990', '#000075','#e6194B', '#f58231', '#ffe119', '#bfef45', '#3cb44b',
  173. '#42d4f4', '#4363d8', '#911eb4', '#f032e6', '#a9a9a9','#fabed4', '#ffd8b1', '#aaffc3', '#dcbeff']
  174. sc.set_figure_params(transparent=True, dpi=80, fontsize=14, dpi_save=300, figsize=(10, 8))
  175. sns.set_palette("deep")
  176. plt.figure(figsize=(10, 8))
  178. sc.pl.umap(adata,
  179. color=['leiden'],
  180. palette=color_s,
  181. use_raw=False,
  182. title=[''],
  183. legend_loc=False,
  184. legend_fontsize=8,
  185. legend_fontoutline=2,
  186. #frameon=False,
  187. size=30,
  188. alpha=1,
  189. linewidth=0.5,
  190. show=False
  191. )
  193. plt.title('Leiden Clusters', fontsize=20)
  194. plt.xlabel('UMAP1', fontsize=20)
  195. plt.ylabel('UMAP2', fontsize=20)
  197. plt.gca().spines['top'].set_visible(False)
  198. plt.gca().spines['right'].set_visible(False)
  200. plt.savefig('final_umap_leiden_res8.png', bbox_inches='tight', dpi=300)
  201. plt.show()
  204. # %%
  205. color_s1 = {'PN5': 'darkred', 'PN21': 'purple', 'PN41': 'darkorange', 'PN59': 'gold'}
  207. # %%
  208. color_s1 = {'PN5': 'darkred', 'PN21': 'purple', 'PN41': 'darkorange', 'PN59': 'gold'}
  209. legend_handles1 = [Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=10, label=label) for label, color in color_s1.items()]
  211. fig, ax = plt.subplots(figsize=(1, 1))
  212. ax.axis('off')
  213. ax.legend(handles=legend_handles1, ncol=1, loc='center')
  215. plt.savefig('legend_clusters_sample.png', bbox_inches='tight', transparent=True, dpi=300)
  216. plt.show()
  218. # %%
  219. # UMAP Samples
  221. sc.set_figure_params(transparent=True, dpi=80, fontsize=14, dpi_save=300, figsize=(10, 8))
  222. sns.set_palette("deep")
  224. plt.figure(figsize=(10, 8))
  225. sc.pl.umap(adata,
  226. color='Sample',
  227. palette=color_s1,
  228. use_raw=False,
  229. title=[''],
  230. legend_loc=False,
  231. size=30,
  232. alpha=0.8,
  233. linewidth=0.5,
  234. show=False
  235. )
  237. #plt.title('Samples', fontsize=20)
  238. plt.xlabel('UMAP1', fontsize=20)
  239. plt.ylabel('UMAP2', fontsize=20)
  240. plt.gca().spines['top'].set_visible(False)
  241. plt.gca().spines['right'].set_visible(False)
  243. plt.savefig('final_umap_sample_res3.png', bbox_inches='tight', dpi=300)
  244. plt.show()
  247. # %%
  248. # Sub UMAP Samples
  250. sc.set_figure_params(transparent=True, dpi=80, fontsize=14, dpi_save=300)
  251. sns.set_palette("deep")
  252. fig, axs = plt.subplots(4, 1, figsize=(2, 5), sharex=True, sharey=True)
  253. sample_groups = ['PN5', 'PN21', 'PN41', 'PN59']
  255. for i, sample_group in enumerate(sample_groups):
  256. ax = axs[i]
  257. sc.pl.umap(adata, color=['Sample'],
  258. groups=[sample_group],
  259. title=[''],
  260. size=2,
  261. alpha=0.8,
  262. linewidth=0.1,
  263. legend_loc=None,
  264. show=False,
  265. ax=ax)
  266. ax.set_xlabel('')
  267. ax.set_ylabel('')
  268. ax.spines['top'].set_visible(True)
  269. ax.spines['right'].set_visible(True)
  270. ax.spines['bottom'].set_visible(True)
  271. ax.spines['left'].set_visible(True)
  272. ax.set_xticks([])
  273. ax.set_yticks([])
  274. ax.set_aspect('equal')
  276. plt.subplots_adjust(hspace=0.05)
  277. plt.savefig('combined_umap_plots_column_with_boxes.png', bbox_inches='tight', dpi=300)
  278. plt.show()
  281. # %%
  282. sample_counts = adata.obs['Sample'].value_counts()
  283. print(sample_counts)
  285. # %%
  286. sample_counts = adata.obs['Sample'].value_counts()
  288. sample_counts = sample_counts.sort_index(ascending=False)
  289. plt.figure(figsize=(4, 4))
  290. ax = plt.axes()
  291. sample_counts.plot(kind='barh', color=[color_s1.get(sample, 'slategray') for sample in sample_counts.index], ax=ax)
  292. ax.grid(False)
  293. plt.xlabel('Number of cells')
  294. alpha=0.8
  296. ax.spines['top'].set_visible(False)
  297. ax.spines['right'].set_visible(False)
  298. #ax.spines['bottom'].set_visible(False)
  299. #ax.spines['left'].set_visible(False)
  301. for i, value in enumerate(sample_counts):
  302. plt.text(value + 0.05, i, str(value), ha='left', va='center')
  304. plt.xlim(0, 5500)
  306. plt.savefig('final_number_of_cells_in_sample.png', bbox_inches='tight', dpi=300, transparent=True)
  307. plt.show()
  310. # %%
  311. genes_to_plot = ['HOXB9', 'HOXB4', 'HOXB2', 'LMX1B', 'TH', 'NEFL', 'SOX1', 'NES', 'SOX2', 'HES5']
  313. num_genes = len(genes_to_plot)
  314. num_cols = min(5, num_genes)
  315. num_rows = -(-num_genes // num_cols)
  316. fig, axs = plt.subplots(num_rows, num_cols, figsize=(3 * num_cols, 3 * num_rows), sharex=True, sharey=True)
  317. sc.set_figure_params(transparent=True, dpi=300, fontsize=20)
  319. for gene, ax in zip(genes_to_plot, axs.flatten()):
  320. sc.pl.umap(adata,
  321. color=gene,
  322. cmap='rocket_r',
  323. frameon=False,
  324. vmin=0,
  325. size=5,
  326. use_raw=True,
  327. alpha=0.8,
  328. linewidth=0.8,
  329. legend_fontsize=1,
  330. show=False,
  331. ax=ax
  332. )
  333. ax.set_aspect('equal')
  334. for ax in axs.flatten()[num_genes:]:
  335. ax.remove()
  337. #plt.savefig('final_genes_on_clusters.png', bbox_inches='tight', dpi=300)
  338. plt.show()
  341. # %%
  342. sc.tl.rank_genes_groups(adata, groupby='Sample',reference='PN59', method='wilcoxon', use_raw=True)
  344. # %%
  345. sc.set_figure_params(dpi=300)
  346. sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, show=False, ncols=3)
  347. plt.savefig('final_wilcoxon ranked genes vs. PN59.png', bbox_inches='tight', dpi=300, transparent=True)
  349. # %%
  350. top_n_genes = 25
  351. marker_genes_df = pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(top_n_genes)
  352. file_path = 'final_wilcoxon_sample_25.csv'
  354. marker_genes_df.to_csv(file_path, index=False)
  355. marker_genes_df
  357. # %%
  358. #subset = adata[adata.obs['Sample'].isin(['PN5', 'PN59'])].copy()
  359. #sc.tl.rank_genes_groups(subset, groupby='Sample', method='t-test', use_raw=True)
  361. # %%
  362. dedf = sc.get.rank_genes_groups_df(adata, group='PN5', key='rank_genes_groups')
  363. upreg_df = dedf[(dedf.pvals < 0.05) & (dedf.logfoldchanges > 1.5)]
  364. downreg_df = dedf[(dedf.pvals < 0.05) & (dedf.logfoldchanges < -1.5)]
  366. upreg_df = upreg_df.sort_values('logfoldchanges', ascending=False).reset_index(drop=True)
  367. downreg_df = downreg_df.sort_values('logfoldchanges', ascending=True).reset_index(drop=True)
  369. upreg_genes = upreg_df.names.tolist()
  370. downreg_genes = downreg_df.names.tolist()
  372. upreg_genes10 = upreg_df[:25].names.tolist()
  373. downreg_genes10 = downreg_df[:25].names.tolist()
  375. # %%
  376. fig, ax1 = plt.subplots(figsize=(12, 10))
  378. ax2 = ax1.twinx()
  379. ax1.barh(range(len(upreg_genes10)), upreg_df[:25]['logfoldchanges'], color='darkgreen', label='Upregulated', alpha=0.7)
  380. ax2.barh(range(len(downreg_genes10)), downreg_df[:25]['logfoldchanges'], color='darkred', label='Downregulated', alpha=0.7)
  382. ax1.axvline(x=0, color='black', linestyle='-')
  383. ax1.set_yticks(range(len(downreg_genes10)))
  384. ax1.set_yticklabels(downreg_genes10)
  385. ax2.set_yticks(range(len(upreg_genes10)))
  386. ax2.set_yticklabels(upreg_genes10)
  388. ax1.set_xlabel('Log-Fold Change')
  389. ax1.set_xlim([-30, 30])
  391. lines, labels = ax1.get_legend_handles_labels()
  392. lines2, labels2 = ax2.get_legend_handles_labels()
  393. ax1.legend(lines + lines2, labels + labels2, loc='upper right', fontsize=14)
  395. ax1.grid(False)
  396. ax2.grid(False)
  398. ax1.spines['top'].set_visible(False)
  399. ax1.spines['right'].set_visible(False)
  400. ax1.spines['left'].set_visible(False)
  401. ax2.spines['top'].set_visible(False)
  402. ax2.spines['right'].set_visible(False)
  403. ax2.spines['left'].set_visible(False)
  405. #plt.suptitle('', fontsize=20, y=0.95)
  406. plt.title('PN5 vs. PN59 (Top 25 DEGs)', fontsize=20)
  407. plt.tight_layout()
  408. plt.savefig('final_top25_DEGs.png', bbox_inches='tight', dpi=300, transparent=True)
  410. plt.show()
  413. # %%
  414. from adjustText import adjust_text
  416. # %%
  417. #subset_genes = pd.DataFrame(subset.uns['rank_genes_groups']['names'])
  418. #subset_genes
  420. # %%
  421. gp.get_library_name()
  423. # %%
  424. upreg_df_new = pd.DataFrame({'SigUpreg': upreg_genes})
  425. downreg_df_new = pd.DataFrame({'SigDownreg': downreg_genes})
  427. # %%
  428. num_clusters = upreg_df_new.shape[1]
  429. background_gene_list = adata.raw.var_names.tolist()
  430. gene_sets = ['PanglaoDB_Augmented_2021', 'GO_Biological_Process_2023', 'GO_Cellular_Component_2023' ]
  432. for cluster_number, cluster_name in enumerate(upreg_df_new.columns):
  433. if cluster_number >= num_clusters:
  434. break
  436. fig, axes = plt.subplots(1, 3, figsize=(40, 10))
  437. plt.subplots_adjust(wspace=1.5, hspace=2)
  438. plt.suptitle('PN5 vs. PN59 upregulated genes ($\it{p}$ < 0.05)', fontsize=36, fontweight='bold', y=1.02)
  440. for gene_set, ax in zip(gene_sets, axes.flatten()):
  441. gene_list = upreg_df_new[cluster_name].dropna().tolist()
  443. if not gene_list:
  444. continue
  446. enr_results = gp.enrichr(gene_list=gene_list,
  447. gene_sets=[gene_set],
  448. organism='human',
  449. outdir=None,
  450. background=background_gene_list)
  452. enr_results_sorted = enr_results.results.sort_values(by='P-value')
  454. top_terms = enr_results_sorted.head(10)['Term'].apply(lambda x: x[:50] + '...' if len(x) > 35 else x)
  455. top_p_values = enr_results_sorted.head(10)['P-value']
  457. ax.barh(top_terms, -1 * np.log10(top_p_values), color='slategray')
  458. ax.set_xlabel('-log10(p-value)', fontsize=18)
  459. ax.set_title(gene_set, fontsize=26)
  460. ax.tick_params(axis='y', labelsize=26)
  461. ax.xaxis.grid(False)
  462. ax.yaxis.grid(False)
  463. ax.spines['top'].set_visible(False)
  464. ax.spines['right'].set_visible(False)
  465. ax.set_aspect('auto')
  466. ax.set_ylim([-0.5, 10])
  468. plt.tight_layout()
  469. plt.savefig(f'final_{cluster_name}_analysis_PN5vsPN59.png', bbox_inches='tight', dpi=300, transparent=True)
  470. plt.show()
  473. # %%
  474. num_clusters = downreg_df_new.shape[1]
  475. background_gene_list = adata.raw.var_names.tolist()
  476. gene_sets = ['PanglaoDB_Augmented_2021', 'GO_Biological_Process_2023', 'GO_Cellular_Component_2023' ]
  478. for cluster_number, cluster_name in enumerate(downreg_df_new.columns):
  479. if cluster_number >= num_clusters:
  480. break
  482. fig, axes = plt.subplots(1, 3, figsize=(40, 10))
  483. plt.subplots_adjust(wspace=1.5, hspace=2)
  484. plt.suptitle('PN5 vs. PN59 downregulated genes ($\it{p}$ < 0.05)', fontsize=36, fontweight='bold', y=1.02)
  486. for gene_set, ax in zip(gene_sets, axes.flatten()):
  487. gene_list = downreg_df_new[cluster_name].dropna().tolist()
  489. if not gene_list:
  490. continue
  492. enr_results = gp.enrichr(gene_list=gene_list,
  493. gene_sets=[gene_set],
  494. organism='human',
  495. outdir=None,
  496. background=background_gene_list)
  498. enr_results_sorted = enr_results.results.sort_values(by='P-value')
  500. top_terms = enr_results_sorted.head(10)['Term'].apply(lambda x: x[:50] + '...' if len(x) > 35 else x)
  501. top_p_values = enr_results_sorted.head(10)['P-value']
  503. ax.barh(top_terms, -1 * np.log10(top_p_values), color='slategray')
  504. ax.set_xlabel('-log10(p-value)', fontsize=18)
  505. ax.set_title(gene_set, fontsize=26)
  506. ax.tick_params(axis='y', labelsize=26)
  507. ax.xaxis.grid(False)
  508. ax.yaxis.grid(False)
  509. ax.spines['top'].set_visible(False)
  510. ax.spines['right'].set_visible(False)
  511. ax.set_aspect('auto')
  512. ax.set_ylim([-0.5, 10])
  514. plt.tight_layout()
  515. plt.savefig(f'final_{cluster_name}_analysis_PN5vsPN59.png', bbox_inches='tight', dpi=300, transparent=True)
  516. plt.show()
  519. # %%
  520. #selected_genes = subset.raw[:, genes_to_show].X.toarray()
  521. #sample_labels = subset.obs['Sample'].tolist()
  522. #df = pd.DataFrame(selected_genes, columns=genes_to_show)
  523. #df['Sample'] = sample_labels
  525. # %%
  526. #df.set_index('Sample', inplace=True)
  527. #df_zscore = df.apply(zscore)
  528. #gene_list
  530. # %%
  531. color_s2 = {'Hindbrain - Differentiated':'#469990',
  532. 'Hindbrain - SCs':'#9A6324',
  533. 'Hindbrain - Serotonergic':'#800000',
  534. 'Midbrain - Dopaminergic':'#808000',
  535. 'Motor neurons':'#000075',
  536. 'Neural - SCs':'#e6194B',
  537. 'Serotonergic':'#f58231',
  538. 'Spinal - GABAergic':'#ffe119',
  539. 'Spinal - SCs':'#bfef45',
  540. 'Spinal - Schwann SCs':'#3cb44b',
  541. }
  543. legend_handles3 = [Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=10, label=label) for label, color in color_s2.items()]
  544. fig, ax = plt.subplots(figsize=(1, 1))
  545. ax.axis('off')
  546. ax.legend(handles=legend_handles3, ncol=1, loc='center', frameon=False)
  547. plt.savefig('legend_celltype.png', bbox_inches='tight', transparent=True, dpi=300)
  548. plt.show()
  550. # %%
  551. color_s2 = ['#469990',
  552. '#9A6324',
  553. '#800000',
  554. '#808000',
  555. '#000075',
  556. '#e6194B',
  557. '#f58231',
  558. '#ffe119',
  559. '#bfef45',
  560. '#3cb44b',
  561. ]
  563. # %%
  565. cluster_names = {
  566. '0': 'Hindbrain - SCs',
  567. '1': 'Hindbrain - SCs',
  568. '2': 'Hindbrain - SCs',
  569. '3': 'Hindbrain - SCs',
  570. '4': 'Hindbrain - SCs',
  571. '5': 'Motor neurons',
  572. '6': 'Spinal - SCs',
  573. '7': 'Hindbrain - SCs',
  574. '8': 'Midbrain - Dopaminergic',
  575. '9': 'Spinal - SCs',
  576. '10': 'Spinal - SCs',
  577. '11': 'Hindbrain - Serotonergic',
  578. '12': 'Spinal - Schwann SCs',
  579. '13': 'Spinal - GABAergic',
  580. '14': 'Neural - SCs',
  581. '15': 'Serotonergic',
  582. '16': 'Hindbrain - Differentiated',
  584. }
  586. adata.obs['cell_type'] = adata.obs['leiden'].map(cluster_names)
  588. # UMAP cell_types
  590. sc.set_figure_params(transparent=True, dpi=80, fontsize=14, dpi_save=300, figsize=(10, 8))
  591. sns.set_palette("deep")
  593. plt.figure(figsize=(10, 8))
  594. sc.pl.umap(adata,
  595. color='cell_type',
  596. palette=color_s2,
  597. use_raw=False,
  598. title=[''],
  599. frameon=False,
  600. legend_loc=False,
  601. legend_fontsize=16,
  602. #legend_fontoutline=1,
  603. size=30,
  604. alpha=0.8,
  605. linewidth=0.5,
  606. show=False
  607. )
  609. #plt.title('Samples', fontsize=20)
  610. #plt.xlabel('UMAP1', fontsize=20)
  611. #plt.ylabel('UMAP2', fontsize=20)
  612. #plt.gca().spines['top'].set_visible(False)
  613. #plt.gca().spines['right'].set_visible(False)
  615. ax.set_aspect('equal')
  616. plt.savefig('final__celltypes_NSCs.png', bbox_inches='tight', dpi=300)
  617. plt.show()
  620. # %%
  621. cell_type_counts = adata.obs['cell_type'].value_counts()
  622. print(cell_type_counts)
  624. # %%
  625. cell_type_counts = adata.obs['cell_type'].value_counts()
  627. plt.figure(figsize=(8, 6))
  628. ax = plt.axes()
  629. bars = cell_type_counts.plot(kind='barh', color=color_t, ax=ax)
  630. ax.grid(False)
  632. ax.spines['top'].set_visible(False)
  633. ax.spines['right'].set_visible(False)
  634. #ax.spines['bottom'].set_visible(False)
  635. #ax.spines['left'].set_visible(False)
  636. plt.xlabel('Number of Cells')
  638. for i, value in enumerate(cell_type_counts):
  639. plt.text(value + 0.05, i, str(value), ha='left', va='center')
  640. plt.savefig('final_number_of_cells_in_celltype.png', bbox_inches='tight', dpi=300, transparent=True)
  641. plt.show()
  644. # %%
  645. sc.tl.rank_genes_groups(adata, groupby='Sample', method='wilcoxon', use_raw=True)
  647. # %%
  648. #sns.set_palette("rocket")
  649. fig, ax = plt.subplots(figsize=(0.2, 2)) # Adjust the figsize as needed for your design
  650. bar_lgd = plt.scatter([0, 1], [0, 1], c=[0, 1], cmap="rocket")
  651. cbar = plt.colorbar(bar_lgd, cax=ax, orientation='vertical', ticks=[0, 1])
  652. cbar.set_label('Scaled expression', rotation=270, labelpad=10, fontsize=10)
  653. cbar.ax.tick_params(labelsize=12)
  654. cbar.ax.set_yticklabels(['0', '1'])
  655. ax.set_xticks([])
  656. ax.set_xticklabels([])
  657. plt.savefig('legend_heatmap.png', bbox_inches='tight', transparent=True, dpi=300)
  658. plt.show()
  662. # %%
  663. sc.set_figure_params(transparent=True, fontsize=14, dpi_save=300)
  664. sc.pl.rank_genes_groups_heatmap(adata,
  665. n_genes=10,
  666. groupby='Sample',
  667. cmap='rocket',
  668. dendrogram=True,
  669. figsize=(6, 10),
  670. use_raw=True,
  671. swap_axes=True,
  672. show_gene_labels=True,
  673. standard_scale='var',
  674. show=False,
  675. save='final_top10genes_heatmap.png'
  676. )
  678. # %%
  679. sc.tl.dendrogram(adata, 'Sample', n_pcs=20)
  680. axes_list = sc.pl.correlation_matrix(adata, 'Sample', show=False)
  682. for ax in axes_list:
  683. ax.yaxis.grid(False)
  684. ax.xaxis.grid(False)
  685. #plt.savefig('correlation_matrix.png', bbox_inches='tight', dpi=300, transparent=True)
  686. plt.show()
  689. # %% [markdown]
  690. # # markers in clusters
  692. # %%
  693. #finding markers in clusters
  694. sc.tl.rank_genes_groups(adata, groupby='leiden', method='t-test')
  695. sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, show=False, ncols=3)
  696. plt.savefig('final_t-test ranked genes leiden.png', bbox_inches='tight', dpi=300, transparent=True)
  699. # %%
  700. color_s2 = ['#0fbabd', '#efd39b', '#fbbd3c', '#bc131f', '#876c56']
  702. # %%
  703. sc.set_figure_params(transparent=True, dpi=80, fontsize=14, dpi_save=300, figsize=(10, 8))
  704. selected_genes = {
  705. 'Forebrain' :['SIX3', 'FOXG1', 'SP8', 'EOMES'],
  706. 'Midbrain':['CLSTN2', 'PTPRO','LMX1B', 'FGF8', 'WNT5A', 'EN2'],
  707. 'Hindbrain':['HOXB4', 'HOXB2', 'HOXA2', 'ERBB4', 'GBX2', 'PAX3'],
  708. 'Spinal\ncord':['HOXB9', 'HOXB8', 'HOXB7', 'HOXB6', 'CADPS2'],
  709. 'Seroton-\nergic' :['FEV', 'TPH2', 'SLC6A4','HTR7', 'HTR2A'],
  710. 'GABA-\nergic' :['GAD1', 'GAD2', 'SLC6A1'],
  711. 'Dopamin-\nergic' :['TH', 'SLC6A3'],
  712. 'Neural stem cells':['SOX1', 'SOX2', 'DLL1', 'FABP7', 'HES5', 'NES','PAX6', 'NOTCH1'],
  713. 'Differentiated ':['L1CAM','MAP2', 'TUBB3', 'NEFM', 'NEFL'],
  714. 'Schwann\ncells' :['MPZ', 'SOX10', 'S100B'],
  715. 'Motor\nneurons' :['LHX3', 'NKX6-1', 'AQP4'],
  716. }
  719. sc.pl.matrixplot(adata,
  720. selected_genes,
  721. #use_raw = False,
  722. groupby='leiden',
  723. dendrogram=True,
  724. cmap='Reds',
  725. #swap_axes=True,
  726. standard_scale='var',
  727. var_group_rotation=0,
  728. colorbar_title='column scaled\nexpression',
  729. #save='final_matrixplot.png'
  730. )
  732. # %%
  733. num_tot_cells = adata.obs.groupby(['Sample']).count()
  734. num_tot_cells = dict(zip(num_tot_cells.index, num_tot_cells.n_genes))
  735. cell_type_counts = adata.obs.groupby(['Sample', 'cell_type']).count()
  736. cell_type_counts = cell_type_counts[cell_type_counts.sum(axis=1)>0].reset_index()
  737. cell_type_counts = cell_type_counts.rename(columns={'Sample': 'sample_column_name'})
  738. cell_type_counts['total_cells'] = cell_type_counts['sample_column_name'].map(num_tot_cells).astype(int)
  739. cell_type_counts['frequency'] = cell_type_counts.leiden/cell_type_counts.total_cells
  741. # %%
  742. leiden_counts = adata.obs['leiden'].value_counts()
  743. total_cells = len(adata.obs)
  745. percentage_cells = (leiden_counts / total_cells) * 100
  746. plt.figure(figsize=(10, 6))
  747. ax = plt.axes() # Get the current axes
  748. percentage_cells.plot(kind='bar', color='skyblue', ax=ax)
  749. ax.grid(False) # Turn off the grid
  750. plt.title('Percentage of Cells in Each Cell Type Cluster')
  751. plt.xlabel('leiden cluster')
  752. plt.ylabel('Percentage of Cells')
  753. plt.xticks(rotation=45, ha='right')
  754. plt.savefig('final_percentage_of_cells.png', bbox_inches='tight', dpi=300)
  755. plt.show()
  758. # %%
  759. agg_data = adata.obs.groupby(['cell_type', 'Sample']).size().unstack().fillna(0)
  760. agg_data_normalized = agg_data.div(agg_data.sum(axis=1), axis=0) * 100
  761. colors = [color_s1[sample] for sample in agg_data_normalized.columns]
  763. plt.figure(figsize=(12, 8))
  764. plt.gca().set_facecolor('none')
  766. for i, sample in enumerate(agg_data_normalized.columns):
  767. plt.bar(
  768. agg_data_normalized.index,
  769. agg_data_normalized[sample],
  770. label=sample,
  771. color=colors[i],
  772. width=0.8,
  773. bottom=agg_data_normalized.iloc[:, :i].sum(axis=1),
  774. alpha=0.8
  775. )
  777. # Customization
  778. plt.ylabel('Percentage of cells', fontsize=20)
  779. plt.legend(bbox_to_anchor=(0.99, 1.024), loc='upper left', fontsize=30)
  780. plt.xticks(rotation=45, ha='right', fontsize=36)
  781. plt.yticks(fontsize=20)
  782. plt.ylim(0, 100)
  783. plt.grid(False)
  784. plt.savefig('final_Cell_Proportion.png', bbox_inches='tight', transparent=True, dpi=300)
  785. plt.show()
  787. # %%
  788. agg_data_normalized
  790. # %%
  792. my_order = ['6', '9', '10', '2', '12', '3', '7', '4', '0', '1', '14', '5', '16', '15', '11', '8', '13']
  793. #my_order = my_order[::-1]
  794. percentage_cells = percentage_cells.loc[my_order]
  796. fig, ax = plt.subplots(figsize=(2, 6))
  797. bars = ax.barh(percentage_cells.index, percentage_cells, color='black')
  798. ax.invert_yaxis()
  799. ax.yaxis.tick_right()
  800. ax.yaxis.set_label_position("right")
  801. ax.set_xlim(0, 25)
  802. ax.invert_xaxis()
  803. ax.yaxis.set_label_position("left")
  805. plt.ylabel('cluster')
  806. plt.xlabel('Percentage of Cells')
  807. plt.grid(False)
  809. for bar in bars:
  810. plt.text(bar.get_width() +8, bar.get_y() + bar.get_height() / 2, f'{bar.get_width():.1f}%',
  811. va='center', ha='left', fontsize=10, color='black')
  813. plt.show()
  816. # %%
  817. gene_index = np.where(adata.raw.var_names == 'HOXB9')[0]
  818. expressing_cells = np.sum(adata.raw.X[:, gene_index] > 0)
  819. print(f"The number of cells expressing SIX3 is: {expressing_cells}")
  822. # %%
  823. cell_genes = [gene for gene_list in selected_genes.values() for gene in gene_list]
  825. # Filter out non-gene entries
  826. gene_names = [gene for gene in cell_genes if isinstance(gene, str)]
  828. # Calculate the number of expressing cells for each gene
  829. expressing_cells_counts = []
  831. for gene in gene_names:
  832. gene_index = np.where(adata.raw.var_names == gene)[0]
  833. expressing_cells = np.sum(adata.raw.X[:, gene_index] > 0)
  834. expressing_cells_counts.append(expressing_cells)
  836. # Create a bar plot
  837. plt.figure(figsize=(22, 1))
  838. bars = plt.bar(gene_names, expressing_cells_counts, color='slategray')
  839. plt.xticks(rotation=90, ha='center')
  840. plt.grid(False)
  842. plt.gca().spines['top'].set_visible(False)
  843. plt.gca().spines['right'].set_visible(False)
  844. plt.gca().spines['bottom'].set_visible(False)
  845. plt.gca().spines['left'].set_visible(False)
  846. plt.yticks([])
  847. plt.tick_params(axis='y', which='both', left=False)
  848. plt.ylabel('Number of cells', fontsize=16, ha='center', va='center')
  851. # Add text labels on top of each bar
  852. for bar, count in zip(bars, expressing_cells_counts):
  853. plt.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 20, str(count),fontsize=14, ha='center', va='bottom', rotation=90)
  855. plt.show()
  858. # %%
  859. sc.set_figure_params(fontsize=20)
  862. fig = plt.figure(figsize=(22, 8))
  863. gs = gridspec.GridSpec(1, 2, width_ratios=[0.1, 1], height_ratios=[0.01]) # 1 row, 2 columns, with different widths
  866. ax = plt.subplot(gs[0])
  867. bars = ax.barh(percentage_cells.index, percentage_cells, color='black')
  868. ax.invert_yaxis()
  869. ax.yaxis.tick_right()
  870. ax.yaxis.set_label_position("right")
  871. ax.set_xlim(0, 20)
  872. ax.invert_xaxis()
  873. plt.xlabel('Percentage of cells')
  874. ax.yaxis.set_label_position("left")
  875. plt.ylabel('Clusters')
  876. plt.grid(False)
  879. for bar in bars:
  880. plt.text(bar.get_width() + 5, bar.get_y() + bar.get_height() / 2, f'{bar.get_width():.1f}%',
  881. va='center', ha='left', fontsize=14, color='black')
  884. ax.set_yticks([])
  885. ax.set_yticklabels([])
  888. ax.spines['top'].set_visible(False)
  889. ax.spines['right'].set_visible(False)
  890. ax.spines['bottom'].set_visible(True)
  891. ax.spines['left'].set_visible(False)
  894. ax_matrixplot = plt.subplot(gs[1])
  897. gs.update(right=0.90)
  898. plt.subplots_adjust(left=0.01, right=0.1, wspace=0.05)
  899. bar_plot_height = 0.705
  900. ax.set_position([ax.get_position().x0, ax.get_position().y0, ax.get_position().width, bar_plot_height])
  901. bar_plot_position = ax.get_position()
  902. new_bar_plot_position = [bar_plot_position.x0, bar_plot_position.y0 - 0.025, bar_plot_position.width, bar_plot_position.height]
  903. ax.set_position(new_bar_plot_position)
  904. sc.pl.matrixplot(adata, selected_genes, groupby='leiden', dendrogram=True, cmap='Reds',
  905. standard_scale='var', var_group_rotation=90, colorbar_title='column scaled\nexpression',
  906. ax=ax_matrixplot, show = False)
  910. plt.savefig('final_matrix_with_bar_plot.png', bbox_inches='tight', transparent=True, dpi=300)
  911. plt.show()
  914. # %%
  915. #pie charts
  916. sample_names = adata.obs['Sample'].unique()
  917. color_s1 = {'PN5': 'darkred', 'PN21': 'purple', 'PN41': 'darkorange', 'PN59': 'gold'}
  919. num_genes = len(cell_genes)
  920. fig, axes = plt.subplots(1, num_genes, figsize=(3 * num_genes, 3))
  923. for gene, ax in zip(cell_genes, axes):
  924. if gene not in adata.raw.var_names:
  925. print(f"Skipping gene '{gene}' as it is not present in the raw data.")
  926. continue
  928. expressing_cells_counts = []
  930. for sample in sample_names:
  931. sample_cells = adata.obs['Sample'] == sample
  932. gene_index = np.where(adata.raw.var_names == gene)[0]
  933. expressing_cells = np.sum((adata.raw.X[sample_cells, gene_index] > 0).A)
  934. expressing_cells_counts.append(expressing_cells)
  936. total_cells_per_sample = np.sum(adata.obs['Sample'].isin(sample_names))
  938. # Check if the total number of expressing cells is less than 10
  939. if np.sum(expressing_cells_counts) < 10:
  940. ax.pie([1], colors=['lightgrey'], autopct='', startangle=90)
  941. ax.set_title(f'{gene} Expression (< 10 cells)')
  942. continue
  944. percentage_expression = [count / total_cells_per_sample * 100 for count in expressing_cells_counts]
  945. colors = [color_s1.get(sample, 'gray') for sample in sample_names]
  946. ax.pie(percentage_expression, autopct='', startangle=90, colors=colors)
  947. ax.set_title(f'{gene}')
  949. plt.tight_layout()
  950. plt.show()
  953. # %%
  954. cell_genes = [gene for gene_list in selected_genes.values() for gene in gene_list]
  956. def plot_gene_pie(ax, gene, adata, color_s1):
  957. if gene not in adata.raw.var_names:
  958. print(f"Skipping gene '{gene}' as it is not present in the raw data.")
  959. return
  960. expressing_cells_counts = []
  962. for sample in sample_names:
  963. sample_cells = adata.obs['Sample'] == sample
  964. gene_index = np.where(adata.raw.var_names == gene)[0]
  965. expressing_cells = np.sum((adata.raw.X[sample_cells, gene_index] > 0).A)
  966. expressing_cells_counts.append(expressing_cells)
  968. total_cells_per_sample = np.sum(adata.obs['Sample'].isin(sample_names))
  970. if np.sum(expressing_cells_counts) < 10:
  971. ax.pie([1], colors=['lightgrey'], autopct='', startangle=90)
  972. ax.set_title('')
  973. ax.axis('equal')
  974. ax.set_xticks([])
  975. ax.set_yticks([])
  976. ax.set_facecolor('None')
  977. return
  979. percentage_expression = [count / total_cells_per_sample * 100 for count in expressing_cells_counts]
  981. colors = [color_s1.get(sample, 'gray') for sample in sample_names]
  983. ax.pie(percentage_expression, autopct='', startangle=90, colors=colors)
  984. ax.set_title('')
  985. ax.axis('equal')
  986. ax.set_xticks([])
  987. ax.set_yticks([])
  988. ax.set_facecolor('None')
  990. sample_names = adata.obs['Sample'].unique()
  991. color_s1 = {'PN5': 'darkred', 'PN21': 'purple', 'PN41': 'darkorange', 'PN59': 'gold'}
  992. num_genes = len(cell_genes)
  993. fig, axes = plt.subplots(1, num_genes, figsize=(3 * num_genes, 3))
  995. for i, gene in enumerate(cell_genes):
  996. ax_pie = axes[i]
  997. plot_gene_pie(ax_pie, gene, adata, color_s1)
  999. plt.subplots_adjust(wspace=0)
  1000. plt.show()
  1003. # %%
  1005. sc.set_figure_params(fontsize=20)
  1007. fig = plt.figure(figsize=(40, 14))
  1008. gs = gridspec.GridSpec(2, 2, width_ratios=[0.1, 1], height_ratios=[1, 1])
  1010. # Bar plot on the left
  1011. ax_left = plt.subplot(gs[1, 0])
  1012. bars_left = ax_left.barh(percentage_cells.index, percentage_cells, color='black')
  1013. ax_left.invert_yaxis()
  1014. ax_left.yaxis.tick_right()
  1015. ax_left.yaxis.set_label_position("right")
  1016. ax_left.set_xlim(0, 20)
  1017. ax_left.invert_xaxis()
  1018. ax_left.set_xlabel('Percentage of cells')
  1019. ax_left.yaxis.set_label_position("left")
  1020. ax_left.set_ylabel('Clusters')
  1021. ax_left.grid(False)
  1023. for bar in bars_left:
  1024. plt.text(bar.get_width() + 5, bar.get_y() + bar.get_height() / 2, f'{bar.get_width():.1f}%',
  1025. va='center', ha='left', fontsize=14, color='black')
  1027. ax_left.set_yticks([])
  1028. ax_left.set_yticklabels([])
  1029. ax_left.spines['top'].set_visible(False)
  1030. ax_left.spines['right'].set_visible(False)
  1031. ax_left.spines['bottom'].set_visible(True)
  1032. ax_left.spines['left'].set_visible(False)
  1034. # Matrix plot in the center
  1035. ax_matrixplot = plt.subplot(gs[1, 1])
  1036. sc.pl.matrixplot(adata, selected_genes, groupby='leiden', dendrogram=True, cmap='Reds',
  1037. standard_scale='var', var_group_rotation=90, colorbar_title='column scaled\nexpression',
  1038. ax=ax_matrixplot, show=False)
  1040. # Bar plot on top
  1041. ax_top = plt.subplot(gs[0, 1])
  1042. bars_top = ax_top.bar(gene_names, expressing_cells_counts, color='slategray')
  1043. ax_top.set_xticks([])
  1044. ax_top.grid(False)
  1045. ax_top.spines['top'].set_visible(False)
  1046. ax_top.spines['right'].set_visible(False)
  1047. ax_top.spines['bottom'].set_visible(False)
  1048. ax_top.spines['left'].set_visible(False)
  1049. ax_top.tick_params(axis='y', which='both', left=False, labelleft=False)
  1053. for bar, count in zip(bars_top, expressing_cells_counts):
  1054. ax_top.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 20, str(count),
  1055. fontsize=14, ha='center', va='bottom', rotation=90)
  1057. ax_top.set_ylabel('Number of cells', fontsize=16, ha='right', va='center', rotation=0)
  1060. gs.update(right=0.90)
  1061. gs.update(bottom=0.1, top=1, wspace=0.05, hspace=0.1)
  1063. ax_left_pos = ax_left.get_position()
  1064. ax_left.set_position([ax_left_pos.x0 - 0.03, ax_left_pos.y0 - 0.02, ax_left_pos.width * 1.5, ax_left_pos.height*0.94])
  1066. ax_top_pos = ax_top.get_position()
  1067. ax_top.set_position([ax_top_pos.x0 - 0.028, ax_top_pos.y0 - -0.13, ax_top_pos.width* 0.965, ax_top_pos.height* 0.2])
  1069. ax_matrixplot_pos = ax_matrixplot.get_position()
  1070. ax_matrixplot.set_position([ax_matrixplot_pos.x0, ax_matrixplot_pos.y0, ax_matrixplot_pos.width, ax_matrixplot_pos.height])
  1072. plt.show()
  1076. # %%
  1078. sc.set_figure_params(fontsize=20)
  1081. fig = plt.figure(figsize=(40, 18))
  1082. gs = gridspec.GridSpec(3, 2, width_ratios=[0.1, 1], height_ratios=[1, 1, 1])
  1084. ax_left = plt.subplot(gs[2, 0])
  1085. bars_left = ax_left.barh(percentage_cells.index, percentage_cells, color='black')
  1086. ax_left.invert_yaxis()
  1087. ax_left.yaxis.tick_right()
  1088. ax_left.yaxis.set_label_position("right")
  1089. ax_left.set_xlim(0, 20)
  1090. ax_left.invert_xaxis()
  1091. ax_left.set_xlabel('Percentage of cells')
  1092. ax_left.yaxis.set_label_position("left")
  1093. ax_left.set_ylabel('Clusters')
  1094. ax_left.grid(False)
  1096. for bar in bars_left:
  1097. plt.text(bar.get_width() + 5, bar.get_y() + bar.get_height() / 2, f'{bar.get_width():.1f}%',
  1098. va='center', ha='left', fontsize=14, color='black')
  1100. ax_left.set_yticks([])
  1101. ax_left.set_yticklabels([])
  1102. ax_left.spines['top'].set_visible(False)
  1103. ax_left.spines['right'].set_visible(False)
  1104. ax_left.spines['bottom'].set_visible(True)
  1105. ax_left.spines['left'].set_visible(False)
  1106. ax_left.set_zorder(2)
  1108. ax_matrixplot = plt.subplot(gs[2, 1])
  1109. sc.pl.matrixplot(adata, selected_genes, groupby='leiden', dendrogram=True, cmap='Reds',
  1110. standard_scale='var', var_group_rotation=0, colorbar_title='column scaled\nexpression',
  1111. ax=ax_matrixplot, show=False)
  1112. ax_matrixplot.set_rasterization_zorder(2)
  1115. ax_top = plt.subplot(gs[0, 1])
  1116. bars_top = ax_top.bar(gene_names, expressing_cells_counts, color='slategray')
  1117. ax_top.set_xticks([])
  1118. ax_top.grid(False)
  1119. ax_top.spines['top'].set_visible(False)
  1120. ax_top.spines['right'].set_visible(False)
  1121. ax_top.spines['bottom'].set_visible(False)
  1122. ax_top.spines['left'].set_visible(False)
  1123. ax_top.tick_params(axis='y', which='both', left=False, labelleft=False)
  1127. for bar, count in zip(bars_top, expressing_cells_counts):
  1128. ax_top.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 20, str(count),
  1129. fontsize=14, ha='center', va='bottom', rotation=90)
  1131. ax_top.set_ylabel('Number of cells', fontsize=16, ha='right', va='center', rotation=0)
  1135. sample_names = adata.obs['Sample'].unique()
  1136. color_s1 = {'PN5': 'darkred', 'PN21': 'purple', 'PN41': 'darkorange', 'PN59': 'gold'}
  1139. num_genes = len(cell_genes)
  1140. ax_pie = plt.subplot(gs[1, 0:2])
  1142. total_width = 0.8
  1143. width_per_pie = total_width / num_genes
  1146. for i, gene in enumerate(cell_genes):
  1147. ax_pie_sub = ax_pie.inset_axes([i * width_per_pie, 0, width_per_pie, 1])
  1148. plot_gene_pie(ax_pie_sub, gene, adata, color_s1)
  1149. ax_pie.grid(False)
  1150. ax_pie.set_frame_on(False)
  1151. ax_pie.set_xticks([])
  1152. ax_pie.set_yticks([])
  1153. ax_pie.spines['top'].set_visible(False)
  1154. ax_pie.spines['right'].set_visible(False)
  1155. ax_pie.spines['bottom'].set_visible(False)
  1156. ax_pie.spines['left'].set_visible(False)
  1159. gs.update(right=0.90)
  1160. gs.update(bottom=0.1, top=1, wspace=0.05, hspace=0.1)
  1162. ax_left_pos = ax_left.get_position()
  1163. ax_left.set_position([ax_left_pos.x0 -0.009, ax_left_pos.y0 - 0.011, ax_left_pos.width * 1.2, ax_left_pos.height*0.94])
  1165. ax_top_pos = ax_top.get_position()
  1166. ax_top.set_position([ax_top_pos.x0 - 0.0278, ax_top_pos.y0 - 0.28, ax_top_pos.width* 0.965, ax_top_pos.height* 0.2])
  1168. ax_matrixplot_pos = ax_matrixplot.get_position()
  1169. ax_matrixplot.set_position([ax_matrixplot_pos.x0, ax_matrixplot_pos.y0, ax_matrixplot_pos.width, ax_matrixplot_pos.height])
  1171. ax_pie_pos = ax_pie.get_position()
  1172. ax_pie.set_position([ax_pie_pos.x0 + 0.081, ax_pie_pos.y0 -0.14, ax_pie_pos.width*0.975, ax_pie_pos.height])
  1175. plt.savefig('final_heatmap.png', bbox_inches='tight', transparent=True, dpi=300)
  1176. plt.show()
  1180. # %%
  1181. legend_handles2 = [Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=10, label=label) for label, color in color_s1.items()]
  1182. legend_handles2.append(Patch(color='lightgrey', label='< 10 cells'))
  1183. fig, ax = plt.subplots(figsize=(1, 1))
  1184. ax.axis('off')
  1185. ax.legend(handles=legend_handles2, ncol=1, loc='center')
  1186. plt.savefig('legend_heatmap.png', bbox_inches='tight', transparent=True, dpi=300)
  1187. plt.show()

scRNAseq_HbL-iNSCs.ipynb, under CC-BY-4.0 · at the source

Overview

Authors: Ziadoon Al-Akashi1, Denise Zujur1, Nicholas Boyd-Gibbins1, Nathalie Eileen Wiguna1, Masato Nakagawa2, Tetsuhiro Kikuchi1, Asuka Morizane1,3, Jun Takahashi1, Makoto Ikeya1
ORCID iDs: Ziadoon Al-Akashi, Masato Nakagawa, Makoto Ikeya
  1. Department of Clinical Application, Center for iPS Cell Research and Application (CiRA), Kyoto University, 53 Kawahara-cho, Shogoin, Sakyo-ku, Kyoto 606-8507, Japan
  2. Department of Life Science Frontiers, Center for iPS Cell Research and Application (CiRA), Kyoto University, 53 Kawahara-cho, Shogoin, Sakyo-ku, Kyoto 606-8507, Japan
  3. Department of Regenerative Medicine, Center for Clinical Research and Innovation, Kobe City Medical Center General Hospital, Hyogo, Japan
Institutions: Kyoto University (Japan); Center for iPS Cell Research and Application, Kyoto University (Japan); Kobe City Medical Center General Hospital (Japan)
Journal: Cell reports methods, volume 6, issue 4, article 101372
Dates: received 23 May 2025; accepted 27 February 2026; published online 30 March 2026; in print April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.crmeth.2026.101372 · PMID 41916305 · PMCID PMC13107045 · OpenAlex W7143485316
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism)
Methods: Spectral & time-frequency
Keywords: hindbrain, neural stem cell, neuron, pluripotent stem cell
MeSH: Induced Pluripotent Stem Cells*, Neural Stem Cells*, Rhombencephalon*, Smad Proteins*, Wnt Proteins*, Wnt Signaling Pathway*, Animals, Cell Differentiation, Humans, Pyrazoles, Pyridines, Pyrimidines (* major topic)
Topic: Neurogenesis and neuroplasticity mechanisms (Developmental Neuroscience, Neuroscience), according to OpenAlex
Funding: Japan Agency for Medical Research and Development (JP23bm1323001, JP15bm0104001); Otsuka Toshimi Scholarship Foundation; Japan Society for the Promotion of Science (16H05447); Japan Science and Technology Agency (JPMJSP2110)
Citations: not cited yet (Europe PMC); 22 references in the paper
Research resources: RRID:CVCL_LJ39

Abstract

Neurons exhibit region-specific identities corresponding to functional distinctions across different brain areas. Region-restricted neural stem cells (NSCs) have previously been generated from pluripotent stem cells; however, maintaining their regional identity over extended passages remains challenging. Here, we report the generation of hindbrain-like induced NSCs (Hb-LiNSCs) with upregulated hindbrain-specific markers and downregulated forebrain, midbrain, and spinal cord markers under xeno-free and basic fibroblast growth factor (bFGF)-free conditions using three chemicals—CHIR99021 (at a high concentration), a potent activator of the Wnt pathway; A-83-01, a potent inhibitor of the TGF-β/Activin/Nodal pathway; and LDN193189, a potent inhibitor of the bone morphogenetic protein pathway. Hb-LiNSCs maintained their chromosomal integrity, multipotency, and differentiation capacity even after long-term culture for more than 60 weeks. This approach enhances our understanding of neurodevelopmental and neurodegenerative processes in the hindbrain region and paves the way for developing targeted cell-based therapy as well as disease modeling for drug discovery.

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

Repository

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

Zenodo 18653633

License: CC-BY-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: R (4), Jupyter (1)
Size: 5 files, 5 scripts
Software Heritage: not checked
Found in: “Data and code availability”
Holds: 5 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: anndata (1 file), Matplotlib (1 file), NumPy (1 file), pandas (1 file), Scanpy (1 file), scikit-learn (1 file), seaborn (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
5 files
At the source: doi.org/10.5281/zenodo.18653633

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:

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

Single-cell RNAseq data and RNAseq data generated in this study have been deposited to GEO: GSE317025 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE317025) and GSE317319 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE317319). We used RNAseq data of human iPSCs deposited to GEO: GSE206128 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE206128). They are publicly available as of the date of publication. Other data and images that support the findings of this study are available on request from the lead contact.

The code for generating most of the figures can be found at https://doi.org/10.5281/zenodo.18653633.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 4 keywords, 12 MeSH terms, 4 funders, 22 references, 1 RRID.

Cite

This paper

Al-Akashi, Z., Zujur, D., Boyd-Gibbins, N., Wiguna, N. E., Nakagawa, M., Kikuchi, T., Morizane, A., Takahashi, J., & Ikeya, M. (2026). Wnt activation and dual SMAD inhibition for induction and maintenance of hindbrain-like neural stem cell from hiPSCs. Cell reports methods, 6(4), 101372. https://doi.org/10.1016/j.crmeth.2026.101372

BibTeX

@article{alakashi2026wnt,
author = {Al-Akashi, Ziadoon and Zujur, Denise and Boyd-Gibbins, Nicholas and Wiguna, Nathalie Eileen and Nakagawa, Masato and Kikuchi, Tetsuhiro and Morizane, Asuka and Takahashi, Jun and Ikeya, Makoto},
title = {{Wnt activation and dual SMAD inhibition for induction and maintenance of hindbrain-like neural stem cell from hiPSCs}},
journal = {Cell reports methods},
year = {2026},
month = mar,
volume = {6},
number = {4},
pages = {101372},
publisher = {Elsevier},
issn = {2667-2375},
doi = {10.1016/j.crmeth.2026.101372},
url = {https://doi.org/10.1016/j.crmeth.2026.101372},
pmid = {41916305},
pmcid = {PMC13107045}
}

RIS

TY - JOUR
AU - Al-Akashi, Ziadoon
AU - Zujur, Denise
AU - Boyd-Gibbins, Nicholas
AU - Wiguna, Nathalie Eileen
AU - Nakagawa, Masato
AU - Kikuchi, Tetsuhiro
AU - Morizane, Asuka
AU - Takahashi, Jun
AU - Ikeya, Makoto
TI - Wnt activation and dual SMAD inhibition for induction and maintenance of hindbrain-like neural stem cell from hiPSCs
T2 - Cell reports methods
J2 - Cell Rep Methods
PY - 2026
DA - 2026/03/30
VL - 6
IS - 4
SP - 101372
SN - 2667-2375
PB - Elsevier
DO - 10.1016/j.crmeth.2026.101372
UR - https://doi.org/10.1016/j.crmeth.2026.101372
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.crmeth.2026.101372",
"type": "article-journal",
"title": "Wnt activation and dual SMAD inhibition for induction and maintenance of hindbrain-like neural stem cell from hiPSCs",
"container-title": "Cell reports methods",
"author": [
{
"family": "Al-Akashi",
"given": "Ziadoon"
},
{
"family": "Zujur",
"given": "Denise"
},
{
"family": "Boyd-Gibbins",
"given": "Nicholas"
},
{
"family": "Wiguna",
"given": "Nathalie Eileen"
},
{
"family": "Nakagawa",
"given": "Masato"
},
{
"family": "Kikuchi",
"given": "Tetsuhiro"
},
{
"family": "Morizane",
"given": "Asuka"
},
{
"family": "Takahashi",
"given": "Jun"
},
{
"family": "Ikeya",
"given": "Makoto"
}
],
"container-title-short": "Cell Rep Methods",
"volume": "6",
"issue": "4",
"page": "101372",
"DOI": "10.1016/j.crmeth.2026.101372",
"PMID": "41916305",
"PMCID": "PMC13107045",
"ISSN": "2667-2375",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.crmeth.2026.101372",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
30
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41592-024-02412-5 [code]
A polarized FGF8 source specifies frontotemporal signatures in spatially oriented cell populations of cortical assembloids.
Journal: Nature methods
In common: anndata, Scanpy, seaborn, 4 other tools, 3 references
[2] doi:10.1038/s41586-025-09151-3 [code]
Morphodynamics of human early brain organoid development.
Journal: Nature
In common: anndata, Scanpy, seaborn, 4 other tools, 2 references
[3] doi:10.1038/s41593-026-02208-0 [code]
Single-nucleus epigenomic profiling of the adult human central nervous system unveils epigenetic memory of developmental programs.
Journal: Nature neuroscience
In common: anndata, Scanpy, seaborn, 4 other tools, 2 references
[4] doi:10.1126/sciadv.adq4731 [code]
An ancient apical patterning system sets the position of the forebrain in chordates.
Journal: Science advances
In common: anndata, Scanpy, seaborn, 3 other tools, 3 references
[5] doi:10.1242/dev.202545 [code]
Exploiting spatiotemporal regulation of FZD5 during neural patterning for efficient ventral midbrain specification.
Journal: Development (Cambridge, England)
In common: anndata, Scanpy, seaborn, 4 other tools, 2 references
[6] doi:10.1016/j.celrep.2025.116498 [code]
In vivo differentiation of embryonic cells devoid of key reprogramming factors.
Journal: Cell reports
In common: anndata, Scanpy, seaborn, 4 other tools, 2 references
[7] doi:10.1038/s41593-026-02200-8 [code]
Distinct radial glia subtypes regulate midbrain dopaminergic neuron development.
Journal: Nature neuroscience
In common: Scanpy, seaborn, pandas, 2 other tools, 3 references
[8] doi:10.1016/j.isci.2025.114482
Embryonic lineage-specific iPSC-derived mesenchymal stem/stromal cells exhibit different morphologies and intrinsic functions.
Journal: iScience
In common: 2 references, author Makoto Ikeya
[9] doi:10.7554/elife.102819 [code]
A single-cell atlas of spatial and temporal gene expression in the mouse cranial neural plate.
Journal: eLife
In common: Scanpy, seaborn, scikit-learn, 3 other tools, 2 references
[10] doi:10.1093/nargab/lqae100 [code]
SIRV: spatial inference of RNA velocity at the single-cell resolution.
Journal: NAR genomics and bioinformatics
In common: Scanpy, seaborn, scikit-learn, 3 other tools, 2 references

Contribute

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

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

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1016/j.crmeth.2026.101372.svg)](https://openscicode.org/badge/doi/10.1016/j.crmeth.2026.101372)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1016/j.crmeth.2026.101372.svg
   :target: https://openscicode.org/badge/doi/10.1016/j.crmeth.2026.101372
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1016/j.crmeth.2026.101372"><img src="https://openscicode.org/badge/doi/10.1016/j.crmeth.2026.101372.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

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