python单细胞分析之Scanpy单细胞

单细胞转录组数据分析|| scanpy教程:可视化套件

2020-04-01  本文已影响0人  周运来就是我

在数据分析中离不开结果的呈现,像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')
上一篇 下一篇

猜你喜欢

热点阅读