美文网首页
空间转录组数据分析之空间邻域特征评分

空间转录组数据分析之空间邻域特征评分

作者: 单细胞空间交响乐 | 来源:发表于2024-05-19 11:32 被阅读0次

    作者,Evil Genius

    昨天,2024年5月19日,发生了两件大事。

    一是在昨天lol在成都举办的2024 MSI决赛中,Gen以3:1获得了冠军,而就在前不久,Gen俱乐部公开表示台湾是一个国家,遭到了国内粉丝的抵制。至于比赛内容,就想问一句,为什么执着于偷家?

    二是昨天相亲回来路过迎泽大桥,亲眼目睹了一个姑娘一跃而下,跳进了汾河,等到救护车赶到的时候已经过了很长时间,估计又一个鲜活的生命里离开了人世,而这个跳桥事件,已经是我听到太原今年发生的第四起跳桥事件,之前在长风大桥,晋阳湖等均有发生。

    与之鲜明对比的是,太原刚被评为全国最具幸福感的城市。

    最后,来分享一下邻域富集

    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(f".h5ad")
    
    Niche_NMF_palette = dict(zip(adata_vis.obs.Niche_NMF.cat.categories, adata_vis.uns["Niche_NMF_colors"]))
    
    sc.pl.umap(adata_vis, color="Niche_NMF", palette=Niche_NMF_palette,)
    
    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
    sc.pl.umap(adata_vis, color = "senescence_hs_score", cmap = "bwr", ax = mysize(6, 5, 120), size = 30, vcenter=0,)# save="homeo_hs_score.pdf")
    
    
    #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")  
    

    生活很好,有你更好

    相关文章

      网友评论

          本文标题:空间转录组数据分析之空间邻域特征评分

          本文链接:https://www.haomeiwen.com/subject/yaggfjtx.html