python： 使用tangram进行空间转录组映射表达量分析

import scanpy as sc
import tangram as tg

xenium = sc.read_10x_mtx("data/Xenium_FFPE_Human_Cancer_Rep1_outs/cell_feature_matrix")
adata  = sc.read_10x_h5('data/SC3pv3_GEX_Breast_Cancer_DTC_Aggr_count_filtered_feature_bc_matrix.h5')

xenium.var_names_make_unique()
adata.var_names_make_unique()

# filter
sc.pp.filter_cells(adata, min_genes=3)
sc.pp.filter_genes(adata, min_cells=3)

sc.pp.filter_cells(xenium, min_genes=3)
sc.pp.filter_genes(xenium, min_cells=3)

# normalization
adata_raw = adata
sc.pp.normalize_total(adata, target_sum = 1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata)
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value  =10)
sc.tl.pca(adata)
sc.pp.neighbors(adata)
sc.tl.leiden(adata, resolution = 0.5)
adata_raw.obs['leiden'] = adata.obs.leiden

tg.pp_adatas(adata_raw, xenium, genes=None) # genes = None, using all genes
ad_map = tg.map_cells_to_space(adata_raw,xenium)
ad_ge = tg.project_genes(ad_map, adata_raw)

tg.pp_adatas(adata_raw, xenium, genes=None) # genes = None, using all genes
ad_map = tg.map_cells_to_space(
  adata_raw,
  xenium, 
  device = 'cpu', 
  mode   = 'clusters', # 'cell', 'clusters', 'constrained'
  cluster_label = 'leiden', 
  density_prior = 'rna_count_based' # 'rna_count_based' or 'uniform', can also be a ndarray
)
ad_ge = tg.project_genes(ad_map, adata_raw, cluster_label = 'leiden')

import scanpy as sc
import tangram as tg
import pandas as pd

# script to run tangram
# according to paper of QUNKUNlab, raw-raw is recommended
def run_tangram(
    singlecell: sc.AnnData,
    sp: sc.AnnData,
    trainning_gene=None,
    predicted_gene=None,
    do_norm: bool = False,
    do_log: bool = False,
    device: str = "cpu",
    mode: str = 'clusters',
    cluster_label: str = 'leiden',
    resolution: int = 0.5,
    density_prior: str = 'rna_count_based'
):
    """
    Function to run tangram.

    Args:
        do_norm: whether to do normalize, default is False
        do_log:  whether to do log1p, default is False
        trainning_gene: genes using to trainning model
        predicted_gene: genes to predict, if not none,  only return dataset of predicted_gene

    Returns:
        a pandas.DataFrame, row is cells and column is gene
    """

    singlecell = singlecell.copy()
    sp = sp.copy()

    # norm
    if do_norm:
        sc.pp.normalize_total(singlecell, target_sum=1e4)
        sc.pp.normalize_total(sp, target_sum=1e4)
    if do_log:
        sc.pp.log1p(singlecell)
        sc.pp.log1p(sp)

    # saving gene name, gene will changed to lower-case after running pp_adatas
    gene_name = singlecell.var_names

    # genes = None, using all genes
    tg.pp_adatas(singlecell, sp, genes=trainning_gene)

    if mode == 'clusters':
        sc_raw = singlecell
        sc.pp.normalize_total(singlecell, target_sum=1e4)
        sc.pp.log1p(singlecell)
        sc.pp.highly_variable_genes(singlecell)
        singlecell = singlecell[:, singlecell.var.highly_variable]
        sc.pp.scale(singlecell, max_value=10)
        sc.tl.pca(singlecell)
        sc.pp.neighbors(singlecell)
        sc.tl.leiden(singlecell, resolution=resolution)
        sc_raw.obs[cluster_label] = singlecell.obs.leiden

        ad_map = tg.map_cells_to_space(
            sc_raw,
            sp,
            device=device,
            mode=mode,  # 'cell', 'clusters', 'constrained'
            cluster_label=cluster_label,
            density_prior=density_prior  # 'rna_count_based' or 'uniform', can also be a ndarray
        )
        ad_ge = tg.project_genes(ad_map, sc_raw, cluster_label=cluster_label)
    else:
        ad_map = tg.map_cells_to_space(singlecell, sp, device=device)
        ad_ge = tg.project_genes(ad_map, singlecell)

    assert all(ad_ge.var_names.str.lower() == gene_name.str.lower())

    predicted_dat = pd.DataFrame(
        ad_ge.X, index=ad_ge.obs_names, columns=gene_name)


    if not predicted_gene is None:
        predicted_dat = predicted_dat.loc[:, predicted_gene]

    return(predicted_dat)

# running analysis
import os
import random
import numpy as np
import pandas as pd
import scanpy as sc


dir_name = "geneImputation_output"

assert os.path.exists(dir_name)
assert os.path.exists('data')


xenium = sc.read_10x_mtx("data/Xenium_FFPE_Human_Cancer_Rep1_outs/cell_feature_matrix")
adata  = sc.read_10x_h5('data/SC3pv3_GEX_Breast_Cancer_DTC_Aggr_count_filtered_feature_bc_matrix.h5')

xenium.var_names_make_unique()
adata.var_names_make_unique()


# filter
sc.pp.filter_cells(adata, min_genes=10)
sc.pp.filter_genes(adata, min_cells=100)

sc.pp.filter_cells(xenium, min_genes=10)  # cell count bigger than 250M
sc.pp.filter_genes(xenium, min_cells=10)  # only 280+ genes


# trainning gene
all_common_genes = sorted(list(set(xenium.var_names) & set(adata.var_names)))
assert len(all_common_genes) > 0

train_ratio = 2/3
train_genes_count = int(train_ratio * len(all_common_genes))
assert train_genes_count > 0

random.seed(1234)
trainning_gene = random.sample(all_common_genes, train_genes_count)
testing_gene = list(set(all_common_genes) - set(trainning_gene))


# subset xenium
if not all(xenium.var_names.isin(all_common_genes)):
    print(f"Some genes ({len(set(xenium.var_names) - set(all_common_genes))}) not existed in singlecell, woule be discared!", flush=True)
    xenium = xenium[:, xenium.var_names.isin(all_common_genes)]

# analysis
for norm in ['raw', 'norm']:
    # 是否norm
    if norm == "raw":
        norm_used = False
    else:
        norm_used = True

    # testing_gene
    if not norm_used:
        testing_gene_used = None
    else:
        testing_gene_used = testing_gene

    # run tangram
    dat_tg = run_tangram(
        adata,
        xenium, 
        do_norm=norm_used,
        do_log=norm_used,
        mode='clusters',
        trainning_gene=trainning_gene, 
        predicted_gene=testing_gene_used
    )
    if not norm_used:
        # 直接导出数据
        save_csv(dat_tg.loc[: , testing_gene] , 'tangram_in_' + str(norm) + '.csv', dir_name=dir_name)
   else:
        save_csv(dat_tg, 'tangram_in_' + str(norm) + '.csv', dir_name=dir_name)