空间转录组细胞niche、细胞频率、通路niche、空间互作、空间hotspot全系列脚本更新（python版本）

作者，Evil Genius

清明了，最近一直在相亲，每周都会见一个女人，注意，是女人，不是女孩子，因为大家都30多岁了，或者快30岁了，法定16岁就是妇女了，不得不说真的被相亲恶心到了，感觉相亲就是贴标签把自己往外卖，女人要房要车是常态，房车还是入场卷，女人们会说自己慢热、看眼缘，然后一般就没有然后了，最恶心的是一般男的联系几次女人不理，停上几天不联系，过几天回头就跟别人说这男的根本不上心，要求李老师（山西太原介绍相亲的）再介绍，而且女人们心机很深啊，看不上对方的条件绝对不明说，就说没眼缘，没感觉，问了一些身边相亲的男的，都差不多一样的经历，女人明明自己不愿意，但是把自己包装成弱势方，避免负责，不得不说，最毒妇人心啊😄。。

工作上也是遇到各种奇葩事，以前看过一个笑话，说人花了好几千检查身体，啥病没有，自己感觉太亏了，当初真当笑话听了，结果现在真的遇到了，一个病人做NGS报告显示都是良性的突变位点，身体状态也显示良好，结果把公司给投诉了，说我花了16000多做检查结果没有想要的结果，自己应该是能拿到那种显示有害突变位点需要药物治疗的结果，我听了真的是笑不出来。

还有以前上课，拔苗助长，直到前几天也觉得世界上哪有这么傻的人，结果是领导们都在拔苗助长，完全忽略其中的问题，也没有充分的测试，“长高了”就行了，真的是不知道世界怎么了。

空间已经有了很大的发展了，4种精度的平台都有了。

• 亚细胞级：10X Visium HD、BGI、百迈客的百创S1000
• 单细胞级：寻因、Akoya CODEX、Nanostring CosMx、10X Xenium
• 接近单细胞级：slide seq。
• 多细胞级：10X Visium、德运康瑞

但是大家要明白，都是空间平台，实力亦有差距：

高精度空间平台

这一篇我们来更新全系列的空间分析脚本，其中我们实现4个目标，全部更新python版本，R版本都有了。

1、细胞Niche，这个niche相对于之前的我们需要量化

R版本在空转第10课共定位内容补充（通路 && 细胞类型）

10X空间转录组数据分析之细胞niche

10X空间转录组数据分析之细胞的空间依赖性

2、细胞频率

3、通路niche

4、空间临近互作（直接空间角度）

5、空间hotspot

我们先来第一个，细胞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 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

加载数据，注意要做了单细胞空间联合

##load
adata_vis = sc.read(f"/.h5ad")
#####细胞类型，根据自己的情况调整
cts = ['AT0', 'AT1', 'AT2', 'Aberrant basaloid', 'Adventitial fibroblast', 'Alveolar fibroblast',
       'Artery', 'B/Plasma', 'Basal', 'Basophil/Mast', 'Bronchial Vessel', 'Capillary', 'Capillary Aerocyte',
       'Ciliated', 'Ciliated SFTPB+/SCGB1A1+', 'Dendritic', 'Ionocyte', 'Lymphatic', 'Macrophage C1Q hi',
       'Macrophage CHI3L1+/CD9 hi/', 'Macrophage FABP4+', 'Macrophage IL1B+', 'Macrophage LYVE1+', 
       'Macrophage RETN+/VCAN+', 'Mesothelial', 'Monocyte', 'Mucous', 'Myofibroblast', 'NK', 
       'Peribronchial fibroblast', 'Pericyte', 'Smooth Muscle', 'T cell', 'TB-SC', 'Vein', 'pDC', 'preTB-SC/RAS']

计算细胞niche

tab = sc.get.obs_df(adata_vis, keys= ["sampleID","Niche_NMF","treatment"]+ cts)

ct_label = "cell_type"
group_by = "sampleID"
xlabel = "Niche_NMF"
ctss = cts
#cm = "sampleID"

tab = tab.loc[:, ctss + [xlabel]].copy()
tab.head()

spot_id	AT0	AT1	AT2	...	NK	Peribronchial fibroblast	Niche_NMF
AAACAAGTATCTCCCA-1	0.001162	0.012275	0.000589	0.002108	0.203605 1.176494	0.027367	Fibroblast
AAACACCAATAACTGC-1	0.284728	0.797208	0.466162	0.047545	0.064358	0.030167	Alveolar

## Calculate the mean per cell type and compartment
test_tab = test_tab2.groupby([xlabel]).mean().reset_index()
print(test_tab.shape)
test_tab.head()

Niche_NMF	AT0	AT1	AT2	Aberrant basaloid	Adventitial	fibroblast	Alveolar fibroblast	Artery	B/Plasma	Basal	...	Myofibroblast	NK	Peribronchial fibroblast	Pericyte	Smooth Muscle	T cell	TB-SC	Vein	pDC	preTB-SC/RAS
Airway	0.009889	0.022928	0.012908	0.130147	0.033568	0.079144	0.028104	0.151954	0.603255	...	0.135885	0.047271	0.100981	0.074290	0.156967	0.094300	0.107136	0.031213	0.035241	0.343569	0.00004
Alveolar	0.249349	0.674018	0.292087	0.063480	0.055637	0.279309	0.157352	0.085331	0.024479	...	0.070915	0.133856	0.049070	0.318479	0.115077	0.189918	0.050761	0.083650	0.041931	0.021899	0.00004

test_tab.index = test_tab["Niche_NMF"].tolist()
test_tab = test_tab.drop("Niche_NMF", axis=1)
test_tab2 = test_tab.div(test_tab.sum(axis=0), axis=1)
test_tab3 = test_tab2.drop(['Monocyte', 'NK','Basophil/Mast','Ionocyte'], axis=1)
cts2 = ['AT0', 'AT1', 'AT2','Capillary', 'Capillary Aerocyte', #'Monocyte', 'NK',
        'Artery','Smooth Muscle', 'Adventitial fibroblast','Mesothelial','Lymphatic', 
        'Vein','Macrophage LYVE1+', 'Peribronchial fibroblast','Alveolar fibroblast',  'Pericyte',  #'Basophil/Mast',
        'Dendritic','pDC','Bronchial Vessel',  'B/Plasma','T cell',  
          'Basal', 'Ciliated', 'Ciliated SFTPB+/SCGB1A1+','Mucous',   'preTB-SC/RAS','TB-SC',
        'Aberrant basaloid','Myofibroblast','Macrophage IL1B+',
        'Macrophage FABP4+', 'Macrophage RETN+/VCAN+', #'Ionocyte',
        'Macrophage C1Q hi','Macrophage CHI3L1+/CD9 hi/',]
test_tab3 = test_tab3[cts2]
test_tab3 = test_tab3.iloc[[1,7,2,4,0,3,5,6]]
test_tab4 = test_tab3*100
sb.reset_defaults()
plt.figure(figsize=(15, 2.25))
sb.heatmap(test_tab4, annot=True, cmap="Reds", fmt=".0f",linewidths=0.25, linecolor='grey', )
plt.savefig("./figures/NMF_niches_celltype_percentages.pdf")