Inferring histology-associated gene expression gradients in spatial transcriptomic studies

Abstract

Spatially resolved transcriptomics has revolutionized RNA studies by aligning RNA abundance with tissue structure, enabling direct comparisons between histology and gene expression. Traditional approaches to identifying signature genes often involve preliminary data grouping, which can overlook subtle expression patterns in complex tissues. We present Spatial Gradient Screening, an algorithm which facilitates the supervised detection of histology-associated gene expression patterns without prior data grouping. Utilizing spatial transcriptomic data along with single-cell deconvolution from injured mouse cortex, and TCR-seq data from brain tumors, we compare our methodology to standard differential gene expression analysis. Our findings illustrate both the advantages and limitations of cluster-free detection of gene expression, offering more profound insights into the spatial architecture of transcriptomes. The algorithm is embedded in SPATA2, an open-source framework written in R, which provides a comprehensive set of tools for investigating gene expression within tissue.

Fig. 1. Histologically defined group-based analysis.

a–c SPATA’s manual annotation tool is employed to delineate borders, facilitating the grouping of spots based on histology. d Gene expression-driven UMAP projection of spots. e A dot plot showcases the 18 most statistically significant unique marker genes, ranked by their average log2-fold change, in accordance with histological areas. f and g Surface plots are color-coded to highlight inferred copy number alterations that are characteristic of glioblastoma. h and i Statistical analyses examine copy number alterations across histological areas (two-sided test, no adjustments for multiple comparisons, Tumor (n = 1307), Transition (n = 478), Infiltrated Cortex (n = 1428). The minima represent the smallest and the maxima represent the largest value within 1.5 times the interquartile range (IQR) below or above the first or third quartile (Q1, Q3), respectively. The median is shown as a line inside the box. The box bounds are the first quartile (Q1, 25th percentile) and the third quartile (Q3, 75th percentile). The whiskers extend from Q1 and Q3 to the minima and maxima. j A heatmap provides a comprehensive view of alterations across all chromosomes in relation to histological areas, corroborating the tumor area’s classification with multiple alterations, the nearly unaltered cortex infiltrated by the tumor, and the transition zone exhibiting intermediate levels of alterations. k A volcano plot from the DEA analysis across histological areas highlights marker genes for the Tumor and Transition areas, characterized by an adjusted p-value of 0 (infinite −log10).

Fig. 2. Spatial expression analysis of genes identified by differential expression analysis (DEA) along a linear trajectory.

a A surface plot provides a visual representation of the trajectory’s course (indicated by an arrow), complemented by a barplot illustrating the proportion of histological groups along this trajectory. b and c Surface plots offer insights into the gene expression patterns of genes predominantly localized within their respective assigned areas. Error bands of line plots indicate the confidence interval (level: 0.95). d–g Heatmaps present the expression profiles of the top 18 differentially expressed genes (DEGs) per histological group (as referenced in Fig. 1e). Notably, while these genes have been identified as unique DEGs, those from the tumor and transition regions do not consistently exhibit confined expression patterns to their designated areas, in contrast to DEGs from the (infiltrated) cortex. e–h Line plots and surface plots further elucidate the expression patterns of selected genes, demonstrating how they may traverse area boundaries in various ways. Error bands of line plots indicate the confidence interval (level: 0.95).

Fig. 3. Integration of visium mouse brain dataset for gene co-expression analysis and eq data examination.

a An H&E image provides an overview of the analyzed sample, with two enlarged windows highlighting the stab wounds inflicted on the brain (n = 1 mouse with bilateral injury for Visium, n = 3 mice for scRNA-seq experiments). b Surface plots reveal clustering results obtained through the Scanpy pipeline, emphasizing significant barcode spot clusters within the injury area. c and d Surface plots and gradient ridge plots illustrate the gene expression patterns of two marker genes for the injury area, Hmox1 and Lcn2. e and f Ridgeplots visualize the expression gradients of genes sharing similar patterns with either of the two example genes, as identified through spatial annotation screening. g–j Surface plots showcase the expression profiles of co-expressed genes identified via spatial annotation screening. k UMAP of the scRNA-seq dataset, from the same mouse model as the Visium sample. l and m Visualization of Tangram results, featuring 2D Density plots that highlight an elevated density of monocytes, macrophages, and microglia in the vicinity of the injury zone, as compared to a control area.

Fig. 4. Comprehensive annotation of histology and gene expression relative to necrosis in glioblastoma.

a Presents the H&E image of the glioblastoma sample UKF313T, emphasizing key areas within the sample that are detailed in the following figure (b). Additionally, a surface plot visualizes the count distribution, supporting our necrosis (dead tissue) annotation. b Provides a visualization of distance values from the necrotic areas, with lines illustrating the assumed gradient direction and orientation based on proximity to necrotic regions. c–e Showcase representative genes exhibiting a pattern reminiscent of association with necrosis, displaying decreasing expression levels with increasing distance from necrotic regions. f–h Feature representative genes showing a pattern resembling the recovery of expression levels with increasing distance from necrotic areas.

Fig. 5. Comprehensive and integrated analysis of spatial relationships between hypoxia and T-cell abundance in glioblastoma.

a Overview of the workflow, encompassing spatial transcriptomic (ST) data, single-cell deconvolution, and spatial T-cell receptor sequencing (SPTCR-seq), showcasing the horizontal integration of six glioblastoma samples featuring prominent hypoxic spatial niches. b A representative Visium glioblastoma sample (UKF260T) with multiple hypoxic areas, annotated using SPATA2’s automatic annotation tool. Lines indicate the screening direction guided by the hypoxic gradient. c Presentation of single-cell deconvolution results for sample UKF260T, highlighting proximity to annotated hypoxic regions. d The inferred gradient of hypoxic gene signatures merged from data obtained from six samples. Error bands of line plots indicate the confidence interval (level: 0.95) e and f Evaluation of T-cell and anti-inflammatory bone-derived macrophage abundance as a function of distance from the hypoxic areas. Error bands of line plots indicate the confidence interval (level: 0.95). f Abundance of cell types described by Neftel et al. (2019), revealing distinct abundance patterns. g Visualization of MES-like and NPC-like cell abundance in sample UKF269T. h An illustration of a glioblastoma example integrated into the comprehensive screening, featuring a solitary hypoxic area and a separate area of T-cell abundance approximately 1 mm distant from the hypoxic region. i Gradient representation of T-cell abundance as a function of distance from hypoxia. j Gradient plot depicting various T-cell subtypes as a function of distance from hypoxia, revealing an inverse correlation between T-cell cytotoxicity and distance to hypoxia, peaking at ~1 mm.

Fig. 6. Statistical comparison of the test statistic from SPATA2’s spatial gradient screening with those of analogous methods PseudotimeDE and tradeSeq.

a–c Scatterplots display the relationship between different test statistics used by different algorithms to determine the noise ratio of simulations. Colors indicate the type of simulated pattern hidden by noise (also see Supplementary Figs. 9b, 10b).