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

本文介绍scanpy的一些Core plotting functions，核心画图功能。

使用的数据为软件自带的数据

1.导包及参数设置

import scanpy as sc
import pandas as pd
from matplotlib.pyplot import rc_context

sc.set_figure_params(dpi=100, color_map = 'viridis_r')
sc.settings.verbosity = 1
sc.logging.print_header()

2.读取数据

自带数据，只有700个细胞，765个高变基因，已经有很多数据计算过了，如pca，umap等

最重要的是有细胞类型信息

pbmc = sc.datasets.pbmc68k_reduced()

3.散点图，scatter plot

①

rc_context控制图的大小

umap算法对数据降维之后，上色可以用很多种数据。一个点代表一个细胞，该细胞有什么数据都可以用来上色，比如基因的表达量，count数，细胞类型等等

with rc_context({'figure.figsize': (4, 4)}):
    sc.pl.umap(pbmc, color=['CD79A', 'MS4A1', 'IGJ', 'CD3D', 'FCER1A', 'FCGR3A', 'n_counts', 'bulk_labels'],frameon=True)

②

给每个簇上一个颜色

首先使用leiden或者louvain算法计算cluster，然后保存下来聚类结果分别为clusters1和clusters2

作图时color='clusters'1或2即可

sc.tl.leiden(pbmc, key_added='clusters1', resolution=0.5)
sc.tl.louvain(pbmc, key_added='clusters2', resolution=0.5)
with rc_context({'figure.figsize': (5, 5)}):
    sc.pl.umap(pbmc, color='clusters1', add_outline=True, legend_loc='on data',
               legend_fontsize=12, legend_fontoutline=2,frameon=False,
               title='leiden clustering of cells', palette='Set1')
with rc_context({'figure.figsize': (5, 5)}):
    sc.pl.umap(pbmc, color='clusters2', add_outline=True, legend_loc='on data',
               legend_fontsize=12, legend_fontoutline=2,frameon=False,
               title='louvain clustering of cells', palette='Set1')

这里有一个疑问，umap和leiden及louvain是什么关系？leiden或者louvain聚类之后，也是用sc.pl.umap函数画图，而该函数是画umap结果的；其次leiden和louvain算法对点的位置没有影响，区别仅在于某些点划分到哪个簇，所以是umap算法对数据进行降维，并作图到二维坐标上，然后聚类算法如leiden及louvain决定哪些点划分到一个簇？嗯！粗浅地这么判断，有大佬知道请留言！

4.点图，dotplot

①

Identification of clusters based on known marker genes

基于已知的maker gene识别簇

#不同细胞的maker基因
marker_genes_dict = {
    'B-cell': ['CD79A', 'MS4A1'],
    'Dendritic': ['FCER1A', 'CST3'],
    'Monocytes': ['FCGR3A'],
    'NK': ['GNLY', 'NKG7'],
    'Other': ['IGLL1'],
    'Plasma': ['IGJ'],
    'T-cell': ['CD3D'],
}
#点的大小代表在该cluster中有多少细胞表达了某基因，点的颜色代表在该cluster中某基因的平均表达量
#dendrogram参数表示将图聚类
sc.pl.dotplot(pbmc, marker_genes_dict, groupby='clusters1', dendrogram=True)

点的大小代表某个簇中，有多少的细胞表达了该基因

点的颜色代表某个簇中，该基因的平均表达量

②

由上图可以推断出各个cluster分别是什么细胞，将序号0-8改成相对应的细胞名字、

在obs属性的clusters1变量是由leiden算法计算的聚类结果，用0-8数字表示，将数字表示的细胞类型在obs属性中新生成一个‘cell type’变量，如下图

cluster2annotation = {
     '0': 'Monocytes',
     '1': 'Dendritic',
     '2': 'T-cell',
     '3': 'NK',
     '4': 'B-cell',
     '5': 'Dendritic',
     '6': 'Plasma',
     '7': 'Other',
     '8': 'Dendritic',
}
# add a new `.obs` column called `cell type` by mapping clusters to annotation using pandas `map` function
pbmc.obs['cell type'] = pbmc.obs['clusters1'].map(cluster2annotation).astype('category')
#将数字换做对应的细胞类型
sc.pl.dotplot(pbmc, marker_genes_dict, groupby='cell type', dendrogram=True)

同样在聚类散点图中，也可以将数字换做对应的细胞类型

sc.pl.umap(pbmc, color='cell type', legend_loc='on data',
           frameon=False, legend_fontsize=10, legend_fontoutline=2)

5.violin plot，小提琴图

①

小提琴图展示不同cluster中基因表达量

with rc_context({'figure.figsize': (4.5, 3)}):
    sc.pl.violin(pbmc, ['CD79A', 'MS4A1'], groupby='clusters1' )

小提琴图可以画obs属性中任何数值类型的数据，如画出各个cluster表达的基因数量，线粒体基因百分比

②

堆叠小提琴图stacked-violin plot

sc.pl.stacked_violin(pbmc, marker_genes_dict, groupby='clusters1', swap_axes=False, dendrogram=True)

③

matrixplot，矩阵图

#表达量标准化为0-1
sc.pl.matrixplot(pbmc, marker_genes_dict, 'clusters1', dendrogram=True, cmap='Blues', standard_scale='var', colorbar_title='column scaled\nexpression')

还可以标准化gene的表达量，自定义设置范围，这里设置的是-2到2，并且新生成一个数据层layer,包含了标准化后的的数据

利用标准化后的数据作图，cmap是表示颜色的字符串，这里使用use a diverging(发散的) color map (in this case RdBu_r where _r means reversed).

pbmc.layers['scaled'] = sc.pp.scale(pbmc, copy=True).X
sc.pl.matrixplot(pbmc, marker_genes_dict, 'clusters1', dendrogram=True,
                 colorbar_title='mean z-score', layer='scaled', vmin=-2, vmax=2, cmap='RdBu_r')

6.Combining plots in subplots，合并若干图

import matplotlib.pyplot as plt
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20,4), gridspec_kw={'wspace':0.9})
ax1_dict = sc.pl.dotplot(pbmc, marker_genes_dict, groupby='bulk_labels', ax=ax1, show=False)
ax2_dict = sc.pl.stacked_violin(pbmc, marker_genes_dict, groupby='bulk_labels', ax=ax2, show=False)
ax3_dict = sc.pl.matrixplot(pbmc, marker_genes_dict, groupby='bulk_labels', ax=ax3, show=False, cmap='viridis')

7.heatmap热图，tracksmap轨迹图

sc.pl.heatmap(pbmc, marker_genes_dict, groupby='clusters1', cmap='viridis', dendrogram=True)
sc.pl.tracksplot(pbmc, marker_genes_dict, groupby='clusters1', dendrogram=True)

8.针对maker gene可视化

①

画出每个簇前四个高变基因

sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=4)

②

还可以将表达量换成foldchange

但是按照官网的代码会报错，报错信息#报错信息Please run `sc.tl.rank_genes_groups` with 'n_genes=adata.shape[1]' to save all gene scores. Currently, only 100 are found

应该是还没有按照groups计算rank gene，就是说在每个簇中基因的顺序还能没有计算，要先计算

sc.tl.rank_genes_groups默认调用method='t-test'，计算后跑出来的图与官网不一致，将其改成wilcoxon后，与官网一致

#默认调用method='t-test'，计算后跑出来的图与官网不一致，将其改成wilcoxon后，与官网一致
sc.tl.rank_genes_groups(pbmc,groupby='clusters1',method='wilcoxon')
##上图是表达量，下图是foldchange，数据是与剩下(rest)相比较计算得到的foldchange，fc能更好显示数据间的差异
#报错信息Please run `sc.tl.rank_genes_groups` with 'n_genes=adata.shape[1]' to save all gene scores. Currently, only 100 are found
sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=4, values_to_plot='logfoldchanges', min_logfoldchange=3, vmax=7, vmin=-7, cmap='bwr')

③

也可以指定某些groups画出来，画出1，5两个簇前30个基因，并且最小fc为4

sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=30, groups=['1', '5'], values_to_plot='logfoldchanges', min_logfoldchange=4, vmax=7, vmin=-7, cmap='bwr')

发现没有那么多基因超过4

这里fold change是跟剩余的簇相比计算的，如1簇跟剩余的簇的合集相比计算fc值

④

其它图可视化maker gene

#Visualize marker genes using matrixplot
sc.pl.rank_genes_groups_matrixplot(pbmc, n_genes=3, use_raw=False, vmin=-3, vmax=3, cmap='bwr', layer='scaled')
#Visualize marker genes using stacked violin plots
sc.pl.rank_genes_groups_stacked_violin(pbmc, n_genes=3, cmap='viridis_r')
#Visualize marker genes using heatmap
sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=3, use_raw=False, swap_axes=True, vmin=-3, vmax=3, cmap='bwr', layer='scaled', figsize=(10,7), show=False);
#Visualize marker genes using tracksplot
sc.pl.rank_genes_groups_tracksplot(pbmc, n_genes=3)

⑤

Comparison of marker genes using split violin plots，分隔小提琴图比较maker gene

这婀娜的图让我想起篮球中的小提琴变向

每个簇与剩余簇比较，这里只展示0簇

with rc_context({'figure.figsize': (9, 1.5)}):
    sc.pl.rank_genes_groups_violin(pbmc, n_genes=20, jitter=False)

9.dendrogram功能

上面很多图都有一个参数dendrogram=True，这代表将图聚类生成一个树状的分支

# compute hierarchical clustering using PCs (several distance metrics and linkage methods are available).
sc.tl.dendrogram(pbmc, 'bulk_labels')
ax = sc.pl.dendrogram(pbmc, 'bulk_labels')

10.相关性correlation

不同细胞之间的相关性

sc.pl.correlation_matrix(pbmc, 'bulk_labels', figsize=(5,3.5))

obs属性里‘cell type’相关性

sc.pl.correlation_matrix(pbmc, 'cell type', figsize=(5,3.5))

obs属性里‘clusters1’相关性

sc.pl.correlation_matrix(pbmc, 'clusters1', figsize=(5,3.5))