SpaGFT Tutorial：Visium Human Lymph Node

Outline

1. Import packages

2. Load data

3. QC and preprocessing

4. Function: identify spatially variable genes

5. Function: gene expression enhancement

6. Function: characterize functional tissue units

The installation steps can be found here. SpaGFT is a python package to analyze spatial transcriptomics data via graph Fourier transform. To install SpaGFT, the python version is required to be >= 3.7. You can check your python version by:

import platform
platform.python_version()

'3.8.18'

1. Import packages

import SpaGFT as spg
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
import os


sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')

scanpy==1.9.1 anndata==0.9.2 umap==0.5.5 numpy==1.21.5 scipy==1.8.0 pandas==1.4.2 scikit-learn==1.0.2 statsmodels==0.14.1 python-igraph==0.9.10 louvain==0.7.1 pynndescent==0.5.11

To begin with, define the folder to save the results which are generated from the following analysis:|

results_folder = './results/lymph_nodes_analysis/'
if not os.path.exists(results_folder):
    os.makedirs(results_folder)

2. Load data

This tutorial uses a publicly available Visium dataset of the human lymph node. The dataset can be downloaded easily through Scanpy. Two key elements are essential for running SpaGFT: the raw gene counts matrix and spatial coordinates of spots. In addition, these two elements should be integrated into an anndata object, which is a frequently-used python object.

# Load data
adata = sc.datasets.visium_sge(sample_id="V1_Human_Lymph_Node")
adata.var_names_make_unique()
# Add raw to anndata object
adata.raw = adata

You can also create your datasets which are anndata objects. Note that the raw count matrix should be found in adata.X, and the spatial coordinates of all spots should be found in adata.obs or adata.obsm. You can use diverse functions provided by scanpy to create the anndata object. For example, you can read Visium datasets by

sc.read_visium() # recommend

Alternatively, users can also create the anndata object using the raw count matrix and spatial coordinates of spots in a direct way:

count_mtx = pd.read_csv(PATH_TO_COUNT_MATRIX)
coord_mtx = pd.read_csv(PATH_TO_COORDINATE_MATRIX)
count_mtx.iloc[:5, :5]
coord_mtx.iloc[:5, :5]
# create the anndata object
adata = sc.AnnData(count_mtx)
adata.obs.loc[:, ['x', 'y']] = coord_mtx
adata.var_names_make_unique()
adata.raw = adata

More approaches to load datasets can be found at Scanpy

3. QC and preprocessing

Before the formal analysis, we perform some basic filtering of genes and normalize Visium counts data with the normalize_total method from Scanpy followed by log-transform.

# QC
sc.pp.filter_genes(adata, min_cells=10)
# Normalization
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)

4. Function: identify spatially variable genes

4.1 Determine the Fourier modes and identify SVGs

To begin with, determine the number of Fourier modes for detecting SVGs. SpaGFT provides the dermine_frequency_ratio function based on the kneedle algorithm to determine how many FMs used automatically. Note that ratio_low \(\sqrt{n}\)* and ratio_high \(\sqrt{n}\)* are the recommended low-frequency and high-frequency FMs, where :math:`sqrt{n}` is the number of spots in the original dataset.

# determine the number of low-frequency FMs and high-frequency FMs
(ratio_low, ratio_high) = spg.gft.determine_frequency_ratio(adata,
                                                            ratio_neighbors=1)

Obatain the Laplacian matrix

# calculation
gene_df = spg.detect_svg(adata,
                         spatial_info=['array_row', 'array_col'],
                         ratio_low_freq=ratio_low,
                         ratio_high_freq=ratio_high,
                         ratio_neighbors=1,
                         filter_peaks=True,
                         S=6)
# S determines the  sensitivity of kneedle algorithm
# extract spaitally variable genes
svg_list = gene_df[gene_df.cutoff_gft_score][gene_df.fdr < 0.05].index.tolist()
print("The number of SVGs: ", len(svg_list))
# the top 20 SVGs
print(svg_list[:20])

The precalculated low-frequency FMs are USED
The precalculated high-frequency FMs are USED
Graph Fourier Transform finished!
svg ranking could be found in adata.var['svg_rank']
The spatially variable genes judged by gft_score could be found
          in adata.var['cutoff_gft_score']
Gene signals in frequency domain when detect svgs could be found
          in adata.varm['freq_domain_svg']
The number of SVGs:  1346
['IFI44L', 'ISG15', 'MT-CO2', 'IFI6', 'MX1', 'MT-CO1', 'EEF1A1', 'CXCL12', 'TMSB4X', 'ACTG1', 'IFI44', 'IFIT1', 'CCL2', 'GAPDH', 'MT-CO3', 'RPL32', 'OASL', 'RPS2', 'SFRP2', 'ACTB']

4.2 Visualize the identified SVGs

We will visualize several detected genes which exhibit diverse spatial patterns. This step can be achieved by sc.pl.spatial function provided by scanpy or our built-in scatter_gene function to visualize gene expression patterns.

plot_svgs = ['CD3E', 'IL7R', 'CCR7', 'PCNA', 'CDK1', 'CDC20', 'CD19', 'CD79B']
sc.pl.spatial(adata, color=plot_svgs, size=1.6, cmap='magma', use_raw=False)

4.3 Plot frequency signals

One of the characteristics of SpaGFT is that it can obtain new representations for genes in the frequency domain. In the above genes’ normalized frequency signals can be found in adata.varm['freq_domain_svg'] and visualized by

spg.plot.gene_freq_signal(adata, gene=plot_svgs[0])
spg.plot.gene_freq_signal(adata, gene=plot_svgs[4])

where the red signals correspond to low-frequency signals, while the blue signals correspond to high-frequency signals.

4.4 UMAP visualization of gene frequency signals

In addition, we can also separate SVGs and non-SVGs in the frequency domain. By projecting the frequency domain to 2-dimensional space by UMAP, we can verify this easily:

# gene umap
spg.plot.gene_signal_umap(adata, svg_list, size=30)
svg_umap_df = pd.DataFrame(adata.varm['freq_umap_svg'],
                           index=adata.var_names,
                           columns=['UMAP_1', 'UMAP_2'])

The umap coordinates of genes when identify svgs could be found in
          adata.varm['freq_umap_svg']

5 Function: gene expression enhancement

SRT expression data suffers from high noise, which will influence the credibility of analytical results. Therefore, SpaGFT adapts a low-pass filter to process signals in the frequency domain and performs inverse GFT to obtain enhanced signals. Note that this step is time-consuming and will take several minutes for Visium datasets.

# copy from the original dataset
new_adata = adata.copy()
new_adata = spg.low_pass_enhancement(new_adata,
                                     inplace=True,
                                     c=0.001,
                                     ratio_low_freq=15)

where, parameter c controls the balance of the smoothness and similarity with the original dataset. For the fine-grained tissue, we recommend a low c to preserve more details.

# Before enhancement
sc.pl.spatial(adata, color=plot_svgs[:4], img_key=None, size=1.6, cmap='magma', use_raw=False)

# After enhancement
sc.pl.spatial(new_adata, color=plot_svgs[:4], img_key=None, size=1.6, cmap='magma', use_raw=False)

In the following, we will use the enhanced gene expression matrix for analysis.

6. Function: identify functional tissue units

6.1 Group SVGs based on frequency signals

The core of SpaGFT is that it can utilize SVGs to characterize functional tissue units by groups of SVGs with coincident spatial patterns. To achieve this, SpaGFT will analyze and process the signal’s low-frequency parts. In SpaGFT, the Louvain algorithm is used to cluster frequency signals of genes whose parameter resolution will be selected based on the designed objective function. Users can input the range of resolution by giving (start, end, step).

# identify functional tissue units (FTU)
gene_df, new_adata = spg.identify_ftu(new_adata,
                                     svg_list=svg_list,
                                     ratio_fms=ratio_low,
                                     resolution=(0.5, 1.5, 0.1),
                                     ratio_neighbors=2,
                                     n_neighbors=15)

resolution: 0.500;  score: 0.1675
resolution: 0.600;  score: 0.1675
resolution: 0.700;  score: 0.1671
resolution: 0.800;  score: 0.1676
resolution: 0.900;  score: 0.1662
resolution: 1.000;  score: 0.1651
resolution: 1.100;  score: 0.1327
resolution: 1.200;  score: 0.1429
resolution: 1.300;  score: 0.1440
resolution: 1.400;  score: 0.1525

The above step will identify gene groups based on their spatial patterns. and the results can be found at adata.var['tissue_module']. The genes in a group support a common spatial pattern, which is called a functional tissue unit. Besides, the functional tissue unit information can be found at adata.obsm['ftu_binary'].

# the genes which support corresponding functional tissue units
gene_df.to_csv(os.path.join(results_folder,
                           'Lymphnode_gene_information_results.csv'))
gene_df.iloc[:5, :]

	gene_ids	feature_types	genome	n_cells	gft_score	svg_rank	cutoff_gft_score	pvalue	fdr	ftu
IFI44L	ENSG00000137959	Gene Expression	GRCh38	1882	4.973869	1	True	5.813716e-80	6.207737e-78	8
ISG15	ENSG00000187608	Gene Expression	GRCh38	3054	4.877024	2	True	2.888858e-55	1.733219e-53	8
MT-CO2	ENSG00000198712	Gene Expression	GRCh38	4030	4.615628	3	True	2.816433e-100	5.385450e-98	6
IFI6	ENSG00000126709	Gene Expression	GRCh38	2681	4.561138	4	True	9.051660e-80	9.630435e-78	8
MX1	ENSG00000157601	Gene Expression	GRCh38	3326	4.508339	5	True	1.832998e-49	9.583907e-48	8

The last column indicates the supported functional tissue units.

6.2 Visualize clustering results

We can visualize the clustering results by UMAP in this way:

spg.plot.scatter_umap_clustering(new_adata, svg_list)

There are nine modules via our computation. In the following, we will visualize the results.

6.3 Extract functional tissue unit information

# The functional tissue unit information can be obtained by
ftu_binary_df = new_adata.obsm['ftu_binary']
ftu_binary_df.to_csv(os.path.join(results_folder,
                           'Lymphnode_tissue_module_information_results.csv'))
ftu_binary_df.iloc[:5, :]

	ftu_1	ftu_2	ftu_3	ftu_4	ftu_5	ftu_6	ftu_7	ftu_8	ftu_9
AAACAAGTATCTCCCA-1	0	1	0	1	0	0	0	0	0
AAACAATCTACTAGCA-1	1	1	1	0	0	1	0	0	0
AAACACCAATAACTGC-1	1	1	0	1	1	1	1	0	0
AAACAGAGCGACTCCT-1	1	1	0	0	1	1	1	1	0
AAACAGCTTTCAGAAG-1	0	1	0	0	0	1	1	0	0

The rows of the above dataframe indicate the spots, and the columns indicate the functional tissue unit. Among them, ‘1’ reflects the structure of a functional tissue unit.

6.4 Visualize interested FTUs

We can visualize the interested functional tissue units by built-in scatter_ftu function:

# plot ftus
spg.plot.scatter_ftu(new_adata,
                    ftu=['ftu_3', 'ftu_5', 'ftu_7'],
                    spatial_info='spatial',
                    size=15)

Specifically, we can also visualize the functional tissue unit and genes that support this functional tissue unit. In the following, we will explore a specific functional tissue unit, tm_3

# plot tm's spatial map and corresponding genes
spg.plot.scatter_ftu_gene(new_adata,
                         ftu="ftu_3",
                         ftu_color='#9aabe1',
                         gene=['CD3E','IL7R','CCR7'],
                         spatial_info='spatial',
                         size=13)

6.5 Enrichment analysis of biological process

A group of genes that support a functional tissue unit has a similar spatial pattern, and we can perform functional enrichment analysis to explore the underlying biological process under this pattern.

# GO enrichment analysis
import gseapy as gp
ftu_gene_list = gene_df[gene_df.ftu == '3'].index.tolist()
enr = gp.enrichr(gene_list=ftu_gene_list,
                  gene_sets=['BioPlanet_2019','GO_Biological_Process_2021',
                            'ChEA_2016'],
                  organism='Human',
                  description='Tissue_module',
                  outdir='./results/lymph_nodes_analysis/enrichr_kegg',
                  no_plot=False,
                  cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr_results = enr.results
from gseapy.plot import barplot, dotplot
barplot(enr.results[enr.results.Gene_set=='GO_Biological_Process_2021'],
        column='P-value',
        color='#FF6879',
        top_term=5,
        title='GO_Biological_Process_2021')