scanpy单细胞分析官网示例全解析（三）全网最详细及细节

本文介绍scanpy的空间转录组分析功能

1.导包及参数设置

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

sc.logging.print_versions()
sc.set_figure_params(facecolor="white", figsize=(8, 8))
sc.settings.verbosity = 3

2.下载数据并读取

代码自动下载保存读取

adata = sc.datasets.visium_sge(sample_id="V1_Human_Lymph_Node")
adata.var_names_make_unique()
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

3.质控和预处理

QC and preprocessing

官网的图字都叠一起了，这是谁写的教程，扣他鸡腿！所以使用subplots_adjust调节子图之间的距离

subplot功能是matplotlib包的，这里设置1行，4列，共4幅小图，图的大小为15、4

用seaborn包画一些简单的图，对count等数据统计

fig, axs = plt.subplots(1, 4, figsize=(15, 4))
fig.subplots_adjust(hspace=0.5, wspace=0.5)
sns.histplot(adata.obs["total_counts"], kde=False, ax=axs[0])
sns.histplot(adata.obs["total_counts"][adata.obs["total_counts"] < 10000], kde=False, bins=40, ax=axs[1])
sns.histplot(adata.obs["n_genes_by_counts"], kde=False, bins=60, ax=axs[2])
sns.histplot(adata.obs["n_genes_by_counts"][adata.obs["n_genes_by_counts"] < 4000], kde=False, bins=60, ax=axs[3])

4.过滤数据

细胞的count最少5000，最多35000；线粒体基因占比大于20%删除，基因至少要在10个细胞表达否则删除

sc.pp.filter_cells(adata, min_counts=5000)
sc.pp.filter_cells(adata, max_counts=35000)
adata = adata[adata.obs["pct_counts_mt"] < 20]
print(f"#cells after MT filter: {adata.n_obs}")
sc.pp.filter_genes(adata, min_cells=10)

5.高变基因

先标准化再对数化

sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)

6.聚类图

计算pca，neighbors，umap，再用leiden聚类

sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="clusters")
#控制图的大小
plt.rcParams["figure.figsize"] = (4, 4)
#点代表着细胞，用不同的数据上色
sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)

7.空间坐标可视化

Visualization in spatial coordinates

图像数据存储在adata.uns['spatial']中，也就是参数img_key="hires"，还有一些其他信息，比如'software_version': 'spaceranger-1.1.0'，应该是spaceranger记录染色之后的图像数据

plt.rcParams["figure.figsize"] = (8, 8)
sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts","clusters"],size=1.5)
#局部的图，group参数控制只显示哪些簇的点，crop_coord控制坐标
sc.pl.spatial(adata, img_key="hires", color="clusters", groups=["0", "3"], crop_coord=[7000, 10000, 0, 6000], alpha=0.5, size=1.3)

可能与官网的配色不一致，不过不需要担心，只是簇上成不一样的颜色罢了

8.maker基因

t-test计算 rank gene

sc.tl.rank_genes_groups(adata, "clusters", method="t-test")
#展示group3的前1个rank gene，
sc.pl.rank_genes_groups_heatmap(adata, groups="3", n_genes=10, groupby="clusters")
#在染色图片上对某些基因可视化
sc.pl.spatial(adata, img_key="hires", color=["clusters", "CR2","COL1A2", "SYPL1"])

官网教程没有他说的CR2，我找到了！

这是基因在不同细胞的表达

9.fish数据

通过荧光原位杂交得到的转录组数据

读取数据

官网下载失败可以用百度网盘：链接：https://pan.baidu.com/s/1ycj2HdGZk3we08al4f44Cw 
提取码：xysm

#coordinates是细胞的坐标数据，counts是细胞counts数
coordinates = pd.read_excel("./data/pnas.1912459116.sd15.xlsx", index_col=0)
counts = sc.read_csv("./data/pnas.1912459116.sd12.csv").transpose()
#选取数据，根据coordinates的index，也就是细胞名，获取counts数据
adata_merfish = counts[coordinates.index, :]
adata_merfish.obsm["spatial"] = coordinates.to_numpy()

聚类

sc.pp.normalize_per_cell(adata_merfish, counts_per_cell_after=1e6)
sc.pp.log1p(adata_merfish)
sc.pp.pca(adata_merfish, n_comps=15)
sc.pp.neighbors(adata_merfish)
sc.tl.umap(adata_merfish)
sc.tl.leiden(adata_merfish, key_added="clusters", resolution=0.5)

sc.pl.umap(adata_merfish, color="clusters")

这是leiden算法聚类结果

basis就是你想使用什么作图，可以换为pca，umap，spatial

看来umap在单细胞数据中能有效将数据分类！可能使用了什么高级的classification机器学习算法吧！pca不行了是吗？

sc.pl.embedding(adata_merfish, basis="pca", color="clusters")