Spatial transcriptome (Xenium)

We demonstrate a noise reduction with RECODE for spatial transcriptome data (FISH based). We use spatial transcriptome data of Xenium 10X, FFPE Human Pancreatic Ductal Adenocarcinoma Data with Human Immuno-Oncology Profiling Panel. The dataset is available from 10X datasets.

We use scanpy to read/write data. Import numpy and scanpy in addlition to screcode.

import scanpy as sc
import numpy as np
import screcode
import warnings
warnings.simplefilter('ignore')
import matplotlib.pyplot as plt
import pandas as pd

Read in the count matrix into an AnnData object.

INPUT_DIR = 'data/Xenium_V1_Human_Ductal_Adenocarcinoma_FFPE_outs'
INPUT_FILE = "cell_feature_matrix.h5"
Raw_key = "count"
adata = sc.read_10x_h5("%s/%s" % (INPUT_DIR,INPUT_FILE))
adata.obs = pd.read_csv("%s/cells.csv.gz" % INPUT_DIR)
adata.var_names_make_unique()
adata = adata[:,np.sum(adata.X,axis=0)>0]
adata = adata[np.sum(adata.X,axis=1)>0]
adata.layers["Raw"] = adata.X.toarray()
adata

AnnData object with n_obs × n_vars = 234857 × 380
    obs: 'cell_id', 'x_centroid', 'y_centroid', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', 'deprecated_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area'
    var: 'gene_ids', 'feature_types', 'genome'
    layers: 'Raw'

Apply RECODE

Apply RECODE to the count matrix (without using spatial coordinates).

import screcode
recode = screcode.RECODE(seq_target='RNA',version=2)
adata = recode.fit_transform(adata)

start RECODE for scRNA-seq data
end RECODE for scRNA-seq
log: {'seq_target': 'RNA', '#significant genes': 380, '#non-significant genes': 0, '#silent genes': 0, 'ell': 128, 'Elapsed time': '0h 1m 15s 913ms', 'solver': 'randomized', '#test_data': 46971}

Performance check

recode.report()

Log normalizaation

target_sum = np.median(np.sum(adata.layers["RECODE"],axis=1))
adata = recode.lognormalize(adata,target_sum=target_sum)
print(np.median(np.sum(adata.layers["RECODE"],axis=1)))

Normalized data are stored in "RECODE_norm" and "RECODE_log"
46.68497

adata.layers["Raw_norm"] = target_sum*adata.layers["Raw"]/np.sum(adata.layers["Raw"],axis=1)[:,np.newaxis]
adata.layers["Raw_log"] = np.log(adata.layers["Raw_norm"]+1)

Plot spatial gene expression

def spatial_gex(
        genes,
        sp_x = adata.obs["x_centroid"],
        sp_y = -adata.obs["y_centroid"],
        psize = 1,
        figsize=(6,3),
        dpi=100,
        percentiles = [10,90],
        fs_title = 20,
        fs_label = 20,
    ):
    fig,ax = plt.subplots(2,len(genes),figsize=(figsize[0]*len(genes),figsize[1]*2),tight_layout=True)
    for i in range(len(genes)):
        idx_gene = adata.var.index == genes[i]
        if sum(idx_gene) == 0: continue
        exp = adata.layers["RECODE_log"][:,idx_gene]
        vmin,vmax = np.percentile(exp,percentiles[0]),np.percentile(exp,percentiles[1])
        ax_ = ax[1,i]
        ax_.scatter(sp_x, sp_y,c=exp,s=psize,marker="H",vmin=vmin,vmax=vmax)
        if i== 0:
            ax_.set_ylabel("RECODE",fontsize=fs_label)
            ax_.tick_params(bottom=False, left=False, right=False, top=False,
                            labelbottom=False, labelleft=False, labelright=False, labeltop=False)
            [ax_.spines[c_].set_visible(False) for c_ in ['right','top','bottom','left']]
        else:
            ax_.axis('off')


        exp = adata.layers["Raw_log"][:,idx_gene]
        ax_ = ax[0,i]
        ax_.scatter(sp_x, sp_y,c=exp,s=psize,marker="H",vmin=vmin,vmax=vmax)
        ax_.set_title("$\it{%s}$" % genes[i],fontsize=fs_title)
        if i== 0:
            ax_.set_ylabel("Raw",fontsize=fs_label)
            ax_.tick_params(bottom=False, left=False, right=False, top=False,
                            labelbottom=False, labelleft=False, labelright=False, labeltop=False)
            [ax_.spines[c_].set_visible(False) for c_ in ['right','top','bottom','left']]
        else:
            ax_.axis('off')


GENES = ["A2M", "ACTA2", "ACTB", "ANPEP", "APOE", "C1R", "C1S", "CD163", "CD44", "CD74", "CEACAM1", "CEACAM6", "CEACAM8", "CTSD", "CXCL1",
         "CXCL12", "CXCL13", "CXCL5", "DCN", "DMBT1", "EPCAM", "FCMR", "FLT1", "FN1", "G0S2", "IFITM3", "IFNGR1", "IGFBP7", "IGHG1", "IGHG2",
         "IGHG3", "IGHG4", "IGHGP", "IGHM", "IGKC", "IL1B", "IL1R1", "IL7R", "ITGAX", "JCHAIN", "LUM", "MS4A1", "MUC5AC", "PLVAP", "REG4", "RGS5",
         "S100A9", "SPARC", "SPARCL1", "TFF3", "TNFRSF13C", "VCAN", "XBP1"]

n_plots = 4
for i in range(int(len(GENES)/n_plots+0.9)):
    if (i+1)*n_plots < len(GENES):
        genes = GENES[i*n_plots:(i+1)*n_plots]
    else:
        genes = GENES[i*n_plots:len(GENES)]
    spatial_gex(genes = genes)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[44], line 53
     51 else:
     52     genes = GENES[i*n_plots:len(GENES)]
---> 53 spatial_gex(genes = genes)

Cell In[44], line 18, in spatial_gex(genes, sp_x, sp_y, psize, figsize, dpi, percentiles, fs_title, fs_label)
     16 exp = adata.layers["RECODE_log"][:,idx_gene]
     17 vmin,vmax = np.percentile(exp,percentiles[0]),np.percentile(exp,percentiles[1])
---> 18 ax_ = ax[1,i]
     19 ax_.scatter(sp_x, sp_y,c=exp,s=psize,marker="H",vmin=vmin,vmax=vmax)
     20 if i== 0:

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

['XBP1']