作者，Evil Genius
在昨天的文章细节分享--关于高精度（Xenium、CosMx）细胞分割外基因表达的分析讨论中表明uRNA可能携带有生物学意义的信号，值得进一步探讨。
而我们今天分享的troutpy，是一个用于空间转录组学数据中定量探索uRNA的Python包。
示例代码分享一下，现在空间分析的代码几乎全是python，也是对大家的要求。
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')Segmentation-free analysis
首先将组织中的每个转录本按指定微米级像素尺寸（binsize）进行空间分箱。同时根据background_filter参数对背景信号进行过滤。由于空间数据集通常在特定分辨率下较为稀疏，还会通过gaussian_kernel_key参数调整高斯模糊对信号进行平滑处理。
在完成分箱、过滤和平滑处理后，计算每个像素的基因表达特征与预定义细胞类型特征之间的余弦相似度。细胞类型特征通过sdata['table'].obs中指定列（由celltype_key参数定义）的细胞类型平均表达值来确定。
针对每个像素，将获得其与所有细胞类型的余弦相似度值。为简化分析，仅保留每个像素与任意细胞类型的最大余弦相似度值及其对应的细胞类型信息。最终，这些计算结果将以新列的形式添加到sdata['transcripts']的每个转录本信息中。
tp.pp.segmentation_free_sainsc(sdata,binsize=5,celltype_key="leiden",background_filter=0.2,gaussian_kernel_key=2.0)由于分析包含了所有转录本，因此可以对比已分割至细胞的转录本与未分配转录本的余弦相似度。这将使我们能够将未分配的RNA分为两类：
第一类：在局部特征上与已分割细胞相似的RNA（类细胞RNA）
第二类：具有与已分割细胞不同局部组成的RNA（其他未分配RNA）
# 3) Histogram of scores tp.pl.histogram(sdata, x="cosine_similarity", hue="overlaps_cell", palette="troutpy")结合分割分析和免分割分析的结果，对非类细胞型转录本进行处理。
细胞外转录本将通过以下两个条件来定义：
位于已分割细胞区域之外；
与细胞类型特征图谱的余弦相似度低于设定的百分位数阈值（percentile_threshold），该阈值仅基于细胞像素计算得出。
# 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))Gene-centric exploration of uRNA compopsition
探索未分配RNA的丰度与比例
首先，将量化各类uRNA（即各基因）在未分配RNA库中的丰度，并将其与背景噪音水平进行比较。为此，实现了troutpy.tl.quantify_overexpression方法，用于识别相对于噪音阈值过表达的基因。
该方法的核心原理是：基于指定对照特征的计数计算阈值，并将各基因计数与该阈值比较以判定过表达情况。该函数会计算每个基因的对数倍数变化，将结果标注于元数据中，最终返回更新后的空间数据以及各基因的得分与计算阈值。
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")探究哪些基因相较于噪音水平具有最高的富集程度。为此，我们通过散点图进行可视化展示：图中每个点代表一个基因，横轴为未分配RNA（uRNA）计数，纵轴为相较于噪音水平的对数比值（logFoldRatio）。通过该可视化分析，我们识别出特定基因在未分配RNA库中呈现高丰度特征。
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)可以并排可视化两个指标——uRNA丰度与uRNA占比，并观察不同基因呈现的差异模式：
低uRNA/细胞外占比的基因（如Aldoa）：这类基因仅表现为高丰度基因
高uRNA占比但低uRNA丰度的基因：提示存在未通过细胞分割充分捕获的uRNA表达模式
同时具有高uRNA丰度与高占比的基因：代表在已分割细胞外显著富集的高表达基因，呈现清晰的特征模式
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 )探索空间变异性与未分配RNA分布模式
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)未分配RNA基因筛选
仅保留明确高于噪音水平且非对照探针的基因。
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)Assesing diffusion
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')Cell-specific uRNAs source score
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')