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

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

scanpy的空间转录组学数据的分析和可视化
空间转录组的分析思路
NBIS系列单细胞转录组数据分析实战（七）：空间转录组分析
NBIS系列单细胞转录组数据分析实战（七）：空间转录组分析
10X单细胞空间联合分析之十（RCTD）
NC：利用SpaceFlow分析组织中细胞的时空模式 | 空转工
科普讲堂|空间转录组测序知多少
Seurat 分析Visium空间转录组
空间转录组教程||用BayesSpace提高分群和基因分辨率
Giotto|| 空间表达数据分析工具箱

作者，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")

生活很好，有你更好

scanpy的空间转录组学数据的分析和可视化
本教程演示了如何使用Scanpy处理空间转录组学数据。scanpy分析空间转录组数据[https://scanpy...
空间转录组的分析思路
空间转录组数据探索空间转录组的一般分析思路如下：1.Characterize: 结合scRNA-seq数据反卷积...
NBIS系列单细胞转录组数据分析实战（七）：空间转录组分析
第七节：空间转录组分析在本节教程中，我们将学习空间转录组的一些分析方法。10X Visium平台的空间转录组数据...
NBIS系列单细胞转录组数据分析实战（七）：空间转录组分析
第七节：空间转录组分析在本节教程中，我们将学习空间转录组的一些分析方法。10X Visium平台的空间转录组数据...
10X单细胞空间联合分析之十（RCTD）
hello,上一篇文章中，10X空间转录组数据分析之空间注释（解卷积，STdeconvolve）[https://...
NC：利用SpaceFlow分析组织中细胞的时空模式 | 空转工
分析空间转录组数据集的一个主要挑战是同时合并细胞转录组相似性及其空间位置。近日《Nature Communicat...
科普讲堂|空间转录组测序知多少
一、空间转录组芯片空间转录组数据分析的核心是根据每个芯片上每个spot的基因表达信息进行聚类，然后将spot根据...
Seurat 分析Visium空间转录组
一、Seurat v3.2 对空间转录组Visium的结果分析中实现的功能：归一化降维和聚类检测空间可变特征...
空间转录组教程||用BayesSpace提高分群和基因分辨率
最近开发的空间基因表达技术，如空间转录组学可以在保留空间背景的同时，全面测量转录组谱。然而，现有的空间基因表达数据...
Giotto|| 空间表达数据分析工具箱
Seurat 新版教程：分析空间转录组数据[https://www.jianshu.com/p/f6da86489...