An emergent disease-associated motor neuron state precedes cell death in ALS.
The 40 matches
- [1] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Cholinergic Neuron Peak Calling with Downsampled Datasets ↔ 07. Multiome 3 - input preparation for Fig. 3D-F, S4E-H.ipynb, lines 358–384 · score 0.97 · bgdGroups, maxCells, useGroups, PeakMatrix, testMethod, getMarkerFeatures
- [2] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Alpha Motor Neuron Differentially Accessible Peaks with Disease ↔ 07. Multiome 3 - input preparation for Fig. 3D-F, S4E-H.ipynb, lines 358–384 · score 0.97 · bgdGroups, maxCells, useGroups, PeakMatrix, testMethod, getMarkerFeatures
- [3] § STAR★METHODS › METHOD DETAILS › snRNA-seq/Multiome Sequencing: Label Transfer from Multiome to snRNA-seq Alpha Motor Neurons ↔ 11. Multiome & snRNA-seq 1 - input preparation (cross-modal label transfer, multiome RNA differential expression).ipynb, lines 169–174 · score 0.96 · FindTransferAnchors, TransferData, query.assay, reference.assay, AddMetaData, multiome RNA
- [4] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Cholinergic Neuron and Alpha Motor Neuron Subclustering ↔ 05. Multiome 1 - initial processing and input preparation for Fig. 3B-C, 3G-J, S4A-D, S5A-B.ipynb, lines 383–400 · score 0.94 · Gex_nUMI, firstSelection, varFeatures, depthCol, clusterParams, addIterativeLSI
- [5] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Clustering and Doublet Removal ↔ 05. Multiome 1 - initial processing and input preparation for Fig. 3B-C, 3G-J, S4A-D, S5A-B.ipynb, lines 177–194 · score 0.94 · Gex_nUMI, firstSelection, varFeatures, depthCol, clusterParams, addIterativeLSI
- [6] § RESULTS › A DM signature in alpha motor neurons ↔ 04. snRNA-seq 4 - Fig. 2A-J, S2E, S3A-H.ipynb, lines 433–531 · score 0.93 · synaptic transmission, axon guidance, proteasomal protein catabolism, endoplasmic reticulum, unfolded protein, amino acid
- [7] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Cholinergic Neuron and Alpha Motor Neuron Subclustering ↔ 01. snRNA-seq 1 - initial processing and input preparation.ipynb, lines 539–550 · score 0.92 · RunUMAP, FindClusters, FindNeighbors, RunPCA, ScaleData, DefaultAssay
- [8] § STAR★METHODS › METHOD DETAILS › MERFISH: Generation of a Custom Motor Neuron Segmentation Model ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 388–487 · score 0.92 · Aldh1l1, Cx3cr1, Slc6a5, Slc5a7, Slc17a6, Trem2
- [9] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Identification of Positive Transcription Factor Regulators ↔ 09. Multiome 5 - input preparation for Fig. 4A-D, S6A-G.ipynb, lines 209–216 · score 0.92 · correlateMatrices, useMatrix1, useMatrix2, LSI_Combined, MotifMatrix, reducedDims
- [10] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Data Pre-Processing and Ambient RNA Removal ↔ 01. snRNA-seq 1 - initial processing and input preparation.ipynb, lines 76–89 · score 0.92 · NormalizeData, FindClusters, FindNeighbors, RunPCA, ScaleData, FindVariableFeatures
- [11] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Cholinergic Neuron Peak Calling with Downsampled Datasets ↔ 07. Multiome 3 - input preparation for Fig. 3D-F, S4E-H.ipynb, lines 168–192 · score 0.90 · minReplicates, minCells, addGroupCoverages, addPeakMatrix, addReproduciblePeakSet, groupBy
- [12] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Cholinergic Neuron Peak Calling with Downsampled Datasets ↔ 07. Multiome 3 - input preparation for Fig. 3D-F, S4E-H.ipynb, lines 168–192 · score 0.90 · minReplicates, minCells, addGroupCoverages, addPeakMatrix, addReproduciblePeakSet, peak matrices
- [13] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Integration, Quality Control, and Clustering ↔ 01. snRNA-seq 1 - initial processing and input preparation.ipynb, lines 215–223 · score 0.86 · RunUMAP, FindClusters, FindNeighbors, RunPCA, ScaleData, DefaultAssay
- [14] § RESULTS › TFs associated with the DM state transition ↔ 10. Multiome 6 - Fig. 4A-D, S6A-G.ipynb, lines 64–126 · score 0.86 · E4f1, ps sum, perturbation score, CellOracle, Jund, Klf6
- [15] § RESULTS › A DM signature in alpha motor neurons ↔ 04. snRNA-seq 4 - Fig. 2A-J, S2E, S3A-H.ipynb, lines 880–949 · score 0.84 · extracellular matrix organization, potassium ion transport, integrated stress response, amino acid, negative regulation, signaling
- [16] § STAR★METHODS › METHOD DETAILS › In Vitro Motor Neuron Differentiation, Lentiviral Transduction, Western Blot Analysis, Bulk RNA Sequencing, and Gene Set Enrichment Analysis ↔ 13. iMN in vitro TF OE - Fig. 5.ipynb, lines 131–205 · score 0.83 · fgseaMultilevel, minSize, maxSize, DESeq2, NES, ranked
- [17] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Differential Expression and Gene Ontology (GO) Enrichment Analyses ↔ 02. snRNA-seq 2 - DESeq2 differential expression analysis.ipynb, lines 136–163 · score 0.82 · DESeqDataSetFromMatrix, minReplicatesForReplace, ex, DESeq2, padj, frame
- [18] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Integration, Quality Control, and Clustering ↔ 01. snRNA-seq 1 - initial processing and input preparation.ipynb, lines 360–362 · score 0.82 · Read10X, endpoint female, nextseqs_11_18, CreateSeuratObject, cellranger, intron
- [19] § STAR★METHODS › METHOD DETAILS › MERFISH: Cholinergic Neuron and Alpha Motor Neuron Subclustering ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 734–748 · score 0.80 · n_iterations, min_dist, igraph, flavor, BBKNN, external
- [20] § STAR★METHODS › METHOD DETAILS › MERFISH: Alpha Motor Neuron Morphological Quantification ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 5867–5891 · score 0.79 · Otsu threshold, DAPI_high_pass, log10 transformed, MERFISH
- [21] § RESULTS › A DM signature in alpha motor neurons ↔ 13. iMN in vitro TF OE - Fig. 5.ipynb, lines 591–660 · score 0.75 · axon guidance, proteasomal protein catabolism, endoplasmic reticulum, axonogenesis, synaptic, transport
- [22] § RESULTS › TFs associated with the DM state transition ↔ 10. Multiome 6 - Fig. 4A-D, S6A-G.ipynb, lines 64–126 · score 0.72 · Nfe2l1, Nfe2l2, Nfe2l3, Xbp1, Nfil3, CREB3
- [23] § STAR★METHODS › METHOD DETAILS › MERFISH: Data Pre-Processing and Initial Clustering ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 734–748 · score 0.71 · n_iterations, igraph, flavor, BBKNN, external, UMAP
- [24] § RESULTS › Gene expression changes in spinal motor neurons during neurodegeneration ↔ 03. snRNA-seq 3 - Fig. 2B, 2D-H, S1B-C, S2A-D.ipynb, lines 663–732 · score 0.70 · sciatic nerve crush, nerve injury, SOD1 G93A, gamma, disease, gene
- [25] § STAR★METHODS › METHOD DETAILS › MERFISH: Reactive Glial Cell Classification ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 3759–3813 · score 0.69 · control astrocytes, control microglia, Apoe, Gfap, thresholds, reactive
- [26] § STAR★METHODS › METHOD DETAILS › snRNA-seq/Multiome Sequencing: In Silico Transcription Factor Perturbation with CellOracle ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 5194–5271 · score 0.69 · fast firing, slow firing, late DM, early DM, UMAP, Scanpy
- [27] § STAR★METHODS › METHOD DETAILS › MERFISH: Label Transfer from snRNA-seq to MERFISH Alpha Motor Neurons ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 5051–5192 · score 0.69 · Fast Firing, Slow Firing, Late DM, Early DM, Predicted, Intermediate
- [28] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Integration, Quality Control, and Clustering ↔ 01. snRNA-seq 1 - initial processing and input preparation.ipynb, lines 141–150 · score 0.65 · nova_CZI, sod1_mn_nuclei_2_9, seq
- [29] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Integration, Quality Control, and Clustering ↔ 11. Multiome & snRNA-seq 1 - input preparation (cross-modal label transfer, multiome RNA differential expression).ipynb, lines 157–160 · score 0.65 · NormalizeData, FindVariableFeatures, vst, selection, seq
- [30] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Cholinergic Neuron Peak Calling with Downsampled Datasets ↔ 08. Multiome 4 - Fig. 3D-F, S4E-H.ipynb, lines 268–334 · score 0.63 · downsampled control, Gamma MNs, Alpha MNs, ArchR, Peak, cholinergic
- [31] § STAR★METHODS › METHOD DETAILS › Multiome Sequencing: Cholinergic Neuron Peak Calling with Downsampled Datasets ↔ 07. Multiome 3 - input preparation for Fig. 3D-F, S4E-H.ipynb, lines 243–250 · score 0.62 · Gamma MNs, metadata column, Alpha MNs, ArchR, downsampled, neurons
- [32] § STAR★METHODS › METHOD DETAILS › Human Spinal Cord snRNA-seq/Fragment-seq: Cross-Species Wilcoxon Rank-Based Gene Set Enrichment and DM Scoring ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 1332–1420 · score 0.61 · log2 fold change, downregulated genes, upregulated genes, snRNA, seq
- [33] § STAR★METHODS › METHOD DETAILS › MERFISH: Data Pre-Processing and Initial Clustering ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 168–222 · score 0.61 · detected genes, filtered AnnData, blank, MERFISH, slide, transcripts
- [34] § STAR★METHODS › METHOD DETAILS › Human Spinal Cord snRNA-seq/Fragment-seq: Data Pre-Processing and Integration ↔ 11. Multiome & snRNA-seq 1 - input preparation (cross-modal label transfer, multiome RNA differential expression).ipynb, lines 157–160 · score 0.60 · NormalizeData, FindVariableFeatures, nfeature, vst, RNA seq
- [35] § STAR★METHODS › METHOD DETAILS › MERFISH: Alpha Motor Neuron Differential Expression Analysis with Disease ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 1250–1324 · score 0.59 · Benjamini Hochberg, Fold changes, raw, log2, MERFISH, mid
- [36] § STAR★METHODS › METHOD DETAILS › MERFISH: Cell Segmentation, Transcript Partitioning, and Cell Metadata Calculation ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 69–166 · score 0.59 · sum signals, cell metadata, DAPI, MERFISH, transcript
- [37] § STAR★METHODS › METHOD DETAILS › MERFISH: Generation of a Custom Motor Neuron Segmentation Model ↔ 05. Multiome 1 - initial processing and input preparation for Fig. 3B-C, 3G-J, S4A-D, S5A-B.ipynb, lines 227–238 · score 0.57 · Slc5a7, Slc17a6, Aqp4, Mog, Atf3, cholinergic
- [38] § STAR★METHODS › METHOD DETAILS › snRNA-seq: Differential Expression and Gene Ontology (GO) Enrichment Analyses ↔ 04. snRNA-seq 4 - Fig. 2A-J, S2E, S3A-H.ipynb, lines 41–61 · score 0.57 · GO Biological Process, DESeq2, downregulated, upregulated, cholinergic, gene
- [39] § STAR★METHODS › METHOD DETAILS › Human Spinal Cord snRNA-seq/Fragment-seq: Differential Gene Expression Analysis ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 1250–1324 · score 0.56 · Benjamini Hochberg, fold change, FDR, raw, Gene
- [40] § STAR★METHODS › METHOD DETAILS › Human Spinal Cord snRNA-seq/Fragment-seq: Cross-Species Wilcoxon Rank-Based Gene Set Enrichment and DM Scoring ↔ 14. MERFISH spatial transcriptomics analysis.ipynb, lines 1332–1420 · score 0.54 · log2 fold change, downregulated genes, v1, upregulated, snRNA, seq
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 · 6,326 lines · 198 KB · no license · 13 matches
- # %%
- import anndata as ad
- import scanpy as sc
- import pandas as pd
- import os
- import numpy as np
- import squidpy as sq
- import bbknn
- import seaborn as sns
- import matplotlib.pyplot as plt
- from matplotlib.colors import ListedColormap
- import shutil
- import scvi
- import scipy.sparse
- from scipy.stats import pearsonr, spearmanr
- # %%
- working_dir = "/home/users/ogautier/oak/Shared/SOD1_Paper/Vizgen/mn_nonmn_segmentation"
- os.chdir(working_dir)
- figures_dir = os.path.join(working_dir, "figures")
- # %% [markdown]
- # ### Create anndata objects with metadata
- # %%
- def is_point_inside_rotated_rect(px, py, rx, ry, width, height, angle):
- """
- Check if a point (px, py) is inside a rotated rectangle.
- - (rx, ry): Bottom-left corner of the rectangle
- - width, height: Dimensions of the rectangle
- - angle: Rotation angle in degrees
- """
- # Convert angle to radians
- theta = np.radians(-angle) # Negative to rotate in the opposite direction
- # Translate point and rectangle to origin
- px, py = px - rx, py - ry
- # Apply inverse rotation to the point
- rotated_x = px * np.cos(theta) - py * np.sin(theta)
- rotated_y = px * np.sin(theta) + py * np.cos(theta)
- # Check if rotated point is within rectangle bounds
- return 0 <= rotated_x <= width and 0 <= rotated_y <= height
- # %%
- slide_list = ['202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8',
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801',
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801',
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101',
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8',
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701',
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8',
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302',
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis',
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802',
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802',
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011',
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802',
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8',
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011',
- #'202405281422_MsSpinalCord-VS223-Lumbar-CE3-SE5-S1_VMSC12002', Exclude due to batch effect
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201',
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701',
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802']
- # %%
- csv_file_path = 'Merfish_Metadata.csv'
- rect_df_all = pd.read_csv(csv_file_path)
- data_fn = "analysis_outputs_7000_900/cell_by_gene.csv"
- meta_fn = "analysis_outputs_7000_900/cell_metadata.csv"
- sumsig_fn = "analysis_outputs_7000_900/sum_signals.csv"
- for slide_dir in slide_list:
- print(slide_dir)
- # read in and make anndata object
- adata = sc.read_csv(os.path.join(slide_dir, data_fn), first_column_names=True)
- print(f'adata initial shape is {adata.shape}')
- metadata = pd.read_csv(os.path.join(slide_dir,meta_fn), index_col=0)
- metadata.index = metadata.index.map(str)
- ### make sumsigs consistent start
- sumsig = pd.read_csv(os.path.join(slide_dir,sumsig_fn), index_col=0)
- sumsig.index = sumsig.index.map(str)
- if 'hSOD1G93A_raw' not in sumsig:
- sumsig['hSOD1G93A_raw'] = [float('nan')] * sumsig['DAPI_raw'].shape[0]
- sumsig['hSOD1G93A_high_pass'] = [float('nan')] * sumsig['DAPI_raw'].shape[0]
- ### end
- result = pd.merge(metadata, sumsig, left_index=True, right_index=True)
- result.index = adata.obs_names
- adata.obs = result
- adata.obsm['spatial'] = adata.obs[["center_x", "center_y"]].to_numpy()
- # get additional critical metadata from rect_df_all
- rect_df = rect_df_all[rect_df_all['slide'] == slide_dir]
- x, y = adata.obsm['spatial'][:, 0], adata.obsm['spatial'][:, 1]
- # initialize metadata columns
- slide = ['na'] * len(x)
- section = ['na'] * len(x)
- region = ['na'] * len(x)
- stage = ['na'] * len(x)
- # loop through rect_df to populate metadata of appropriate cells
- # rect_df allows us to map a coordinate to section and its
- # corresponding conditions
- for index, row in rect_df.iterrows():
- rect_x, rect_y = row['x']-150, row['y']-150
- rect_width, rect_height = int(row['w']*1.2), int(row['h']*1.2)
- rect_angle = row['angle']
- for i in range(len(x)):
- point_x, point_y = adata.obsm['spatial'][i, 0], adata.obsm['spatial'][i, 1]
- inside = is_point_inside_rotated_rect(point_x,
- point_y,
- rect_x,
- rect_y,
- rect_width,
- rect_height,
- rect_angle)
- if inside == True:
- section[i] = row['section']
- slide[i] = row['slide']
- region[i] = row['region']
- stage[i] = row['stage']
- # append metadata lists to adata.obs
- adata.obs["section"] = section
- adata.obs["slide"] = slide
- adata.obs["region"] = region
- adata.obs["stage"] = stage
- # remove cells that are not assigned section metadata
- # Ensure 'section' is treated as a string
- adata.obs["section"] = adata.obs["section"].astype(str)
- adata = adata[adata.obs["section"] != "na"].copy()
- print(f'adata shape after filtering is {adata.shape}')
- # make VS metadata column
- VS = []
- slide = adata.obs['slide'].to_list()
- for s in slide:
- if "VS119" in s:
- VS.append("VS119")
- elif "VS223" in s:
- VS.append("VS223")
- else:
- print("uh oh")
- adata.obs['VS'] = VS
- list(set(adata.obs['stage'].to_list()))
- print(f"adata.obs.shape is {adata.obs.shape}")
- # save anndata
- save_name = f"{slide_dir}.h5ad"
- adata.write_h5ad(os.path.join(working_dir,
- slide_dir,
- "analysis_outputs_7000_900",
- save_name))
- # %% [markdown]
- # ### Filter and combine anndata objects
- # %%
- def filter_anndata(slide_list, min_count, min_genes):
- """
- Filters anndata objects for each slide in slide_list based on transcript count and detected genes.
- Parameters:
- slide_list (list): List of slide names.
- min_count (int): Minimum transcript count threshold.
- min_genes (int): Minimum number of detected genes.
- Returns:
- None
- """
- for slide in slide_list:
- # Construct file paths
- input_path = f"{slide}/analysis_outputs_7000_900/{slide}.h5ad"
- output_dir = f"{slide}/analysis_outputs_7000_900/"
- output_path = os.path.join(output_dir, f"{slide}_filtered.h5ad")
- # Check if the file exists
- if not os.path.exists(input_path):
- print(f"File not found: {input_path}, skipping...")
- continue
- print(f"Processing {slide}...")
- # Load the AnnData object
- adata = sc.read_h5ad(input_path)
- # Remove blank genes
- non_blank_genes = [gene for gene in adata.var_names if "blank" not in gene.lower()]
- adata = adata[:, non_blank_genes].copy()
- # Compute number of detected genes and total transcript count
- adata.obs["num_genes"] = np.count_nonzero(adata.X, axis=1)
- adata.obs["num_transcripts"] = np.sum(adata.X, axis=1)
- # Filter cells based on criteria
- prev_len = adata.shape[0]
- adata = adata[
- (adata.obs["num_transcripts"] > min_count) &
- (adata.obs["num_genes"] > min_genes), :
- ].copy()
- print(f"Filtered from {prev_len} to {adata.shape[0]} cells.")
- # Ensure output directory exists
- os.makedirs(output_dir, exist_ok=True)
- # Save the processed AnnData object
- adata.write_h5ad(output_path)
- print(f"Saved filtered data to {output_path}\n")
- # %%
- # Filter anndata objects
- min_count = 20
- min_genes = 5
- filter_anndata(slide_list, min_count, min_genes)
- # %%
- # Combine filtered anndata objects
- idx = 0
- for slide_dir in slide_list:
- print(slide_dir)
- # read in and make anndata object
- load_dir = os.path.join(slide_dir,
- "analysis_outputs_7000_900",
- (slide_dir+'_filtered.h5ad'))
- adata = sc.read_h5ad(load_dir)
- if idx == 0:
- adata_combined = adata.copy()
- idx += 1
- continue
- adata_combined = ad.concat([adata_combined, adata]).copy()
- adata_combined.obs_names_make_unique()
- print(f'adata_combined shape is {adata_combined.shape}')
- idx += 1
- # %%
- # Saving count data
- adata_combined.layers["counts"] = adata_combined.X.copy()
- # %%
- # save anndata
- save_name = f"adata_objects/combined_filtered_anndata.h5ad"
- adata_combined.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # # All cells BBKNN integration
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(f"adata_objects/combined_filtered_anndata.h5ad")
- # %%
- # Add metadata columns
- # Slide + section and slide + stage
- adata.obs["slide_section"] = adata.obs["slide"].astype(str) + "_" + adata.obs["section"].astype(str)
- adata.obs["slide_stage"] = adata.obs["slide"].astype(str) + "_" + adata.obs["stage"].astype(str)
- # Auxillary channels normalized by volume
- columns_to_normalize = [
- "Gfap_raw", "Gfap_high_pass",
- "Apoe_raw", "Apoe_high_pass",
- "hSOD1G93A_raw", "hSOD1G93A_high_pass"
- ]
- for col in columns_to_normalize:
- norm_col = col + "_norm"
- adata.obs[norm_col] = adata.obs[col] / adata.obs["volume"]
- # %%
- # Remove low-quality tissue sections
- # List of slide_section values to remove
- to_remove = [
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R2C2',
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R2C1',
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C1',
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C1'
- ]
- # Keep cells NOT in the to_remove list
- adata = adata[~adata.obs['slide_section'].isin(to_remove)].copy()
- # %%
- # Normalization of data by volume
- adata.X = adata.X/np.array(adata.obs['volume'])[:,None]
- # Normalize total to 250
- sc.pp.normalize_total(adata, target_sum=250)
- # Log transform
- sc.pp.log1p(adata)
- # Z-score (need to do)
- sc.pp.scale(adata, max_value=10)
- # Run pca
- sc.tl.pca(adata)
- # %%
- # Neighbors and UMAP
- sc.external.pp.bbknn(adata, batch_key='VS') # use inplace of sc.pp.neighbors()
- sc.tl.umap(adata)
- # %%
- # Clustering
- sc.tl.leiden(adata,
- key_added="leiden_0.5",
- resolution=0.5,
- flavor="igraph",
- n_iterations=2)
- # %%
- sc.pl.umap(adata, color="leiden_0.5")
- # %% [markdown]
- # ## Fig. S1D
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- # 0) Set all font sizes to 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "figure.titlesize": 7
- })
- # 1) Create a high‐DPI figure (this sets the base resolution for rasterized artists)
- fig, ax = plt.subplots(dpi=300) # <-- 300 DPI here
- # 2) Draw UMAP off‐screen, no title, into that axes
- sc.pl.umap(
- adata,
- color="leiden_0.5",
- legend_loc='on data',
- legend_fontsize=7, # redundant now, but explicit
- title="",
- show=False,
- ax=ax
- )
- # 3) Rasterize only the scatter collections
- for coll in ax.collections:
- coll.set_rasterized(True)
- # 4) Fix the output size in inches
- fig.set_size_inches(2.15, 2.15)
- # 5) Save with an even higher DPI if you like (affects only raster parts in a vector file)
- fig.savefig(
- "figures/umap/UMAP_all_leiden.svg",
- format="svg",
- bbox_inches="tight",
- dpi=600 # <-- rasterized points will now be 600 DPI
- )
- # %%
- sc.pl.umap(adata, color="stage")
- # %%
- sc.pl.umap(adata, color="slide")
- # %% [markdown]
- # ## Fig. S1E
- # %%
- import matplotlib.pyplot as plt
- import seaborn as sns
- import os
- genes = ["Chat", "Prph", "Slc5a7", "Rbfox3", "Slc17a6", "Gad1", "Slc6a5", "Tnfrsf12a", "Atf3", "Cx3cr1", "Trem2",
- "Mog", "Stmn4", "Aqp4", "Aldh1l1", "Slc2a1", "Slc7a1", "Il34", "Anxa2", "Atp1a1", "Hsp90ab1", "Gria3",
- "Grik4"]
- # Define the desired cluster order
- desired_order = ["6", "10", "9", "12", "1", "2", "7", "8", "0", "11", "5", "4", "3"]
- adata_counts = adata[:, genes].copy()
- adata_counts.X = adata_counts.layers['counts']
- # Convert to a DataFrame
- counts_df = pd.DataFrame(adata_counts.X, index=adata_counts.obs_names, columns=adata_counts.var_names)
- # Normalize total counts to 250
- sc.pp.normalize_total(adata_counts, target_sum=250)
- # Add cluster information
- cluster_key = "leiden_0.5" # Update this if using a different cluster key
- counts_df["cluster"] = adata_counts.obs[cluster_key].values
- # Compute average expression per cluster
- cluster_avg = counts_df.groupby("cluster")[genes].mean().T # Transpose to have genes on y-axis
- # Convert row_min and row_max to NumPy arrays
- row_min = cluster_avg.min(axis=1).to_numpy()
- row_max = cluster_avg.max(axis=1).to_numpy()
- # Avoid division by zero by checking where min == max
- constant_rows = row_max == row_min
- # Perform row normalization (Min-Max Scaling)
- cluster_avg_norm = (cluster_avg - row_min[:, np.newaxis]) / (row_max - row_min)[:, np.newaxis]
- # Set constant rows to 0 (or another value like NaN if needed)
- cluster_avg_norm[constant_rows] = 0
- # Reorder the columns of the DataFrame
- cluster_avg_norm = cluster_avg_norm[desired_order]
- # Plot the heatmap
- fig, ax = plt.subplots(figsize=(5.6, 3.5))
- # draw heatmap with only two ticks at 0 and 1
- heatmap = sns.heatmap(
- cluster_avg_norm,
- cmap="viridis",
- annot=False,
- fmt=".2f",
- linewidths=0.5,
- ax=ax,
- cbar_kws={
- "ticks": [0.0, 1.0], # only these two positions
- "shrink": 0.9,
- "pad": 0.13
- }
- )
- # 1) Remove all tick-marks (but keep the labels)
- ax.tick_params(axis="both", which="both", length=0)
- # 2) Rotate the x-axis labels 90°
- ax.set_xticklabels(ax.get_xticklabels(), rotation=-90)
- # 3) Move the gene/row names (y-tick labels) to the right & make them horizontal
- ax.tick_params(axis="y", labelleft=False, labelright=True, left=False, right=True)
- ax.set_yticklabels(ax.get_yticklabels(), rotation=0)#, ha="right")
- # 4) Clear axis titles and labels
- ax.set_xlabel("")
- ax.set_ylabel("")
- ax.set_title("")
- # 5) Tweak the colorbar
- cbar = heatmap.collections[0].colorbar
- # a) remove the little tick‐lines
- cbar.ax.tick_params(length=0)
- # b) replace the two tick labels
- cbar.ax.set_yticklabels(["Min", "Max"])
- # c) remove any other labels (there won’t be any beyond your two ticks)
- # (no extra step needed since we only set two ticks)
- # d) add a vertically-oriented label “By Row” centered along the bar
- cbar.ax.set_ylabel(
- "By Row",
- rotation=-90, # text runs horizontally
- va="center", # center along the bar
- labelpad=-6 # push it out away from the bar
- )
- # 6) Save
- heatmap_path = os.path.join(figures_dir, "heatmap_all.svg")
- fig.savefig(heatmap_path, dpi=300, bbox_inches="tight", transparent=True)
- plt.show()
- # %%
- # Add cell class labels
- type_map = {'0': 'Astrocytes',
- '1': 'Microglia/Macrophages',
- '2': 'Oligodendrocytes',
- '3': 'Other',
- '4': 'Putative Ependymal Cells',
- '5': 'Putative Perivascular/Meningeal Cells',
- '6': 'Cholinergic Neurons',
- '7': 'Oligodendrocytes',
- '8': 'Oligodendrocytes',
- '9': 'Non-Cholinergic Interneurons',
- '10': 'Non-Cholinergic Interneurons',
- '11': 'Putative Vascular Cells',
- '12': 'Disease-Associated Interneurons'}
- cell_class = []
- leiden_cluster = adata.obs['leiden_0.5'].to_list()
- for l in leiden_cluster:
- cell_class.append(type_map[l])
- adata.obs['cell_class'] = cell_class
- # %%
- # Plot the new sub-clusters
- sc.pl.umap(adata, color="cell_class")
- # %% [markdown]
- # ## Fig. S1F
- # %%
- # 0) Set all font sizes to 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "figure.titlesize": 7
- })
- # 1) Create a high‐DPI figure (this sets the base resolution for rasterized artists)
- fig, ax = plt.subplots(dpi=300) # <-- 300 DPI here
- # your abbr dict
- abbr = {
- "Astrocytes": "Ast",
- "Cholinergic Neurons": "CholN",
- "Disease-Associated Interneurons": "DAI",
- "Microglia/Macrophages": "MG",
- "Non-Cholinergic Interneurons": "NCI",
- "Oligodendrocytes": "Oligo",
- "Other": "Other",
- "Putative Ependymal Cells": "Epen",
- "Putative Perivascular/Meningeal Cells": "PVM",
- "Putative Vascular Cells": "Vasc"
- }
- # 2) Draw UMAP off‐screen, no title, into that axes
- sc.pl.umap(
- adata,
- color="cell_class",
- legend_loc='on data',
- legend_fontsize=7,
- title="",
- show=False,
- ax=ax
- )
- # 2b) Replace each on‐data label with its abbr
- for txt in ax.texts:
- orig = txt.get_text()
- if orig in abbr:
- txt.set_text(abbr[orig])
- # now continue with your rasterization + sizing + saving…
- # 3) Rasterize only the scatter collections
- for coll in ax.collections:
- coll.set_rasterized(True)
- # 4) Fix the output size in inches
- fig.set_size_inches(2.15, 2.15)
- # 5) Save...
- fig.savefig(
- "figures/umap/UMAP_cell_class.svg",
- format="svg",
- bbox_inches="tight",
- dpi=600
- )
- # %%
- # save anndata
- save_name = f"adata_objects/all_anndata.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata.h5ad"))
- # %% [markdown]
- # ## Subcluster cholinergic neurons
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata.h5ad"))
- # %%
- # Select cholinergic neurons
- adata_chol = adata[adata.obs["cell_class"] == "Cholinergic Neurons"].copy()
- # %%
- # Neighbors and UMAP
- sc.external.pp.bbknn(adata_chol, batch_key='VS') # use inplace of sc.pp.neighbors()
- sc.tl.umap(adata_chol,
- min_dist=0.1,
- spread=1.0)
- # %%
- # Clustering
- sc.tl.leiden(adata_chol,
- key_added="chol_leiden_1",
- resolution=1,
- flavor="igraph",
- n_iterations=2)
- # %%
- # Plot the sub-clusters
- sc.pl.umap(adata_chol, color="chol_leiden_1")
- # %%
- # Remove Mog+ oligodendrocyte doublets
- sc.pl.umap(adata_chol, color=["Mog"])
- # %% [markdown]
- # #### Get barcodes for peri-motor neurons oligodendrocytes
- # %%
- # Assess volume distribution of exisiting Oligodendrocytes
- import matplotlib.pyplot as plt
- import seaborn as sns
- # Subset to Oligodendrocytes
- adata_oligo = adata[adata.obs['cell_class'] == 'Oligodendrocytes']
- # Extract the volume values
- vals = adata_oligo.obs['volume'].dropna()
- # Plot
- fig, ax = plt.subplots(figsize=(4, 6))
- sns.violinplot(y=vals, inner="box", ax=ax)
- ax.set_ylabel('volume')
- ax.set_title('Oligodendrocytes — Volume Distribution')
- plt.tight_layout()
- plt.show()
- # Compute quartiles
- q1 = vals.quantile(0.25)
- median = vals.quantile(0.50)
- q3 = vals.quantile(0.75)
- iqr = q3 - q1
- # Compute whisker positions (1.5 × IQR rule)
- lower_whisker = vals[vals >= (q1 - 1.5 * iqr)].min()
- upper_whisker = vals[vals <= (q3 + 1.5 * iqr)].max()
- # Print them out
- print(f"Count: {len(vals)}")
- print(f"Q1 (25th pct): {q1:.3f}")
- print(f"Median (50th): {median:.3f}")
- print(f"Q3 (75th pct): {q3:.3f}")
- print(f"IQR: {iqr:.3f}")
- print(f"Lower whisker: {lower_whisker:.3f}")
- print(f"Upper whisker: {upper_whisker:.3f}")
- # %%
- # Assess volume distribution of putative Oligodendrocytes/Oligo-MN doublets
- import matplotlib.pyplot as plt
- import seaborn as sns
- # Subset to putative Oligodendrocytes/Oligo-MN doublets
- adata_periMN_oligo = adata_chol[adata_chol.obs["chol_leiden_1"].isin(["12", "13", "14"])]
- # Extract the volume values
- vals = adata_periMN_oligo.obs['volume'].dropna()
- # Plot
- fig, ax = plt.subplots(figsize=(4, 6))
- sns.violinplot(y=vals, inner="box", ax=ax)
- ax.set_ylabel('volume')
- ax.set_title('Volume Distribution')
- plt.tight_layout()
- plt.show()
- # Compute quartiles
- q1 = vals.quantile(0.25)
- median = vals.quantile(0.50)
- q3 = vals.quantile(0.75)
- iqr = q3 - q1
- # Compute whisker positions (1.5 × IQR rule)
- lower_whisker = vals[vals >= (q1 - 1.5 * iqr)].min()
- upper_whisker = vals[vals <= (q3 + 1.5 * iqr)].max()
- # Print them out
- print(f"Count: {len(vals)}")
- print(f"Q1 (25th pct): {q1:.3f}")
- print(f"Median (50th): {median:.3f}")
- print(f"Q3 (75th pct): {q3:.3f}")
- print(f"IQR: {iqr:.3f}")
- print(f"Lower whisker: {lower_whisker:.3f}")
- print(f"Upper whisker: {upper_whisker:.3f}")
- # %%
- # Save metadata for cells within the lower and upper whiskers
- # Build the boolean mask
- mask = (
- (adata_periMN_oligo.obs['volume'] >= lower_whisker) &
- (adata_periMN_oligo.obs['volume'] <= upper_whisker)
- )
- # Get the filtered obs DataFrame
- filtered_obs = adata_periMN_oligo.obs.loc[mask]
- # Change cell_class to Oligodendrocytes
- filtered_obs['cell_class'] = 'Oligodendrocytes'
- # Save with the index as a column
- filtered_obs.to_csv(os.path.join(working_dir, 'periMN_oligo_metadata.csv'), index=True)
- # %% [markdown]
- # #### Back to cholinergic neurons
- # %%
- clusters_to_remove = ["12", "13", "14"]
- adata_chol_filtered = adata_chol[~adata_chol.obs["chol_leiden_1"].isin(clusters_to_remove)].copy()
- sc.pl.umap(adata_chol_filtered, color="chol_leiden_1")
- # %%
- # Re-run BBKNN, UMAP, and clustering
- # Neighbors and UMAP
- sc.external.pp.bbknn(adata_chol_filtered, batch_key='VS') # use inplace of sc.pp.neighbors()
- sc.tl.umap(adata_chol_filtered,
- min_dist=0.5,
- spread=1.0)
- # Clustering
- sc.tl.leiden(adata_chol_filtered,
- key_added="chol_leiden_1.2",
- resolution=1.2,
- flavor="igraph",
- n_iterations=2)
- # %%
- # Plot the new sub-clusters
- sc.pl.umap(adata_chol_filtered, color="chol_leiden_1.2")
- sc.pl.umap(adata_chol_filtered, color="stage")
- sc.pl.umap(adata_chol_filtered, color="slide")
- # %% [markdown]
- # ## Fig. S2F
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- # 0) Set all font sizes to 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "figure.titlesize": 7
- })
- # 1) Create a high‐DPI figure (this sets the base resolution for rasterized artists)
- fig, ax = plt.subplots(dpi=300) # <-- 300 DPI here
- # 2) Draw UMAP off‐screen, no title, into that axes
- sc.pl.umap(
- adata_chol_filtered,
- color="chol_leiden_1.2",
- size=5,
- legend_loc='on data',
- legend_fontsize=7, # redundant now, but explicit
- title="",
- show=False,
- ax=ax
- )
- # 3) Rasterize only the scatter collections
- for coll in ax.collections:
- coll.set_rasterized(True)
- # 4) Fix the output size in inches
- fig.set_size_inches(2.15, 2.15)
- # 5) Save with an even higher DPI if you like (affects only raster parts in a vector file)
- fig.savefig(
- "figures/umap/UMAP_chol_leiden.svg",
- format="svg",
- bbox_inches="tight",
- dpi=600 # <-- rasterized points will now be 600 DPI
- )
- # %% [markdown]
- # ## Fig. S2G
- # %%
- import matplotlib.pyplot as plt
- import seaborn as sns
- import os
- genes = ["Bcl6", "Stk32a", "Vipr2", "Npas1", "Creb5", "Plch1", "Gap43", "Chat", "Slc5a7", "Prph"]
- # Define the desired cluster order
- desired_order = ['12', '4', '10', '18', '2', '7', '8', '15', '14', '19',
- '11', '17', '5', '3', '6', '16', '9', '13', '1', '0']
- adata_chol_filtered_counts = adata_chol_filtered[:, genes].copy()
- adata_chol_filtered_counts.X = adata_chol_filtered_counts.layers['counts']
- # Convert to a DataFrame
- counts_df = pd.DataFrame(adata_chol_filtered_counts.X, index=adata_chol_filtered_counts.obs_names, columns=adata_chol_filtered_counts.var_names)
- # Normalize total counts to 250
- sc.pp.normalize_total(adata_chol_filtered_counts, target_sum=250)
- # Add cluster information
- cluster_key = "chol_leiden_1.2" # Update this if using a different cluster key
- counts_df["cluster"] = adata_chol_filtered_counts.obs[cluster_key].values
- # Compute average expression per cluster
- cluster_avg = counts_df.groupby("cluster")[genes].mean().T # Transpose to have genes on y-axis
- # Convert row_min and row_max to NumPy arrays
- row_min = cluster_avg.min(axis=1).to_numpy()
- row_max = cluster_avg.max(axis=1).to_numpy()
- # Avoid division by zero by checking where min == max
- constant_rows = row_max == row_min
- # Perform row normalization (Min-Max Scaling)
- cluster_avg_norm = (cluster_avg - row_min[:, np.newaxis]) / (row_max - row_min)[:, np.newaxis]
- # Set constant rows to 0 (or another value like NaN if needed)
- cluster_avg_norm[constant_rows] = 0
- # Reorder the columns of the DataFrame
- cluster_avg_norm = cluster_avg_norm[desired_order]
- # Plot the heatmap
- fig, ax = plt.subplots(figsize=(5.7, 2.4))
- # draw heatmap with only two ticks at 0 and 1
- heatmap = sns.heatmap(
- cluster_avg_norm,
- cmap="viridis",
- annot=False,
- fmt=".2f",
- linewidths=0.5,
- ax=ax,
- cbar_kws={
- "ticks": [0.0, 1.0], # only these two positions
- "shrink": 0.9,
- "pad": 0.13
- }
- )
- # 1) Remove all tick-marks (but keep the labels)
- ax.tick_params(axis="both", which="both", length=0)
- # 2) Rotate the x-axis labels 90°
- ax.set_xticklabels(ax.get_xticklabels(), rotation=-90)
- # 3) Move the gene/row names (y-tick labels) to the right & make them horizontal
- ax.tick_params(axis="y", labelleft=False, labelright=True, left=False, right=True)
- ax.set_yticklabels(ax.get_yticklabels(), rotation=0)#, ha="right")
- # 4) Clear axis titles and labels
- ax.set_xlabel("")
- ax.set_ylabel("")
- ax.set_title("")
- # 5) Tweak the colorbar
- cbar = heatmap.collections[0].colorbar
- # a) remove the little tick‐lines
- cbar.ax.tick_params(length=0)
- # b) replace the two tick labels
- cbar.ax.set_yticklabels(["Min", "Max"])
- # c) remove any other labels (there won’t be any beyond your two ticks)
- # (no extra step needed since we only set two ticks)
- # d) add a vertically-oriented label “By Row” centered along the bar
- cbar.ax.set_ylabel(
- "By Row",
- rotation=-90, # text runs horizontally
- va="center", # center along the bar
- labelpad=-6 # push it out away from the bar
- )
- # 6) Save
- heatmap_path = os.path.join(figures_dir, "heatmap_chol.svg")
- fig.savefig(heatmap_path, dpi=300, bbox_inches="tight", transparent=True)
- plt.show()
- # %%
- # Add cholinergic type labels
- type_map = {'0': 'Cholinergic Interneurons',
- '1': 'Visceral MNs',
- '2': 'Alpha MNs',
- '3': 'Alpha MNs',
- '4': 'Alpha MNs',
- '5': 'Alpha MNs',
- '6': 'Alpha MNs',
- '7': 'Alpha MNs',
- '8': 'Alpha MNs',
- '9': 'Gamma MNs',
- '10': 'Alpha MNs',
- '11': 'Alpha MNs',
- '12': 'Alpha MNs',
- '13': 'Gamma* MNs',
- '14': 'Alpha MNs',
- '15': 'Alpha MNs',
- '16': 'Alpha MNs',
- '17': 'Alpha MNs',
- '18': 'Alpha MNs',
- '19': 'Alpha MNs'}
- cholinergic_type = []
- leiden_cluster = adata_chol_filtered.obs['chol_leiden_1.2'].to_list()
- for l in leiden_cluster:
- cholinergic_type.append(type_map[l])
- adata_chol_filtered.obs['cholinergic_type'] = cholinergic_type
- # %%
- # Plot the annotated UMAP
- sc.pl.umap(adata_chol_filtered, color="cholinergic_type", save="/UMAP_cholinergic_type.png")
- # %% [markdown]
- # ## Fig. S2H
- # %%
- # 0) Set all font sizes to 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "figure.titlesize": 7
- })
- # 1) Create a high‐DPI figure (this sets the base resolution for rasterized artists)
- fig, ax = plt.subplots(dpi=300) # <-- 300 DPI here
- # your abbr dict
- abbr = {
- "Alpha MNs": "Alpha MNs",
- "Gamma* MNs": "Gamma* MNs",
- "Gamma MNs": "Gamma MNs",
- "Visceral MNs": "Visceral \nMNs",
- "Cholinergic Interneurons": "Chol. \nInt."
- }
- # 2) Draw UMAP off‐screen, no title, into that axes
- sc.pl.umap(
- adata_chol_filtered,
- size=5,
- color="cholinergic_type",
- legend_loc='on data',
- legend_fontsize=7,
- title="",
- show=False,
- ax=ax
- )
- # 2b) Replace each on‐data label with its abbr
- for txt in ax.texts:
- orig = txt.get_text()
- if orig in abbr:
- txt.set_text(abbr[orig])
- # now continue with your rasterization + sizing + saving…
- # 3) Rasterize only the scatter collections
- for coll in ax.collections:
- coll.set_rasterized(True)
- # 4) Fix the output size in inches
- fig.set_size_inches(2.15, 2.15)
- # 5) Save...
- fig.savefig(
- "figures/umap/UMAP_chol_type.svg",
- format="svg",
- bbox_inches="tight",
- dpi=600
- )
- # %%
- sc.pl.violin(adata_chol_filtered, keys=['volume'], groupby='cholinergic_type', rotation=90)
- # %%
- # Build the mask for cells to keep
- keep_mask = adata_chol_filtered.obs['volume'] >= 1000
- # Make a new AnnData
- adata_chol_filtered_highvol = adata_chol_filtered[keep_mask, :].copy()
- # %%
- # Plot the new sub-clusters
- sc.pl.umap(adata_chol_filtered_highvol, color="chol_leiden_1.2")
- # %%
- # Plot the annotated UMAP
- sc.pl.umap(adata_chol_filtered_highvol, color="cholinergic_type")
- # %%
- sc.pl.violin(adata_chol_filtered_highvol, keys=['volume'], groupby='cholinergic_type', rotation=90)
- # %%
- # Save anndata
- save_name = f"adata_objects/cholinergic_anndata.h5ad"
- adata_chol_filtered.write_h5ad(os.path.join(working_dir, save_name))
- save_name = f"adata_objects/cholinergic_anndata_highvol.h5ad"
- adata_chol_filtered_highvol.write_h5ad(os.path.join(working_dir, save_name))
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/cholinergic_anndata.h5ad"))
- # %%
- adata_chol_filtered = adata
- # %% [markdown]
- # ## Subcluster alpha MNs
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/cholinergic_anndata_highvol.h5ad"))
- # %%
- # Select alpha MNs
- adata_alpha = adata[adata.obs["cholinergic_type"] == "Alpha MNs"].copy()
- # %%
- # Neighbors and UMAP
- sc.external.pp.bbknn(adata_alpha, batch_key='VS') # use inplace of sc.pp.neighbors()
- sc.tl.umap(adata_alpha,
- min_dist=0.5,
- spread=1.0)
- # %%
- # Clustering
- sc.tl.leiden(adata_alpha,
- key_added="alpha_leiden_0.45",
- resolution=0.45,
- flavor="igraph",
- n_iterations=2)
- # %%
- # Plot the sub-clusters
- sc.pl.umap(adata_alpha, color="alpha_leiden_0.45", save="/UMAP_alpha_leiden.png")
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- # 0) Set all font sizes to 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "figure.titlesize": 7
- })
- # 1) Create a high‐DPI figure (this sets the base resolution for rasterized artists)
- fig, ax = plt.subplots(dpi=300) # <-- 300 DPI here
- # 2) Draw UMAP off‐screen, no title, into that axes
- sc.pl.umap(
- adata_alpha,
- color="alpha_leiden_0.45",
- size=10,
- legend_loc='on data',
- legend_fontsize=7, # redundant now, but explicit
- title="",
- show=False,
- ax=ax
- )
- # 3) Rasterize only the scatter collections
- for coll in ax.collections:
- coll.set_rasterized(True)
- # 4) Fix the output size in inches
- fig.set_size_inches(2.15, 2.15)
- # 5) Save with an even higher DPI if you like (affects only raster parts in a vector file)
- fig.savefig(
- "figures/umap/UMAP_alpha_leiden.svg",
- format="svg",
- bbox_inches="tight",
- dpi=600 # <-- rasterized points will now be 600 DPI
- )
- # %%
- sc.pl.violin(adata_alpha, keys=['volume'], groupby='alpha_leiden_0.45')
- # %%
- # Plot the sub-clusters
- sc.pl.umap(adata_alpha, color="stage")
- sc.pl.umap(adata_alpha, color="slide")
- # %%
- sc.pl.umap(adata_alpha, color=["Prkcd", "Chodl", "Atf3", "Gap43"])
- # %%
- # Save anndata
- save_name = f"adata_objects/alpha_anndata.h5ad"
- adata_alpha.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ## Alpha MN UMAP labeled by stage
- # %%
- adata_alpha = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata.h5ad"))
- # %% [markdown]
- # ## Fig. S3I
- # %%
- import scanpy as sc
- import matplotlib.pyplot as plt
- import pandas as pd
- # ---- Enforce stage order ----
- desired_order = ["Control", "Early", "Mid", "End"]
- adata_alpha.obs["stage"] = pd.Categorical(
- adata_alpha.obs["stage"],
- categories=desired_order,
- ordered=True
- )
- # Global font size = 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 6
- })
- # Create UMAP and return figure
- fig = sc.pl.umap(
- adata_alpha,
- color="stage",
- title="",
- size=5,
- show=False,
- return_fig=True
- )
- # Set exact figure size AFTER creation
- fig.set_size_inches(1.85, 1.3306)
- plt.tight_layout()
- plt.show()
- fig.savefig(
- "figures/umap/adata_alpha_stage_umap.svg",
- format="svg",
- bbox_inches="tight"
- )
- # %% [markdown]
- # ## Differentially expressed genes in alpha MNs
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata_label_transfer.h5ad"))
- # %%
- # Get all gene names from the AnnData object
- gene_names = adata.var_names
- # Convert to DataFrame with a column name
- genes_df = pd.DataFrame(gene_names, columns=["Gene"])
- # Save to CSV (no row index, just gene names)
- genes_df.to_csv("MERFISH_genes.csv", index=False)
- # %%
- # Read the CSV file of differentially expressed genes with disease in Alpha MNs from snRNA-seq
- alpha_snrna_end = pd.read_csv("/home/users/ogautier/oak/Shared/SOD1_Paper/RNA/DESeq2/Cholinergic_Type/Alpha MNs_sod.end_vs_ctl.csv", index_col=0)
- # %%
- # Filter alpha_snrna_end for padj < 0.01 & gene/row name in adata
- alpha_snrna_end_sig = alpha_snrna_end[
- (alpha_snrna_end["padj"] < 0.01) & alpha_snrna_end.index.isin(adata.var_names)
- ]
- # %%
- # Get upregulated and downregulated genes
- alpha_snrna_upregulated = alpha_snrna_end_sig[
- (alpha_snrna_end_sig["log2FoldChange"] > 0)
- ]
- alpha_snrna_downregulated = alpha_snrna_end_sig[
- (alpha_snrna_end_sig["log2FoldChange"] < 0)
- ]
- # %%
- alpha_snrna_end_sig
- # %%
- alpha_snrna_upregulated
- # %%
- alpha_snrna_downregulated
- # %%
- # Filter adata to keep genes in alpha_snrna_end_sig
- mask = adata.var_names.isin(alpha_snrna_end_sig.index)
- adata_filtered = adata[:, mask]
- # %%
- adata_filtered = adata_filtered[:, :].copy()
- adata_filtered.X = adata_filtered.layers['counts']
- # Normalization of data by volume
- adata_filtered.X = adata_filtered.X/np.array(adata_filtered.obs['volume'])[:,None]
- # Normalize total to 250
- sc.pp.normalize_total(adata_filtered, target_sum=250)
- # Log transform
- sc.pp.log1p(adata_filtered)
- # %%
- import pandas as pd
- import numpy as np
- from scipy.stats import mannwhitneyu
- from statsmodels.stats.multitest import multipletests
- # Initialize lists to store the results
- genes = []
- p_values_comparison1 = []
- log2foldchange_comparison1 = []
- p_values_comparison2 = []
- log2foldchange_comparison2 = []
- # Iterate over all genes
- for gene in adata_filtered.var_names:
- # Extract expression data (dense format) for the gene
- control_expr = np.asarray(adata_filtered[adata_filtered.obs["stage"] == "Control", gene].X).flatten()
- mid_expr = np.asarray(adata_filtered[adata_filtered.obs["stage"] == "Mid", gene].X).flatten()
- end_expr = np.asarray(adata_filtered[adata_filtered.obs["stage"] == "End", gene].X).flatten()
- # Control vs Mid
- p_mid = float(mannwhitneyu(control_expr, mid_expr, alternative='two-sided')[1])
- fc_mid = np.nanmean(np.exp(mid_expr) - 1) / np.nanmean(np.exp(control_expr) - 1)
- log2_fc_mid = np.log2(fc_mid)
- # Control vs End
- p_end = float(mannwhitneyu(control_expr, end_expr, alternative='two-sided')[1])
- fc_end = np.nanmean(np.exp(end_expr) - 1) / np.nanmean(np.exp(control_expr) - 1)
- log2_fc_end = np.log2(fc_end)
- # Store results
- genes.append(gene)
- p_values_comparison1.append(p_mid)
- log2foldchange_comparison1.append(log2_fc_mid)
- p_values_comparison2.append(p_end)
- log2foldchange_comparison2.append(log2_fc_end)
- # Create a DataFrame with the raw p-values and log2 fold changes
- results_df = pd.DataFrame({
- "Gene": genes,
- "p-value_Control_vs_Mid": p_values_comparison1,
- "log2FoldChange_Control_vs_Mid": log2foldchange_comparison1,
- "p-value_Control_vs_End": p_values_comparison2,
- "log2FoldChange_Control_vs_End": log2foldchange_comparison2
- })
- # Adjust p-values using Benjamini-Hochberg (FDR) correction
- # For Control vs Mid:
- adj_mid = multipletests(results_df["p-value_Control_vs_Mid"], method="fdr_bh")
- results_df["adj-p-value_Control_vs_Mid"] = adj_mid[1]
- # For Control vs End:
- adj_end = multipletests(results_df["p-value_Control_vs_End"], method="fdr_bh")
- results_df["adj-p-value_Control_vs_End"] = adj_end[1]
- # Sort by adjusted p-value for one of the comparisons:
- results_df = results_df.sort_values("adj-p-value_Control_vs_End")
- # Specify the desired column order
- desired_order = [
- "Gene",
- "log2FoldChange_Control_vs_Mid",
- "p-value_Control_vs_Mid",
- "adj-p-value_Control_vs_Mid",
- "log2FoldChange_Control_vs_End",
- "p-value_Control_vs_End",
- "adj-p-value_Control_vs_End"
- ]
- # Reorder the DataFrame
- results_df = results_df[desired_order]
- # Display the reordered DataFrame
- results_df
- # %%
- # Save to CSV in the figures directory
- output_path = os.path.join(figures_dir, "alpha_differential_gene_expression_results.csv")
- results_df.to_csv(output_path, index=False)
- # %% [markdown]
- # ## Fig. S3J
- # %%
- import matplotlib.pyplot as plt
- from matplotlib_venn import venn3
- # 1) Prepare gene sets as before
- # Upregulated Genes (adj‐p < 0.01, log2FC > 0)
- mid_upregulated = results_df[
- (results_df["adj-p-value_Control_vs_Mid"] < 0.01) &
- (results_df["log2FoldChange_Control_vs_Mid"] > 0)
- ]
- end_upregulated = results_df[
- (results_df["adj-p-value_Control_vs_End"] < 0.01) &
- (results_df["log2FoldChange_Control_vs_End"] > 0)
- ]
- # Downregulated Genes (adj‐p < 0.01, log2FC < 0)
- mid_downregulated = results_df[
- (results_df["adj-p-value_Control_vs_Mid"] < 0.01) &
- (results_df["log2FoldChange_Control_vs_Mid"] < 0)
- ]
- end_downregulated = results_df[
- (results_df["adj-p-value_Control_vs_End"] < 0.01) &
- (results_df["log2FoldChange_Control_vs_End"] < 0)
- ]
- mid_up_set = set(mid_upregulated["Gene"])
- end_up_set = set(end_upregulated["Gene"])
- mid_down_set = set(mid_downregulated["Gene"])
- end_down_set = set(end_downregulated["Gene"])
- snrna_up_set = set(alpha_snrna_upregulated.index)
- snrna_down_set = set(alpha_snrna_downregulated.index)
- # 2) Create one figure with 2 rows, 1 column
- # Size specified in inches: width=1.075", height=2.15"
- fig, axes = plt.subplots(
- nrows=2, ncols=1,
- figsize=(1.075, 1.95),
- dpi=300
- )
- # 3) Top subplot: Upregulated venn3
- v1 = venn3(
- subsets=[mid_up_set, end_up_set, snrna_up_set],
- set_labels=("MERFISH (Mid)", "MERFISH (End)", "snRNA-seq"),
- ax=axes[0]
- )
- axes[0].set_title("Upregulated Genes", fontsize=7)
- # Adjust all label fonts in the first Venn
- for text in v1.set_labels:
- text.set_fontsize(7)
- for text in v1.subset_labels:
- if text: # some subset regions may be empty (None)
- text.set_fontsize(7)
- # 4) Bottom subplot: Downregulated venn3
- v2 = venn3(
- subsets=[mid_down_set, end_down_set, snrna_down_set],
- set_labels=("MERFISH (Mid)", "MERFISH (End)", "snRNA-seq"),
- ax=axes[1]
- )
- axes[1].set_title("Downregulated Genes", fontsize=7)
- # Adjust all label fonts in the second Venn
- for text in v2.set_labels:
- text.set_fontsize(7)
- for text in v2.subset_labels:
- if text:
- text.set_fontsize(7)
- # 5) Manually adjust margins so no labels get cut off
- fig.subplots_adjust(
- left=0.12,
- right=0.98,
- top=0.98,
- bottom=0.02,
- hspace=0.3 # space between the two plots
- )
- # 6) Save the combined figure
- fig.savefig("figures/alpha_venn_up_down.svg", dpi=300)
- # 7) Display
- plt.show()
- # %% [markdown]
- # ## Apoptosis/Stress genes
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata_label_transfer.h5ad"))
- # %%
- adata_violin = adata[:, :].copy()
- adata_violin.X = adata_violin.layers['counts']
- # Normalization of data by volume
- adata_violin.X = adata_violin.X/np.array(adata_violin.obs['volume'])[:,None]
- # Normalize total to 250
- sc.pp.normalize_total(adata_violin, target_sum=250)
- # Log transform
- sc.pp.log1p(adata_violin)
- # %% [markdown]
- # ## Fig. S3K
- # %%
- import os
- import numpy as np
- import matplotlib as mpl
- import matplotlib.pyplot as plt
- import scanpy as sc
- from matplotlib.collections import PolyCollection
- # global font size
- mpl.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "axes.linewidth": 0.5,
- "xtick.major.width": 0.5,
- "ytick.major.width": 0.5,
- "xtick.minor.width": 0.5,
- "ytick.minor.width": 0.5
- })
- # genes and stage order
- keys = ['Hrk', 'Ddit3', 'Trib3']
- order = ['Control', 'Early', 'Mid', 'End']
- # prepare figure with independent y-axes
- fig, axes = plt.subplots(
- nrows=1, ncols=len(keys),
- figsize=(2, 0.9),
- sharey=False
- )
- # draw one violins
- # then compute and plot a single median dot
- for idx, (ax, gene) in enumerate(zip(axes, keys)):
- sc.pl.violin(
- adata_violin,
- keys=[gene],
- groupby='stage',
- order=order,
- rotation=90,
- multi_panel=False,
- stripplot=False, # no individual points
- jitter=False, # no jitter
- inner=None, # no quartile/median lines
- ax=ax,
- show=False
- )
- # thin the violin outlines
- for coll in ax.collections:
- if isinstance(coll, PolyCollection):
- coll.set_linewidth(0.75)
- # thin the axes spines
- for spine in ax.spines.values():
- spine.set_linewidth(0.5)
- # title & axis labels
- ax.set_title(gene)
- if idx == 0:
- ax.set_ylabel("Expression")
- else:
- ax.set_ylabel("")
- ax.set_xlabel("")
- ax.tick_params(axis='y', labelleft=True)
- ax.grid(False)
- # get expression values and overlay medians
- expr = adata_violin.obs_vector(gene)
- for xi, stage in enumerate(order):
- mask = adata_violin.obs['stage'] == stage
- med = np.median(expr[mask.values])
- ax.scatter(xi, med,
- color='black',
- marker='o',
- s=1,
- zorder=10)
- # Reduce horizontal space between panels
- fig.subplots_adjust(wspace=0.8) # try 0.05–0.2
- plt.show()
- fig.savefig(
- "figures/violin/apoptosis_stress_genes_violin_stage.svg",
- format="svg",
- bbox_inches="tight",
- pad_inches=0.02 # reduces outer whitespace
- )
- # %% [markdown]
- # ## Alpha MN label transfer from snRNA-seq data
- # %% [markdown]
- # #### Create h5ad file for all data
- # %%
- # read the matrix
- adata = sc.read_10x_mtx(
- "/oak/stanford/groups/agitler/Shared/SOD1_Paper/RNA/files_to_make_h5ad/final_obj/",
- var_names='gene_symbols',
- make_unique=True
- )
- # %%
- # read the metadata
- md = pd.read_csv(
- "/oak/stanford/groups/agitler/Shared/SOD1_Paper/RNA/files_to_make_h5ad/final_obj/cell_metadata.tsv",
- sep="\t",
- index_col=0
- )
- # %%
- # make sure the index matches adata.obs_names
- md = md.reindex(adata.obs_names)
- # %%
- # assign
- adata.obs = md
- # %%
- # Save AnnData
- save_name = f"adata_objects/all_snrna_data.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # #### Create h5ad file from alpha MN snRNA-seq data
- # %%
- # read the matrix
- adata = sc.read_10x_mtx(
- "/oak/stanford/groups/agitler/Shared/SOD1_Paper/RNA/files_to_make_h5ad/alpha_label_transfer_50/",
- var_names='gene_symbols',
- make_unique=True
- )
- # %%
- # read the metadata
- md = pd.read_csv(
- "/oak/stanford/groups/agitler/Shared/SOD1_Paper/RNA/files_to_make_h5ad/alpha_label_transfer_50/cell_metadata.tsv",
- sep="\t",
- index_col=0
- )
- # %%
- # make sure the index matches adata.obs_names
- md = md.reindex(adata.obs_names)
- # %%
- # assign
- adata.obs = md
- # %%
- # Save AnnData
- save_name = f"adata_objects/alpha_snrna_data.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # #### Label transfer
- # %%
- # Load MERFISH and snRNA-seq datasets
- merfish = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata.h5ad"))
- snrna = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_snrna_data.h5ad"))
- # %%
- merfish.X = merfish.layers["counts"]
- # %%
- # Find common genes between MERFISH and snRNA-seq
- common_genes = list(set(merfish.var_names) & set(snrna.var_names))
- # %%
- # Filter both datasets to keep only the common genes
- merfish = merfish[:, common_genes]
- snrna = snrna[:, common_genes]
- # %%
- merfish.obs["batch"] = "MERFISH"
- snrna.obs["batch"] = "snRNA"
- # Concatenate datasets
- adata_combined = ad.concat([merfish, snrna], join="outer", label="batch", keys=["MERFISH", "snRNA"])
- # %%
- scvi.model.SCVI.setup_anndata(adata_combined, batch_key="batch")
- model = scvi.model.SCVI(adata_combined)
- # %%
- model.train()
- # %%
- adata_combined.obs['predicted.id'] = adata_combined.obs['predicted.id'].cat.add_categories('Unknown')
- adata_combined.obs = adata_combined.obs.fillna(value = {'predicted.id': 'Unknown'})
- # %%
- model2 = scvi.model.SCANVI.from_scvi_model(model, adata = adata_combined, unlabeled_category = 'Unknown',
- labels_key = 'predicted.id')
- # %%
- model2.train(max_epochs = 400)
- # %%
- adata_combined.obs['predicted'] = model2.predict(adata_combined)
- # %% [markdown]
- # #### Add the predicted labels to the alpha MN anndata object
- # %%
- cell_mapper = dict(zip(adata_combined.obs.index, adata_combined.obs.predicted))
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata.h5ad"))
- # %%
- adata.obs['predicted.id'] = adata.obs.index.map(cell_mapper)
- # %%
- # Plot the transferred labels
- sc.pl.umap(adata, color="predicted.id")
- # %%
- sc.pl.violin(
- adata,
- keys=['Prkcd', 'Chodl', 'Atf3', 'Gap43'],
- groupby='predicted.id',
- order=["Slow-Firing", "Intermediate", "Fast-Firing", "Early DAMN", "Late DAMN"],
- rotation=90
- )
- # %%
- # Save anndata
- save_name = f"adata_objects/alpha_anndata_label_transfer.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ## Alpha MN violin plots
- # %%
- # Read in AnnData
- adata_alpha = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata_label_transfer.h5ad"))
- # %%
- adata_violin = adata_alpha[:, :].copy()
- adata_violin.X = adata_violin.layers['counts']
- # Normalization of data by volume
- adata_violin.X = adata_violin.X/np.array(adata_violin.obs['volume'])[:,None]
- # Normalize total to 250
- sc.pp.normalize_total(adata_violin, target_sum=250)
- # Log transform
- sc.pp.log1p(adata_violin)
- # %% [markdown]
- # ## Fig. 2K
- # %%
- import os
- import numpy as np
- import matplotlib as mpl
- import matplotlib.pyplot as plt
- import scanpy as sc
- from matplotlib.collections import PolyCollection
- # ensure output folder exists
- os.makedirs("figures/violin", exist_ok=True)
- # global font size and thin axis spines
- mpl.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "axes.linewidth": 0.5,
- "xtick.major.width": 0.5,
- "ytick.major.width": 0.5,
- "xtick.minor.width": 0.5,
- "ytick.minor.width": 0.5
- })
- # genes and stage order
- keys = [
- 'Tenm1','Kcng1','Kcnj14','Nalcn','Grm8','Grik4','Gria3',
- 'Gabra5','Nefh','Ina','Slc7a11','Slc7a1','Sqstm1','Taf15'
- ]
- order = ['Control','Early','Mid','End']
- # figure: 2×7 panels
- ncol = len(keys)//2
- fig, axes = plt.subplots(2, ncol, figsize=(6.5936, 2.5296), sharey=False)
- axes_flat = axes.flatten()
- for idx, (ax, gene) in enumerate(zip(axes_flat, keys)):
- sc.pl.violin(
- adata_violin,
- keys=[gene],
- groupby='stage',
- order=order,
- rotation=90,
- multi_panel=False,
- stripplot=False,
- jitter=False,
- inner=None,
- ax=ax,
- show=False
- )
- # thin the violin outlines
- for coll in ax.collections:
- if isinstance(coll, PolyCollection):
- coll.set_linewidth(0.75)
- # thin the axes spines
- for spine in ax.spines.values():
- spine.set_linewidth(0.5)
- # title & labels
- ax.set_title(gene)
- if idx % ncol == 0:
- ax.set_ylabel("Expression")
- else:
- ax.set_ylabel("")
- ax.set_xlabel("")
- ax.tick_params(axis='y', labelleft=True)
- # overlay median dot
- expr = adata_violin.obs_vector(gene)
- for xi, stage in enumerate(order):
- med = np.median(expr[adata_violin.obs['stage'] == stage])
- ax.scatter(xi, med, color='black', s=1, zorder=10)
- for idx, ax in enumerate(axes_flat):
- # first ncol are top row
- if idx < ncol:
- ax.tick_params(axis='x', labelbottom=False)
- else:
- ax.tick_params(axis='x', labelbottom=True)
- plt.tight_layout()
- plt.show()
- fig.savefig(
- "figures/violin/disease_genes_violin.svg",
- format="svg",
- bbox_inches="tight"
- )
- # %% [markdown]
- # ## Fig. 4E
- # %%
- import os
- import numpy as np
- import matplotlib as mpl
- import matplotlib.pyplot as plt
- import scanpy as sc
- from matplotlib.collections import PolyCollection
- # global font size
- mpl.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "axes.linewidth": 0.5,
- "xtick.major.width": 0.5,
- "ytick.major.width": 0.5,
- "xtick.minor.width": 0.5,
- "ytick.minor.width": 0.5
- })
- # genes and stage order
- keys = ['Nfia', 'Nfil3', 'Atf5', 'Atf3', 'Jun']
- order = ['Control', 'Early', 'Mid', 'End']
- # prepare figure with independent y-axes
- fig, axes = plt.subplots(
- nrows=1, ncols=len(keys),
- figsize=(6.5, 1.4712),
- sharey=False
- )
- # draw one violins
- # then compute and plot a single median dot
- for idx, (ax, gene) in enumerate(zip(axes, keys)):
- sc.pl.violin(
- adata_violin,
- keys=[gene],
- groupby='stage',
- order=order,
- rotation=90,
- multi_panel=False,
- stripplot=False, # no individual points
- jitter=False, # no jitter
- inner=None, # no quartile/median lines
- ax=ax,
- show=False
- )
- # thin the violin outlines
- for coll in ax.collections:
- if isinstance(coll, PolyCollection):
- coll.set_linewidth(0.75)
- # thin the axes spines
- for spine in ax.spines.values():
- spine.set_linewidth(0.5)
- # title & axis labels
- ax.set_title(gene)
- if idx == 0:
- ax.set_ylabel("Expression")
- else:
- ax.set_ylabel("")
- ax.set_xlabel("")
- ax.tick_params(axis='y', labelleft=True)
- ax.grid(False)
- # get expression values and overlay medians
- expr = adata_violin.obs_vector(gene)
- for xi, stage in enumerate(order):
- mask = adata_violin.obs['stage'] == stage
- med = np.median(expr[mask.values])
- ax.scatter(xi, med,
- color='black',
- marker='o',
- s=1,
- zorder=10)
- # tidy up, display & save
- plt.tight_layout()
- plt.show()
- fig.savefig(
- "figures/violin/disease_TFs_violin.svg",
- format="svg",
- bbox_inches="tight"
- )
- # %% [markdown]
- # ## Create a final anndata object where initial cholinergic neurons are replaced
- # %%
- # Read in AnnData
- adata_all = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata.h5ad"))
- adata_chol = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/cholinergic_anndata_highvol.h5ad"))
- adata_alpha = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata_label_transfer.h5ad"))
- # %%
- # Change cell class of peri-MN oligodendrocytes
- df = pd.read_csv(os.path.join(working_dir, 'periMN_oligo_metadata.csv'), index_col=0)
- common = adata_all.obs.index.intersection(df.index)
- adata_all.obs.loc[common, 'cell_class'] = 'Oligodendrocytes'
- # %%
- # Remove cells where cell_class == "Cholinergic Neurons" from adata_all
- adata_all_filtered = adata_all[adata_all.obs["cell_class"] != "Cholinergic Neurons"].copy()
- # Concatenate adata_all_filtered with adata_chol
- adata_all_chol_cat = adata_all_filtered.concatenate(adata_chol, index_unique=None)
- # Remove cells where cholinergic_type == "Alpha MNs" from adata_all_chol_cat
- adata_all_chol_filtered = adata_all_chol_cat[adata_all_chol_cat.obs["cholinergic_type"] != "Alpha MNs"].copy()
- # Concatenate adata_all_chol_cat with adata_alpha
- adata_final = adata_all_chol_filtered.concatenate(adata_alpha, index_unique=None)
- # %%
- # Neighbors and UMAP
- sc.external.pp.bbknn(adata_final, batch_key='VS') # use inplace of sc.pp.neighbors()
- sc.tl.umap(adata_final)
- sc.pl.umap(adata_final, color="cell_class")
- # %% [markdown]
- # #### Add DAMN_status metadata
- # %%
- # Create DAMN_status column
- adata_final.obs['DAMN_status'] = pd.NA # initialize with NA
- # Mark DAMN cells
- adata_final.obs.loc[
- adata_final.obs['predicted.id'].isin(['Early DAMN', 'Late DAMN']),
- 'DAMN_status'
- ] = 'DAMN'
- # Mark Non-DAMN cells (non-NA and not DAMN)
- adata_final.obs.loc[
- (~adata_final.obs['predicted.id'].isin(['Early DAMN', 'Late DAMN'])) &
- (adata_final.obs['predicted.id'].notna()),
- 'DAMN_status'
- ] = 'Non-DAMN'
- # %%
- # Save anndata
- save_name = f"adata_objects/all_anndata_final.h5ad"
- adata_final.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ## Rotate tissue sections for visualization
- # %% [markdown]
- # #### Adapted from Sun et al., Nature (2025)
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final.h5ad"))
- # %%
- # Check initial visualization
- fig, axarr = plt.subplots(20, 6, sharex=False, sharey=False, figsize=(20, 40)) # 20 rows, 6 columns
- mids = np.sort(np.unique(adata.obs['slide_section']))
- for i in range(20):
- for j in range(6):
- index = i * 6 + j # Compute the 1D index from 2D indices
- if index < len(mids): # Ensure we don't go out of bounds
- mid = mids[index]
- sub_adata = adata[adata.obs['slide_section'] == mid].copy()
- sc.pl.embedding(sub_adata, 'spatial', color="cell_class", size=50, show=False,
- ax=axarr[i, j], vmin=-10, vmax=10, title=None)
- axarr[i, j].legend_.remove() if axarr[i, j].get_legend() else None # Remove legend if it exists
- plt.tight_layout()
- plt.show()
- # %%
- # Function for rotating about the origin
- def rotate(p, origin=(0, 0), degrees=0):
- # Rigid rotation by degrees around origin
- angle = np.deg2rad(degrees)
- R = np.array([[np.cos(angle), -np.sin(angle)],
- [np.sin(angle), np.cos(angle)]])
- o = np.atleast_2d(origin)
- p = np.atleast_2d(p)
- return np.squeeze((R @ (p.T-o.T) + o.T).T)
- # %%
- # Center spatial coordinates
- new_spatial = np.zeros(adata.obsm['spatial'].shape)
- mids = np.sort(np.unique(adata.obs['slide_section']))
- for mid in mids:
- X = np.array(adata[adata.obs.slide_section==mid].obsm["spatial"].copy())
- new_spatial[adata.obs.slide_section==mid,:] = X-np.mean(X,axis=0)
- adata.obsm['spatial'] = new_spatial
- # %%
- # Rotate by manually determined angles
- rotation_dict = {
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C1': 180,
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C2': 176,
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C3': -177,
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R1C1': 168,
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R1C2': 168,
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R1C3': 173,
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R2C1': -12,
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R2C3': 123,
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C1': -175,
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C2': -167,
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C3': 170,
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C1': 7,
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C2': 7,
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C3': 8,
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C1': 170,
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C2': 166,
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C3': 158,
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R2C2': -18,
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R2C3': -175, #19
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C1': 180,
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C2': -166,
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C3': -172,
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C4': 179,
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R2C1': 0,
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R2C2': 0,
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R2C3': 3,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C1': 174,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C2': 180,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C3': 174,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C1': 0,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C2': 3,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C3': 12,
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C4': 25,
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C1': -175,
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C2': -178,
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C3': 165,
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R2C1': 6,
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R2C2': 7,
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R2C3': -4,
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R1C1': -30,
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R1C2': -19,
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R1C3': -16,
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C2': 155,
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C3': 152,
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C4': -28,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C1': 35,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C2': 34,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C3': 32,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C1': -144,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C2': -146,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C3': -144,
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C4': -144,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R1C1': 41,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R1C2': 51,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R1C3': 36,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C1': -146,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C2': -133,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C3': -145,
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C4': -155,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C1': -43,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C2': -39,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C3': -44,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C4': -41,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C1': 138,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C2': 135,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C3': 127,
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C4': 145,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C1': -27,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C2': -29,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C3': -30,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C4': -28,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C1': 145,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C2': 155,
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C3': 139,
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R1C1': 49,
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R1C2': 42,
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R1C3': 37,
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C2': -138,
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C3': -143,
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C4': -128,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R1C1': -17,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R1C2': -4,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R1C3': -7,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C1': 174,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C2': 177,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C3': -172,
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C4': 178,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C1': 0,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C2': -16,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C3': 2,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C4': 8,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C1': -171,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C2': -175,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C3': -173,
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C4': -165,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C1': -171,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C2': -12,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C3': -1,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C4': 12,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C1': 180,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C2': -174,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C3': -176,
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C4': 166,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C1': -57,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C2': -45,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C3': -38,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C4': -47,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C1': 149,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C2': 147,
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C3': 149,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C1': -13,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C2': -17,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C3': -16,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C1': 170,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C2': 166,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C3': 169,
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C4': 163}
- # %%
- # Rotational alignment
- new_spatial = np.zeros(adata.obsm['spatial'].shape)
- mids = np.sort(np.unique(adata.obs['slide_section']))
- for mid in mids:
- X = np.array(adata[adata.obs.slide_section==mid].obsm["spatial"].copy())
- new_spatial[adata.obs.slide_section==mid,:] = rotate(X, degrees=360-rotation_dict[mid])
- adata.obsm['spatial'] = new_spatial
- # %%
- # Check alignment
- fig, axarr = plt.subplots(20, 6, sharex=False, sharey=False, figsize=(20, 40)) # 20 rows, 6 columns
- mids = np.sort(np.unique(adata.obs['slide_section']))
- for i in range(20):
- for j in range(6):
- index = i * 6 + j # Compute the 1D index from 2D indices
- if index < len(mids): # Ensure we don't go out of bounds
- mid = mids[index]
- sub_adata = adata[adata.obs['slide_section'] == mid].copy()
- sc.pl.embedding(sub_adata, 'spatial', color="cell_class", size=50, show=False,
- ax=axarr[i, j], vmin=-10, vmax=10, title=None)
- axarr[i, j].legend_.remove() if axarr[i, j].get_legend() else None # Remove legend if it exists
- plt.tight_layout()
- plt.show()
- # %%
- # Save anndata
- save_name = f"adata_objects/all_anndata_final_rotated.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ## Get spatial plot
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated.h5ad"))
- # %% [markdown]
- # ## Fig. 1C
- # %%
- import matplotlib.pyplot as plt
- import squidpy as sq
- import numpy as np
- import pandas as pd
- from matplotlib.colors import ListedColormap
- # 1. Subset the data
- slide_id = '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C2'
- sub_adata = adata[adata.obs['slide_section'] == slide_id].copy()
- # 2. Assign custom cell group
- def assign_group(row):
- if row["cell_class"] == "Cholinergic Neurons":
- return "CholN"
- elif row["cell_class"] == "Non-Cholinergic Interneurons":
- return "NCI"
- elif row["cell_class"] in ["Astrocytes", "Microglia/Macrophages", "Oligodendrocytes"]:
- return "Glia"
- else:
- return "Other"
- sub_adata.obs["cell_group"] = sub_adata.obs.apply(assign_group, axis=1)
- # 3. Remove "Other" and set categorical order
- sub_adata = sub_adata[sub_adata.obs["cell_group"] != "Other"].copy()
- sub_adata.obs["cell_group"] = pd.Categorical(
- sub_adata.obs["cell_group"],
- categories=["CholN", "NCI", "Glia"],
- ordered=True
- )
- # 4. Define colormap matching the order above
- colormap = ListedColormap(["#FF0000", "#0000FF", "#C0C0C0"]) # red, blue, light gray
- # 5. Plot
- fig, ax = plt.subplots(figsize=(10, 10), dpi=600)
- sq.pl.spatial_scatter(
- sub_adata,
- shape=None,
- color="cell_group",
- size=150,
- library_id="spatial",
- palette=colormap,
- ax=ax,
- legend_loc="upper right",
- title="",
- axis_label=""
- )
- # Rasterize the scatter points
- for coll in ax.collections:
- coll.set_rasterized(True)
- # Adjust legend font size manually
- legend = ax.get_legend()
- if legend is not None:
- for text in legend.get_texts():
- text.set_fontsize(10)
- # Remove ticks and box
- ax.set_xticks([])
- ax.set_yticks([])
- ax.axis("off")
- # Rotate 180°
- ax.invert_xaxis()
- ax.invert_yaxis()
- plt.show()
- # 6. Save high-resolution SVG
- fig_path = os.path.join(figures_dir, "spatial_plot_clean.svg")
- fig.savefig(fig_path, format="svg", bbox_inches="tight", dpi=600)
- # %% [markdown]
- # ## Get a subset of high-quality tissue sections for cell type ratio and spatial analyses
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated.h5ad"))
- # %%
- # Specify which sections/section parts to keep or remove
- sections_to_keep_dict = {
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C1': "remove",
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C2': "remove",
- '202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8_R1C3': "keep",
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R1C1': "keep",
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R1C2': "keep",
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R1C3': "keep",
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R2C1': "remove",
- '202401150953_MsSpinalCord-VS119-Lumbar-Ctrl5-SOD11_Beta8_R2C3': "remove",
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C1': "left",
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C2': "keep",
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C3': "remove",
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C1': "remove",
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C2': "keep",
- '202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C3': "right",
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C1': "remove",
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C2': "left",
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C3': "remove",
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R2C2': "keep",
- '202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R2C3': "remove",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C1': "remove",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C2': "keep",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C3': "keep",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C4': "right",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R2C1': "remove",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R2C2': "keep",
- '202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R2C3': "keep",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C1': "keep",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C2': "right",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C3': "right",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C1': "remove",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C2': "keep",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C3': "right",
- '202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C4': "remove",
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C1': "right",
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C2': "keep",
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C3': "right",
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R2C1': "remove",
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R2C2': "keep",
- '202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R2C3': "remove",
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R1C1': "keep",
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R1C2': "keep",
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R1C3': "keep",
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C2': "left",
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C3': "keep",
- '202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C4': "keep",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C1': "keep",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C2': "keep",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C3': "right",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C1': "keep",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C2': "keep",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C3': "keep",
- '202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R2C4': "keep",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R1C1': "keep",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R1C2': "remove",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R1C3': "keep",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C1': "keep",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C2': "keep",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C3': "right",
- '202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C4': "remove",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C1': "right",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C2': "remove",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C3': "remove",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C4': "remove",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C1': "remove",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C2': "remove",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C3': "keep",
- '202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R2C4': "remove",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C1': "keep",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C2': "keep",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C3': "right",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C4': "keep",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C1': "left",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C2': "right",
- '202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C3': "remove",
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R1C1': "keep",
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R1C2': "remove",
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R1C3': "keep",
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C2': "remove",
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C3': "keep",
- '202406031050_MsSpinalCord-VS223-Cervical-CM5-SM5-S1_VMSC07201_R2C4': "remove",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R1C1': "keep",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R1C2': "keep",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R1C3': "keep",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C1': "right",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C2': "left",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C3': "left",
- '202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C4': "remove",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C1': "remove",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C2': "keep",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C3': "right",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C4': "keep",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C1': "keep",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C2': "left",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C3': "keep",
- '202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C4': "right",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C1': "remove",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C2': "keep",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C3': "remove",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R1C4': "keep",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C1': "remove",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C2': "left",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C3': "left",
- '202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C4': "right",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C1': "right",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C2': "keep",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C3': "keep",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C4': "keep",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C1': "remove",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C2': "remove",
- '202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C3': "left",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C1': "right",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C2': "keep",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C3': "keep",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C1': "right",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C2': "left",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C3': "remove",
- '202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C4': "right"}
- # %%
- # Create a new column in adata.obs that maps each cell's slide_section to its label in sections_to_keep_dict
- adata.obs["section_label"] = adata.obs["slide_section"].map(sections_to_keep_dict)
- # Subset the AnnData object based on the section label
- adata_keep = adata[adata.obs["section_label"] == "keep"].copy()
- adata_remove = adata[adata.obs["section_label"] == "remove"].copy()
- adata_left = adata[adata.obs["section_label"] == "left"].copy()
- adata_right = adata[adata.obs["section_label"] == "right"].copy()
- # %%
- highlight_classes = {
- "Non-Cholinergic Interneurons": "#1f77b4", # blue
- "Cholinergic Neurons": "#ff7f0e", # orange
- "Putative Ependymal Cells": "#2ca02c" # green
- }
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- import numpy as np
- # Make sure cell_class is categorical
- adata_keep.obs['cell_class'] = adata_keep.obs['cell_class'].astype('category')
- # Unique slide sections
- keep_mids = np.sort(adata_keep.obs['slide_section'].unique())
- # Loop through each slide_section
- for mid in keep_mids:
- print(f"Plotting slide_section: {mid}")
- # Subset data
- sub_keep = adata_keep[adata_keep.obs['slide_section'] == mid].copy()
- # Create color list for categories
- unique_classes = sub_keep.obs['cell_class'].cat.categories
- class_colors = [highlight_classes.get(cls, "#d3d3d3") for cls in unique_classes]
- sub_keep.uns['cell_class_colors'] = class_colors
- # Plot spatial embedding
- sc.pl.embedding(
- sub_keep,
- basis='spatial',
- color='cell_class',
- size=40,
- title=f"Slide Section: {mid}",
- show=True
- )
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- import numpy as np
- # Make sure cell_class is categorical
- adata_remove.obs['cell_class'] = adata_remove.obs['cell_class'].astype('category')
- # Unique slide sections
- remove_mids = np.sort(adata_remove.obs['slide_section'].unique())
- # Loop through each slide_section
- for mid in remove_mids:
- print(f"Plotting slide_section: {mid}")
- # Subset data
- sub_remove = adata_remove[adata_remove.obs['slide_section'] == mid].copy()
- # Create color list for categories
- unique_classes = sub_remove.obs['cell_class'].cat.categories
- class_colors = [highlight_classes.get(cls, "#d3d3d3") for cls in unique_classes]
- sub_remove.uns['cell_class_colors'] = class_colors
- # Plot spatial embedding
- sc.pl.embedding(
- sub_remove,
- basis='spatial',
- color='cell_class',
- size=40,
- title=f"Slide Section: {mid}",
- show=True
- )
- # %%
- import scanpy as sc
- import matplotlib.pyplot as plt
- import numpy as np
- from matplotlib.colors import to_hex
- # Degrees of rotation and x-shift as fraction of section width
- line_config = {
- "202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R1C1": {"angle": -6, "x_shift_frac": 0.0375},
- "202401151110_MsSpinalCord-VS119-Lumbar-Ctrl1-SOD8_VMSC07101_R1C2": {"angle": -7, "x_shift_frac": 0.0325},
- "202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302_R2C2": {"angle": 1, "x_shift_frac": 0.02},
- "202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C1": {"angle": -2, "x_shift_frac": -0.015},
- "202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C2": {"angle": -5, "x_shift_frac": 0.07},
- "202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C3": {"angle": -3, "x_shift_frac": 0.035},
- "202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C2": {"angle": 0, "x_shift_frac": 0.02},
- "202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C2": {"angle": -1, "x_shift_frac": 0.01},
- "202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C3": {"angle": 0, "x_shift_frac": 0.01},
- "202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R2C3": {"angle": -2, "x_shift_frac": 0.025},
- "202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C2": {"angle": -2, "x_shift_frac": 0.025}
- }
- # Store filtered subsets
- filtered_subsets = []
- mids = np.sort(np.unique(adata_left.obs['slide_section']))
- for mid in mids:
- print(f"Plotting and filtering slide_section: {mid}")
- sub_adata = adata_left[adata_left.obs['slide_section'] == mid].copy()
- sub_adata.obs['cell_class'] = sub_adata.obs['cell_class'].astype('category')
- coords = sub_adata.obsm['spatial']
- x = coords[:, 0]
- y = coords[:, 1]
- config = line_config.get(mid, {"angle": 0, "x_shift_frac": 0})
- angle_deg = config["angle"]
- x_shift_frac = config["x_shift_frac"]
- section_width = x.max() - x.min()
- x_shift = x_shift_frac * section_width
- x_center = np.median(x) + x_shift
- y_min, y_max = y.min(), y.max()
- print(f"Slide: {mid}")
- print(f" x range: {x.min()} to {x.max()} (width = {section_width:.2f})")
- print(f" x_shift_frac: {x_shift_frac} => x_shift: {x_shift:.2f}")
- print(f" x_center: {x_center:.2f}\n")
- # === Filter for left side ===
- theta = np.deg2rad(-angle_deg)
- rot_matrix = np.array([
- [np.cos(theta), -np.sin(theta)],
- [np.sin(theta), np.cos(theta)]
- ])
- coords_shifted = coords.copy()
- coords_shifted[:, 0] -= x_center # translate
- coords_rotated = coords_shifted @ rot_matrix.T # rotate
- mask_left = coords_rotated[:, 0] < 0
- sub_filtered = sub_adata[mask_left].copy()
- filtered_subsets.append(sub_filtered)
- # === Plotting ===
- unique_classes = sub_adata.obs['cell_class'].cat.categories
- class_colors = [highlight_classes.get(cls, "#d3d3d3") for cls in unique_classes]
- sub_adata.uns['cell_class_colors'] = class_colors
- line_x = np.array([0, 0])
- line_y = np.array([y_min - 50, y_max + 50])
- theta_plot = np.deg2rad(angle_deg)
- rot_matrix_plot = np.array([
- [np.cos(theta_plot), -np.sin(theta_plot)],
- [np.sin(theta_plot), np.cos(theta_plot)]
- ])
- line_coords = np.stack([line_x, line_y])
- rotated = rot_matrix_plot @ line_coords
- rotated[0, :] += x_center
- fig, ax = plt.subplots(1, 1, figsize=(5, 5))
- sc.pl.embedding(
- sub_adata,
- basis='spatial',
- color='cell_class',
- size=50,
- show=False,
- ax=ax,
- title=f"Slide Section: {mid}"
- )
- ax.plot(rotated[0], rotated[1], color="black", linestyle="--", linewidth=2)
- ax.axvline(x.min(), color='gray', linestyle=':', linewidth=1)
- ax.axvline(x.max(), color='gray', linestyle=':', linewidth=1)
- if ax.get_legend():
- ax.legend_.remove()
- plt.tight_layout()
- plt.show()
- # === Combine filtered subsets into new AnnData ===
- adata_left_filtered = filtered_subsets[0].concatenate(
- *filtered_subsets[1:],
- batch_key=None
- )
- print(f"Original shape: {adata_left.shape}")
- print(f"Filtered (left side) shape: {adata_left_filtered.shape}")
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- import numpy as np
- # Make sure cell_class is categorical
- adata_left_filtered.obs['cell_class'] = adata_left_filtered.obs['cell_class'].astype('category')
- # Get sorted slide_section values
- filtered_mids = np.sort(adata_left_filtered.obs['slide_section'].unique())
- for mid in filtered_mids:
- print(f"Plotting filtered (left-side) cells for slide_section: {mid}")
- # Subset filtered AnnData
- sub_filtered = adata_left_filtered[adata_left_filtered.obs['slide_section'] == mid].copy()
- # Use consistent coloring as before
- unique_classes = sub_filtered.obs['cell_class'].cat.categories
- class_colors = [highlight_classes.get(cls, "#d3d3d3") for cls in unique_classes]
- sub_filtered.uns['cell_class_colors'] = class_colors
- # Plot
- sc.pl.embedding(
- sub_filtered,
- basis='spatial',
- color='cell_class',
- size=50,
- title=f"Filtered Left Cells: {mid}",
- show=True
- )
- # %%
- import scanpy as sc
- import matplotlib.pyplot as plt
- import numpy as np
- from matplotlib.colors import to_hex
- # Degrees of rotation and x-shift as fraction of section width
- line_config = {
- "202401150956_MsSpinalCord-VS119-Cervical-Ctrl2-SOD10_VMSC10802_R2C3": {"angle": -2, "x_shift_frac": -0.05},
- "202401241119_HuSpinalCord-VS119-Cervical-Ctrl3-SOD11_VMSC07101_reanalysis_R1C4": {"angle": 6, "x_shift_frac": 0},
- "202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C2": {"angle": -4, "x_shift_frac": -0.02},
- "202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R1C3": {"angle": 0, "x_shift_frac": -0.065},
- "202401241136_HuSpinalCord-VS119-Lumbar-Ctrl2-SOD9_Beta8_R2C3": {"angle": 0, "x_shift_frac": -0.005},
- "202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C1": {"angle": -2, "x_shift_frac": -0.01},
- "202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802_R1C3": {"angle": 1, "x_shift_frac": -0.18},
- "202405281419_MsSpinalCord-VS223-Lumbar-CM5-SM8-S1_VMSC01801_R1C3": {"angle": -3, "x_shift_frac": -0.02},
- "202405281527_MsSpinalCord-VS223-Lubmar-CM4-SM6-S1_VMSC02011_R2C3": {"angle": 6, "x_shift_frac": -0.045},
- "202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701_R1C1": {"angle": 2, "x_shift_frac": -0.015},
- "202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R1C3": {"angle": 4, "x_shift_frac": -0.015},
- "202406031049_MsSpinalCord-VS223-Cervical-CE2-SE3-S1_VMSC10802_R2C2": {"angle": 2, "x_shift_frac": 0},
- "202406071118_MsSpinalCord-VS223-Cervical-CM2-SM4-S1_VMSC10802_R2C1": {"angle": 2, "x_shift_frac": -0.01},
- "202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R1C3": {"angle": -5, "x_shift_frac": -0.0075},
- "202406071119_MsSpinalCord-VS223-Lumbar-CE1-SE3-S2_Beta8_R2C4": {"angle": 3, "x_shift_frac": -0.0525},
- "202406071120_MsSpinalCord-VS223-Lumbar-CE5-SE6-S1_VMSC01801_R2C4": {"angle": 4, "x_shift_frac": -0.01},
- "202406071222_MsSpinalCord-VS223-Cervical-CE5-SE5-S2_VMSC02011_R1C1": {"angle": 6, "x_shift_frac": -0.035},
- "202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R1C1": {"angle": -1, "x_shift_frac": -0.005},
- "202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C1": {"angle": 0, "x_shift_frac": -0.008},
- "202406071239_MsSpinalCord-VS223-Lumbar-CM3-SM5-S2_VMSC02701_R2C4": {"angle": 5, "x_shift_frac": -.03}
- }
- # Store filtered subsets
- filtered_subsets = []
- mids = np.sort(np.unique(adata_right.obs['slide_section']))
- for mid in mids:
- print(f"Plotting and filtering slide_section: {mid}")
- sub_adata = adata_right[adata_right.obs['slide_section'] == mid].copy()
- sub_adata.obs['cell_class'] = sub_adata.obs['cell_class'].astype('category')
- coords = sub_adata.obsm['spatial']
- x = coords[:, 0]
- y = coords[:, 1]
- config = line_config.get(mid, {"angle": 0, "x_shift_frac": 0})
- angle_deg = config["angle"]
- x_shift_frac = config["x_shift_frac"]
- section_width = x.max() - x.min()
- x_shift = x_shift_frac * section_width
- x_center = np.median(x) + x_shift
- y_min, y_max = y.min(), y.max()
- print(f"Slide: {mid}")
- print(f" x range: {x.min()} to {x.max()} (width = {section_width:.2f})")
- print(f" x_shift_frac: {x_shift_frac} => x_shift: {x_shift:.2f}")
- print(f" x_center: {x_center:.2f}\n")
- # === Filter for right side ===
- theta = np.deg2rad(-angle_deg)
- rot_matrix = np.array([
- [np.cos(theta), -np.sin(theta)],
- [np.sin(theta), np.cos(theta)]
- ])
- coords_shifted = coords.copy()
- coords_shifted[:, 0] -= x_center # translate
- coords_rotated = coords_shifted @ rot_matrix.T # rotate
- mask_right = coords_rotated[:, 0] > 0
- sub_filtered = sub_adata[mask_right].copy()
- filtered_subsets.append(sub_filtered)
- # === Plotting ===
- unique_classes = sub_adata.obs['cell_class'].cat.categories
- class_colors = [highlight_classes.get(cls, "#d3d3d3") for cls in unique_classes]
- sub_adata.uns['cell_class_colors'] = class_colors
- line_x = np.array([0, 0])
- line_y = np.array([y_min - 50, y_max + 50])
- theta_plot = np.deg2rad(angle_deg)
- rot_matrix_plot = np.array([
- [np.cos(theta_plot), -np.sin(theta_plot)],
- [np.sin(theta_plot), np.cos(theta_plot)]
- ])
- line_coords = np.stack([line_x, line_y])
- rotated = rot_matrix_plot @ line_coords
- rotated[0, :] += x_center
- fig, ax = plt.subplots(1, 1, figsize=(5, 5))
- sc.pl.embedding(
- sub_adata,
- basis='spatial',
- color='cell_class',
- size=50,
- show=False,
- ax=ax,
- title=f"Slide Section: {mid}"
- )
- ax.plot(rotated[0], rotated[1], color="black", linestyle="--", linewidth=2)
- ax.axvline(x.min(), color='gray', linestyle=':', linewidth=1)
- ax.axvline(x.max(), color='gray', linestyle=':', linewidth=1)
- if ax.get_legend():
- ax.legend_.remove()
- plt.tight_layout()
- plt.show()
- # === Combine filtered subsets into new AnnData ===
- adata_right_filtered = filtered_subsets[0].concatenate(
- *filtered_subsets[1:],
- batch_key=None
- )
- print(f"Original shape: {adata_right.shape}")
- print(f"Filtered (right side) shape: {adata_right_filtered.shape}")
- # %%
- import matplotlib.pyplot as plt
- import scanpy as sc
- import numpy as np
- # Make sure cell_class is categorical
- adata_right_filtered.obs['cell_class'] = adata_right_filtered.obs['cell_class'].astype('category')
- # Get sorted slide_section values
- filtered_mids = np.sort(adata_right_filtered.obs['slide_section'].unique())
- for mid in filtered_mids:
- print(f"Plotting filtered (right-side) cells for slide_section: {mid}")
- # Subset filtered AnnData
- sub_filtered = adata_right_filtered[adata_right_filtered.obs['slide_section'] == mid].copy()
- # Use consistent coloring as before
- unique_classes = sub_filtered.obs['cell_class'].cat.categories
- class_colors = [highlight_classes.get(cls, "#d3d3d3") for cls in unique_classes]
- sub_filtered.uns['cell_class_colors'] = class_colors
- # Plot
- sc.pl.embedding(
- sub_filtered,
- basis='spatial',
- color='cell_class',
- size=50,
- title=f"Filtered right Cells: {mid}",
- show=True
- )
- # %%
- adata_filtered = adata_keep.concatenate(
- adata_left_filtered,
- adata_right_filtered,
- index_unique=None
- )
- # %%
- # Save anndata
- save_name = f"adata_objects/all_anndata_final_rotated_filtered.h5ad"
- adata_filtered.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ## Cell types neighboring motor neurons in control and disease
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated_filtered.h5ad"))
- # %% [markdown]
- # ## Fig. 6A
- # %%
- import os
- import pandas as pd
- import numpy as np
- import seaborn as sns
- import matplotlib.pyplot as plt
- from collections import defaultdict
- import squidpy as sq
- # Ensure all text in the figure is at least size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Parameters
- groups = ["Non-DAMN", "DAMN"]
- title_map = {"Non-DAMN": "Non-DM", "DAMN": "DM"}
- # Abbreviations for cell class names
- abbr = {
- "Astrocytes": "Ast",
- "Cholinergic Neurons": "CholN",
- "Disease-Associated Interneurons": "DAI",
- "Microglia/Macrophages": "MG",
- "Non-Cholinergic Interneurons": "NCI",
- "Oligodendrocytes": "Oligo",
- "Other": "Other",
- "Putative Ependymal Cells": "Epen",
- "Putative Perivascular/Meningeal Cells": "PVM",
- "Putative Vascular Cells": "Vasc"
- }
- # Prepare accumulators
- neighbor_counts = {g: [defaultdict(int) for _ in range(6)] for g in groups}
- cell_totals = {g: 0 for g in groups}
- # Loop over each slide_section
- for sec in adata.obs["slide_section"].unique():
- sub = adata[adata.obs["slide_section"] == sec].copy()
- # compute 6-NN on this section
- sq.gr.spatial_neighbors(sub, coord_type="generic", n_neighs=6, key_added="spatial_sec")
- conn = sub.obsp["spatial_sec_connectivities"]
- classes = sub.obs["cell_class"]
- status = sub.obs["DAMN_status"]
- for g in groups:
- mask = (status == g)
- idx = np.where(mask)[0]
- cell_totals[g] += len(idx)
- for i in idx:
- row = conn[i].tocoo()
- top6 = sorted(zip(row.col, row.data), key=lambda x: -x[1])[:6]
- for rank, (nbr, _) in enumerate(top6):
- neighbor_counts[g][rank][classes.iloc[nbr]] += 1
- # Build DataFrame of proportions: count / total cells
- rows = []
- for g in groups:
- total_cells = cell_totals[g]
- for rank, ctr in enumerate(neighbor_counts[g], start=1):
- for cell_type, count in ctr.items():
- rows.append({
- "Neighbor Rank": str(rank),
- "Cell Class": cell_type,
- "Proportion": count / total_cells,
- "Group": g
- })
- df_combined = pd.DataFrame(rows)
- # Pivot for plotting
- df_pivot = df_combined.pivot_table(
- index=["Neighbor Rank", "Group"],
- columns="Cell Class",
- values="Proportion",
- aggfunc="sum"
- ).fillna(0)
- # Plot: Non-DM vs DM with the desired size and save as SVG
- fig, axes = plt.subplots(1, 2, figsize=(3.25, 2.15), sharey=True)
- for ax, g in zip(axes, groups):
- data = df_pivot.xs(g, level="Group")
- bottom = np.zeros(len(data))
- for ct in data.columns:
- ax.bar(
- data.index,
- data[ct],
- bottom=bottom,
- label=abbr.get(ct, ct)
- )
- bottom += data[ct]
- ax.set_title(title_map.get(g, g))
- ax.set_xlabel("Neighbor Rank")
- ax.set_ylabel("Proportion")
- # legend on the right of the second subplot, with abbreviations
- axes[1].legend(
- title="Cell Class",
- bbox_to_anchor=(1.02, 0.5),
- loc="center left",
- ncol=1
- )
- plt.tight_layout()
- # Define the filename and directory
- fig_path = os.path.join(figures_dir, "neighbor_barplots_dm.svg")
- # Save high-quality SVG
- fig.savefig(
- fig_path,
- format="svg",
- bbox_inches="tight"
- )
- # %%
- from scipy.stats import fisher_exact
- from statsmodels.stats.multitest import multipletests
- # Create a pivot table from df_combined for easier comparison;
- # here, the index will be (Neighbor Rank, Group)
- df_pivot = df_combined.pivot_table(
- index=["Neighbor Rank", "Group"],
- columns="Cell Class",
- values="Proportion",
- aggfunc="sum"
- ).fillna(0)
- # Get indices of DAMN and Non-DAMN cells
- damn_indices = np.where(adata.obs["DAMN_status"] == "DAMN")[0]
- nondamn_indices = np.where(adata.obs["DAMN_status"] == "Non-DAMN")[0]
- # Precompute total numbers in each group using the same criteria from Block 1.
- damn_total = len(damn_indices) # DAMN group total
- nondamn_total = len(nondamn_indices) # Non-DAMN group total
- fisher_results = []
- # Loop over each unique neighbor rank found in df_combined
- for rank in df_combined["Neighbor Rank"].unique():
- # For each cell class, compare proportions between DAMN and Non-DAMN at this rank.
- for cell_class in df_pivot.columns:
- # Use the pivot table to get proportions:
- try:
- prop_damn = df_pivot.loc[(rank, "DAMN"), cell_class]
- except KeyError:
- prop_damn = 0
- try:
- prop_nondamn = df_pivot.loc[(rank, "Non-DAMN"), cell_class]
- except KeyError:
- prop_nondamn = 0
- # Convert proportions to counts using precomputed totals:
- count_damn = prop_damn * damn_total
- count_nondamn = prop_nondamn * nondamn_total
- # Build the 2x2 contingency table:
- contingency_table = np.array([
- [count_damn, damn_total - count_damn],
- [count_nondamn, nondamn_total - count_nondamn]
- ])
- # Perform Fisher's Exact Test:
- try:
- odds_ratio, p_value = fisher_exact(contingency_table)
- except Exception:
- odds_ratio, p_value = np.nan, np.nan
- fisher_results.append({
- "Neighbor Rank": rank,
- "Cell Class": cell_class,
- "Odds Ratio": odds_ratio,
- "p-value": p_value
- })
- # Convert the results list to a DataFrame
- fisher_results_df = pd.DataFrame(fisher_results)
- # Apply multiple hypothesis corrections:
- fisher_results_df["padj"] = multipletests(fisher_results_df["p-value"], method="bonferroni")[1]
- # Display the Fisher's test results DataFrame
- fisher_results_df
- # %%
- # Save as CSV without row names
- fisher_results_df.to_csv("DAMN_nearest_neighbors_fisher_results.csv", index=False)
- # %% [markdown]
- # ## Fig. S7A
- # %%
- import os
- import pandas as pd
- import numpy as np
- import seaborn as sns
- import matplotlib.pyplot as plt
- from collections import defaultdict
- import squidpy as sq
- # Ensure all text in the figure is at least size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Abbreviations for cell class names
- abbr = {
- "Astrocytes": "Ast",
- "Cholinergic Neurons": "CholN",
- "Disease-Associated Interneurons": "DAI",
- "Microglia/Macrophages": "MG",
- "Non-Cholinergic Interneurons": "NCI",
- "Oligodendrocytes": "Oligo",
- "Other": "Other",
- "Putative Ependymal Cells": "Epen",
- "Putative Perivascular/Meningeal Cells": "PVM",
- "Putative Vascular Cells": "Vasc"
- }
- # Parameters
- cholinergic_types = ["Alpha MNs", "Gamma MNs", "Gamma* MNs"]
- conditions = {
- "Control": lambda st: st == "Control",
- "Disease": lambda st: st.isin(["Early", "Mid", "End"])
- }
- # Prepare accumulators
- neighbor_counts = {}
- cell_totals = {}
- for sec in adata.obs["slide_section"].unique():
- sub = adata[adata.obs["slide_section"] == sec].copy()
- sq.gr.spatial_neighbors(sub, coord_type="generic", n_neighs=6, key_added="spatial_sec")
- conn = sub.obsp["spatial_sec_connectivities"]
- classes = sub.obs["cell_class"]
- stages = sub.obs["stage"]
- types = sub.obs["cholinergic_type"]
- for chol in cholinergic_types:
- for cond_label, cond_fn in conditions.items():
- mask = (types == chol) & cond_fn(stages)
- idx = np.where(mask)[0]
- if len(idx) == 0:
- continue
- key = (chol, cond_label)
- if key not in neighbor_counts:
- neighbor_counts[key] = [defaultdict(int) for _ in range(6)]
- cell_totals[key] = 0
- cell_totals[key] += len(idx)
- for i in idx:
- row = conn[i].tocoo()
- top6 = sorted(zip(row.col, row.data), key=lambda x: -x[1])[:6]
- for rank, (nbr, _) in enumerate(top6, start=1):
- neighbor_counts[key][rank-1][classes.iloc[nbr]] += 1
- # Build DataFrame of proportions
- rows = []
- for (chol, cond_label), rank_dicts in neighbor_counts.items():
- total_cells = cell_totals[(chol, cond_label)]
- for rank, ctr in enumerate(rank_dicts, start=1):
- for cell_type, count in ctr.items():
- rows.append({
- "Neighbor Rank": str(rank),
- "Cell Class": cell_type,
- "Proportion": count / total_cells,
- "Cholinergic Type": chol,
- "Stage": cond_label
- })
- df_combined = pd.DataFrame(rows)
- # Pivot for plotting
- df_pivot = df_combined.pivot_table(
- index=["Neighbor Rank", "Cholinergic Type", "Stage"],
- columns="Cell Class",
- values="Proportion",
- aggfunc="sum"
- ).fillna(0)
- # Plot: 3 rows = skeletal MN types, 2 cols = conditions
- fig, axes = plt.subplots(
- nrows=3,
- ncols=2,
- figsize=(3.25, 4.3),
- sharey=True
- )
- for i, chol in enumerate(cholinergic_types):
- for j, cond_label in enumerate(conditions):
- ax = axes[i, j]
- try:
- data = df_pivot.xs(
- (chol, cond_label),
- level=("Cholinergic Type", "Stage")
- )
- except KeyError:
- data = pd.DataFrame()
- bottom = np.zeros(len(data))
- for ct in data.columns:
- ax.bar(
- data.index,
- data[ct],
- bottom=bottom,
- label=abbr.get(ct, ct)
- )
- bottom += data[ct]
- # custom titles with extra padding for Gamma MNs rows
- if chol == "Gamma MNs" and cond_label == "Control":
- title_str = "Gamma MNs - Control "
- elif chol == "Gamma MNs" and cond_label == "Disease":
- title_str = " Gamma MNs - Disease"
- elif chol == "Gamma* MNs" and cond_label == "Control":
- title_str = "Gamma* MNs - Control "
- elif chol == "Gamma* MNs" and cond_label == "Disease":
- title_str = " Gamma* MNs - Disease"
- else:
- title_str = f"{chol} – {cond_label}"
- ax.set_title(title_str)
- ax.set_xlabel("Neighbor Rank")
- if j == 0:
- ax.set_ylabel("Proportion")
- else:
- ax.set_ylabel("")
- # only add legend on the top‐right subplot
- if i == 0 and j == 1:
- handles, labels = ax.get_legend_handles_labels()
- ax.legend(
- handles,
- labels,
- title="Cell Class",
- bbox_to_anchor=(1.02, 0.5),
- loc="center left",
- labelspacing=0.2,
- handletextpad=0.3,
- columnspacing=0.5
- )
- plt.tight_layout()
- fig_path = os.path.join(figures_dir, "neighbor_barplots_skeletal_disease.svg")
- fig.savefig(fig_path, format="svg", bbox_inches="tight")
- # %%
- from scipy.stats import fisher_exact
- from statsmodels.stats.multitest import multipletests
- # Extract needed data
- stage = adata.obs["stage"]
- # Pre-compute totals for each condition
- control_total = len(np.where(stage == "Control")[0])
- disease_total = len(np.where(stage.isin(["Early", "Mid", "End"]))[0])
- fisher_results = []
- # Loop over each cholinergic type, neighbor rank (1-6), and each cell class
- for cholinergic_group in cholinergic_types:
- for rank in range(1, 7):
- rank_str = str(rank)
- # Get proportions for Control and Disease for this (cholinergic_group, rank)
- try:
- control_data = df_pivot.loc[(rank_str, cholinergic_group, "Control")]
- except KeyError:
- control_data = pd.Series(0, index=df_pivot.columns)
- try:
- disease_data = df_pivot.loc[(rank_str, cholinergic_group, "Disease")]
- except KeyError:
- disease_data = pd.Series(0, index=df_pivot.columns)
- # Loop over each cell class
- for cell_class in df_pivot.columns:
- control_prop = control_data[cell_class]
- disease_prop = disease_data[cell_class]
- control_count = control_prop * control_total
- disease_count = disease_prop * disease_total
- contingency_table = np.array([
- [disease_count, disease_total - disease_count], # Disease group as the first row
- [control_count, control_total - control_count] # Control group as the second row
- ])
- # Perform Fisher's Exact Test
- try:
- odds_ratio, p_value = fisher_exact(contingency_table)
- except Exception:
- odds_ratio, p_value = np.nan, np.nan
- fisher_results.append({
- "Cholinergic Type": cholinergic_group,
- "Neighbor Rank": rank,
- "Cell Class": cell_class,
- "Odds Ratio": odds_ratio,
- "p-value": p_value
- })
- # Convert results to DataFrame
- fisher_results_df = pd.DataFrame(fisher_results)
- # Apply Bonferroni correction to the p-values
- fisher_results_df["padj"] = multipletests(fisher_results_df["p-value"], method="bonferroni")[1]
- # Display the final results
- fisher_results_df
- # %%
- # Save as CSV without row names
- fisher_results_df.to_csv("Skeletal_nearest_neighbors_fisher_results.csv", index=False)
- # %% [markdown]
- # ## Cell type abundance changes with disease
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated_filtered.h5ad"))
- # %% [markdown]
- # ## Fig. 6B and S7B
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.stats import ttest_ind
- from statannotations.Annotator import Annotator
- from statsmodels.stats.multitest import multipletests
- # ——————————————————————————————————
- # 1) Ensure all text in the figures is size 7
- # ——————————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # 2) Color mapping
- point_colors = {
- "Microglia/Macrophages": "#d62728",
- "Astrocytes": "#1f77b4",
- "Oligodendrocytes": "#8c564b",
- "Alpha MNs": "#1f77b4",
- "Gamma MNs": "#2ca02c",
- "Gamma* MNs": "#d62728",
- }
- stage_order = ["Control", "Early", "Mid", "End"]
- # 3) Ratio function (denominator = all interneurons)
- def compute_ratio(x, column, value):
- num = (x[column] == value).sum()
- den = (
- (x["cell_class"] == "Non-Cholinergic Interneurons").sum() +
- (x["cholinergic_type"] == "Cholinergic Interneurons").sum() +
- (x["cholinergic_type"] == "Disease-Associated Interneurons").sum()
- )
- return num / den if den > 0 else np.nan
- # =============================================================================
- # First plot: 1×2 for Microglia/Macrophages & Astrocytes
- # =============================================================================
- cell_types_g1 = [
- ("cell_class", "Microglia/Macrophages", "Microglia/Macrophages"),
- ("cell_class", "Astrocytes", "Astrocytes"),
- ]
- fig1, axs1 = plt.subplots(1, 2, figsize=(3.25, 2.15), sharey=False)
- stat_tests1 = {}
- for ax, (col, val, title) in zip(axs1, cell_types_g1):
- # a) per-slide_stage raw ratios
- df = (
- adata.obs
- .groupby(["slide_section", "slide_stage", "stage"], observed=True)
- .apply(lambda x: compute_ratio(x, col, val))
- .reset_index(name="celltype_ratio")
- )
- # b) average per slide_stage
- df_avg = (
- df
- .groupby("slide_stage", as_index=False)
- .agg(celltype_ratio=("celltype_ratio", "mean"),
- stage=("stage", "first"))
- )
- df_avg["stage"] = pd.Categorical(df_avg["stage"], categories=stage_order, ordered=True)
- # c) scale by Control mean
- ctrl_mean = df_avg.loc[df_avg.stage == "Control", "celltype_ratio"].mean()
- df_avg["scaled"] = df_avg["celltype_ratio"] / ctrl_mean
- # d) plotting with dot size = 4
- sns.stripplot(
- data=df_avg, x="stage", y="scaled",
- jitter=True, alpha=0.7, size=4,
- color="black", ax=ax
- )
- sns.pointplot(
- data=df_avg, x="stage", y="scaled",
- estimator="mean", errorbar=("ci", 95),
- markers="o", color=point_colors[title],
- ax=ax
- )
- ax.set_title(title)
- ax.set_xlabel("")
- ax.set_ylabel("Relative Abundance")
- # e) two-sided, unequal-variance t-tests vs Control
- ctrl_vals = df_avg.loc[df_avg.stage == "Control", "scaled"]
- raw_p = []
- for s in ["Early", "Mid", "End"]:
- grp = df_avg.loc[df_avg.stage == s, "scaled"]
- raw_p.append(
- ttest_ind(ctrl_vals, grp, equal_var=False).pvalue
- if len(grp) > 0 else np.nan
- )
- raw_p = np.array(raw_p)
- mask = ~np.isnan(raw_p)
- corr = np.full_like(raw_p, np.nan, dtype=float)
- if mask.sum() > 0:
- _, p_corr, _, _ = multipletests(raw_p[mask], method="bonferroni")
- corr[mask] = p_corr
- pvals_corr = corr.tolist()
- # f) annotate
- annot = Annotator(
- ax,
- [("Control", "Early"), ("Control", "Mid"), ("Control", "End")],
- data=df_avg, x="stage", y="scaled",
- order=stage_order, perform_stat_test=False
- )
- annot.configure(test=None, text_format="star", loc="inside")
- annot.set_pvalues_and_annotate(pvals_corr)
- # Nudge tick labels on fig1
- for ax in axs1:
- labels = [t.get_text() for t in ax.get_xticklabels()]
- new_labels = [
- "Control " if lbl == "Control" else
- " Early" if lbl == "Early" else
- lbl
- for lbl in labels
- ]
- ax.set_xticklabels(new_labels)
- fig1.tight_layout()
- fig1.savefig(
- os.path.join(figures_dir, "celltype_group1.svg"),
- format="svg", bbox_inches="tight"
- )
- plt.show()
- # =============================================================================
- # Second plot: 2×2 for Oligodendrocytes & Motor Neuron subtypes
- # =============================================================================
- cell_types_g2 = [
- ("cell_class", "Oligodendrocytes", "Oligodendrocytes"),
- ("cholinergic_type", "Alpha MNs", "Alpha MNs"),
- ("cholinergic_type", "Gamma MNs", "Gamma MNs"),
- ("cholinergic_type", "Gamma* MNs", "Gamma* MNs"),
- ]
- fig2, axs2 = plt.subplots(2, 2, figsize=(3.25, 4.3), sharey=False)
- axs2_flat = axs2.flatten()
- stat_tests2 = {}
- for ax, (col, val, title) in zip(axs2_flat, cell_types_g2):
- # a) per-slide_stage raw ratios
- df = (
- adata.obs
- .groupby(["slide_section", "slide_stage", "stage"], observed=True)
- .apply(lambda x: compute_ratio(x, col, val))
- .reset_index(name="celltype_ratio")
- )
- # b) average per slide_stage
- df_avg = (
- df
- .groupby("slide_stage", as_index=False)
- .agg(celltype_ratio=("celltype_ratio", "mean"),
- stage=("stage", "first"))
- )
- df_avg["stage"] = pd.Categorical(df_avg["stage"], categories=stage_order, ordered=True)
- # c) scale by Control mean
- ctrl_mean = df_avg.loc[df_avg.stage == "Control", "celltype_ratio"].mean()
- df_avg["scaled"] = df_avg["celltype_ratio"] / ctrl_mean
- # d) plotting
- sns.stripplot(
- data=df_avg, x="stage", y="scaled",
- jitter=True, alpha=0.7, size=4,
- color="black", ax=ax
- )
- sns.pointplot(
- data=df_avg, x="stage", y="scaled",
- estimator="mean", errorbar=("ci", 95), markers="o",
- color=point_colors[title], ax=ax
- )
- ax.set_title(title)
- ax.set_xlabel("")
- # only leftmost column gets y-label
- if ax in (axs2[0, 0], axs2[1, 0]):
- ax.set_ylabel("Relative Abundance")
- else:
- ax.set_ylabel("")
- # e) one-sided t-tests vs Control
- ctrl_vals = df_avg.loc[df_avg.stage == "Control", "scaled"]
- raw_p = []
- for s in ["Early", "Mid", "End"]:
- grp = df_avg.loc[df_avg.stage == s, "scaled"]
- raw_p.append(
- ttest_ind(ctrl_vals, grp, equal_var=False).pvalue
- if len(grp) > 0 else np.nan
- )
- raw_p = np.array(raw_p)
- mask = ~np.isnan(raw_p)
- corr = np.full_like(raw_p, np.nan, dtype=float)
- if mask.sum() > 0:
- _, p_corr, _, _ = multipletests(raw_p[mask], method="bonferroni")
- corr[mask] = p_corr
- pvals_corr = corr.tolist()
- # f) annotate
- annot = Annotator(
- ax,
- [("Control", "Early"), ("Control", "Mid"), ("Control", "End")],
- data=df_avg, x="stage", y="scaled",
- order=stage_order, perform_stat_test=False
- )
- annot.configure(test=None, text_format="star", loc="inside")
- annot.set_pvalues_and_annotate(pvals_corr)
- # Nudge tick labels on fig2
- for ax in axs2_flat:
- labels = [t.get_text() for t in ax.get_xticklabels()]
- new_labels = [
- "Control " if lbl == "Control" else
- " Early" if lbl == "Early" else
- lbl
- for lbl in labels
- ]
- ax.set_xticklabels(new_labels)
- fig2.tight_layout()
- fig2.savefig(
- os.path.join(figures_dir, "celltype_group2.svg"),
- format="svg", bbox_inches="tight"
- )
- plt.show()
- # %% [markdown]
- # ## Fig. 6H
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.stats import ttest_ind
- from statannotations.Annotator import Annotator
- from statsmodels.stats.multitest import multipletests
- # ——————————————————————————————————
- # Ensure all text is at least size 7
- # ——————————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Define cell class of interest
- cell_class_value = "Disease-Associated Interneurons"
- # Order for stage categories
- stage_order = ["Control", "Early", "Mid", "End"]
- # 1) Compute ratios per slide_section using NCI + Cholinergic as denominator
- df = (
- adata.obs
- .groupby(["slide_section", "slide_stage", "stage", "region"], observed=True)
- .apply(lambda x: (x["cell_class"] == cell_class_value).sum() /
- ((x["cell_class"] == "Non-Cholinergic Interneurons") |
- (x["cell_class"] == "Cholinergic Interneurons")).sum()
- if ((x["cell_class"] == "Non-Cholinergic Interneurons") |
- (x["cell_class"] == "Cholinergic Interneurons")).sum() > 0
- else np.nan
- )
- .reset_index(name="celltype_ratio")
- )
- # 2) Average ratio per slide_stage
- df_stage_avg = df.groupby("slide_stage", as_index=False).agg(
- celltype_ratio=("celltype_ratio", "mean"),
- stage=("stage", "first"),
- region=("region", "first")
- )
- # 3) Categorical ordering
- df_stage_avg["stage"] = pd.Categorical(df_stage_avg["stage"],
- categories=stage_order,
- ordered=True)
- # 4) Scale by Control mean
- control_mean = df_stage_avg.loc[df_stage_avg.stage == "Control", "celltype_ratio"].mean()
- df_stage_avg["celltype_ratio_scaled"] = df_stage_avg["celltype_ratio"] / control_mean
- # 5) Plot
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- sns.stripplot(
- data=df_stage_avg,
- x="stage", y="celltype_ratio_scaled",
- jitter=True, alpha=0.7, color="black",size=4,
- ax=ax
- )
- sns.pointplot(
- data=df_stage_avg,
- x="stage", y="celltype_ratio_scaled",
- estimator="mean", errorbar=("ci", 95),
- markers="o", color="#55A868",
- ax=ax
- )
- ax.set_title("DAI")
- ax.set_xlabel(None)
- ax.set_ylabel("Relative Abundance")
- # 6) Comparisons
- box_pairs = [("Control", "Early"), ("Control", "Mid"), ("Control", "End")]
- control_vals = df_stage_avg.loc[df_stage_avg.stage == "Control", "celltype_ratio_scaled"]
- early_vals = df_stage_avg.loc[df_stage_avg.stage == "Early", "celltype_ratio_scaled"]
- mid_vals = df_stage_avg.loc[df_stage_avg.stage == "Mid", "celltype_ratio_scaled"]
- end_vals = df_stage_avg.loc[df_stage_avg.stage == "End", "celltype_ratio_scaled"]
- # 7) Unequal‑variance t‑tests
- p_ce = ttest_ind(control_vals, early_vals, equal_var=False).pvalue if len(early_vals)>0 else np.nan
- p_cm = ttest_ind(control_vals, mid_vals, equal_var=False).pvalue if len(mid_vals)>0 else np.nan
- p_cE = ttest_ind(control_vals, end_vals, equal_var=False).pvalue if len(end_vals)>0 else np.nan
- pvals = [p_ce, p_cm, p_cE]
- # 8) Bonferroni correction
- _, pvals_corr, _, _ = multipletests(pvals, method="bonferroni")
- # 9) Annotate
- annotator = Annotator(
- ax, box_pairs,
- data=df_stage_avg,
- x="stage", y="celltype_ratio_scaled",
- order=stage_order,
- perform_stat_test=False
- )
- annotator.configure(test=None, text_format="star", loc="inside")
- annotator.set_pvalues_and_annotate(pvals_corr)
- # 9b) Nudge “Control” and “Early” apart via padded labels
- padded = []
- for l in stage_order:
- if l == "Control":
- padded.append("Control ")
- elif l == "Early":
- padded.append(" Early")
- else:
- padded.append(l)
- ax.set_xticklabels(padded)
- # 10) Save & show
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "DAI_ratio_small.svg")
- fig.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # 11) Print corrected p‑values
- print("Bonferroni‑corrected p‑values:")
- for cmp, p in zip(["Control vs Early","Control vs Mid","Control vs End"], pvals_corr):
- print(f" {cmp}: p = {p:.3e}")
- # %% [markdown]
- # ## Microglia/Macrophage follow-up
- # %% [markdown]
- # ### Check if control microglia/macrophages are near gene & transcript cutoff
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated.h5ad"))
- # %% [markdown]
- # ## Fig. S7C
- # %%
- import os
- import seaborn as sns
- import matplotlib.pyplot as plt
- import pandas as pd
- import numpy as np
- # Filter to only Microglia/Macrophages
- microglia = adata[adata.obs['cell_class'] == "Microglia/Macrophages"].copy()
- # Split into control and disease
- control = microglia[microglia.obs['stage'] == "Control"].copy()
- disease = microglia[microglia.obs['stage'] != "Control"].copy()
- # Create a combined DataFrame for plotting
- control_df = pd.DataFrame({
- 'num_transcripts': control.obs['num_transcripts'],
- 'num_genes': control.obs['num_genes'],
- 'condition': 'Control'
- })
- disease_df = pd.DataFrame({
- 'num_transcripts': disease.obs['num_transcripts'],
- 'num_genes': disease.obs['num_genes'],
- 'condition': 'Disease'
- })
- plot_df = pd.concat([control_df, disease_df])
- # Plotting
- fig, axs = plt.subplots(1, 2, figsize=(1.625, 2.15))
- for ax, ycol, ylabel in zip(
- axs,
- ['num_transcripts', 'num_genes'],
- ['# Transcripts', '# Genes']
- ):
- # draw violin, no inner box
- sns.violinplot(
- data=plot_df,
- x='condition',
- y=ycol,
- color="#d62728",
- inner=None,
- ax=ax
- )
- # drop x-label, set y-label
- ax.set_xlabel("")
- ax.set_ylabel(ylabel)
- # adjust y-axis spacing
- ax.tick_params(axis='y', pad = 0.5) # space between ticks and tick labels
- ax.yaxis.labelpad = 0.5 # space between label and tick labels
- # rotate x-tick labels
- ax.tick_params(axis='x', rotation=90)
- # compute & plot medians
- medians = plot_df.groupby('condition')[ycol].median()
- for i, median_val in enumerate(medians):
- ax.scatter(i, median_val, color='black', s=8, zorder=3)
- # add a single title above both plots
- fig.suptitle("Microglia/Macrophages", x=0.6, y=0.93, fontsize=7)
- # adjust so that suptitle and plots don't overlap
- plt.tight_layout(rect=[0, 0, 1, 1])
- # save and show
- out_path = os.path.join(figures_dir, "microglia_transcripts_genes.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # Function to compute median and mean
- def describe_metric(data, metric):
- values = data.obs[metric]
- median = np.median(values)
- mean = np.mean(values)
- return median, mean
- # Compute for num_transcripts and num_genes
- metrics = ["num_transcripts", "num_genes"]
- for condition_name, data in zip(["Control", "Disease"], [control, disease]):
- print(f"\n{condition_name} condition:")
- for metric in metrics:
- median_val, mean_val = describe_metric(data, metric)
- print(f" {metric}: median = {median_val:.2f}, mean = {mean_val:.2f}")
- # %% [markdown]
- # ### Spatial plots of microglia/macrophages with disease
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated.h5ad"))
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- # 1) Ensure all text in the figures is size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # 2) Define (slide, section) pairs
- sections = [
- ("202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701", "R1C2"),
- ("202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302", "R2C3"),
- ("202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802", "R1C2"),
- ("202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701", "R2C3"),
- ("202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302", "R1C3"),
- ("202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802", "R2C2")
- ]
- # 3) Extract subsets
- sub_adata_list = [
- adata[(adata.obs["slide"] == slide_id) & (adata.obs["section"] == section_id)].copy()
- for slide_id, section_id in sections
- ]
- # 4) Compute global bounds with padding
- xmins, xmaxs, ymins, ymaxs = [], [], [], []
- for sub in sub_adata_list:
- coords = sub.obsm["spatial"]
- xmins.append(coords[:, 0].min()); xmaxs.append(coords[:, 0].max())
- ymins.append(coords[:, 1].min()); ymaxs.append(coords[:, 1].max())
- global_xmin, global_xmax = min(xmins), max(xmaxs)
- global_ymin, global_ymax = min(ymins), max(ymaxs)
- dx, dy = global_xmax - global_xmin, global_ymax - global_ymin
- global_xmin -= 0.05 * dx
- global_xmax += 0.00 * dx
- global_ymin -= 0.05 * dy
- global_ymax += 0.075 * dy
- # ------------------ Microglia/Macrophages (2×3) — Swapped Rows & Titles ------------------ #
- fig1, axes1 = plt.subplots(2, 3, figsize=(4.35, 2.2))
- axes1 = axes1.reshape(2, 3)
- axes1 = np.flipud(axes1).flatten() # Flip the subplot rows
- # Correctly reordered titles to match flipped layout
- titles1 = [
- "Early – Disease", "Mid – Disease", "End – Disease",
- "Early – Control", "Mid – Control", "End – Control"
- ]
- for ax, sub, title in zip(axes1, sub_adata_list, titles1):
- coords = sub.obsm["spatial"]
- mask = (sub.obs["cell_class"] == "Microglia/Macrophages").values
- ax.scatter(coords[~mask, 0], coords[~mask, 1],
- s=4, c="lightgray", marker="o", linewidth=0)
- ax.scatter(coords[mask, 0], coords[mask, 1],
- s=4, c="#d62728", marker="o", linewidth=0)
- ax.set_xlim(global_xmin, global_xmax)
- ax.set_ylim(global_ymin, global_ymax)
- ax.set_xticks([]); ax.set_yticks([])
- ax.set_title(title, pad=2)
- fig1.suptitle("Microglia/Macrophages", y=0.93, fontsize=7)
- fig1.subplots_adjust(top=0.88)
- plt.tight_layout()
- fig1_path = os.path.join(figures_dir, "microglia_macrophages.svg")
- fig1.savefig(fig1_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Reactive glia with disease
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated.h5ad"))
- # %% [markdown]
- # ### Define reactive glia
- # %%
- # Step 1: Compute the 95th percentiles for Microglia and Astrocytes in each slide
- control_95th_percentile_microglia = {}
- control_95th_percentile_astrocytes = {}
- # Exclude control only slide
- excluded_slide = "202312251515_MsSpinalCord-L4-S1-VS119-YS_Beta8"
- unique_slides = [s for s in adata.obs["slide"].unique() if s != excluded_slide]
- for slide in unique_slides:
- # Microglia: use Apoe
- control_microglia = adata.obs[
- (adata.obs["slide"] == slide) &
- (adata.obs["stage"] == "Control") &
- (adata.obs["cell_class"] == "Microglia/Macrophages")
- ]["Apoe_high_pass_norm"].dropna()
- control_95th_percentile_microglia[slide] = (
- np.percentile(control_microglia, 95) if not control_microglia.empty else np.nan
- )
- # Astrocytes: use Gfap
- control_astrocytes = adata.obs[
- (adata.obs["slide"] == slide) &
- (adata.obs["stage"] == "Control") &
- (adata.obs["cell_class"] == "Astrocytes")
- ]["Gfap_high_pass_norm"].dropna()
- control_95th_percentile_astrocytes[slide] = (
- np.percentile(control_astrocytes, 95) if not control_astrocytes.empty else np.nan
- )
- # Step 2: Define functions to determine reactivity
- def is_reactive_microglia(row):
- slide = row["slide"]
- if row["cell_class"] == "Microglia/Macrophages" and slide in control_95th_percentile_microglia:
- threshold = control_95th_percentile_microglia[slide]
- if not np.isnan(threshold):
- return row["Apoe_high_pass_norm"] > threshold
- return False
- def is_reactive_or_wm_astrocyte(row):
- slide = row["slide"]
- if row["cell_class"] == "Astrocytes" and slide in control_95th_percentile_astrocytes:
- threshold = control_95th_percentile_astrocytes[slide]
- if not np.isnan(threshold):
- return row["Gfap_high_pass_norm"] > threshold
- return False
- # Step 3: Apply the functions to create new metadata columns
- adata.obs["reactive_microglia"] = adata.obs.apply(is_reactive_microglia, axis=1)
- adata.obs["reactive_or_wm_astrocytes"] = adata.obs.apply(is_reactive_or_wm_astrocyte, axis=1)
- # %%
- # Save anndata
- save_name = f"adata_objects/all_anndata_final_rotated_w_anno.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ### Reactive glia spatial plots
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated_w_anno.h5ad"))
- # %% [markdown]
- # ## Fig. 6C
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.spatial import ConvexHull
- # Ensure all text in the figures is size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Define (slide, section) pairs
- sections = [
- ("202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701", "R1C2"),
- ("202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302", "R2C3"),
- ("202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802", "R1C2"),
- ("202405281702_MsSpinalCord-VS223-Lumbar-CE2-SE4-S1_VMSC02701", "R2C3"),
- ("202405281416_MsSpinalCord-VS223-Lumbar-CM2-SM4-S1_VMSC15302", "R1C3"),
- ("202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802", "R2C2")
- ]
- # Extract subsets
- sub_adata_list = [
- adata[(adata.obs["slide"] == slide_id) & (adata.obs["section"] == section_id)].copy()
- for slide_id, section_id in sections
- ]
- # Compute global bounds
- xmins, xmaxs, ymins, ymaxs = [], [], [], []
- for sub in sub_adata_list:
- coords = sub.obsm["spatial"]
- xmins.append(coords[:, 0].min()); xmaxs.append(coords[:, 0].max())
- ymins.append(coords[:, 1].min()); ymaxs.append(coords[:, 1].max())
- global_xmin, global_xmax = min(xmins), max(xmaxs)
- global_ymin, global_ymax = min(ymins), max(ymaxs)
- dx, dy = global_xmax - global_xmin, global_ymax - global_ymin
- global_xmin -= 0.05 * dx
- global_xmax += 0.00 * dx
- global_ymin -= 0.05 * dy
- global_ymax += 0.075 * dy
- # Subplot titles
- titles = [
- "Early – Disease", "Mid – Disease", "End – Disease",
- "Early – Control", "Mid – Control", "End – Control"
- ]
- # ------------------ Combined 3×3 Plot ------------------ #
- fig, axes = plt.subplots(3, 3, figsize=(3.1132, 1.9816))
- axes = axes.reshape(3, 3)
- # ------------------ Top 2 rows: Reactive Microglia ------------------ #
- # Flip rows vertically
- marker_axes = np.flipud(axes[:2]).flatten()
- for ax, sub, title in zip(marker_axes, sub_adata_list, titles):
- coords = sub.obsm["spatial"]
- mask = sub.obs["reactive_microglia"].astype(bool).values
- ax.scatter(coords[~mask, 0], coords[~mask, 1], s=4, c="lightgray", marker="o", linewidth=0)
- ax.scatter(coords[ mask, 0], coords[ mask, 1], s=4, c="#d62728", marker="o", linewidth=0)
- ax.set_xlim(global_xmin, global_xmax)
- ax.set_ylim(global_ymin, global_ymax)
- ax.set_xticks([]); ax.set_yticks([])
- ax.set_title(title, pad=2)
- # ------------------ Top 2 rows: Rasterize points ------------------ #
- for ax in axes[:2].flatten():
- for coll in ax.collections:
- coll.set_rasterized(True)
- # ------------------ Bottom row: Reactive Microglia Density + Convex Hull ------------------ #
- for ax, sub, title in zip(axes[2], sub_adata_list[:3], titles[:3]):
- coords = sub.obsm["spatial"]
- # Plot reactive microglia density
- mask = (sub.obs["cell_class"] == "Microglia/Macrophages") & (sub.obs["reactive_microglia"] == True)
- coords_mask = coords[mask.values]
- sns.kdeplot(
- x=coords_mask[:, 0], y=coords_mask[:, 1],
- fill=True, thresh=0, levels=100, cmap="Reds", bw_method="scott", ax=ax
- )
- # Add convex hull of all cells
- if coords.shape[0] >= 3:
- try:
- hull = ConvexHull(coords)
- hull_coords = coords[hull.vertices]
- hull_path = np.append(hull_coords, [hull_coords[0]], axis=0) # close the loop
- ax.plot(hull_path[:, 0], hull_path[:, 1], color="black", linewidth=0.5)
- except:
- pass # silently skip if hull fails
- ax.set_xlim(global_xmin, global_xmax)
- ax.set_ylim(global_ymin, global_ymax)
- ax.set_xticks([]); ax.set_yticks([])
- ax.set_title(title, pad=2)
- # Final formatting
- plt.tight_layout(pad=0.325)
- fig.suptitle("Reactive Microglia/Macrophages", y=1, fontsize=7)
- fig.subplots_adjust(top=0.90)
- out_path = os.path.join(figures_dir, "reactive_microglia_combined_with_hull.svg")
- fig.savefig(out_path, format="svg", bbox_inches="tight", dpi=600)
- plt.show()
- # %% [markdown]
- # ### Dorsal/ventral location of reactive glia and disease-associated interneurons
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated_w_anno.h5ad"))
- # %%
- # Normalize y values and add to spatial data
- import pandas as pd
- # 1) pull out the y‑coordinate and slide_section (as string)
- y = pd.Series(adata.obsm["spatial"][:, 1], index=adata.obs.index, name="y")
- slide_sec = adata.obs["slide_section"].astype(str)
- # 2) shift each section so its minimum y becomes zero
- y_shifted = y.groupby(slide_sec).transform(lambda v: v - v.min())
- # 3) compute the median of the shifted values for each section
- y_median = y_shifted.groupby(slide_sec).transform("median")
- # 4) normalize by dividing by that median
- adata.obs["y_normalized"] = (y_shifted / y_median).values
- # %% [markdown]
- # ## Fig. 6D
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.stats import f_oneway
- from statsmodels.stats.multicomp import pairwise_tukeyhsd
- from statannotations.Annotator import Annotator
- # ——————————————————————————————————
- # Ensure all text is at least size 7
- # ——————————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Subset adata to only include "Lumbar" region
- adata_lumbar = adata[adata.obs["region"] == "Lumbar"].copy()
- # 1) Filter reactive microglia & drop Control
- stage_order = ["Early", "Mid", "End"]
- reactive = (
- adata_lumbar.obs
- .loc[
- adata_lumbar.obs["reactive_microglia"] == True,
- ["y_normalized", "stage", "slide_stage", "slide_section"]
- ]
- .copy()
- )
- reactive = reactive[reactive["stage"].isin(stage_order)]
- # 2) Per‑slide_section average of y_normalized
- ss_avg = (
- reactive
- .groupby(["slide_section", "slide_stage", "stage"], observed=True)
- .agg(avg_y=("y_normalized", "mean"))
- .reset_index()
- )
- # 3) Per‑slide_stage average of those slide_section avgs
- stg_avg = (
- ss_avg
- .groupby(["slide_stage", "stage"], observed=True)
- .agg(avg_y=("avg_y", "mean"))
- .reset_index()
- )
- # -------------------------
- # 4) STATISTICS
- # -------------------------
- # ANOVA across stages
- groups = [stg_avg.loc[stg_avg.stage == s, "avg_y"] for s in stage_order]
- F, p_an = f_oneway(*groups)
- print(f"ANOVA F = {F:.3f}, p = {p_an:.3e}")
- # Tukey HSD
- tukey = pairwise_tukeyhsd(
- endog=stg_avg["avg_y"],
- groups=stg_avg["stage"],
- alpha=0.05
- )
- print("\nTukey HSD results:")
- print(tukey.summary())
- # extract p‑values for annotator
- tuk_df = pd.DataFrame(tukey._results_table.data[1:], columns=tukey._results_table.data[0])
- pairs = [("Early", "Mid"), ("Early", "End"), ("Mid", "End")]
- pvals = []
- for a, b in pairs:
- m = tuk_df.query("(group1==@a & group2==@b) or (group1==@b & group2==@a)")
- pvals.append(m["p-adj"].values[0] if not m.empty else np.nan)
- # -------------------------
- # 5) PLOT
- # -------------------------
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- # boxplot with all boxes colored red
- sns.boxplot(
- data=stg_avg,
- x="stage",
- y="avg_y",
- order=stage_order,
- color="#d62728",
- showfliers=False,
- ax=ax
- )
- # stripplot on top
- sns.stripplot(
- data=stg_avg,
- x="stage",
- y="avg_y",
- order=stage_order,
- color="black",
- size=4,
- jitter=True,
- alpha=0.8,
- ax=ax
- )
- # add stat annotations
- annotator = Annotator(
- ax, pairs, data=stg_avg,
- x="stage", y="avg_y",
- order=stage_order
- )
- annotator.configure(test=None, text_format="star", loc="inside")
- annotator.set_pvalues(pvals)
- annotator.annotate()
- ax.set_xlabel("")
- ax.set_ylabel("Avg. Dorsal/Ventral Location")
- ax.set_title("Reactive MG - Lumbar")
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "reactive_microglia_DV_small.svg")
- fig.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Fig. 6J
- # %%
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.stats import ttest_rel
- from statannotations.Annotator import Annotator
- # ——————————————————————————————————
- # Ensure all text in the figure is size 7
- # ——————————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "figure.titlesize": 7
- })
- # Subset to Lumbar region AND exclude Control-stage
- adata_sub = adata[
- (adata.obs["region"] == "Lumbar") &
- (adata.obs["stage"] != "Control")
- ].copy()
- # Define the cell types (for filtering only)
- cell_types = {
- "Disease-Associated Interneurons": "Disease-Associated Interneurons",
- "Motor Neurons": ["Alpha MNs", "Gamma MNs", "Gamma* MNs"]
- }
- # 1) Filter each and tag a simplified label
- dai = (
- adata_sub.obs
- .loc[
- adata_sub.obs["cell_class"] == cell_types["Disease-Associated Interneurons"],
- ["y_normalized", "slide_section", "slide_stage", "stage"]
- ]
- .copy()
- )
- dai["cell_type"] = "DAI"
- mot = (
- adata_sub.obs
- .loc[
- adata_sub.obs["cholinergic_type"].isin(cell_types["Motor Neurons"]),
- ["y_normalized", "slide_section", "slide_stage", "stage"]
- ]
- .copy()
- )
- mot["cell_type"] = "Skeletal MNs"
- combined = pd.concat([dai, mot], ignore_index=True)
- # 2) Per‑slide_section average of y_normalized
- ss_avg = (
- combined
- .groupby(["slide_section", "slide_stage", "stage", "cell_type"], observed=True)
- .agg(avg_y=("y_normalized", "mean"))
- .reset_index()
- )
- # 3) Per‑slide_stage average of those slide_section avgs
- stg_lvl_avg = (
- ss_avg
- .groupby(["slide_stage", "stage", "cell_type"], observed=True)
- .agg(stage_avg_y=("avg_y", "mean"))
- .reset_index()
- )
- # 4) ONE‑SIDED PAIRED t‑TEST (H₁: DAI > Skeletal MNs)
- paired = (
- stg_lvl_avg
- .pivot(index="slide_stage", columns="cell_type", values="stage_avg_y")
- .dropna()
- )
- from scipy.stats import ttest_rel
- t_stat, p_two_sided = ttest_rel(paired["DAI"], paired["Skeletal MNs"])
- if t_stat > 0:
- p_one_sided = p_two_sided / 2
- else:
- p_one_sided = 1 - p_two_sided / 2
- print(f"paired one‑sided t = {t_stat:.3f}, p = {p_one_sided:.3e}")
- pairs = [("DAI", "Skeletal MNs")]
- pvals = [p_one_sided]
- # 5) PLOT using slide_stage averages with custom colors
- plt.figure(figsize=(1.625, 2.15))
- ax = sns.boxplot(
- data=stg_lvl_avg,
- x="cell_type",
- y="stage_avg_y",
- order=["DAI", "Skeletal MNs"],
- palette={"DAI": "#2ca02c", "Skeletal MNs": "#ff7f0e"},
- showfliers=False,
- ax=plt.gca()
- )
- sns.stripplot(
- data=stg_lvl_avg,
- x="cell_type",
- y="stage_avg_y",
- order=["DAI", "Skeletal MNs"],
- color="black",
- size=4,
- jitter=True,
- alpha=0.8,
- ax=ax
- )
- # connect paired points from same slide_stage
- for slide_stage in paired.index:
- ax.plot(
- ["DAI", "Skeletal MNs"],
- [paired.loc[slide_stage, "DAI"], paired.loc[slide_stage, "Skeletal MNs"]],
- color="gray", linewidth=0.8, alpha=0.6
- )
- # annotate
- annotator = Annotator(
- ax, [("DAI", "Skeletal MNs")], data=stg_lvl_avg,
- x="cell_type", y="stage_avg_y",
- order=["DAI", "Skeletal MNs"]
- )
- annotator.configure(test=None, text_format="star", loc="inside")
- annotator.set_pvalues([p_one_sided])
- annotator.annotate()
- ax.set_xlabel("")
- ax.set_ylabel("Avg. Dorsal/Ventral Location")
- ax.set_title("DAI vs. Skeletal MNs")
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "DAI_MN_DV_paired_small.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Reactive and non-reactive MG - distance to nearest MN
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated_w_anno.h5ad"))
- # %% [markdown]
- # ## Fig. 6E
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.spatial.distance import cdist
- from scipy.stats import ttest_rel
- from statannotations.Annotator import Annotator
- # ——————————————————————————————
- # Ensure all text is size 7
- # ——————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # ——————————————————————————————
- # 1) Subset to Lumbar & Early
- # ——————————————————————————————
- adata_sub = adata[
- (adata.obs["region"] == "Lumbar") &
- (adata.obs["stage"] == "Early")
- ].copy()
- # ——————————————————————————————
- # 2) Compute per-slide_stage avg distances
- # ——————————————————————————————
- records = []
- for slide_stage in adata_sub.obs["slide_stage"].unique():
- mask_ss = adata_sub.obs["slide_stage"] == slide_stage
- sections = adata_sub.obs.loc[mask_ss, "slide_section"].unique()
- for reactive_status in [True, False]:
- sec_means = []
- for sec in sections:
- mask_sec = mask_ss & (adata_sub.obs["slide_section"] == sec)
- coords = adata_sub.obsm["spatial"][mask_sec.values]
- is_motor = (
- adata_sub.obs.loc[mask_sec, "cholinergic_type"]
- .isin(["Alpha MNs","Gamma MNs","Gamma* MNs"])
- .to_numpy()
- )
- # restrict to microglia/macrophages first
- is_micro_class = (
- adata_sub.obs.loc[mask_sec, "cell_class"] == "Microglia/Macrophages"
- ).to_numpy()
- is_reactive_flag = (
- adata_sub.obs.loc[mask_sec, "reactive_microglia"]
- .astype(bool)
- .to_numpy()
- )
- is_micro = is_micro_class & (is_reactive_flag == reactive_status)
- pts_micro = coords[is_micro]
- pts_motor = coords[is_motor]
- if pts_micro.size and pts_motor.size:
- D = cdist(pts_micro, pts_motor, "euclidean")
- sec_means.append(D.min(axis=1).mean())
- if sec_means:
- records.append({
- "slide_stage": slide_stage,
- "reactive": reactive_status,
- "avg_dist": np.mean(sec_means)
- })
- dist_df = pd.DataFrame(records)
- # ——————————————————————————————
- # 3) Build plot_df
- # ——————————————————————————————
- plot_df = dist_df.copy()
- plot_df["cell_type"] = plot_df["reactive"].map({
- True: "Reactive",
- False: "Non-Reactive"
- })
- # ——————————————————————————————
- # 4) Compute one-sided p-value (Non-Reactive > Reactive)
- # ——————————————————————————————
- pivot = plot_df.pivot(
- index="slide_stage",
- columns="cell_type",
- values="avg_dist"
- ).dropna(subset=["Reactive", "Non-Reactive"])
- stat, p_two = ttest_rel(pivot["Non-Reactive"], pivot["Reactive"])
- p_one = p_two / 2 if stat > 0 else 1 - p_two / 2
- print(f"One-sided p-value (Non-Reactive > Reactive): {p_one:.3e}")
- # ——————————————————————————————
- # 5) Plot paired box+strip with lines and annotate
- # ——————————————————————————————
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- sns.boxplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive","Non-Reactive"],
- palette={"Reactive":"#d62728","Non-Reactive":"#ea9393"},
- showfliers=False,
- linewidth=1,
- ax=ax, zorder=1
- )
- sns.stripplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive","Non-Reactive"],
- color="black", size=4, jitter=True, alpha=0.8,
- ax=ax, zorder=2
- )
- for ss, grp in plot_df.groupby("slide_stage"):
- y0 = grp.loc[grp.cell_type=="Reactive", "avg_dist"].values
- y1 = grp.loc[grp.cell_type=="Non-Reactive","avg_dist"].values
- if y0.size and y1.size:
- ax.plot(
- ["Reactive","Non-Reactive"],
- [y0[0], y1[0]],
- color="gray", linewidth=0.8, alpha=0.6, zorder=3
- )
- pairs = [("Reactive", "Non-Reactive")]
- annot = Annotator(
- ax, pairs,
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive","Non-Reactive"]
- )
- annot.set_pvalues([p_one])
- annot.configure(text_format="star", loc="inside")
- annot.annotate()
- ax.set_title("Early Disease - Lumbar", fontsize=7)
- ax.set_xlabel("Microglia/Macrophages")
- ax.set_ylabel("Avg. Dist. to Nearest MN")
- # ——————————————————————————————
- # FORCE tick + spine thickness to 1
- # ——————————————————————————————
- ax.tick_params(width=1)
- for spine in ax.spines.values():
- spine.set_linewidth(1)
- plt.tight_layout()
- out_path = os.path.join(
- figures_dir,
- "reactive_microglia_nearest_MN_small_fixed.svg"
- )
- fig.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Reactive/WM and other Astrocytes - distance to nearest reactive MG
- # %%
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.spatial.distance import cdist
- from scipy.stats import ttest_rel
- from statannotations.Annotator import Annotator
- # ——————————————————————————————
- # 1) Subset to Lumbar & Early
- # ——————————————————————————————
- adata_sub = adata[
- (adata.obs["region"] == "Lumbar") &
- (adata.obs["stage"] == "Early")
- ].copy()
- # ——————————————————————————————
- # 2) Compute per-slide_stage avg distances
- # ——————————————————————————————
- records = []
- for slide_stage in adata_sub.obs["slide_stage"].unique():
- mask_ss = adata_sub.obs["slide_stage"] == slide_stage
- sections = adata_sub.obs.loc[mask_ss, "slide_section"].unique()
- for is_reactive_wm, label in [(True, "Reactive/WM"), (False, "Other")]:
- sec_means = []
- for sec in sections:
- mask_sec = mask_ss & (adata_sub.obs["slide_section"] == sec)
- coords = adata_sub.obsm["spatial"][mask_sec.values]
- # Query: astrocytes of this group
- is_astro = (
- (adata_sub.obs.loc[mask_sec, "cell_class"] == "Astrocytes") &
- (adata_sub.obs.loc[mask_sec, "reactive_or_wm_astrocytes"] == is_reactive_wm)
- ).to_numpy()
- # Reference: reactive microglia
- is_reactive_micro = (
- (adata_sub.obs.loc[mask_sec, "cell_class"] == "Microglia/Macrophages") &
- (adata_sub.obs.loc[mask_sec, "reactive_microglia"])
- ).to_numpy()
- pts_astro = coords[is_astro]
- pts_micro = coords[is_reactive_micro]
- if pts_astro.size and pts_micro.size:
- D = cdist(pts_astro, pts_micro, metric="euclidean")
- sec_means.append(D.min(axis=1).mean())
- if sec_means:
- records.append({
- "slide_stage": slide_stage,
- "astro_group": label,
- "avg_dist": np.mean(sec_means)
- })
- # ——————————————————————————————
- # 3) Create plot_df and perform stats
- # ——————————————————————————————
- dist_df = pd.DataFrame(records)
- plot_df = dist_df.copy()
- plot_df["cell_type"] = plot_df["astro_group"]
- # Paired one-sided t-test (Other > Reactive/WM)
- pivot = plot_df.pivot(index="slide_stage", columns="cell_type", values="avg_dist")
- stat, p_two = ttest_rel(pivot["Other"], pivot["Reactive/WM"])
- p_one = p_two / 2 if stat > 0 else 1 - p_two / 2
- print(f"One-sided p-value (Other > Reactive/WM): {p_one:.3e}")
- # ——————————————————————————————
- # 4) Plot
- # ——————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7
- })
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- sns.boxplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"],
- palette={"Reactive/WM": "#1f77b4", "Other": "#8fbbd9"},
- showfliers=False,
- ax=ax, zorder=1
- )
- sns.stripplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"],
- color="black", size=4, jitter=True, alpha=0.8,
- ax=ax, zorder=2
- )
- # Paired lines
- for ss, grp in plot_df.groupby("slide_stage"):
- y_reactive = grp.loc[grp.cell_type == "Reactive/WM", "avg_dist"].values
- y_other = grp.loc[grp.cell_type == "Other", "avg_dist"].values
- if y_reactive.size and y_other.size:
- ax.plot(
- ["Reactive/WM", "Other"],
- [y_reactive[0], y_other[0]],
- color="gray", linewidth=0.8, alpha=0.6, zorder=3
- )
- # Annotate p-value
- pairs = [("Reactive/WM", "Other")]
- annot = Annotator(
- ax, pairs,
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"]
- )
- annot.set_pvalues([p_one])
- annot.configure(text_format="star", loc="inside")
- annot.annotate()
- # Labels
- ax.set_title("Early Disease - Lumbar")
- ax.set_xlabel("Astrocytes")
- ax.set_ylabel("Avg. Dist. to Nearest Reactive MG")
- plt.tight_layout()
- fig.savefig(os.path.join(figures_dir, "astro_dist_to_reactive_microglia_early.svg"), format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Fig. 6F
- # %%
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.spatial.distance import cdist
- from scipy.stats import ttest_rel
- from statannotations.Annotator import Annotator
- # ——————————————————————————————
- # 1) Subset to Lumbar & Mid
- # ——————————————————————————————
- adata_sub = adata[
- (adata.obs["region"] == "Lumbar") &
- (adata.obs["stage"] == "Mid")
- ].copy()
- # ——————————————————————————————
- # 2) Compute per-slide_stage avg distances
- # ——————————————————————————————
- records = []
- for slide_stage in adata_sub.obs["slide_stage"].unique():
- mask_ss = adata_sub.obs["slide_stage"] == slide_stage
- sections = adata_sub.obs.loc[mask_ss, "slide_section"].unique()
- for is_reactive_wm, label in [(True, "Reactive/WM"), (False, "Other")]:
- sec_means = []
- for sec in sections:
- mask_sec = mask_ss & (adata_sub.obs["slide_section"] == sec)
- coords = adata_sub.obsm["spatial"][mask_sec.values]
- # Query: astrocytes of this group
- is_astro = (
- (adata_sub.obs.loc[mask_sec, "cell_class"] == "Astrocytes") &
- (adata_sub.obs.loc[mask_sec, "reactive_or_wm_astrocytes"] == is_reactive_wm)
- ).to_numpy()
- # Reference: reactive microglia
- is_reactive_micro = (
- (adata_sub.obs.loc[mask_sec, "cell_class"] == "Microglia/Macrophages") &
- (adata_sub.obs.loc[mask_sec, "reactive_microglia"])
- ).to_numpy()
- pts_astro = coords[is_astro]
- pts_micro = coords[is_reactive_micro]
- if pts_astro.size and pts_micro.size:
- D = cdist(pts_astro, pts_micro, metric="euclidean")
- sec_means.append(D.min(axis=1).mean())
- if sec_means:
- records.append({
- "slide_stage": slide_stage,
- "astro_group": label,
- "avg_dist": np.mean(sec_means)
- })
- # ——————————————————————————————
- # 3) Create plot_df and perform stats
- # ——————————————————————————————
- dist_df = pd.DataFrame(records)
- plot_df = dist_df.copy()
- plot_df["cell_type"] = plot_df["astro_group"]
- # Paired one-sided t-test (Other > Reactive/WM)
- pivot = plot_df.pivot(index="slide_stage", columns="cell_type", values="avg_dist")
- stat, p_two = ttest_rel(pivot["Other"], pivot["Reactive/WM"])
- p_one = p_two / 2 if stat > 0 else 1 - p_two / 2
- print(f"One-sided p-value (Other > Reactive/WM): {p_one:.3e}")
- # ——————————————————————————————
- # 4) Plot
- # ——————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7
- })
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- sns.boxplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"],
- palette={"Reactive/WM": "#1f77b4", "Other": "#8fbbd9"},
- showfliers=False,
- ax=ax, zorder=1
- )
- sns.stripplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"],
- color="black", size=4, jitter=True, alpha=0.8,
- ax=ax, zorder=2
- )
- # Paired lines
- for ss, grp in plot_df.groupby("slide_stage"):
- y_reactive = grp.loc[grp.cell_type == "Reactive/WM", "avg_dist"].values
- y_other = grp.loc[grp.cell_type == "Other", "avg_dist"].values
- if y_reactive.size and y_other.size:
- ax.plot(
- ["Reactive/WM", "Other"],
- [y_reactive[0], y_other[0]],
- color="gray", linewidth=0.8, alpha=0.6, zorder=3
- )
- # Annotate p-value
- pairs = [("Reactive/WM", "Other")]
- annot = Annotator(
- ax, pairs,
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"]
- )
- annot.set_pvalues([p_one])
- annot.configure(text_format="star", loc="inside")
- annot.annotate()
- # Labels
- ax.set_title("Mid Disease - Lumbar")
- ax.set_xlabel("Astrocytes")
- ax.set_ylabel("Avg. Dist. to Nearest Reactive MG")
- plt.tight_layout()
- fig.savefig(os.path.join(figures_dir, "astro_dist_to_reactive_microglia_mid.svg"), format="svg", bbox_inches="tight")
- plt.show()
- # %%
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.spatial.distance import cdist
- from scipy.stats import ttest_rel
- from statannotations.Annotator import Annotator
- # ——————————————————————————————
- # 1) Subset to Lumbar & End
- # ——————————————————————————————
- adata_sub = adata[
- (adata.obs["region"] == "Lumbar") &
- (adata.obs["stage"] == "End")
- ].copy()
- # ——————————————————————————————
- # 2) Compute per-slide_stage avg distances
- # ——————————————————————————————
- records = []
- for slide_stage in adata_sub.obs["slide_stage"].unique():
- mask_ss = adata_sub.obs["slide_stage"] == slide_stage
- sections = adata_sub.obs.loc[mask_ss, "slide_section"].unique()
- for is_reactive_wm, label in [(True, "Reactive/WM"), (False, "Other")]:
- sec_means = []
- for sec in sections:
- mask_sec = mask_ss & (adata_sub.obs["slide_section"] == sec)
- coords = adata_sub.obsm["spatial"][mask_sec.values]
- # Query: astrocytes of this group
- is_astro = (
- (adata_sub.obs.loc[mask_sec, "cell_class"] == "Astrocytes") &
- (adata_sub.obs.loc[mask_sec, "reactive_or_wm_astrocytes"] == is_reactive_wm)
- ).to_numpy()
- # Reference: reactive microglia
- is_reactive_micro = (
- (adata_sub.obs.loc[mask_sec, "cell_class"] == "Microglia/Macrophages") &
- (adata_sub.obs.loc[mask_sec, "reactive_microglia"])
- ).to_numpy()
- pts_astro = coords[is_astro]
- pts_micro = coords[is_reactive_micro]
- if pts_astro.size and pts_micro.size:
- D = cdist(pts_astro, pts_micro, metric="euclidean")
- sec_means.append(D.min(axis=1).mean())
- if sec_means:
- records.append({
- "slide_stage": slide_stage,
- "astro_group": label,
- "avg_dist": np.mean(sec_means)
- })
- # ——————————————————————————————
- # 3) Create plot_df and perform stats
- # ——————————————————————————————
- dist_df = pd.DataFrame(records)
- plot_df = dist_df.copy()
- plot_df["cell_type"] = plot_df["astro_group"]
- # Paired one-sided t-test (Other > Reactive/WM)
- pivot = plot_df.pivot(index="slide_stage", columns="cell_type", values="avg_dist")
- stat, p_two = ttest_rel(pivot["Other"], pivot["Reactive/WM"])
- p_one = p_two / 2 if stat > 0 else 1 - p_two / 2
- print(f"One-sided p-value (Other > Reactive/WM): {p_one:.3e}")
- # ——————————————————————————————
- # 4) Plot
- # ——————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7
- })
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- sns.boxplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"],
- palette={"Reactive/WM": "#1f77b4", "Other": "#8fbbd9"},
- showfliers=False,
- ax=ax, zorder=1
- )
- sns.stripplot(
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"],
- color="black", size=4, jitter=True, alpha=0.8,
- ax=ax, zorder=2
- )
- # Paired lines
- for ss, grp in plot_df.groupby("slide_stage"):
- y_reactive = grp.loc[grp.cell_type == "Reactive/WM", "avg_dist"].values
- y_other = grp.loc[grp.cell_type == "Other", "avg_dist"].values
- if y_reactive.size and y_other.size:
- ax.plot(
- ["Reactive/WM", "Other"],
- [y_reactive[0], y_other[0]],
- color="gray", linewidth=0.8, alpha=0.6, zorder=3
- )
- # Annotate p-value
- pairs = [("Reactive/WM", "Other")]
- annot = Annotator(
- ax, pairs,
- data=plot_df,
- x="cell_type", y="avg_dist",
- order=["Reactive/WM", "Other"]
- )
- annot.set_pvalues([p_one])
- annot.configure(text_format="star", loc="inside")
- annot.annotate()
- # Labels
- ax.set_title("End Disease - Lumbar")
- ax.set_xlabel("Astrocytes")
- ax.set_ylabel("Avg. Dist. to Nearest Reactive MG")
- # Ensure y-axis ticks are integers only
- ax.yaxis.set_major_locator(plt.MaxNLocator(integer=True))
- plt.tight_layout()
- fig.savefig(os.path.join(figures_dir, "astro_dist_to_reactive_microglia_end.svg"), format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Reactive glia, DAI, and DAMN proportions with disease
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_rotated_w_anno.h5ad"))
- # %% [markdown]
- # ## Fig. S7D
- # %%
- import matplotlib.pyplot as plt
- import seaborn as sns
- # 0) Set all fonts to size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Step 1: Subset to Microglia/Macrophages only
- mg_obs = adata.obs[adata.obs["cell_class"] == "Microglia/Macrophages"].copy()
- # Step 2: Count reactive vs non‑reactive by stage
- stage_reactive_counts = (
- mg_obs
- .groupby(["stage", "reactive_microglia"], observed=True)
- .size()
- .unstack(fill_value=0)
- )
- # Step 3: Convert to proportions within each stage
- stage_reactive_props = stage_reactive_counts.div(stage_reactive_counts.sum(axis=1), axis=0)
- # Step 4: Ensure stage order
- stage_order = ["Control", "Early", "Mid", "End"]
- stage_reactive_props = stage_reactive_props.reindex(stage_order).fillna(0)
- # Step 5: Force True first (bottom) then False (top)
- stage_reactive_props = stage_reactive_props.reindex(columns=[True, False])
- # Step 6: Plot
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- stage_reactive_props.plot(
- kind="bar",
- stacked=True,
- color=["#d62728", "#ea9393"], # True (bottom), False (top)
- edgecolor="black",
- ax=ax
- )
- # Step 7: Labels, legend, title
- ax.set_xlabel("") # no x-axis label
- ax.set_ylabel("Proportion")
- ax.set_xticklabels(stage_order, rotation=0)
- ax.legend(
- title="Reactive",
- labels=["True", "False"],
- loc="upper center",
- bbox_to_anchor=(0.5, -0.12),
- ncol=2,
- handlelength=1, # shorten the legend color box length
- handletextpad=0.4, # reduce space between color box and text
- columnspacing=0.8, # reduce space between columns
- borderpad=0.3 # reduce padding inside the legend box
- )
- ax.set_title("Microglia/Macrophages")
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "reactive_microglia_proportion.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Fig. S7E
- # %%
- import matplotlib.pyplot as plt
- import seaborn as sns
- # 0) Set all fonts to size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Step 1: Subset to Astrocytes only
- ast_obs = adata.obs[adata.obs["cell_class"] == "Astrocytes"].copy()
- # Step 2: Count reactive vs non‑reactive by stage
- stage_reactive_counts = (
- ast_obs
- .groupby(["stage", "reactive_or_wm_astrocytes"], observed=True)
- .size()
- .unstack(fill_value=0)
- )
- # Step 3: Convert to proportions within each stage
- stage_reactive_props = stage_reactive_counts.div(stage_reactive_counts.sum(axis=1), axis=0)
- # Step 4: Ensure stage order
- stage_order = ["Control", "Early", "Mid", "End"]
- stage_reactive_props = stage_reactive_props.reindex(stage_order).fillna(0)
- # Step 5: Force True first (bottom) then False (top)
- stage_reactive_props = stage_reactive_props.reindex(columns=[True, False])
- # Step 6: Plot
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- stage_reactive_props.plot(
- kind="bar",
- stacked=True,
- color=["#1f77b4", "#8fbbd9"], # True (bottom), False (top)
- edgecolor="black",
- ax=ax
- )
- # Step 7: Labels, legend, title
- ax.set_xlabel("") # no x-axis label
- ax.set_ylabel("Proportion")
- ax.set_xticklabels(stage_order, rotation=0)
- ax.legend(
- title="Reactive",
- labels=["True", "False"],
- loc="upper center",
- bbox_to_anchor=(0.5, -0.12),
- ncol=2,
- handlelength=1, # shorten the legend color box length
- handletextpad=0.4, # reduce space between color box and text
- columnspacing=0.8, # reduce space between columns
- borderpad=0.3 # reduce padding inside the legend box
- )
- ax.set_title("Astrocytes")
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "reactive_astrocyte_proportion.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Fig. S7F
- # %%
- import matplotlib.pyplot as plt
- import seaborn as sns
- # 0) Set all fonts to size 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # Step 1: Subset to Alpha MNs only
- alpha_obs = adata.obs[adata.obs["cholinergic_type"] == "Alpha MNs"].copy()
- # Step 2: Count DAMN vs non‑DAMN by stage
- stage_DAMN_counts = (
- alpha_obs
- .groupby(["stage", "DAMN_status"], observed=True)
- .size()
- .unstack(fill_value=0)
- )
- # Step 3: Convert to proportions within each stage
- stage_DAMN_props = stage_DAMN_counts.div(stage_DAMN_counts.sum(axis=1), axis=0)
- # Step 4: Ensure stage order
- stage_order = ["Control", "Early", "Mid", "End"]
- stage_DAMN_props = stage_DAMN_props.reindex(stage_order).fillna(0)
- # 5) reorder to DAMN bottom, Non-DAMN top
- stage_DAMN_props = stage_DAMN_props[["DAMN","Non-DAMN"]].fillna(0)
- # 6) plot
- fig, ax = plt.subplots(figsize=(1.625, 2.15))
- stage_DAMN_props.plot(
- kind="bar",
- stacked=True,
- color=["#1f77b4", "#8fbbd9"], # bottom=“DAMN”, top=“Non-DAMN”
- edgecolor="black",
- ax=ax
- )
- # Step 7: Labels, legend, title
- ax.set_xlabel("") # no x-axis label
- ax.set_ylabel("Proportion")
- ax.set_xticklabels(stage_order, rotation=0)
- ax.legend(
- title="DM",
- labels=["True", "False"],
- loc="upper center",
- bbox_to_anchor=(0.5, -0.12),
- ncol=2,
- handlelength=1, # shorten the legend color box length
- handletextpad=0.4, # reduce space between color box and text
- columnspacing=0.8, # reduce space between columns
- borderpad=0.3 # reduce padding inside the legend box
- )
- ax.set_title("Alpha MNs")
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "DM_proportion.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %%
- print("\nProportion DAMN_status == 'DAMN' by stage (percent):")
- for stage, prop in stage_DAMN_props["DAMN"].items():
- print(f"{stage}: {prop:.1%}")
- # %% [markdown]
- # ## Alpha MN subtype proportion plot
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata_label_transfer.h5ad"))
- # %% [markdown]
- # ## Fig. S7G
- # %%
- import pandas as pd
- import seaborn as sns
- import matplotlib.pyplot as plt
- import os
- # ------------------------------------------------------------------
- # 0) Set all fonts to size 7 (match your other figure)
- # ------------------------------------------------------------------
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # ------------------------------------------------------------------
- # 1) Subset to nuclear cells (optional)
- # ------------------------------------------------------------------
- adata_sub = adata # (or: adata[adata.obs["nuclear"] == True])
- # ------------------------------------------------------------------
- # 2) Extract relevant metadata (now includes region)
- # ------------------------------------------------------------------
- df = adata_sub.obs[["stage", "predicted.id", "region"]].copy()
- df = df.dropna(subset=["stage", "predicted.id", "region"])
- stage_order = ["Control", "Early", "Mid", "End"]
- pred_order = ["Slow-Firing", "Intermediate", "Fast-Firing", "Early DAMN", "Late DAMN"]
- pred_order_stack = list(reversed(pred_order))
- # Display-only legend labels
- label_map = {"Early DAMN": "Early DM", "Late DAMN": "Late DM"}
- df["stage"] = pd.Categorical(df["stage"], categories=stage_order, ordered=True)
- df["predicted.id"] = pd.Categorical(df["predicted.id"], categories=pred_order, ordered=True)
- # ------------------------------------------------------------------
- # 3) Define colors
- # ------------------------------------------------------------------
- alpha_colors = {
- "Slow-Firing": "#F47D2B",
- "Intermediate": "#272E6A",
- "Fast-Firing": "#208A42",
- "Early DAMN": "#89288F",
- "Late DAMN": "#D51F26"
- }
- # ------------------------------------------------------------------
- # 4) Helper: compute proportions for a given region
- # ------------------------------------------------------------------
- def compute_prop_for_region(df_in: pd.DataFrame, region_name: str) -> pd.DataFrame:
- d = df_in[df_in["region"] == region_name].copy()
- counts = d.groupby(["stage", "predicted.id"]).size().to_frame("count")
- counts["proportion"] = counts["count"] / counts.groupby(level=0)["count"].transform("sum")
- prop = counts.reset_index()
- # Ensure all stage x predicted.id combos exist (so bars stack cleanly even if missing)
- full = (
- pd.MultiIndex.from_product([stage_order, pred_order], names=["stage", "predicted.id"])
- .to_frame(index=False)
- )
- prop = full.merge(prop, on=["stage", "predicted.id"], how="left").fillna({"proportion": 0, "count": 0})
- prop["stage"] = pd.Categorical(prop["stage"], categories=stage_order, ordered=True)
- prop["predicted.id"] = pd.Categorical(prop["predicted.id"], categories=pred_order, ordered=True)
- return prop
- # ------------------------------------------------------------------
- # 5) Plot: two stacked barplots (Lumbar + Cervical)
- # ------------------------------------------------------------------
- regions = ["Lumbar", "Cervical"]
- fig, axes = plt.subplots(
- 1, 2,
- figsize=(3.25, 1.625), # ~2x width of your single panel (1.625 -> 3.25)
- sharey=True
- )
- legend_handles = None # we will grab once and use a single legend for the figure
- for ax, region_name in zip(axes, regions):
- prop = compute_prop_for_region(df, region_name)
- bottom = pd.Series([0.0] * len(stage_order), index=stage_order)
- handles = {}
- for pred in pred_order_stack:
- sub = prop[prop["predicted.id"] == pred].sort_values("stage")
- bars = ax.bar(
- sub["stage"].astype(str),
- sub["proportion"].values,
- bottom=bottom[sub["stage"].astype(str)].values if bottom.index.dtype == object else bottom[sub["stage"]].values,
- color=alpha_colors[pred],
- width=0.8,
- edgecolor="none"
- )
- handles[pred] = bars[0]
- # Update bottom by stage order
- bottom.loc[sub["stage"].astype(str).values] = (
- bottom.loc[sub["stage"].astype(str).values].values + sub["proportion"].values
- )
- ax.set_title(region_name)
- ax.set_xlabel("")
- ax.set_ylim(0, 1)
- ax.set_xticklabels(stage_order, rotation=0)
- sns.despine(ax=ax)
- if legend_handles is None:
- legend_handles = handles # save for figure-level legend
- # Y label only on left
- axes[0].set_ylabel("Proportion")
- axes[1].set_ylabel("")
- # One legend for the whole figure (no title)
- fig.legend(
- handles=[legend_handles[p] for p in pred_order],
- labels=[label_map.get(p, p) for p in pred_order],
- bbox_to_anchor=(1.02, 1),
- loc="upper left",
- frameon=False
- )
- plt.tight_layout()
- # ------------------------------------------------------------------
- # 6) Save figure
- # ------------------------------------------------------------------
- out_path = os.path.join(figures_dir, "alpha_mn_subtype_barplot_Lumbar_Cervical_DM.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %% [markdown]
- # ## Fig. S7H
- # %%
- import scanpy as sc
- import matplotlib.pyplot as plt
- import pandas as pd
- import os
- from matplotlib.collections import PathCollection
- # ---- Enforce predicted.id order ----
- pred_order = ["Slow-Firing", "Intermediate", "Fast-Firing", "Early DAMN", "Late DAMN"]
- adata.obs["predicted.id"] = pd.Categorical(
- adata.obs["predicted.id"],
- categories=pred_order,
- ordered=True
- )
- # ---- Display-only labels for plotting ----
- label_map = {"Early DAMN": "Early DM", "Late DAMN": "Late DM"}
- pred_order_plot = [label_map.get(c, c) for c in pred_order]
- adata.obs["predicted.id_plot"] = adata.obs["predicted.id"].cat.rename_categories(label_map)
- adata.obs["predicted.id_plot"] = pd.Categorical(
- adata.obs["predicted.id_plot"],
- categories=pred_order_plot,
- ordered=True
- )
- # ---- Colors ----
- alpha_colors = {
- "Slow-Firing": "#F47D2B",
- "Intermediate": "#272E6A",
- "Fast-Firing": "#208A42",
- "Early DM": "#89288F",
- "Late DM": "#D51F26"
- }
- adata.uns["predicted.id_plot_colors"] = [alpha_colors[c] for c in pred_order_plot]
- # Global font size = 7
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 6
- })
- # Create UMAP and return figure
- fig = sc.pl.umap(
- adata,
- color="predicted.id_plot",
- title="",
- size=5,
- show=False,
- return_fig=True
- )
- # Hard-remove any lingering title
- for ax in fig.axes:
- ax.set_title("")
- # Set exact figure size AFTER creation
- fig.set_size_inches(2.3, 1.525)
- # ---- Force VECTOR points: turn OFF rasterization on the scatter collection(s)
- for ax in fig.axes:
- for coll in ax.collections:
- if isinstance(coll, PathCollection):
- coll.set_rasterized(False)
- plt.tight_layout()
- plt.show()
- out_path = os.path.join(figures_dir, "alpha_mn_subtype_umap_DM.svg")
- fig.savefig(out_path, format="svg", bbox_inches="tight")
- # %% [markdown]
- # ## Fig. 6I
- # %%
- import numpy as np
- import matplotlib.pyplot as plt
- import seaborn as sns
- # ——————————————————————————————————
- # Ensure all text is size 7
- # ——————————————————————————————————
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7
- })
- # 1) Define the three groups in bottom→top stacking order
- group_names = [
- "Disease-Associated Interneurons", # bottom
- "Cholinergic Interneurons", # middle
- "Non-Cholinergic Interneurons", # top
- ]
- # 2) Filter and tag each cell into one of those groups
- mask = (
- (adata.obs["cell_class"] == "Disease-Associated Interneurons") |
- (adata.obs["cholinergic_type"] == "Cholinergic Interneurons") |
- (adata.obs["cell_class"] == "Non-Cholinergic Interneurons")
- )
- df = adata.obs.loc[mask].copy()
- df["cell_group"] = np.select(
- [
- df["cell_class"] == "Disease-Associated Interneurons",
- df["cholinergic_type"] == "Cholinergic Interneurons",
- df["cell_class"] == "Non-Cholinergic Interneurons",
- ],
- group_names
- )
- # 3) Count per (stage, cell_group)
- stage_order = ["Control", "Early", "Mid", "End"]
- counts = (
- df
- .groupby(["stage", "cell_group"], observed=True)
- .size()
- .unstack(fill_value=0)
- .reindex(stage_order)
- .fillna(0)
- )
- # 4) Proportions within each stage
- props = counts.div(counts.sum(axis=1), axis=0)
- # 5) Reorder columns to match stacking bottom→top
- props = props[group_names]
- # 6) Define colors in the same bottom→top order
- colors = [
- "#2ca02c", # DAI (bottom)
- "#ff7f0e", # CI (middle)
- "#9467bd", # NCI (top)
- ]
- # 7) Plot
- fig, ax = plt.subplots(figsize=(1.625, 2.2852))
- props.plot(
- kind="bar",
- stacked=True,
- color=colors,
- edgecolor="black",
- ax=ax
- )
- # 8) Labels & legend
- ax.set_xlabel("")
- ax.set_ylabel("Proportion")
- ax.set_xticklabels(stage_order, rotation=0)
- ax.set_title("Interneurons")
- short_labels = ["DAI", "CI", "NCI"]
- ax.legend(
- title="",
- labels=short_labels,
- loc="lower center",
- bbox_to_anchor=(0.5, -0.35),
- ncol=3,
- columnspacing=0.5, # ↓ Reduce space between columns
- handletextpad=0.3, # ↓ Reduce space between handle and text
- borderpad=0.2, # ↓ Reduce border padding inside legend box
- handlelength=1.0 # ↓ Optional: reduce legend key size
- )
- plt.tight_layout()
- out_path = os.path.join(figures_dir, "interneuron_proportions_small.svg")
- plt.savefig(out_path, format="svg", bbox_inches="tight")
- plt.show()
- # %%
- print("\nProportions of each interneuron type by stage (percent):")
- for stage in stage_order:
- print(f"{stage}:")
- for group in group_names:
- pct = props.loc[stage, group]
- print(f" {group}: {pct:.1%}")
- # %% [markdown]
- # ## Add additional metadata to all_anndata_final.h5ad
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final.h5ad"))
- adata_w_anno = sc.read_h5ad(os.path.join(working_dir, "adata_objects/all_anndata_final_rotated_w_anno.h5ad"))
- # %%
- adata.obs["predicted.id"] = (
- adata.obs["predicted.id"]
- .astype("category")
- .cat.add_categories(["Not_Alpha_MN"])
- .fillna("Not_Alpha_MN")
- )
- adata.obs["DAMN_status"] = (
- adata.obs["DAMN_status"]
- .astype("category")
- .cat.add_categories(["Not_Alpha_MN"])
- .fillna("Not_Alpha_MN")
- )
- # %%
- cols_to_add = ["reactive_microglia", "reactive_or_wm_astrocytes"]
- adata.obs[cols_to_add] = adata_w_anno.obs[cols_to_add].to_numpy()
- # %%
- # Make sure the needed columns exist
- need = ["cell_class", "reactive_microglia", "reactive_or_wm_astrocytes"]
- missing = [c for c in need if c not in adata.obs.columns]
- if missing:
- raise ValueError(f"Missing required columns in adata.obs: {missing}")
- # Start with default
- adata.obs["cell_class_reactivity"] = "Other"
- # Cholinergic neurons
- mask_chol = adata.obs["cell_class"] == "Cholinergic Neurons"
- adata.obs.loc[mask_chol, "cell_class_reactivity"] = "Cholinergic Neurons"
- # Microglia/Macrophages → Reactive vs Non-Reactive
- mask_mg = adata.obs["cell_class"] == "Microglia/Macrophages"
- adata.obs.loc[mask_mg & (adata.obs["reactive_microglia"] == True), "cell_class_reactivity"] = "Reactive MG"
- adata.obs.loc[mask_mg & (adata.obs["reactive_microglia"] == False), "cell_class_reactivity"] = "Non-Reactive MG"
- # Astrocytes → Reactive vs Non-Reactive
- mask_ast = adata.obs["cell_class"] == "Astrocytes"
- adata.obs.loc[mask_ast & (adata.obs["reactive_or_wm_astrocytes"] == True), "cell_class_reactivity"] = "Reactive Astrocytes"
- adata.obs.loc[mask_ast & (adata.obs["reactive_or_wm_astrocytes"] == False), "cell_class_reactivity"] = "Non-Reactive Astrocytes"
- # %%
- # -------------------------------------------------
- # Create combined DAMN / Reactivity annotation
- # -------------------------------------------------
- # Make sure required columns exist
- need = ["DAMN_status", "cell_class_reactivity"]
- missing = [c for c in need if c not in adata.obs.columns]
- if missing:
- raise ValueError(f"Missing required columns in adata.obs: {missing}")
- # Initialize with default
- adata.obs["DAMN_reactivity_class"] = "Other"
- # DAMN motor neurons
- adata.obs.loc[
- adata.obs["DAMN_status"] == "DAMN",
- "DAMN_reactivity_class"
- ] = "DAMN"
- # Non-DAMN motor neurons
- adata.obs.loc[
- adata.obs["DAMN_status"] == "Non-DAMN",
- "DAMN_reactivity_class"
- ] = "Non-DAMN"
- # Reactive Microglia
- adata.obs.loc[
- adata.obs["cell_class_reactivity"] == "Reactive MG",
- "DAMN_reactivity_class"
- ] = "Reactive MG"
- # Non-Reactive Microglia
- adata.obs.loc[
- adata.obs["cell_class_reactivity"] == "Non-Reactive MG",
- "DAMN_reactivity_class"
- ] = "Non-Reactive MG"
- # Convert to ordered categorical (optional but recommended for Vizualizer)
- adata.obs["DAMN_reactivity_class"] = pd.Categorical(
- adata.obs["DAMN_reactivity_class"],
- categories=["DAMN", "Non-DAMN", "Reactive MG", "Non-Reactive MG", "Other"],
- ordered=True
- )
- # %%
- # Save anndata
- save_name = f"adata_objects/all_anndata_final_vizualizer.h5ad"
- adata.write_h5ad(os.path.join(working_dir, save_name))
- # %% [markdown]
- # ## Save files to use with MERSCOPE Vizualizer
- # %%
- # Paths
- input_file = os.path.join(working_dir, "adata_objects/all_anndata_final_vizualizer.h5ad")
- output_dir = os.path.join(working_dir, "adata_objects_vizualizer")
- # load full AnnData
- adata = sc.read_h5ad(input_file)
- # get all unique slide IDs
- slide_ids = adata.obs["slide"].unique()
- # loop through and save one file per slide
- for slide_id in slide_ids:
- # subset
- ad = adata[adata.obs["slide"] == slide_id].copy()
- # construct filename
- # e.g. "202405281416_MsSpinalCord-...-VMSC15302.hdf5"
- fname = f"{slide_id}.hdf5"
- out_path = os.path.join(output_dir, fname)
- # write
- ad.write_h5ad(out_path)
- print(f"Wrote {out_path}")
- # %% [markdown]
- # ## MN segmentations from parquet files
- # %%
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/all_anndata_final_vizualizer.h5ad"))
- # %%
- import numpy as np
- import matplotlib.pyplot as plt
- import seaborn as sns
- from skimage.filters import threshold_otsu
- # -----------------------------
- # Restrict to Alpha MNs
- # -----------------------------
- alpha_mask = adata.obs["cholinergic_type"] == "Alpha MNs"
- # Pull DAPI values for Alpha MNs only
- vals = adata.obs.loc[alpha_mask, "DAPI_high_pass"].values
- log_vals = np.log10(vals + 1)
- # -----------------------------
- # Otsu threshold (Alpha MNs only)
- # -----------------------------
- thresh = threshold_otsu(log_vals)
- print("Otsu log10(DAPI) threshold (Alpha MNs) =", thresh)
- # Initialize column (False by default)
- adata.obs["nuclear_alphaMN"] = False
- # Assign nuclear status ONLY for Alpha MNs
- adata.obs.loc[alpha_mask, "nuclear_alphaMN"] = (
- np.log10(adata.obs.loc[alpha_mask, "DAPI_high_pass"].values + 1) >= thresh
- )
- # -----------------------------
- # QC plot
- # -----------------------------
- plt.figure(figsize=(6,4))
- sns.histplot(log_vals, bins=300, color="steelblue", edgecolor=None)
- plt.axvline(thresh, color="red", linestyle="--", label=f"Otsu={thresh:.2f}")
- plt.xlabel("log10(DAPI_high_pass + 1)")
- plt.ylabel("Alpha MN count")
- plt.legend()
- plt.tight_layout()
- plt.show()
- # %% [markdown]
- # ## Fig. S9A
- # %%
- import os, glob
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- from shapely import wkb
- from scipy.spatial import cKDTree
- from matplotlib.patches import Polygon as MplPolygon
- # -----------------------------
- # Paths / inputs
- # -----------------------------
- base = "/home/users/ogautier/oak/Shared/SOD1_Paper/Vizgen/mn_nonmn_segmentation/202401241200_HuSpinalCord-VS119-Lubmar-Ctrl3-SOD10_VMSC10802/analysis_outputs_7000_900"
- tile_files = sorted(glob.glob(os.path.join(base, "result_tiles", "cell_*.parquet")))
- FIG_DIR = "/home/users/ogautier/oak/Shared/SOD1_Paper/Vizgen/mn_nonmn_segmentation/figures"
- os.makedirs(FIG_DIR, exist_ok=True)
- # -----------------------------
- # Filters / columns
- # -----------------------------
- slide_section_col = "slide_section"
- mn_filter_col = "cholinergic_type"
- mn_filter_value = "Alpha MNs"
- nuclear_col = "nuclear_alphaMN"
- nuclear_value = True
- color_by_col = "predicted.id"
- tol = 1.7
- # -----------------------------
- # Plot settings
- # -----------------------------
- plt.rcParams.update({
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "font.family": "DejaVu Sans",
- "svg.fonttype": "path",
- })
- FIG_W_IN, FIG_H_IN = 4.875, 3.5
- LEFT, RIGHT, BOTTOM, TOP = 0.06, 0.98, 0.18, 0.98
- WSPACE, HSPACE = -0.035, -0.035
- SPINE_LW = 0.9
- BOX_SIZE_UM = 600
- HALF = BOX_SIZE_UM / 2.0
- LEGEND_ORDER = ["Slow-Firing", "Intermediate", "Fast-Firing", "Early DAMN", "Late DAMN"]
- LEGEND_LABEL_MAP = {"Early DAMN": "Early DM", "Late DAMN": "Late DM"}
- COLOR_MAP = {
- "Slow-Firing": "#f47d2b",
- "Intermediate": "#272e6a",
- "Fast-Firing": "#208a42",
- "Early DAMN": "#89288f",
- "Late DAMN": "#d51f26"
- }
- # -----------------------------
- # Checks
- # -----------------------------
- for col in [slide_section_col, mn_filter_col, nuclear_col, color_by_col]:
- if col not in adata.obs.columns:
- raise ValueError(f"adata.obs missing '{col}'")
- if "spatial" not in adata.obsm.keys():
- raise ValueError("adata.obsm missing 'spatial'")
- sp = np.asarray(adata.obsm["spatial"])
- tree = cKDTree(sp)
- # -----------------------------
- # Helpers
- # -----------------------------
- def is_R_prefix(adata_slide_section_vals, prefix):
- suffix = pd.Series(adata_slide_section_vals.astype(str)).str.split("_").str[-1]
- return suffix.str.startswith(prefix).values
- def nuclear_alpha_mask(adata):
- return (
- (adata.obs[mn_filter_col].astype(str) == mn_filter_value) &
- (adata.obs[nuclear_col].astype(bool) == nuclear_value)
- )
- def compute_centroids_from_wkb_bytes(wkb_bytes_arr):
- cx = np.empty(len(wkb_bytes_arr), dtype=float)
- cy = np.empty(len(wkb_bytes_arr), dtype=float)
- ok = np.ones(len(wkb_bytes_arr), dtype=bool)
- for i, b in enumerate(wkb_bytes_arr):
- try:
- g = wkb.loads(b)
- c = g.centroid
- cx[i] = c.x
- cy[i] = c.y
- except Exception:
- ok[i] = False
- cx[i] = np.nan
- cy[i] = np.nan
- return cx, cy, ok
- def refine_center_fixed_window(centroids, half, n_iter=2):
- if len(centroids) == 0:
- return None
- cx, cy = np.median(centroids[:, 0]), np.median(centroids[:, 1])
- for _ in range(n_iter):
- m = (
- (centroids[:, 0] >= cx-half) & (centroids[:, 0] <= cx+half) &
- (centroids[:, 1] >= cy-half) & (centroids[:, 1] <= cy+half)
- )
- if m.sum() == 0:
- break
- xs = centroids[m, 0]
- ys = centroids[m, 1]
- cx = (xs.min() + xs.max()) / 2.0
- cy = (ys.min() + ys.max()) / 2.0
- return cx, cy
- def iter_polygons(g):
- if g.geom_type == "Polygon":
- yield g
- elif g.geom_type == "MultiPolygon":
- for p in g.geoms:
- yield p
- def add_geom_as_patches(ax, geom, facecolor, edgecolor="black", lw=0.35):
- for poly in iter_polygons(geom):
- x, y = poly.exterior.xy
- coords = np.column_stack([x, y])
- ax.add_patch(MplPolygon(
- coords, closed=True,
- facecolor=facecolor, edgecolor=edgecolor,
- linewidth=lw, joinstyle="miter", capstyle="butt",
- antialiased=True
- ))
- def count_nuclear_alphas_in_tile_for_group(tile_path, group_prefix, na_mask_values):
- df = pd.read_parquet(tile_path)
- if df.empty:
- return 0, 0
- df = df[df["Type"] == "cell"]
- if df.empty:
- return 0, 0
- cx, cy, ok = compute_centroids_from_wkb_bytes(df["Geometry"].values)
- cent = np.column_stack([cx, cy])[ok]
- if cent.shape[0] == 0:
- return 0, 0
- dists, idx = tree.query(cent, k=1)
- idx = idx[dists <= tol]
- if idx.size == 0:
- return 0, 0
- ss_vals = adata.obs[slide_section_col].values[idx]
- in_group = is_R_prefix(ss_vals, group_prefix)
- idx_group = idx[in_group]
- n_group = int(idx_group.size)
- if n_group == 0:
- return 0, 0
- n_nuclear_alpha = int(na_mask_values[idx_group].sum())
- return n_nuclear_alpha, n_group
- def get_top_tiles_for_group(group_prefix, top_k=3):
- na_mask_values = nuclear_alpha_mask(adata).values
- rows = []
- for tf in tile_files:
- n_na, n_group = count_nuclear_alphas_in_tile_for_group(tf, group_prefix, na_mask_values)
- rows.append({"tile": tf, "n_nuclear_alpha": n_na, "n_group_cells": n_group})
- summ = (
- pd.DataFrame(rows)
- .sort_values(["n_nuclear_alpha", "n_group_cells"], ascending=False)
- .reset_index(drop=True)
- )
- top = summ.head(top_k).copy()
- if top["n_nuclear_alpha"].max() == 0:
- raise ValueError(f"No nuclear Alpha MNs found for group {group_prefix}.")
- return top
- def load_tile_geoms_cats_centroids(tile_path, group_prefix):
- df = pd.read_parquet(tile_path)
- df = df[df["Type"] == "cell"].copy()
- if df.empty:
- return [], [], np.zeros((0, 2))
- df["geom"] = df["Geometry"].apply(lambda b: wkb.loads(b) if b is not None else None)
- df = df.dropna(subset=["geom"]).copy()
- if df.empty:
- return [], [], np.zeros((0, 2))
- cent = np.column_stack([
- df["geom"].apply(lambda g: g.centroid.x).to_numpy(),
- df["geom"].apply(lambda g: g.centroid.y).to_numpy()
- ])
- dists, idx = tree.query(cent, k=1)
- keep = dists <= tol
- df = df.iloc[keep].copy()
- cent = cent[keep]
- idx = idx[keep]
- if df.empty:
- return [], [], np.zeros((0, 2))
- ss_vals = adata.obs[slide_section_col].values[idx]
- keep = is_R_prefix(ss_vals, group_prefix)
- df = df.iloc[keep].copy()
- cent = cent[keep]
- idx = idx[keep]
- if df.empty:
- return [], [], np.zeros((0, 2))
- na = nuclear_alpha_mask(adata).values[idx]
- df = df.iloc[na].copy()
- cent = cent[na]
- idx = idx[na]
- if df.empty:
- return [], [], np.zeros((0, 2))
- cats = adata.obs[color_by_col].astype(str).values[idx]
- return df["geom"].tolist(), cats.tolist(), cent
- # -----------------------------
- # Select top tiles
- # -----------------------------
- top_end = get_top_tiles_for_group("R1", top_k=3)
- top_ctl = get_top_tiles_for_group("R2", top_k=3)
- control_tiles = top_ctl["tile"].tolist()
- end_tiles = top_end["tile"].tolist()
- # -----------------------------
- # Build panels
- # -----------------------------
- panels = []
- for t in control_tiles:
- geoms, cats, cent = load_tile_geoms_cats_centroids(t, "R2")
- panels.append(("Control", geoms, cats, cent))
- for t in end_tiles:
- geoms, cats, cent = load_tile_geoms_cats_centroids(t, "R1")
- panels.append(("End", geoms, cats, cent))
- # -----------------------------
- # Plot montage (SVG)
- # -----------------------------
- fig = plt.figure(figsize=(FIG_W_IN, FIG_H_IN))
- gs = fig.add_gridspec(
- 2, 3,
- left=LEFT, right=RIGHT,
- bottom=BOTTOM, top=TOP,
- wspace=WSPACE, hspace=HSPACE
- )
- axes = [[fig.add_subplot(gs[r, c]) for c in range(3)] for r in range(2)]
- for i, (row, geoms, cats, cent) in enumerate(panels):
- r = 0 if row == "Control" else 1
- c = i % 3
- ax = axes[r][c]
- for g, cat in zip(geoms, cats):
- fc = COLOR_MAP.get(cat, "#999999")
- add_geom_as_patches(ax, g, facecolor=fc, edgecolor="black", lw=0.35)
- center = refine_center_fixed_window(cent, HALF, n_iter=2)
- if center is not None:
- cx, cy = center
- ax.set_xlim(cx - HALF, cx + HALF)
- ax.set_ylim(cy + HALF, cy - HALF)
- ax.set_aspect("equal")
- ax.set_xticks([]); ax.set_yticks([])
- for s in ax.spines.values():
- s.set_visible(True)
- s.set_linewidth(SPINE_LW)
- s.set_color("black")
- axes[0][0].text(-0.12, 0.5, "Control", transform=axes[0][0].transAxes,
- rotation=90, va="center", ha="right")
- axes[1][0].text(-0.12, 0.5, "End", transform=axes[1][0].transAxes,
- rotation=90, va="center", ha="right")
- handles = [
- plt.Line2D([0], [0], marker='s', linestyle='',
- markersize=6, markerfacecolor=COLOR_MAP[k],
- markeredgecolor='none', label=LEGEND_LABEL_MAP.get(k, k))
- for k in LEGEND_ORDER
- ]
- fig.legend(handles=handles, loc="lower center", ncol=len(handles),
- frameon=False, bbox_to_anchor=(0.5, 0.04))
- out_svg = os.path.join(
- FIG_DIR,
- "Control_End_MN_segmentations_DM.svg"
- )
- fig.savefig(out_svg)
- plt.show()
- # %% [markdown]
- # ## Alpha motor neuron morphological changes
- # %%
- # Read in AnnData
- adata = sc.read_h5ad(os.path.join(working_dir, f"adata_objects/alpha_anndata_label_transfer.h5ad"))
- # %%
- import numpy as np
- from skimage.filters import threshold_otsu
- import matplotlib.pyplot as plt
- import seaborn as sns
- # log10-transform
- vals = adata.obs["DAPI_high_pass"].values
- log_vals = np.log10(vals + 1)
- # Otsu threshold on log-values
- thresh = threshold_otsu(log_vals)
- print("Otsu log10(DAPI) threshold =", thresh)
- adata.obs["nuclear"] = log_vals >= thresh
- # Plot with threshold line
- plt.figure(figsize=(6,4))
- sns.histplot(log_vals, bins=300, color="steelblue", edgecolor=None)
- plt.axvline(thresh, color="red", linestyle="--", label=f"Otsu={thresh:.2f}")
- plt.xlabel("log10(DAPI_high_pass + 1)")
- plt.ylabel("Count")
- plt.legend()
- plt.tight_layout()
- plt.show()
- # %% [markdown]
- # ## Fig. S9B
- # %% [markdown]
- # #### Control FF vs. SF volume
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.stats import wilcoxon
- # -----------------------------
- # 0. Aesthetics
- # -----------------------------
- sns.set_theme(
- style="ticks",
- rc={
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "axes.linewidth": 0.5,
- "xtick.major.width": 0.5,
- "ytick.major.width": 0.5,
- "pdf.fonttype": 42,
- "ps.fonttype": 42,
- "svg.fonttype": "path"
- }
- )
- group_order = ["Fast-Firing", "Slow-Firing"]
- group_colors = {
- "Fast-Firing": "#2ca02c", # green
- "Slow-Firing": "#ff7f0e" # orange
- }
- # -----------------------------
- # 1. Prepare dataframe (Control, nuclear-only)
- # -----------------------------
- need_cols = ["slide_stage", "slide_section", "stage", "predicted.id", "volume", "nuclear"]
- df = adata.obs[need_cols].copy()
- df = df[df["nuclear"] == True].dropna(
- subset=["slide_stage", "slide_section", "stage", "predicted.id", "volume"]
- )
- df = df[df["stage"] == "Control"].copy()
- df = df[df["predicted.id"].isin(group_order)].copy()
- # -----------------------------
- # 2. Option A: section means -> slide means
- # -----------------------------
- sec_mean = (
- df.groupby(["slide_section", "slide_stage", "predicted.id"], observed=True)
- .agg(sec_mean_volume=("volume", "mean"))
- .reset_index()
- )
- slide_mean = (
- sec_mean.groupby(["slide_stage", "predicted.id"], observed=True)
- .agg(mean_volume=("sec_mean_volume", "mean"))
- .reset_index()
- )
- wide = slide_mean.pivot(
- index="slide_stage",
- columns="predicted.id",
- values="mean_volume"
- )
- paired = wide.dropna(subset=group_order).copy()
- # -----------------------------
- # 3. Paired Wilcoxon signed-rank test (one-sided: SF < FF)
- # -----------------------------
- ff = paired["Fast-Firing"].values
- sf = paired["Slow-Firing"].values
- if len(paired) >= 2:
- res = wilcoxon(
- sf,
- ff,
- alternative="less", # Slow-Firing < Fast-Firing
- zero_method="wilcox"
- )
- p = res.pvalue
- else:
- p = np.nan
- print(f"Paired Wilcoxon test (Slow-Firing < Fast-Firing): p = {p:.3e}")
- # convert to star
- if not np.isfinite(p):
- star = "na"
- elif p < 0.001:
- star = "***"
- elif p < 0.01:
- star = "**"
- elif p < 0.05:
- star = "*"
- else:
- star = "ns"
- # -----------------------------
- # 4. Long format for plotting
- # -----------------------------
- long_means = (
- paired.reset_index()
- .melt(
- id_vars="slide_stage",
- value_vars=group_order,
- var_name="group",
- value_name="mean_volume"
- )
- )
- x_positions = {g: i for i, g in enumerate(group_order)}
- # -----------------------------
- # 5. Boxplot + paired lines
- # -----------------------------
- fig, ax = plt.subplots(figsize=(1.625, 1.7))
- sns.boxplot(
- data=long_means,
- x="group",
- y="mean_volume",
- order=group_order,
- palette=[group_colors[g] for g in group_order],
- showfliers=False,
- width=0.6,
- linewidth=0.8,
- ax=ax
- )
- sns.stripplot(
- data=long_means,
- x="group",
- y="mean_volume",
- order=group_order,
- color="black",
- size=4,
- jitter=0.12,
- alpha=0.8,
- ax=ax
- )
- # connect paired points
- for slide, sub in long_means.groupby("slide_stage"):
- if set(sub["group"]) == set(group_order):
- x1 = x_positions["Fast-Firing"]
- y1 = sub.loc[sub["group"] == "Fast-Firing", "mean_volume"].iloc[0]
- x2 = x_positions["Slow-Firing"]
- y2 = sub.loc[sub["group"] == "Slow-Firing", "mean_volume"].iloc[0]
- ax.plot([x1, x2], [y1, y2], lw=1, c="gray", alpha=0.7)
- # aesthetics
- ax.grid(False)
- for spine in ax.spines.values():
- spine.set_color("black")
- # ↓↓↓ Shorter ticks (Option 2 incorporated here) ↓↓↓
- ax.tick_params(axis="both", which="major", length=2, width=0.5)
- ax.set_xlabel("")
- ax.set_ylabel("Volume")
- ax.set_title("Control")
- ax.set_xticklabels(["FF", "SF"])
- # -----------------------------
- # Stats bracket (ORIGINAL geometry)
- # -----------------------------
- data_min = long_means["mean_volume"].min()
- data_max = long_means["mean_volume"].max()
- data_range = max(1e-12, data_max - data_min)
- y_base = data_max * 1.02
- y_step = 0.12 * data_range
- h = 0.25 * y_step
- x1, x2 = x_positions["Fast-Firing"], x_positions["Slow-Firing"]
- ax.plot(
- [x1, x1, x2, x2],
- [y_base, y_base + h, y_base + h, y_base],
- lw=1.2, c="k"
- )
- ax.text(
- (x1 + x2) / 2,
- y_base + h * 1.05,
- star,
- ha="center",
- va="bottom",
- fontsize=7
- )
- ymin, ymax = ax.get_ylim()
- ax.set_ylim(ymin, max(ymax, y_base + y_step * 1.3))
- plt.tight_layout()
- # -----------------------------
- # 6. Save plot
- # -----------------------------
- out_svg = os.path.join(figures_dir, "control_FF_vs_SF_volume_paired_wilcoxon.svg")
- out_png = os.path.join(figures_dir, "control_FF_vs_SF_volume_paired_wilcoxon.png")
- fig.savefig(out_svg, format="svg", dpi=300, bbox_inches="tight", transparent=True)
- fig.savefig(out_png, format="png", dpi=300, bbox_inches="tight", transparent=True)
- plt.show()
- # %% [markdown]
- # ## Fig. S9C
- # %% [markdown]
- # #### Volume changes with disease stage
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
- from scipy.stats import mannwhitneyu
- from statsmodels.stats.multitest import multipletests
- from statannotations.Annotator import Annotator
- # ------------------ Config ------------------
- order = ["Control", "Early", "Mid", "End"]
- comparisons = [("Early", "Control"), ("Mid", "Control"), ("End", "Control")]
- stage_colors = {
- "Control": "#1f77b4",
- "Early": "#ff7f0e",
- "Mid": "#d62728",
- "End": "#2ca02c"
- }
- # ------------------ Global aesthetics ------------------
- sns.set_theme(
- style="ticks",
- rc={
- "font.size": 7,
- "axes.titlesize": 7,
- "axes.labelsize": 7,
- "xtick.labelsize": 7,
- "ytick.labelsize": 7,
- "legend.fontsize": 7,
- "legend.title_fontsize": 7,
- "axes.linewidth": 0.5,
- "xtick.major.width": 0.5,
- "ytick.major.width": 0.5,
- "pdf.fonttype": 42,
- "ps.fonttype": 42,
- "svg.fonttype": "path"
- }
- )
- # ------------------ Load & clean (NUCLEAR ONLY) ------------------
- df = (
- adata.obs[['slide_stage', 'slide_section', 'stage', 'volume', 'nuclear']]
- .copy()
- .dropna(subset=['slide_stage', 'slide_section', 'stage', 'volume', 'nuclear'])
- )
- df = df[df['nuclear'] == True].copy()
- df = df.drop(columns=['nuclear'])
- df['stage'] = df['stage'].astype(str).str.strip()
- df['volume'] = pd.to_numeric(df['volume'], errors='coerce')
- df = df.dropna(subset=['volume'])
- df = df[df['stage'].isin(order)].copy()
- # ------------------ Option A: slide_section means -> slide_stage means ------------------
- sec_mean = (
- df.groupby(['slide_section', 'slide_stage', 'stage'], observed=True)
- .agg(sec_mean_volume=('volume', 'mean'))
- .reset_index()
- )
- stage_per_slide = (
- df.groupby('slide_stage')['stage']
- .agg(['nunique', 'first'])
- .rename(columns={'nunique': 'n_stage_labels', 'first': 'stage'})
- .reset_index()
- )
- bad = stage_per_slide[stage_per_slide['n_stage_labels'] != 1]
- if len(bad):
- print("[WARN] Some slide_stage have multiple stage labels; using the first label:")
- print(bad.head())
- stage_per_slide = stage_per_slide[['slide_stage', 'stage']]
- slide_mean = (
- sec_mean.groupby('slide_stage', observed=True)
- .agg(mean_volume=('sec_mean_volume', 'mean'))
- .reset_index()
- )
- per_slide = slide_mean.merge(stage_per_slide, on='slide_stage', how='left')
- per_slide = per_slide[per_slide['stage'].isin(order)].copy()
- per_slide['stage'] = pd.Categorical(per_slide['stage'], categories=order, ordered=True)
- print("Slides per stage:")
- print(per_slide.groupby('stage')['slide_stage'].nunique(), "\n")
- # ------------------ Build arrays for tests ------------------
- stage_arrays = {
- g: per_slide.loc[per_slide['stage'].eq(g), 'mean_volume'].to_numpy()
- for g in order
- }
- for g in order:
- print(f"{g}: n_slides = {len(stage_arrays[g])}")
- # ------------------ Statistics (MWU one-sided + Bonferroni) ------------------
- rows = []
- for s_label, ctl_label in comparisons:
- x = stage_arrays[s_label]
- y = stage_arrays[ctl_label]
- if len(x) >= 1 and len(y) >= 1:
- try:
- p = mannwhitneyu(x, y, alternative='less', method='asymptotic').pvalue
- except TypeError:
- p = mannwhitneyu(x, y, alternative='less').pvalue
- else:
- p = np.nan
- rows.append((s_label, ctl_label, p))
- res = pd.DataFrame(rows, columns=['stage', 'control', 'p_raw'])
- mask = res['p_raw'].notna()
- res['p_adj'] = np.nan
- if mask.any():
- _, p_adj, _, _ = multipletests(res.loc[mask, 'p_raw'], method='bonferroni')
- res.loc[mask, 'p_adj'] = p_adj
- print(res, "\n")
- # p-values in the same order as `comparisons`
- pvals_adj = [res.loc[i, 'p_adj'] for i in range(len(comparisons))]
- # ------------------ Plot ------------------
- fig, ax = plt.subplots(figsize=(1.625, 1.8))
- sns.boxplot(
- data=per_slide,
- x='stage',
- y='mean_volume',
- order=order,
- palette=[stage_colors[s] for s in order],
- showfliers=False,
- width=0.6,
- linewidth=0.8,
- ax=ax
- )
- sns.stripplot(
- data=per_slide,
- x='stage',
- y='mean_volume',
- order=order,
- color='black',
- size=4,
- jitter=0.15,
- alpha=0.8,
- ax=ax
- )
- ax.grid(False)
- for spine in ax.spines.values():
- spine.set_color("black")
- # shorter ticks (same as your other panels)
- ax.tick_params(axis="both", which="major", length=2, width=0.5)
- ax.set_xlabel("")
- ax.set_ylabel("Volume")
- ax.set_title("All")
- for label in ax.get_xticklabels():
- label.set_rotation(90)
- # ------------------ Stat annotations (automatic bracket geometry) ------------------
- annotator = Annotator(
- ax,
- comparisons,
- data=per_slide,
- x="stage",
- y="mean_volume",
- order=order
- )
- annotator.configure(
- test=None,
- text_format="star",
- loc="inside", # same behavior as your good-looking example
- fontsize=6 # your request
- )
- annotator.set_pvalues(pvals_adj)
- annotator.annotate()
- ax.set_yticks([7500, 10000, 12500, 15000, 17500])
- plt.tight_layout()
- # ------------------ Save ------------------
- out_svg = os.path.join(figures_dir, "stage_volume_per_slide_nuclear_wilcoxon_bonferroni.svg")
- out_png = os.path.join(figures_dir, "stage_volume_per_slide_nuclear_wilcoxon_bonferroni.png")
- fig.savefig(out_svg, dpi=300, bbox_inches='tight', transparent=True)
- fig.savefig(out_png, dpi=300, bbox_inches='tight', transparent=True)
- plt.show()
- print("Saved:", out_svg, "and", out_png)
- # %% [markdown]
- # ## Fig. S9D
- # %%
- import os
- import numpy as np
- import pandas as pd
- import matplotlib.pyplot as plt
- impo
14. MERFISH spatial transcriptomics analysis.ipynb at commit 10e89db, no license · at the source
Overview
and 6 other authors
Kailee Ong13, Don W. Cleveland12, John Ravits13, Jessica E. Rexach5,14, William J. Greenleaf1,15,9, Aaron D. Gitler1,9,16,1717 affiliations
- Department of Genetics, Stanford University School of Medicine, Stanford, CA 94305, USA
- Stanford Neurosciences Graduate Program, Stanford University School of Medicine, Stanford, CA 94305, USA
- These authors contributed equally
- Research, Biogen Inc., Cambridge, MA 02142, USA
- Program in Neurogenetics, Department of Neurology, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
- Herbold Computational Biology Program, Public Health Sciences Division, Fred Hutchinson Cancer Center, Seattle, WA 98195, USA
- Department of Genome Sciences, University of Washington, Seattle, WA 98195, USA
- Brotman Baty Institute, University of Washington, Seattle, WA 98195, USA
- Biohub – San Francisco, San Francisco, CA 94158, USA
- Department of Neurobiology, Stanford University, Stanford, CA 94305, USA
- Vollum Institute, Oregon Health and Sciences University, Portland, OR 97239, USA
- Departments of Cellular and Molecular Medicine, University of California, San Diego, La Jolla, CA 92093, USA
- Department of Neurosciences, University of California, San Diego, La Jolla, CA 92093, USA
- Department of Human Genetics, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
- Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
- The Phil and Penny Knight Initiative for Brain Resilience, Stanford University, Stanford, CA 94305, USA
- Lead contact
Abstract
To define molecular determinants of motor neuron degeneration in amyotrophic lateral sclerosis (ALS), we generated longitudinal single-nucleus transcriptomes and chromatin accessibility profiles of spinal motor neurons together with spatial transcriptomics from the SOD1-G93A mouse model. Vulnerable alpha motor neurons showed thousands of molecular changes, marking a transition into a distinct cell state we named “disease-associated motor neurons” (DMs). We identified transcription factor networks that govern how healthy cells transition into DMs and those associated with motor neuron subtype-selective vulnerability. Upregulation of DM-associated transcription factors in human motor neurons induced key features of DMs, demonstrating an active regulatory component. Human ALS spinal cord single-nucleus RNA sequencing data demonstrated conservation of the DM signature in alpha motor neurons, and human orthologs of regions differentially accessible in SOD1-G93A mouse motor neurons were enriched for ALS genetic risk variants. Together, these findings establish a conserved, genetically linked motor neuron signature in ALS.
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 40 matches between paragraphs and lines of code.
omgautier/Gautier_Blum_2026
10e89dba3fdbd6b3f9d642a161995d913d6bed0e, 23 June 2026Availability: 1 check, the latest on 27 September 2026: the link answers
- 27 September 2026: the link answers
14 files
- 01. snRNA-seq 1 - initial processing and input preparation.ipynb, Jupyter, 667 lines, 5 matches
- 02. snRNA-seq 2 - DESeq2 differential expression analysis.ipynb, Jupyter, 755 lines, 1 match
- 03. snRNA-seq 3 - Fig. 2B, 2D-H, S1B-C, S2A-D.ipynb, Jupyter, 732 lines, 1 match
- 04. snRNA-seq 4 - Fig. 2A-J, S2E, S3A-H.ipynb, Jupyter, 997 lines, 3 matches
- 05. Multiome 1 - initial processing and input preparation for Fig. 3B-C, 3G-J, S4A-D, S5A-B.ipynb, Jupyter, 647 lines, 3 matches
- 06. Multiome 2 - Fig. 3B-C, 3G-J, S4A-D, S5A-B.ipynb, Jupyter, 260 lines
- 07. Multiome 3 - input preparation for Fig. 3D-F, S4E-H.ipynb, Jupyter, 388 lines, 5 matches
- 08. Multiome 4 - Fig. 3D-F, S4E-H.ipynb, Jupyter, 376 lines, 1 match
- 09. Multiome 5 - input preparation for Fig. 4A-D, S6A-G.ipynb, Jupyter, 328 lines, 1 match
- 10. Multiome 6 - Fig. 4A-D, S6A-G.ipynb, Jupyter, 675 lines, 2 matches
- 11. Multiome & snRNA-seq 1 - input preparation (cross-modal label transfer, multiome RNA differential expression).ipynb, Jupyter, 359 lines, 3 matches
- 12. Multiome & snRNA-seq 2 - Fig. S5C-F, S10C-J.ipynb, Jupyter, 847 lines
- 13. iMN in vitro TF OE - Fig. 5.ipynb, Jupyter, 660 lines, 2 matches
- 14. MERFISH spatial transcriptomics analysis.ipynb, Jupyter, 6,326 lines, 13 matches
- repository limit reached (2,000 files or 30 MB): the rest is at the source (3 files)
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;
- 14 scripts, each with its path and the digest of its content;
- 40 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
- zenodo:16938739, at Zenodo; found in “Data and code availability”
- zenodo:5503876, at Zenodo; found in the text, “Partitioned Heritability and H-MAGMA Analyses”
Data and code availability
Raw and processed sequencing data have been deposited in the NCBI Gene Expression Omnibus (GEO) under accession numbers GEO: GSE306676 for the snRNA-seq data and GEO: GSE306675 for the multiome (paired snATAC/
Original code is available on GitHub: https://
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 2, 28 September 2026
- Publisher: n/a → Cell Press
Version 1, 27 September 2026: the first record
Recorded: type, language, journal, volume, issue, pages, dates, 26 authors, 10 keywords, 14 MeSH terms, 5 funders, 137 references, 9 RRIDs.
Cite
This paper
Gautier, O., Blum, J. A., Nguyen, T. P., Cao, S., Klemm, S., Yamakawa, M., Huh, D., Hurt, J. A., Sinnott-Armstrong, N., Zeng, Y., Davis, C.-h. O., Bombosch, J., Liu, C., Encarnacion, L. N., Guttenplan, K. A., Chen, D., Kathiria, A., Zhao, L., Moore, S., . . . Gitler, A. D. (2026). An emergent disease-associated motor neuron state precedes cell death in ALS. Cell, 189(16), 5044-5064.e12. https://
BibTeX
@article{gautier2026emer
author = {Gautier, Olivia and Blum, Jacob A. and Nguyen, Thao P. and Cao, Shaolong and Klemm, Sandy and Yamakawa, Mai and Huh, Dann and Hurt, Jessica A. and Sinnott-Armstrong, Nasa and Zeng, Yi and Davis, Chung-ha O. and Bombosch, Juliane and Liu, Chang and Encarnacion, Lisa N. and Guttenplan, Kevin A. and Chen, Derek and Kathiria, Arwa and Zhao, Luke and Moore, Stephen and Meng, Alex and Ong, Kailee and Cleveland, Don W. and Ravits, John and Rexach, Jessica E. and Greenleaf, William J. and Gitler, Aaron D.},
title = {{An emergent disease-associated motor neuron state precedes cell death in ALS}},
journal = {Cell},
year = {2026},
month = jun,
volume = {189},
number = {16},
pages = {5044--5064.e12},
publisher = {Cell Press},
issn = {0092-8674},
doi = {10.1016/
url = {https://
pmid = {42335888},
pmcid = {PMC13446465}
}
RIS
TY - JOUR
AU - Gautier, Olivia
AU - Blum, Jacob A.
AU - Nguyen, Thao P.
AU - Cao, Shaolong
AU - Klemm, Sandy
AU - Yamakawa, Mai
AU - Huh, Dann
AU - Hurt, Jessica A.
AU - Sinnott-Armstrong, Nasa
AU - Zeng, Yi
AU - Davis, Chung-ha O.
AU - Bombosch, Juliane
AU - Liu, Chang
AU - Encarnacion, Lisa N.
AU - Guttenplan, Kevin A.
AU - Chen, Derek
AU - Kathiria, Arwa
AU - Zhao, Luke
AU - Moore, Stephen
AU - Meng, Alex
AU - Ong, Kailee
AU - Cleveland, Don W.
AU - Ravits, John
AU - Rexach, Jessica E.
AU - Greenleaf, William J.
AU - Gitler, Aaron D.
TI - An emergent disease-associated motor neuron state precedes cell death in ALS
T2 - Cell
J2 - Cell
PY - 2026
DA - 2026/
VL - 189
IS - 16
SP - 5044
EP - 5064.e12
SN - 0092-8674
PB - Cell Press
DO - 10.1016/
UR - https://
LA - en
ER -
CSL-JSON
{
"id": "10.1016/
"type": "article-journal",
"title": "An emergent disease-associated motor neuron state precedes cell death in ALS",
"container-title": "Cell",
"author": [
{
"family": "Gautier",
"given": "Olivia"
},
{
"family": "Blum",
"given": "Jacob A."
},
{
"family": "Nguyen",
"given": "Thao P."
},
{
"family": "Cao",
"given": "Shaolong"
},
{
"family": "Klemm",
"given": "Sandy"
},
{
"family": "Yamakawa",
"given": "Mai"
},
{
"family": "Huh",
"given": "Dann"
},
{
"family": "Hurt",
"given": "Jessica A."
},
{
"family": "Sinnott-Armstrong",
"given": "Nasa"
},
{
"family": "Zeng",
"given": "Yi"
},
{
"family": "Davis",
"given": "Chung-ha O."
},
{
"family": "Bombosch",
"given": "Juliane"
},
{
"family": "Liu",
"given": "Chang"
},
{
"family": "Encarnacion",
"given": "Lisa N."
},
{
"family": "Guttenplan",
"given": "Kevin A."
},
{
"family": "Chen",
"given": "Derek"
},
{
"family": "Kathiria",
"given": "Arwa"
},
{
"family": "Zhao",
"given": "Luke"
},
{
"family": "Moore",
"given": "Stephen"
},
{
"family": "Meng",
"given": "Alex"
},
{
"family": "Ong",
"given": "Kailee"
},
{
"family": "Cleveland",
"given": "Don W."
},
{
"family": "Ravits",
"given": "John"
},
{
"family": "Rexach",
"given": "Jessica E."
},
{
"family": "Greenleaf",
"given": "William J."
},
{
"family": "Gitler",
"given": "Aaron D."
}
],
"container-title-short":
"volume": "189",
"issue": "16",
"page": "5044-5064.e12",
"DOI": "10.1016/
"PMID": "42335888",
"PMCID": "PMC13446465",
"ISSN": "0092-8674",
"publisher": "Cell Press",
"URL": "https://
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
23
]
]
}
}
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/s41593-026-02300-5 [code]
- Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.Journal: Nature neuroscienceIn common: Squidpy, anndata, DESeq2, 11 other tools, genetics / omics, other condition, cellular / molecular, 6 references
- [2] doi:10.1186/s13024-026-00944-2
- TDP-43: [GU]-ardian of the transcriptome.Journal: Molecular neurodegenerationIn common: genetics / omics, other condition, cellular / molecular, 15 references
- [3] doi:10.1016/j.xcrm.2026.102766 [code]
- A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.Journal: Cell reports. MedicineIn common: anndata, DESeq2, Scanpy, 12 other tools, genetics / omics, other condition, cellular / molecular, 5 references
- [4] doi:10.1038/s42003-026-10034-0 [code]
- Region- and cell type-specific changes in gene expression in the cerebellum after classical fear conditioning.Journal: Communications biologyIn common: Squidpy, anndata, DESeq2, 12 other tools, genetics / omics, mouse, cellular / molecular, 3 references
- [5] doi:10.1038/s41467-026-71803-3 [code]
- Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.Journal: Nature communicationsIn common: Squidpy, anndata, DESeq2, 12 other tools, genetics / omics, mouse, cellular / molecular, 2 references
- [6] doi:10.1038/s41586-026-10629-x [code]
- Whole-genome duplication shaped cell-type evolution in the vertebrate brain.Journal: NatureIn common: anndata, DESeq2, Scanpy, 12 other tools, genetics / omics, mouse, cellular / molecular, 3 references
- [7] doi:10.1016/j.isci.2026.116906 [code]
- Evaluating exon skipping in the central nervous system in Duchenne muscular dystrophy using spatial transcriptomics.Journal: iScienceIn common: Squidpy, statannotations, anndata, 10 other tools, genetics / omics, other condition, mouse, 3 references
- [8] doi:10.1002/imt2.70163 [code]
- Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.Journal: iMetaIn common: anndata, DESeq2, Scanpy, 12 other tools, genetics / omics, mouse, cellular / molecular, 2 references
- [9] doi:10.1038/s41514-026-00391-9 [code]
- Region-specific transcriptional signatures of brain aging in the absence of neuropathology at the single-cell level.Journal: npj agingIn common: anndata, broom, Scanpy, 11 other tools, genetics / omics, cellular / molecular, 4 references
- [10] doi:10.1186/s13059-026-04177-w [code]
- Genomic sequence evolution underlying human neocortical interareal diversification.Journal: Genome biologyIn common: Squidpy, anndata, Scanpy, 10 other tools, genetics / omics, mouse, cellular / molecular, 3 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.
Claim this paper
Correct its record
Say what each link of this record is, remove the ones that are not the paper's, add the ones that are missing. The correction becomes a new version of the record, in its Versions section.
Validate its tracing map
You validate the map as this page shows it: 1 repository of the authors' code, each at its verified commit and with its license, 14 scripts, and 40 matches between paragraphs and code (see the Code and Map sections). It then receives a DOI on Zenodo, with you (your ORCID iD) and OSCR as its creators; the code itself is not deposited.
The map's fingerprint: sha256:4f5c7742097771e8…
Add the badge to its README
The badge links the code to this page. Copy one of these into the README of the paper's code: only you decide where it goes, and nothing is changed for you.
Markdown
[, paste the snippet at the top, then “Commit changes…” and, to review it first, “Create a new branch and start a pull request”. You open the pull request; OSCR asks for no permission.
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.
