利用scanpy进行单细胞测序分析（三）Marker基因的可视化

其实这一部分在前面就已经涉及到一些，不过官网既然把这部分拿出来单独作为一大块讲解，可能也是因为这一部分可供选择的可视化方法有很多。对于图片的优化上也有比较详细的介绍。

官网这部分讲解的地址：https://scanpy-tutorials.readthedocs.io/en/latest/visualizing-marker-genes.html

开始之前需要注意一点的是，这部分官网分别使用了marker 基因和高变基因，来进行可视化。所以看上去好像来来回回都是那么几张图，但实际上前半部分展示的是特定细胞群的marker基因的可视化。而后半部分展示的是每一个cluster里的高变基因的可视化。

#调用软件，设置参数
>>> import scanpy as sc
>>> import pandas as pd
>>> from matplotlib import rcParams

>>> sc.set_figure_params(dpi=80, color_map='viridis')
>>> sc.settings.verbosity = 2
>>> sc.logging.print_versions()
scanpy==1.4.6 anndata==0.7.1 umap==0.4.1 numpy==1.18.2 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.2.post1 statsmodels==0.11.1 python-igraph==0.8.0 louvain==0.6.1

这一部分练习的数据来自10x PBMC 68k数据库。这个dataset已经被过滤了，包含700个细胞和765个高变基因。另外，这个dataset也已经经过PCA和UMAP计算。louvain方法聚类和细胞周期检测的结果都存放在了pbmc.obs里:

#你不需要下载，直接导入数据
>>> pbmc = sc.datasets.pbmc68k_reduced()
>>> pbmc
AnnData object with n_obs × n_vars = 700 × 765 
    obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain'
    var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
    uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'

可以看到这个AnnData对象已经是一个预处理相当全面的dataset了，这一章也是主要介绍经过前期处理后，有哪些方法可以让你的data很漂亮的展示出来。
可以看一下UMAP的图：

#用rcParams修改图片的默认大小
>>> rcParams['figure.figsize'] = 4, 4
>>> sc.pl.umap(pbmc, color=['bulk_labels'], s=50)

给定一个marker基因的list:

>>> marker_genes = ['CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ',  'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1','FCGR3A', 'FCER1A', 'CST3']

小提琴图

画在每一个cluster里这些marker基因的表达情况：

>>> ax = sc.pl.stacked_violin(pbmc, marker_genes, groupby='bulk_labels',
                         var_group_positions=[(7, 8)],var_group_labels=['NK'])
#可以看一下上面代码里bulk_labels里是什么：
>>> print(pbmc.obs['bulk_labels'])
index
AAAGCCTGGCTAAC-1     CD14+ Monocyte
AAATTCGATGCACA-1          Dendritic
AACACGTGGTCTTT-1           CD56+ NK
AAGTGCACGTGCTA-1    CD4+/CD25 T Reg
ACACGAACGGAGTG-1          Dendritic
                         ...       
TGGCACCTCCAACA-8          Dendritic
TGTGAGTGCTTTAC-8          Dendritic
TGTTACTGGCGATT-8    CD4+/CD25 T Reg
TTCAGTACCGGGAA-8            CD19+ B
TTGAGGTGGAGAGC-8          Dendritic
Name: bulk_labels, Length: 700, dtype: category
Categories (10, object): [CD4+/CD25 T Reg, CD4+/CD45RA+/CD25- Naive T, CD4+/CD45RO+ Memory, CD8+ Cytotoxic T, ..., CD19+ B, CD34+, CD56+ NK, Dendritic]

你还可以将x轴和y轴交换，这样就是每一个marker基因在所有cluster里的表达情况：

>>> ax = sc.pl.stacked_violin(pbmc, marker_genes, groupby='bulk_labels', swap_axes=True, var_group_positions=[(7, 8)], var_group_labels=['NK'], dendrogram=True)

点图

只有当你的count矩阵里有０的时候才有意义，点图不适用于那些已经被矫正过不存在０的矩阵。

Marker基因可以是列表，也可以是字典。如果存为字典，你画出的图将会被分组并标记:

>>> marker_genes_dict = {'B-cell': ['CD79A', 'MS4A1'],
                     'T-cell': 'CD3D',
                     'T-cell CD8+': ['CD8A', 'CD8B'],
                     'NK': ['GNLY', 'NKG7'],
                     'Myeloid': ['CST3', 'LYZ'],
                     'Monocytes': ['FCGR3A'],
                     'Dendritic': ['FCER1A']}
# use marker genes as dict to group them
>>> ax = sc.pl.dotplot(pbmc, marker_genes_dict, groupby='bulk_labels')

为了让上面的图更加能清楚的表达我们data的内容，可以加点东西进去：
（１）dendrogram=True：在图里加入树状图，也就是分支图。
（２）dot_max=0.5 最大的点设置为0.5，基因表达量大于50%的都用这个型号的点表示
（３）dot_min=0.3 最小的点设置为0.3，基因表达量小于30%的都用这个型号
（４）standard_scale=’var’ 对矩阵进行scale矫正（０－１）

>>> ax = sc.pl.dotplot(pbmc, marker_genes, groupby='bulk_labels', dendrogram=True, dot_max=0.5, dot_min=0.3, standard_scale='var')

上图还可以换颜色，改变图片大小：

>>> ax = sc.pl.dotplot(pbmc, marker_genes_dict, groupby='bulk_labels', dendrogram=True,
                   standard_scale='var', smallest_dot=40, color_map='Blues', figsize=(8,5))

给基因加上group标签：

>>> ax = sc.pl.dotplot(pbmc, marker_genes, groupby='louvain',
              var_group_positions=[(0,1), (11, 12)],
              var_group_labels=['B cells', 'dendritic'],
              figsize=(12,4), var_group_rotation=0, dendrogram='dendrogram_louvain')

矩阵图

矩阵图以热图的形式显示在不同cluster中基因的平均表达。与点图相比，矩阵图可以用于校正的矩阵。默认情况下使用原始count。

>>> gs = sc.pl.matrixplot(pbmc, marker_genes_dict, groupby='bulk_labels')

接下来我们给图里加入系统树，并用standar_scale=var把矩阵进行scale到０－１：

>>> gs = sc.pl.matrixplot(pbmc, marker_genes_dict, groupby='bulk_labels', dendrogram=True, standard_scale='var')

如果我们用use_raw=False画图呢？使用pbmc里存储的已经经过sc.pp.scale处理后的矩阵来画图，并将坐标轴进行交换：

# to use the 'non-raw' data we select marker genes present in this data.
>>> marker_genes_2 = [x for x in marker_genes if x in pbmc.var_names]
>>> gs = sc.pl.matrixplot(pbmc, marker_genes_2, groupby='bulk_labels', dendrogram=True,
                      use_raw=False, vmin=-3, vmax=3, cmap='bwr',  swap_axes=True, figsize=(5,6))

热图

热图是将每一个细胞作为一行或者一列，groupby信息就是cluster的信息

>>> ax = sc.pl.heatmap(pbmc,marker_genes_dict, groupby='louvain')

你也可以将use_raw=False画图，代码如下，图就不放上来了：

>>> ax = sc.pl.heatmap(pbmc, marker_genes, groupby='louvain',
              var_group_positions=[(0,1), (11, 12)], use_raw=False, vmin=-3, vmax=3, cmap='bwr',
              var_group_labels=['B cells', 'dendritic'], var_group_rotation=0, dendrogram='dendrogram_louvain')

Tracksplots

Trackplot和热图是一样的，只不过基因的表达不是用标尺来表示，而是用高度来表示：

# Track plot 用来表示非log的count值会更好
>>> import numpy as np
>>> ad = pbmc.copy()
>>> ad.raw.X.data = np.exp(ad.raw.X.data)
>>> ax = sc.pl.tracksplot(ad,marker_genes, groupby='louvain')

过滤Marker基因

利用sc.tl.filter_rank_genes_groups工具，我们可以根据一些条件来选择性的可视化marker基因，比如说，在一个cluster里，选择那些变化倍数（fold change）至少是相对于其他cluster里的２倍以上，并且这个cluster里的marker基因变化２倍以上的细胞数至少要50%以上。

建立一个新的category（名为broad_type），把所有的T细胞放进去：

>>> t_cell = ['CD4+/CD25 T Reg', 'CD4+/CD45RA+/CD25- Naive T', 'CD4+/CD45RO+ Memory','CD8+ Cytotoxic T', 'CD8+/CD45RA+ Naive Cytotoxic']
pbmc.obs['broad_type'] = pd.Categorical(pbmc.obs.bulk_labels.apply(lambda x: x if x not in t_cell  else 'T-cell'))

在broad_type里寻找新的top基因：

>>> sc.tl.rank_genes_groups(pbmc, 'broad_type', method='wilcoxon', n_genes=20)
#这里官网写错了！！！他们把20个基因写成了200个！还是那句话：不要轻易相信任何人的代码，除非自己跑一遍！！！！
>>> sc.tl.filter_rank_genes_groups(pbmc)
>>> rcParams['figure.figsize'] = 4,4
>>> rcParams['axes.grid'] = True
>>> sc.pl.rank_genes_groups(pbmc, key='rank_genes_groups_filtered', ncols=3)

用过滤完marker基因画热图：

>>> sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=200, key='rank_genes_groups_filtered',
                                swap_axes=True, use_raw=False, vmax=3, vmin=-3, cmap='bwr', dendrogram=True)

可以设置更为严格的marker基因过滤标准：

>>> sc.tl.filter_rank_genes_groups(pbmc,
                               min_in_group_fraction=0.5,
                               max_out_group_fraction=0.3,
                               min_fold_change=1.5)

然后比较一下过滤前和过滤后的点图：

#过滤前
>>> sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=6)

过滤前

#过滤后
>>> sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=6, key='rank_genes_groups_filtered')

过滤后

过滤后的看看矩阵图和小提琴图：

>>> sc.pl.rank_genes_groups_matrixplot(pbmc, n_genes=6, key='rank_genes_groups_filtered', standard_scale='var', cmap='Blues')

# Track plot data and violin is better visualized using the non-log counts
>>> ad = pbmc.copy()
>>> ad.raw.X.data = np.exp(ad.raw.X.data)
>>> sc.pl.rank_genes_groups_stacked_violin(ad, n_genes=6, key='rank_genes_groups_filtered', standard_scale='var')

>>> sc.pl.rank_genes_groups_tracksplot(ad, n_genes=6, key='rank_genes_groups_filtered')

用panels可视化高变基因

其实这个方法在这一系列笔记的第一篇里就提到了。不知道为什么这里又讲了一遍。我认为区别在于上面的图，我们都是根据给定的基因列表（marker gene）来做的，下面的图是根据每一个cluster里的高变基因做的。也就是说每一个cluster里都有自己的基因列表。

>>> import warnings
>>> warnings.simplefilter(action='ignore', category=FutureWarning)
>>> sc.tl.rank_genes_groups(pbmc, groupby='bulk_labels', method='logreg')
ranking genes
    finished (0:00:00)
#visualize marker genes using panels
>>> rcParams['figure.figsize'] = 4,4
>>> rcParams['axes.grid'] = True
>>> sc.pl.rank_genes_groups(pbmc)

把上面的图转换成点图，每个cluster里显示４个高变基因:

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

上面是显示了全部的cluster，现在可以让它只显示其中的两个，每个cluster显示15个基因：

>>> axs = sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=15, groups=['Dendritic', 'CD19+ B'])

用矩阵点图可视化高变基因

>>> axs = sc.pl.rank_genes_groups_matrixplot(pbmc, n_genes=3, standard_scale='var', cmap='Blues')

换颜色：

>>> axs = sc.pl.rank_genes_groups_matrixplot(pbmc, n_genes=3, use_raw=False, vmin=-3, vmax=3, cmap='bwr')

image.png

用小提琴图展示高变基因

# instead of pbmc we use the 'ad' object (created earlier) in which the raw matrix is exp(pbmc.raw.matrix). This highlights better the differences between the markers.
>>> sc.pl.rank_genes_groups_stacked_violin(ad, n_genes=3)

把所有的“小提琴”设置成同样的颜色：

>>> sc.pl.rank_genes_groups_stacked_violin(ad, n_genes=3, row_palette='slateblue')

坐标轴交换：

>>> sc.pl.rank_genes_groups_stacked_violin(ad, n_genes=3, swap_axes=True, figsize=(6, 10), width=0.4)

用热图可视化高变基因

#展示前３个高变基因
>>> sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=3, standard_scale='var')

同样的，你可以交换坐标轴，更改颜色，图不放了，放代码：

>>> #坐标轴交换
sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=3, use_raw=False, swap_axes=True, vmin=-3, vmax=3, cmap='bwr')

下面可以显示每个cluster里的top10高变基因，去掉基因标签，将图旋转90度：

>>> sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=10, use_raw=False, swap_axes=True, show_gene_labels=False,
                                vmin=-3, vmax=3, cmap='bwr')

你也许注意到了这个热图有点奇怪，感觉两边缺两条，实际不是缺，我感觉可能是其他原因。而且下面的cluster也没有和热图对其。但是并不妨碍结果的展示，因为如果你运行出来的图也有这样的问题，不用担心，你可以将图片保存为svg格式的图片，然后放到illustrator里进行编辑（illustrator是什么就不用说了），把下面的bar的长度调整到正确长度，两边的黑线往里面拉一拉就好了。

tracksplot可视化高变基因

>>> sc.pl.rank_genes_groups_tracksplot(ad, n_genes=3)

Plot相关性

>>> sc.tl.dendrogram(pbmc, 'bulk_labels', n_pcs=30)
>>> ax = sc.pl.correlation_matrix(pbmc, 'bulk_labels')