课前准备----空间转录组hotspot

原创

追风少年i

发布于 2024-07-29 11:59:15

840

发布于 2024-07-29 11:59:15

作者，Evil Genius

参考文章Hotspot identifies informative gene modules across modalities of single-cell genomics（cell systems）。

Hotspot是一个基于图的程序来识别信息基因和基因模块

当我们需要在空间上衡量两个指标的关联度的时候，就需要借助hotspot。hotspots是探索组织微环境的一种手段，就像肿瘤分为“冷”肿瘤和“热”肿瘤，其微环境必然存在很大的差异，在比例肿瘤转移，并不是所有的肿瘤细胞都有转移能力，具有转移能力的肿瘤细胞通常和巨噬细胞“狼狈为奸”，这个例子中具有转移性的肿瘤细胞就属于hotspot，其“hot”的特征就是周围富集了促进其转移的巨噬细胞。

空间转录组学在描述具有不同细胞组成的主要组织结构域、癌症特征、免疫抑制中心、具有不同临床结果的整个肿瘤生态型或特定药物对抑制肿瘤进展的影响方面的优势。（在WES的研究中，基因分为了免疫正相关、免疫负相关、肿瘤易感基因等，大家需要了解）。定义:特定细胞类型密集"居住"或由自定义的基因特征标记的区域(热点)和细胞类型或基因特征耗尽的区域(冷点)，并统计评估这些区域与其他预定义热点或冷点的接近程度。此外，使用热点方法检测到的空间关系与仅观察单个点的邻近区域时观察到的空间关系进行了比较。评估这些变量之间的关系如何变化，hotspot大小或周围的一个spot的大小变化。这里的hotspot不仅指细胞类型，还有通路等因素。分析模块如下： 邻域富集分析（Neighbourhood enrichment analysis）:用于检查单个空间转录组学spot及其邻近区域内细胞状态、细胞类型或过程之间的相关性。在这种情况下，一个“neighbourhood”被描绘成一个环，围绕一个指定的中心spot包含六个Visium spot，通过将空间点视为一个网络来计算。 Hotspot identification：分析揭示了统计上显著的“hotspots”或“coldspots”。hotspots代表特定细胞类型或特征高度集中的区域，表明这种聚集不太可能是偶然的。相反，coldspots表示目标细胞或特征稀少的区域，也超出了随机分布的预期。 距离度量：测量和解释识别具体特征之间的距离，例如肿瘤和免疫热点之间的距离。这方面的主要方法是计算到hotspots的最短路径，定义为从定义的hotspots内的任何spot到指定比较热点内最近点的最小距离。 敏感性分析:通过改变邻域或感兴趣的hotspot的区域大小，系统地评估细胞-细胞关系如何在组织内进化的能力。邻域富集方法，这可以通过改变中心点周围同心圆的数量来评估。热点方法，计算不同邻域大小的特征统计量，从而能够在不同的空间尺度上识别niche；以及空间hotspot的距离变化关系。

而我们最终需要拿到的结果

其生物学意义我们课上讨论，完整的代码如下，当然了，大家可以参考官网代码，文章有提供，不过更多需要对代码更进一步：

import scanpy as sc
import squidpy as sq
import numpy as np
import pandas as pd
import os
import sys
import seaborn as sb
import matplotlib.pyplot as plt
from matplotlib import colors

#import scvi
import anndata as ad

import warnings
warnings.filterwarnings("ignore")

from collections import Counter

import ipywidgets as widgets
from ipywidgets import interact, interact_manual

plt.rcParams['figure.figsize'] = (6, 6)

from IPython.core.display import display, HTML
import random

#Define a colour map for gene expression
colors2 = plt.cm.Reds(np.linspace(0, 1, 128))
colors3 = plt.cm.Greys_r(np.linspace(0.7,0.8,20))
#colorsComb = np.vstack([colors3, colors2])
#mymap = colors.LinearSegmentedColormap.from_list('my_colormap', colorsComb)
from matplotlib import colors
colorsComb = np.vstack([plt.cm.Reds(np.linspace(0, 1, 128)), plt.cm.Greys_r(np.linspace(0.7, 0.8, 0))])
mymap = colors.LinearSegmentedColormap.from_list('my_colormap', colorsComb)

# Helper function to split list in chunks
def chunks(lista, n):
    for i in range(0, len(lista), n):
        yield lista[i:i + n]
        
        plt.rcParams['figure.figsize'] = (6, 5)
sc.set_figure_params(dpi=100, vector_friendly=True)
def mysize(w, h, d):
    fig, ax = plt.subplots(figsize = (w, h), dpi = d)
    return(fig.gca())
plt.rcParams['figure.figsize'] = (6, 5)
sc.set_figure_params(dpi=100, vector_friendly=True)
sc.settings.figdir = "./figures/"

import scvi

## frequently used variables
from matplotlib import colors
import matplotlib.pyplot as plt
colorsComb = np.vstack([plt.cm.Reds(np.linspace(0, 1, 128)), plt.cm.Greys_r(np.linspace(0.7, 0.8, 0))])
mymap = colors.LinearSegmentedColormap.from_list('my_colormap', colorsComb)

## Along these Lines, a colourmap diverging from gray to red
gray_red = colors.LinearSegmentedColormap.from_list("grouping", ["lightgray", "red", "darkred"], N = 128)

## Some more Colour Maps
gray_violet = colors.LinearSegmentedColormap.from_list("grouping", ["lightgray", "mediumvioletred", "indigo"], N = 128)
gray_blue = colors.LinearSegmentedColormap.from_list("grouping", ["lightgray", "cornflowerblue", "darkblue"], N = 128)


def mysize(w, h, d):
    fig, ax = plt.subplots(figsize = (w, h), dpi = d)
    return(fig.gca())
#plt.rcParams['figure.figsize'] = (6, 5)
#sc.set_figure_params(dpi=120, vector_friendly=True)

import matplotlib.colors as colors
c_low = colors.colorConverter.to_rgba('orange', alpha = 0)
c_high = colors.colorConverter.to_rgba('red',alpha = 1)
cmap_transparent = colors.LinearSegmentedColormap.from_list('rb_cmap',[c_low, c_high], 512)

import matplotlib.colors as colors
c_low2 = colors.colorConverter.to_rgba('green', alpha = 0)
c_high2 = colors.colorConverter.to_rgba('darkblue',alpha = 1)
cmap_transparent2 = colors.LinearSegmentedColormap.from_list('rb_cmap',[c_low2, c_high2], 512)

import cell2location as c2l
from cell2location.utils import select_slide

##load
adata_vis = sc.read("test.h5ad")

Niche_NMF_palette = dict(zip(adata_vis.obs.Niche_NMF.cat.categories, adata_vis.uns["Niche_NMF_colors"]))

import hotspot

def compute_scores_hotspot(adata,
                             genes,
                             layer=None,
                             n_neighbors = 30,
                             neighborhood_factor = 3,
                             gene_symbols = None,
                             use_rep = "X_scVI",
                             model = 'danb'
                            ):
    
    
    
    if use_rep is None:
        if layer is None:
            index = pynndescent.NNDescent(adata.X, n_neighbors=n_neighbors+1)
        else:
            index = pynndescent.NNDescent(adata.layers[layer], n_neighbors=n_neighbors+1)
    else:
        index = pynndescent.NNDescent(adata.obsm[use_rep], n_neighbors=n_neighbors+1)
    
    ind, dist = index.neighbor_graph
    ind, dist = ind[:, 1:], dist[:, 1:]
    
    ind = pd.DataFrame(ind, index=list(range(adata.n_obs)))
    neighbors = ind
    
    weights = compute_weights(dist, neighborhood_factor=neighborhood_factor)
    weights = pd.DataFrame(weights, index=neighbors.index,
                           columns=neighbors.columns)
    
       
    if layer is None:
        if scipy.sparse.issparse(adata.X):
            umi_counts = adata.X.sum(axis=1).A1
        else:
            umi_counts = adata.X.sum(axis=1)

    else:
        if scipy.sparse.issparse(adata.layers[layer]):
            umi_counts = adata.layers[layer].sum(axis=1).A1
        else:
            umi_counts = adata.layers[layer].sum(axis=1)
        
    if gene_symbols is None:
        counts_dense = Hotspot._counts_from_anndata(
                    adata[:, adata.var_names.isin(genes)], layer, dense=True)
    else: 
        counts_dense = Hotspot._counts_from_anndata(
                    adata[:, adata.var[gene_symbols].isin(genes)], layer, dense=True)
        

    scores = modules.compute_scores(
        counts_dense,
        model,
        umi_counts,
        neighbors.values,
        weights.values,
    )

    return scores

import pynndescent
import scipy.sparse 

import hotspot
from hotspot import modules
from hotspot.knn import compute_weights
from hotspot.hotspot import Hotspot

senescence =["GLB1","TP53","SERPINE1","CDKN1A","CDKN1B","CDKN2B"]

scores = compute_scores_hotspot(adata_vis, senescence,layer="log1p",n_neighbors = 30,neighborhood_factor = 3,gene_symbols = None,use_rep = "X_scVI",model = 'danb')
adata_vis.obs["senescence_hs_score"] = scores

#for i in fibs:
    gene = "senescence_hs_score"
    vmins = None
    vcenters= None
    vmaxs = 2
    img_keys='lowres'
    cmaps= "bwr"
    palettes=Niche_NMF_palette
    group= None #"Fibrotic"
    import matplotlib.pyplot as plt
    from cell2location.utils import select_slide
    fig, (ax1, ax2, ax3, ax4, ax5,ax6) = plt.subplots(1, 6, figsize=(26,4), )
    plt.suptitle("                                                                          control                           ", y=1.05)
    with plt.rc_context({'axes.facecolor':  'white','figure.figsize': [4, 4]}):
        slide = select_slide(adata_vis, "90_C1_RO-730_Healthy_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax1, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "91_A1_RO-727_Healthy_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax2, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "91_B1_RO-728_Healthy_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax3, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "1217_0002_processed_aligned")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax4, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "92_A1_RO-3203_Healthy_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax5, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "1217_0004_processed_aligned")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax6, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
    plt.savefig("./figures/spatial_plot_senescence_hs_score_healthy_vmax2.pdf")   
    fig, (ax7, ax8, ax9, ax10, ax11,) = plt.subplots(1, 5, figsize=(22,4), )
    plt.suptitle("                                                                          IPF                           ", y=1.05)
    with plt.rc_context({'axes.facecolor':  'white','figure.figsize': [4, 4]}):
        slide = select_slide(adata_vis, "90_A1_H237762_IPF_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax7, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "1217_0001_processed_aligned")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax8, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "91_D1_24513-17_IPF_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax9, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "1217_0003_processed_aligned")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax10, show=False,vcenter= vcenters,palette=palettes, legend_loc=None)
        slide = select_slide(adata_vis, "92_D1_RO-3736_IPF_processed_CM")
        sc.pl.spatial(slide, cmap=cmaps,color=gene,size=1.3,img_key=img_keys,vmin=vmins, use_raw=False,layer="log1p",groups=group, vmax=vmaxs, ax=ax11, show=False,vcenter= vcenters,palette=palettes)
    plt.savefig("./figures/spatial_plot_senescence_hs_score_IPF_vmax2.pdf")