- 單細胞轉(zhuǎn)錄組數(shù)據(jù)分析|| scanpy教程：預處理與聚類
- 單細胞轉(zhuǎn)錄組數(shù)據(jù)分析|| scanpy教程：PAGA軌跡推斷
- 單細胞轉(zhuǎn)錄組數(shù)據(jù)分析|| scanpy教程：使用ingest和BBKNN整合多樣本
在數(shù)據(jù)分析中離不開結(jié)果的呈現(xiàn)净刮，像seurat一樣，scanpy也提供了大量的可視化的函數(shù)印蔗。在python生態(tài)中，繪圖主要由matplotlib和seaborn來完成。數(shù)據(jù)可視化是一門藝術(shù)游昼，每一種可視化的呈現(xiàn)都給我們一個解讀數(shù)據(jù)的視角仪缸，也為我們下一步數(shù)據(jù)分析提供了視覺上的依據(jù)。
今天痴荐，我們就來看看scanpy的可視化結(jié)果吧血柳。
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
在前面的學習中我們發(fā)現(xiàn)scanpy的可視化函數(shù)都是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其實是一些函數(shù)的接口生兆，所有的繪圖函數(shù)匯總难捌，預處理函數(shù)匯總，計算工具匯總鸦难，這樣便于代碼閱讀：
# 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文件的話击敌，就不能這樣寫，而應(yīng)該寫成from .A import XXX辈末，A是指當前文件夾下你想導入的函數(shù)(或者其他的)的python程序名愚争。所以有時候我們也可能看到 from .. import XXX(即上一個文件夾中的__init_.py)映皆，或者from ..A import XXX(即上一個文件夾中的文件A)。
如果細看的話轰枝，可能今天我們一張圖也看不到了捅彻。
讀入之前我們處理過的數(shù)據(jù)：
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))
強調(diào)展示某些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')
應(yīng)用在自己的差異基因中
不同的計算差異基因方法：
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()
看到這個結(jié)果诚撵，作何感想缭裆？
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')
