DOI: 10.5061/dryad.kprr4xhg7
title: Single nucleus sequencing of mouse hippocampi four days after pilocarpine induced status epilepticus

description: All current drug treatments for epilepsy, a neurological disorder
 affecting over 50 million people merely treat symptoms, and a third of
 patients with epilepsy do not respond to medication. There are no disease
 modifying treatments that may be administered briefly to patients to
 enduringly eliminate spontaneous seizures and reverse cognitive deficits.
 Applying network approaches to whole tissue and single-nuclei
 transcriptomic data collected from mouse models of temporal lobe epilepsy
 and publicly available transcriptomic data from human temporal lobectomy
 samples, we confirmed a previously described pattern of rapid and
 transient induction of the Janus kinase/signal transducers and activators
 of transcription (JAK/STAT) pathway within days of epileptogenic insult.
 This was followed by a resurgent activation weeks to months later with the
 onset of spontaneous seizures. Targeting the first wave of JAK/STAT
 activation after epileptic insult did not prevent seizures. However,
 inhibition of the second wave with CP690550 (Tofacitinib) over a 2-week
 period enduringly suppressed seizures, rescued deficits in spatial memory,
 and alleviated epilepsy-associated histopathological alterations. Seizure
 suppression lasted for at least 2 months after the final dose. These
 results indicate that reignition of inflammatory JAK/STAT3 signaling in
 chronic epilepsy opens a window for disease modification with the
 FDA-approved, orally available drug CP690550.
description: <strong><em>Pilocarpine model of status
 epilepticus.</em></strong><strong> </strong>Status
 epilepticus was induced in male and female C57BL/6CRN mice with
 pilocarpine (8-12 weeks old). Mice were injected with terbutaline (4 mg/kg
 i.p., Sigma T2528) and methylscopolamine (4 mg/kg ip, Sigma S8502) to
 alleviate respiratory and cardiovascular effects of pilocarpine, then 30
 min later received saline or pilocarpine (280 mg/kg i.p. as free base,
 Sigma P6503). After 1 hour of SE, seizures were interrupted by diazepam
 (10mg/kg i.p.). Mice were weighed and scored for behavioral deficits
 (Irwin score) daily, with injections of lactated Ringer’s solution as
 needed. <strong><em>Tissue Harvest and
 Sequencing. </em></strong>To select which epileptic mice would
 be sent for sequencing, we determined a composite score of seizure
 severity during pilocarpine SE, weight loss 24 hours after SE, ability to
 build a nest and neurobehavioral recovery scores on day 4 after SE (See
 Fig. S14 in linked manuscript Hoffman et al). The mean number of standard
 deviations from the mean of each of these 4 measures was calculated across
 all 6 saline-treated and 14 pilocarpine-treated mice, and mice were
 selected based on wide separation between saline and pilo groups. Both
 hippocampi and a small portion of neocortex overlaying one hippocampus
 were snap frozen and shipped to Novagene.
 <strong><em>Sequencing, Preprocessing, Expression
 Quantification and Quality
 Control.</em></strong><strong> </strong>Each
 sample was individually subjected to library preparation and Illumina
 sequencing. Reads were processed and mapped with Cell Ranger software (10x
 Genomics) and the Seurat R package at Novagene.  Counts matrices for each
 sample were combined into an adata object in the Roopra lab and subjected
 to quality control, cell and gene filtering using the Scanpy package.
 Genes were kept that were expressed in at least 100 cells across all
 samples. Cells with genes numbering between 1000 and 6500 and with fewer
 than 5% mitochondrial genes were kept. This resulted in 64,647 cells and
 21,423 genes. For PCA analysis of samples, reads were
 summed across all cells for each sample to produce a pseudo-bulked sample
 (n=10) x gene (21,423) array.  4000 most variant genes were selected to
 subject samples to PCA analysis.  All 5 pilo samples clustered together
 when projecting either the 1<sup>st</sup> 2 or the
 2<sup>nd</sup> 2  PCA components. Four saline samples
 clustered together away from the pilo samples.  One saline sample did not
 cluster with either group and was discarded (Fig. S14 in linked manuscript
 Hoffman et al).  This resulted in 5 pilo samples, 4 saline samples, and
 58,365 cells. <strong><em>Normalization and
 Clustering.</em></strong> Data was normalized using the scanpy
 preprocessing normalization tool.  Highly expressed genes were omitted,
 and the target sum was set to 10,000 (transcripts per 10k or TPTK). Counts
 were converted to log(1+tptk) and highly variable genes identified with
 min_mean=0.0125, max_mean=3, min_disp=0.5. Genes were scaled to unit
 variance and values&gt;10 standard deviations were clipped. PCA with
 50 components was performed and used to generate a neighborhood graph with
 n_neighbors=15. Leiden clustering with resolution=0.25 yielded 27
 clusters.
 <strong><em>Annotations. </em></strong> <em>Class annotation:</em><strong> </strong>Clusters were allocated to classes using the following markers: Inhibitory Neurons: GAD1, Excitatory Neurons: SLC17A7, Astrocytes: AQP4, Oligodendrocytes: SOX10, Endothelial cells: RGS5, Myeloid: PTPRC.  Each Class were broken up into Cell Classes as follows: <em>Excitatory Neurons</em>: The single cell dataset from Yao et al. (ref <em>5 in linked manuscript Hoffman et al</em>) was converted to an anndata object and trimmed to glutamatergic neurons (CA1-ProS, CA2-IG-FC, CA3, CR, DG). DEGs/markers were identified for each subtype using the rank_genes_groups scanpy function with method=’wilcoxon’. The same was done for our Excitatory adata object using groupby=’leiden’.  The top 100 Yao markers per cell type and our leiden clusters were used to perform pairwise fisher exact tests to look for overlaps between Yao cell types and our leiden clusters. The Yao cell type with highest overlap with Leiden pairs as judged by highest score (log(odds)*-log(p)) was assigned to the cluster. This yielded the following subclasses: CA1_ProS, CA1_ve, CA1_do, CA3, CR, DG. Using Hipposeq (<em>6 in linked manuscript Hoffman et al</em>) the CA3 subclass (markers: [SPOCK1, MGAT4C, ELAVL2]) was broken up into CA3_ve, CA3_do and mossy cells using the following markers: CA3_ve: [NECAB1,COCH,ADGRA1,CPNE7], CA3_do: [CHGB,RIMBP2], mossy: [CALCRL,CALB2]. Yao et al CR cell markers were RELN, CACNA2D2, STMN1. DG markers were PROX1, PCP4, ADARB2, STXBP6. <em>Inhibitory Neurons</em>: As with excitatory neurons, the Yao dataset was trimmed to GABAergic neurons. Rank genes groups was run on the inhibitory adata object and the Yao adata object. Fisher exact tests with the top 100 genes per group was performed for each Yao subtype and Inhibitory leiden cluster. This approach defined 10 inhibitory subtypes in our dataset: Lamp5 Lhx6, Lamp5, Meis2, Ntng1, Pvalb Vipr2, Pvalb, Sst Ntng1, Sst, Vip, Vip Igfbp6. <em>Astrocytes</em>: AQP4 was used to define astrocytes. GFAP was used to define activated and resting states. <em>Myeloid cells:</em> The myeloid cell class was first defined by PTPRC. The myeloid adata object was then reclustered with n_neighbors=15 and leiden resolution=0.25. This gave 12 clusters. GFAP defined activated microglia. SIGLCECH positive and PTPRC negative marked resting microglia. Cells with high TOP2A, PCNA and MKI67 defined proliferating microglia. CCR2 and CD44 positive denoted invading monocytes. MRC1 and VCAM marked Perivascular Macrophages. <em>Oligodendrocytes:</em> The oligodendrocyte cluster was reclustered with n_neighbors=15 and leiden resolution=0.25 to give 9 clusters. Mature oligodendrocytes were marked by MOG. Progenitors (OPC in Fig. 2 in linked manuscript Hoffman et al) were marked by VCAN.
description: # Single nucleus sequencing of mouse hippocampi four days after
 pilocarpine induced status epilepticus
 [https://doi.org/10.5061/dryad.kprr4xhg7](https://doi.org/10.5061/dryad.kprr4xhg7) ## Description of the data and file structure All current drug treatments for epilepsy, a neurological disorder affecting over 50 million people merely treat symptoms, and a third of patients with epilepsy do not respond to medication. There are no disease modifying treatments that may be administered briefly to patients to enduringly eliminate spontaneous seizures and reverse cognitive deficits. Applying network approaches to mouse models of temporal lobe epilepsy and human temporal lobectomy samples at both whole tissue and single-nuclei resolutions, we observe the well-characterized pattern of rapid induction and subsequent quenching exhibited of the JAK/STAT pathway within days of epileptogenic insult. This is followed by a resurgent activation weeks to months later with the onset of spontaneous seizures. Targeting the first wave of activation after epileptic insult does not prevent disease.  However, brief inhibition of the second wave with CP690550 (Tofacitinib) enduringly suppresses seizures, rescues deficits in spatial memory, and alleviates epilepsy-associated histopathological alterations. Seizure suppression lasts for at least 2 months after the final dose.  We find that reignition of inflammatory JAK/STAT3 signaling in chronic epilepsy opens a window for disease modification with the FDA-approved, orally available drug CP690550. ### Files and variables The h5ad file is an anndata object containing single nucleus Seq profiling of 5 naive and 4 epileptic male and female C57BL/6CRN mouse hippocampi harvested 4 days post Status Epilepticus induced with pilocarpine. #### File: annotated\_adata\_in\_Hoffman\_et\_al\_2025.h5ad.gz **Description:**  This is an anndata object with 58,365 cells, 21,423 genes and the following structure: AnnData object with n_obs × n_vars = 58365 × 2142  obs: 'sample', 'condition', 'leiden', 'class', 'new_leiden', 'subclass', 'cell_type', 'cell_class', 'DIM1', 'DIM2' var: 'highly_variable'  uns: 'cell_class_colors', 'cell_class_sizes', 'cell_type_color_dict', 'cell_type_colors', 'dendrogram_cell_type', 'dendrogram_subclass', 'neighbors', 'paga', 'rank_genes_groups', 'sample_condition', 'umap' obsm: 'X_pca', 'X_umap' layers: 'raw_counts' obsp: 'connectivities', 'distances' Values are log2(1+tptk) where tptk is the transcripts per ten thousand for each gene in a cell. Cell attributes are kept in the ‘obs’ slot. ‘sample’: {'A231', 'A234', 'A235', 'A237', 'A238', 'A239', 'A240', 'A241', 'A246'} ‘condition’: {‘sal’, ‘pilo’} ‘leiden’: original leiden cluster cell was assigned ‘class’: Level 1 (broadest) cell category: {'Astrocyte', 'Endothelial', 'Excitatory', 'Inhibitory', 'Myeloid', 'Oligodendrocyte'} ‘new_leiden’: modified leiden designation after re-combining individual class adata objects during annotation. ‘subclass’: Level 2 cell classification: {'Astrocyte', 'CA1-ProS', 'CA2-IG-FC', 'CA3', 'DG', 'Endothelial', 'Lamp5', 'Meis2', 'Mossy', 'Myeloid', 'Ntng1', 'Oligodendrocyte', 'Pvalb', 'Sst', 'Vip'} ‘cell_class’: Level 3 cell classification: {'Astrocyte', 'CA1', 'CA2', 'CA3', 'CR', 'DG', 'Endothelial', 'Invading_Monocyte', 'Lamp5', 'Lamp5 Lhx6', 'Meis2', 'Microglia', 'Mossy', 'Ntng1', 'Oligo', 'Perivascular_Macrophage',  'Pvalb','Pvalb Vipr2', 'Sst', 'Sst Ntng1', 'Subiculum', 'Vip', 'Vip Igfbp6'} ‘cell_type’: Level 4 (most granular) cell classification:  {'361_DG', '362_DG', '363_DG', '364_DG', 'Activated_Astrocyte', 'Activated_Microglia', 'CA1_do', 'CA1_ve', 'CA2', 'CA3_do', 'CA3_ve', 'CR', 'Endothelial', 'Invading_Monocyte', 'Lamp5', 'Lamp5 Lhx6', 'Mature_Oligo', 'Meis2', 'Mossy', 'Ntng1', 'OPC', 'Perivascular_Macrophage', 'Pvalb', 'Pvalb Vipr2', 'Replicating_Microglia', 'Resting_Astrocyte', 'Resting_Microglia', 'Sst', 'Sst Ntng1', 'Subiculum', 'Vip', 'Vip Igfbp6'} 'DIM1','DIM2': X Y coordinates of cell in UMAP Dictionaries containing preprocessing and other information are stored in ‘uns’ slot. Gene attributes are kept in the ‘var’ slot.                      ‘highly_variable’: genes flagged as highly variable using default scanpy settings Dimension reduced arrays are kept in the obsm slot:                      ‘X_pca’: 1st 50 dimensions of Principle Component Analysis (rows = cells, cols = dimensions).                      ‘X_umap’: 1st 2 dimensions of UMAP projection. (Same as DIM1 and DIM2 in obs slot) Raw counts for genes in each cell are kept in a numpy array in:                      layers: ‘raw_counts’ #### File: Deseq2\_output.zip **Description:**  This is a folder containing Deseq2 outputs for all 27 cell classes in the h5ad file. Deseq2 was used for differential gene analysis between saline and epileptic samples for each cell class (Astrocyte, CA1, CA2, CA3, CR, DG, Endothelial, Invading_Monocyte, Lamp5, Lamp5 Lhx6, Meis2, Microglia, Mossy, Ntng1, Oligo, Perivascular_Macrophage, Pvalb, Pvalb Vipr2, Sst, Sst Ntng1, Subiculum, Vip, Vip Igfbp6).  For each cell class, reads for each gene were summed within a sample such that each sample contributed a single integer value (total reads) for each gene.  Thus, a sample x gene matrix was generated for every cell class.  Each cell class sample x gene matrix was passed to Deseq2 and Wald tests were performed to find DEGs between saline and pilo samples. **Deseq2 output was then filtered for expressed genes in each cell class.  Expressed genes were defined as those genes that had at least 1 read in every sample of either condition.  After filtering Deseq2 output, padj was recalculated using the lower gene number.** DEGs were genes with an absolute log fold change of &gt;1.2 and a padj&lt;0.05 for further analysis. ## Code/software ### Ontomancer. ###### Drag Ontomancer.py into terminal or enter path to file to start and follow instructions.  This script accepts a file containing gene lists and performs ontological analysis using gene sets from [https://www.gsea-msigdb.org/gsea/index.jsp](https://www.gsea-msigdb.org/gsea/index.jsp).  It also analyses 1Kb_gene_2000.GMT which contains the top 2000 target genes for each transcription factor or cofactor in the MAGIC Matrix file 1Kb_gene.mtx (Roopra, Plos.Comp.Biol (2020)). The following databases are used: KEGG.gmt, Reactome.gmt, Biocarta.gmt, Molecular_Function.gmt, Protein_Interaction_Database.gmt, Hallmark.gmt, Biological_Function.gmt, Cellular_Component.gmt, Expression_Signatures.gmt, Oncogenic_Signatures.gmt, Immunologic_Signatures.gmt, Transcription_Factor_Motifs.gmt, MicroRNA_Motifs.gmt, 1Kb_gene_2000.GMT Ontomancer produces one Excel workbook per gene list.  Each workbook has 14 tabs: one per database.  A Fisher Exact Test is performed between the input gene list and gene sets associated with each ontological term. The following columns are provided for each database analysis: Term:                           The ontological term name URL:                            Hyperlink to a description of the term Pathway Genes:         Number of genes associates with pathway or term Query Genes:             Number of genes in input list Overlap:                      Number of genes in common between term and input Odds Ratio:                 Fisher Exact Odds ratio Raw p:                         Fisher Exact p value Adjusted p:                  Benjamini Hochberg correct p value for number of tests in database Score:                         Pythagorean distance between adjusted p and log(odds). (See Hoffman et al 2025) Matching Genes:        List of overlapping genes between gene set and input list. Only terms with odds ratio&gt; 1 and adjusted p &lt;0.25 are reported in Excel. A folder named ‘Raw_Results’ is also generated with a tab delimited file per input list containing all results regardless of odds ratio or p value. The Python script accepts a lists file which is a tab (.txt) or comma (.csv) delimited text file. An example of a lists file is provided (gene_lists_example.TSV). The first column contains a list of all genes expressed in the experiment (background list)*.  Any number of other columns are then added containing query lists.  For example, a query list may be all genes that go up under some criterion and another may be all genes that go down.  The first row is the header and must have unique names for each column.  These names will be used to generate file names so do not use special characters, spaces, period etc in column names.  \* The background list is used to define the size of the Universal set of all genes.   If a background list of expressed genes is not available, enter !! in cell A2 and Ontomancer will use N=22400 for size of background. ### Ontomancer\_network\_maker.py This script will produce a network of ontological terms that encompasses terms from all databases in an Ontomancer-generated workbook. ###### Drag ontomancer\_network.py into terminal or enter path to file to start and follow instructions.  As input, it requires an Ontomancer Excel workbook and 3 user supplied parameters: List of databases:       Which of the ontomancer databases in workbook to use in the analysis. Resolution:                  The resolution parameter for leiden clustering of nodes in the generated network.  The higher the resolution, the more communities/clusters. Beta:                           The power factor to scale edge strengths.  The higher the beta value, the more dispersed the network will be. A similarity matrix is produced for every term in the list of databases (i) versus every term in the list of databases (j).  The values are fisher exact log(odds) for the overlap of genes in the Matching Genes column of term i and j divided by the maximum log(odds) in the matrix to scale values from 0 to 1.  These values are transformed to weights by raising to power Beta (default=6).   The Louvain clustering algorithm with resolution = Resolution (default=1) and weight from Python’s Networkx is used to define communities in the network.  For each community, the term with the highest Score in the List of databases donates its name to the community.  Three networks are generated where node size is set by adjusted p, log(odds) or Score in the orginal ontomancer list of databases. Along with each network, 3 html bar graphs are generated showing original ontomancer log(odds), adjusted p and Score for the highest scoring term in each community.  Finally, 3 tab delimited text files are generated with node, edge and network information. Resolution and beta values are repeatedly asked for in the terminal so that the most pleasing network can be generated.  An example ontomancer workbook is provided (ontomancer_workbook_example.xls).  Try using databases 1,2,3,6 as List of databases to start with.