课后补充----空间转录组hotspot分析

作者，Evil Genius

今日参考文献

其中提到了一个分析hotspot，其中文章的描述是

We therefore calculated a custom senescence gene score for every spot in the spatial data using the hotspot gene-module scoring, visualized as dot plot across the niches, and on one tissue section

那么分析方法呢？

其中提供了分析方法，软件hotspot，函数compute_scores。将每个基因的计数拟合到 DANB 模型(深度调整的负二项式)中，然后对每个基因的值进行归一化(中心化)，然后根据每个细胞的邻域对其进行平滑处理。最后，使用主成分分析对所有基因特征中的这些值进行降维。文章使用log1p layers、30 个neighbors和 X_scVI 嵌入（降维后的坐标）来计算得分。

这其中基因集哪里来的？文章采用的NMF方法分析得到的空间微环境（niche）特征基因。

来看看具体的方法

基础加载

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 seaborn as sns
#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)

print(f"squidpy=={sq.__version__}")
print(f"scanpy=={sc.__version__}")

import cell2location as c2l
from cell2location.utils import select_slide

import hotspot
import pynndescent
import scipy.sparse 

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

定义hotspot计算函数，我们分析空间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 

数据分析

##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"]))

###基因列表要根据之前的分析定义
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")  

sns.reset_defaults()
sc.pl.dotplot(adata_vis, var_names=["senescence_hs_score"], groupby="Niche_NMF", layer="log1p", vcenter=0, cmap = "bwr",
              return_fig=True, save="dotplot_senescence_hs_score_bwr.pdf").add_totals().legend(width=2).show()

生活很好，有你更好