- 单细胞转录组数据分析|| scanpy教程：预处理与聚类
- 单细胞转录组数据分析|| scanpy教程：PAGA轨迹推断
- 单细胞转录组数据分析|| scanpy教程：使用ingest和BBKNN整合多样本
在数据分析中离不开结果的呈现，像seurat一样，scanpy也提供了大量的可视化的函数。在python生态中，绘图主要由matplotlib和seaborn来完成。数据可视化是一门艺术，每一种可视化的呈现都给我们一个解读数据的视角，也为我们下一步数据分析提供了视觉上的依据。
今天，我们就来看看scanpy的可视化结果吧。
import scanpy as sc
import pandas as pd
from matplotlib import rcParams
import numpy as np
sc.set_figure_params(dpi=80, color_map='viridis')
sc.settings.verbosity = 2
sc.logging.print_versions()
scanpy==1.4.5.1
anndata==0.7.1
umap==0.3.10
numpy==1.18.2
scipy==1.3.1
pandas==0.25.1
scikit-learn==0.21.3
statsmodels==0.10.1
python-igraph==0.8.0
在前面的学习中我们发现scanpy的可视化函数都是pl来指引的如：sc.pl.violin，我们就来看看这个pl是如何写的：
from ._anndata import scatter, violin, ranking, clustermap, stacked_violin, heatmap, dotplot, matrixplot, tracksplot, dendrogram, correlation_matrix
from ._preprocessing import filter_genes_dispersion, highly_variable_genes
from ._tools.scatterplots import embedding, pca, diffmap, draw_graph, tsne, umap
from ._tools import pca_loadings, pca_scatter, pca_overview, pca_variance_ratio
from ._tools.paga import paga, paga_adjacency, paga_compare, paga_path
from ._tools import dpt_timeseries, dpt_groups_pseudotime
from ._tools import rank_genes_groups, rank_genes_groups_violin
from ._tools import rank_genes_groups_dotplot, rank_genes_groups_heatmap, rank_genes_groups_stacked_violin, rank_genes_groups_matrixplot, rank_genes_groups_tracksplot
from ._tools import sim
from ._tools import embedding_density
from ._rcmod import set_rcParams_scanpy, set_rcParams_defaults
from . import palettes
from ._utils import matrix
from ._utils import timeseries, timeseries_subplot, timeseries_as_heatmap
from ._qc import highest_expr_genes
我们看到pl其实是一些函数的接口，所有的绘图函数汇总，预处理函数汇总，计算工具汇总，这样便于代码阅读：
# some technical stuff
import sys
from ._utils import pkg_version, check_versions, annotate_doc_types
__author__ = ', '.join([
'Alex Wolf',
'Philipp Angerer',
'Fidel Ramirez',
'Isaac Virshup',
'Sergei Rybakov',
'Gokcen Eraslan',
'Tom White',
'Malte Luecken',
'Davide Cittaro',
'Tobias Callies',
'Marius Lange',
'Andrés R. Muñoz-Rojas',
])
__email__ = ', '.join([
'f.alex.wolf@gmx.de',
'philipp.angerer@helmholtz-muenchen.de',
# We don’t need all, the main authors are sufficient.
])
try:
from setuptools_scm import get_version
__version__ = get_version(root='..', relative_to=__file__)
del get_version
except (LookupError, ImportError):
__version__ = str(pkg_version(__name__))
check_versions()
del pkg_version, check_versions
# the actual API
from ._settings import settings, Verbosity # start with settings as several tools are using it
from . import tools as tl
from . import preprocessing as pp
from . import plotting as pl
from . import datasets, logging, queries, external, get
from anndata import AnnData
from anndata import read_h5ad, read_csv, read_excel, read_hdf, read_loom, read_mtx, read_text, read_umi_tools
from .readwrite import read, read_10x_h5, read_10x_mtx, write
from .neighbors import Neighbors
set_figure_params = settings.set_figure_params
# has to be done at the end, after everything has been imported
annotate_doc_types(sys.modules[__name__], 'scanpy')
del sys, annotate_doc_types
我们看到一种新的语法： from . import XXX，这是什么意思呢？
在当前程序所在文件夹里__init_-.py程序中导入XXX，如果当前程序所在文件夹里没有_init_.py文件的话，就不能这样写，而应该写成from .A import XXX，A是指当前文件夹下你想导入的函数(或者其他的)的python程序名。所以有时候我们也可能看到 from .. import XXX(即上一个文件夹中的__init_.py)，或者from ..A import XXX(即上一个文件夹中的文件A)。
如果细看的话，可能今天我们一张图也看不到了。
读入之前我们处理过的数据：
results_file = 'E:\learnscanpy\write\pbmc3k.h5ad'
sdata = sc.read_h5ad(results_file)
sdata
AnnData object with n_obs × n_vars = 2223 × 2208
obs: 'n_genes', 'percent_mito', 'n_counts', 'percent_mito1', 'leiden'
var: 'gene_ids', 'feature_types', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
uns: 'leiden', 'leiden_colors', 'leiden_sizes', 'neighbors', 'paga', 'pca', 'rank_genes_groups', 'umap'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
堆积小提琴图
核心的可视化基本都是针对marker 基因的，所以我们先读进来一套：
marker_genes = ['CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
'FCGR3A', 'FCER1A', 'CST3']
ax = sc.pl.stacked_violin(sdata, marker_genes, groupby='leiden',
var_group_positions=[(7, 8)], var_group_labels=['6'])
可以看出所给基因在每个分组中的表达分布，也可以对分群聚类：
ax = sc.pl.stacked_violin(sdata, marker_genes, groupby='leiden', swap_axes=True, dendrogram=True)
点图
给定marker基因列表，可以看出他们在指定分组中的表达占比情况。
marker_genes_dict = {'1': ['CD79A', 'MS4A1'],
'2': ['PDXK'],
'3': ['CD8A', 'CD8B'],
'4': ['GNLY', 'NKG7'],
'5': ['CST3', 'LYZ'],
'6': ['FCGR3A'],
'7': ['FCER1A']}
ax = sc.pl.dotplot(sdata, marker_genes_dict, groupby='leiden')
ax = sc.pl.dotplot(sdata, marker_genes, groupby='leiden', dendrogram=True, dot_max=0.5, dot_min=0.3, standard_scale='var')
自定义颜色主题：
ax = sc.pl.dotplot(sdata, marker_genes_dict, groupby='leiden', dendrogram=True,
standard_scale='var', smallest_dot=40, color_map='Blues', figsize=(8,5))
强调展示某些gene列表：
ax = sc.pl.dotplot(sdata, marker_genes, groupby='leiden',
var_group_positions=[(0,1), (11, 12)],
var_group_labels=['1', '2'],
figsize=(12,4), var_group_rotation=0, dendrogram='dendrogram_leiden')
热图
gs = sc.pl.matrixplot(sdata, marker_genes_dict, groupby='leiden')
gs = sc.pl.matrixplot(sdata, marker_genes_dict, groupby='leiden', dendrogram=True, standard_scale='var')
marker_genes_2 = [x for x in marker_genes if x in sdata.var_names]
marker_genes_2
Out[30]:
['CD79A',
'MS4A1',
'CD8A',
'CD8B',
'LYZ',
'LGALS3',
'S100A8',
'GNLY',
'NKG7',
'KLRB1',
'FCGR3A',
'FCER1A',
'CST3']
gs = sc.pl.matrixplot(sdata, marker_genes_2, groupby='leiden', dendrogram=True,
use_raw=False, vmin=-3, vmax=3, cmap='bwr', swap_axes=True, figsize=(5,6))
ax = sc.pl.heatmap(sdata,marker_genes_dict, groupby='leiden')
ax = sc.pl.heatmap(sdata, marker_genes, groupby='leiden', figsize=(5, 8),
var_group_positions=[(0,1), (11, 12)], use_raw=False, vmin=-3, vmax=3, cmap='bwr',
var_group_labels=['1', '2'], var_group_rotation=0, dendrogram='dendrogram_leiden')
tracksplot
import numpy as np
ad = sdata.copy()
ad.X.data = np.exp(ad.X.data)
ax = sc.pl.tracksplot(ad,marker_genes, groupby='leiden')
ax = sc.pl.tracksplot(sdata,marker_genes, groupby='leiden')
应用在自己的差异基因中
不同的计算差异基因方法：
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
sc.tl.rank_genes_groups(sdata, groupby='leiden', method='logreg',key_added = "logreg")
sc.tl.rank_genes_groups(sdata, groupby='leiden', method='wilcoxon',key_added ='wilcoxon')
sc.tl.rank_genes_groups(sdata, groupby='leiden', method='t-test',key_added = 't-test')
venn图比较不同方法的异同：
from matplotlib_venn import venn3
wc = sdata.uns['wilcoxon']['names']['NK']
tt = sdata.uns['t-test']['names']['NK']
lr = sdata.uns['logreg']['names']['NK']
venn3([set(wc),set(tt),set(lr)], ('Wilcox','T-test','logreg') )
plt.show()
看到这个结果，作何感想？
rcParams['figure.figsize'] = 4,4
rcParams['axes.grid'] = True
sc.pl.rank_genes_groups(sdata)
sc.pl.rank_genes_groups_dotplot(sdata, n_genes=4)
axs = sc.pl.rank_genes_groups_dotplot(sdata, groupby='leiden', n_genes=10, dendrogram='dendrogram_leiden')
axs = sc.pl.rank_genes_groups_matrixplot(sdata, n_genes=10, standard_scale='var', cmap='Blues')
axs = sc.pl.rank_genes_groups_matrixplot(sdata, n_genes=10, use_raw=False, vmin=-3, vmax=3, cmap='bwr')
sc.pl.rank_genes_groups_stacked_violin(sdata, n_genes=3)
sc.pl.rank_genes_groups_stacked_violin(sdata, 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(sdata, n_genes=3, standard_scale='var')
sc.pl.rank_genes_groups_heatmap(sdata, n_genes=10, use_raw=False, swap_axes=True, show_gene_labels=False,
vmin=-3, vmax=3, cmap='bwr')
sc.pl.rank_genes_groups_tracksplot(ad, n_genes=20)
ax = sc.pl.dendrogram(sdata, 'leiden')
sc.tl.dendrogram(sdata, 'leiden', var_names=marker_genes, use_raw=False)
ax = sc.pl.dendrogram(sdata, 'leiden', orientation='left')
ax = sc.pl.correlation_matrix(sdata, 'leiden')
