代码更新--高精度空间（Xenium、CosMx）细胞外基因表达的数据分析

import spatialdata as sd
import troutpy as tp
import scanpy as sc
import numpy as np
###Xenium数据
sdata=sd.read_zarr('./data.zarr') #please modify the path to your conveniece

###Processing cellular information & cell type definition
adata = sdata["table"]
adata.raw = adata.copy()

all_cells = adata.shape[0]
sc.pp.filter_cells(adata, min_genes=10)
sc.pp.filter_cells(adata, min_counts=30)
sc.pp.filter_genes(adata, min_cells=3)

print(f"Proportion of cells retained after filtering: {adata.shape[0] / all_cells:.2%}")

adata.layers["raw"] = adata.X.copy()
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)

sc.tl.pca(adata)
sc.pp.neighbors(adata, n_neighbors=10)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="leiden", resolution=0.6)

sc.tl.rank_genes_groups(adata, groupby="leiden")
sc.tl.dendrogram(adata,groupby='leiden')
sc.pl.rank_genes_groups_dotplot(adata, groupby="leiden",n_genes=4,cmap='Blues')

dict_anno={0:'DG',1:'Oligodendrocytes',2:'CA1-2',3:'Astrocytes',4:'Endothelial',5:'EXC neurons Thal.',6:'SMC/Pericytes',7:'CA3',8:'EXC. neurons CTX',
9:'Ependymal cells',10:'Microglia',11:'INH neurons',12:'OPC',13:'Astrocytes',14:'SGZ Neuroblasts',15:'Cajal-Retzius',16:'Choroid plexus',17:'INH neurons CP',
18:'Oligodendrocytes'}
adata.obs['cell type']=adata.obs['leiden'].astype(int).map(dict_anno)
sc.pl.umap(adata,color='cell type',palette='tab20')
sc.pl.spatial(adata,spot_size=30,color='cell type',palette='tab20')

tp.pp.segmentation_free_sainsc(sdata,binsize=5,celltype_key="leiden",background_filter=0.2,gaussian_kernel_key=2.0)

# 3) Histogram of scores
tp.pl.histogram(sdata, x="cosine_similarity", hue="overlaps_cell", palette="troutpy")

# 4) Define extracellular
tp.pp.define_urna(sdata, layer="transcripts", method="sainsc", percentile_threshold=5)

  # 5) Crosstab / pie / heatmap
tp.pl.crosstab(sdata,yvar="extracellular",xvar="overlaps_cell",normalize=False,cmap="troutpy",kind="barh",stacked=True,figsize=(5, 2))

control_codewords = ["negative_control_probe", "unassigned_codeword", "deprecated_codeword", "genomic_control_probe", "negative_control_codeword"]

tp.tl.quantify_overexpression(sdata,layer="transcripts",codeword_key="codeword_category",control_codewords=control_codewords,gene_key="gene",
percentile_threshold=99.5)
tp.pl.logfoldratio_over_noise(sdata, test_method="auto")

tp.pl.metric_scatter(
            sdata,
            x='count',
            y='logfoldratio_over_noise',
            label_top_n_x=7,
            label_top_n_y=5,
            label_bottom_n_x=0,
            label_bottom_n_y=0, size=20
        )

tp.tl.extracellular_enrichment(sdata)

tp.pl.metric_scatter(
            sdata,
            x='extracellular_proportion',
            y='logfoldratio_over_noise',
            label_top_n_x=7,
            label_top_n_y=5,
            label_bottom_n_x=0,
            label_bottom_n_y=0, size=5,linewidth=0.1
        )

tp.pp.aggregate_urna(sdata, gene_key="gene", square_size=50)
tp.tl.spatial_variability(sdata,gene_key="gene",n_neighbors=10)
tp.pl.spatial_transcripts(sdata,gene_list=['Adcy1'],scatter_size=2, boundary_linewidth=0.1)

sdata['transcripts']['gene']=sdata['transcripts']['gene'].astype('category')
tp.pp.aggregate_urna(sdata, gene_key="gene", square_size=50)
tp.tl.in_out_correlation(sdata, n_neighbors=20)
 tp.pl.metric_scatter(
            sdata,
            x='moran_I',
            y='in_out_spearmanR',
            label_top_n_x=0,
            label_top_n_y=3,
            label_bottom_n_x=0,
            label_bottom_n_y=3, size=5,linewidth=0.1
        )

tp.pl.spatial_transcripts(sdata,gene_list=['Ptk2b'],scatter_size=2, boundary_linewidth=0.1)

tp.pp.filter_urna(sdata, min_logfoldratio_over_noise=7, gene_key="gene") # filter out genes with a logfoldratio over noise below 7.
tp.pp.filter_urna(sdata,gene_key="gene",control_probe=True) # filter control probes

#### Identifying uRNAs source
 tp.tl.compute_source_score(sdata,layer="transcripts",gene_key="gene",
        lambda_decay=0.1,celltype_key="cell type",n_jobs=1)

tp.tl.assess_diffusion(sdata, gene_key="gene", distance_key="distance")
tp.pl.diffusion_results(sdata,label_top_n_y=10,y_logscale=True)

tp.pl.gene_distribution_from_source(sdata,['Camk2b'])

tp.pl.source_score_by_celltype(sdata,figsize=(6,6),cmap='terrain')

sc.pl.umap(sdata['table'],color=['urna_source_score','cell type'],cmap='terrain',title='uRNA source score',vmax='p99.995')

sdata.write('./data_filt.zarr')