- 學(xué)習(xí)資料來源：
- scanpy主頁：https://scanpy.readthedocs.io/en/stable/
- 官網(wǎng)：https://scanpy-tutorials.readthedocs.io/en/latest/plotting/core.html【注意教程有兩個版本，這里是latest版本的學(xué)習(xí)筆記】
本教程將探索 Scanpy 的可視化可能性，并將其分為三個部分：
- Scatter plots for embeddings (eg. UMAP, t-SNE)
- Identification of clusters using known marker genes
- Visualization of differentially expressed genes
在本教程中占调，我們將使用來自10x 的數(shù)據(jù)集捆等，其中包含來自 PBMC 的68k 個細(xì)胞。
scanpy包中封裝了這個數(shù)據(jù)集的少部分?jǐn)?shù)據(jù)：700 cells and 765 highly variable genes漫玄，經(jīng)過了預(yù)處理以及UMAP降維柳骄。
本教程使用的marker如下：
- B-cell: CD79A, MS4A1
- Plasma: IGJ (JCHAIN)
- T-cell: CD3D
- NK: GNLY, NKG7
- Myeloid: CST3, LYZ
- Monocytes: FCGR3A
- Dendritic: FCER1A
01 tSNE, UMAP散點(diǎn)圖繪制
可以使用sc.pl.tsne, sc.pl.umap等函數(shù)繪制散點(diǎn)圖。這些函數(shù)訪問存儲在 adata.obms 中的數(shù)據(jù)。
- sc.pl.umap使用adata.obsm['X_umap']
import scanpy as sc
import pandas as pd
from matplotlib.pyplot import rc_context
import matplotlib.pyplot as pl
sc.set_figure_params(dpi=1000, color_map = 'viridis_r')
sc.settings.verbosity = 1
sc.logging.print_header()
## 加載數(shù)據(jù)
pbmc = sc.datasets.pbmc68k_reduced()
# inspect pbmc contents
pbmc
數(shù)據(jù)情況如下：
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'
obsp: 'distances', 'connectivities'
基因表達(dá)可視化：
繪制的參數(shù)可以是.obs中的任意列名辙浑，如基因或者其他參數(shù)狸捅，.obs是一個數(shù)據(jù)框衷蜓，每一行為一個細(xì)胞，有點(diǎn)類似Seurat數(shù)據(jù)結(jié)構(gòu)中的metadata
outdir = '/Pub/Users/zhangjuan/project/scanpy/Plot/'
# rc_context is used for the figure size, in this case 4x4
with rc_context({'figure.figsize': (4, 4)}):
sc.pl.umap(pbmc, color='CD79A')
pl.savefig(outdir + "./01-UMAP_CD79A.png")
CD79A基因表達(dá)：
可以繪制多個基因或者變量：
- ncols：控制每列繪制幾幅圖
- vmax：控制圖中最大值
with rc_context({'figure.figsize': (3, 3)}):
sc.pl.umap(pbmc, color=['CD79A', 'MS4A1', 'IGJ', 'CD3D', 'FCER1A', 'FCGR3A', 'n_counts', 'bulk_labels'], s=50, frameon=False, ncols=4, vmax='p99')
pl.savefig(outdir + "./01-UMAP_Gene.png")
結(jié)果：可以看見marker基因在特定群中特異性高表達(dá)
聚類圖修改：
# compute clusters using the leiden method and store the results with the name `clusters`
sc.tl.leiden(pbmc, key_added='clusters', resolution=0.5)
with rc_context({'figure.figsize': (5, 5)}):
sc.pl.umap(pbmc, color='clusters', add_outline=True, legend_loc='on data',
legend_fontsize=12, legend_fontoutline=2,frameon=False,
title='clustering of cells', palette='Set1')
pl.savefig(outdir + "./01-UMAP_clusters.png")
修改之后如下：比之前的好看一些
02 基于已知markers的細(xì)胞類型鑒定
通常薪贫，細(xì)胞cluster需要使用已知的標(biāo)記基因進(jìn)行標(biāo)記恍箭。利用散點(diǎn)圖，我們可以看到一個基因的表達(dá)瞧省，也許可以把它與一個cluster聯(lián)系起來扯夭。
在這里鳍贾，我們將展示其他可視化的方法，使用點(diǎn)圖交洗，小提琴圖骑科，熱圖和我們稱之為“軌跡圖tracksplot”的東西，將標(biāo)記基因關(guān)聯(lián)到cluster.
所有這些可視化總結(jié)相同的信息构拳，在不同的cluster中的表達(dá)情況咆爽，和最佳結(jié)果的選擇是留給研究者做決定。
首先置森，我們建立了一個標(biāo)記基因字典斗埂，因為這將允許 Scanpy 自動標(biāo)記基因組：
marker_genes_dict = {
'B-cell': ['CD79A', 'MS4A1'],
'Dendritic': ['FCER1A', 'CST3'],
'Monocytes': ['FCGR3A'],
'NK': ['GNLY', 'NKG7'],
'Other': ['IGLL1'],
'Plasma': ['IGJ'],
'T-cell': ['CD3D'],
}
03 dotplotk可視化
這種類型的圖總結(jié)了兩種類型的信息: 顏色表示每個類別內(nèi)的平均表達(dá)(在這種情況下是每個簇) ，點(diǎn)大小表示表達(dá)基因的類別中的細(xì)胞比例
此外凫海，向圖中添加一個樹狀圖也很有用呛凶，可以將類似的集群聚集在一起。利用聚類之間 PCA 成分的相關(guān)性自動計算層次聚類行贪。
sc.pl.dotplot(pbmc, marker_genes_dict, 'clusters', dendrogram=True)
pl.savefig(outdir + "./02-Dotplot_markers.png")
使用這個圖漾稀，我們可以看到第4組對應(yīng)于 B 細(xì)胞，第2組是 T 細(xì)胞等建瘫。此信息可用于手動注釋單元格崭捍，如下所示:
# create a dictionary to map cluster to annotation label
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['clusters'].map(cluster2annotation).astype('category')
sc.pl.dotplot(pbmc, marker_genes_dict, 'cell type', dendrogram=True)
pl.savefig(outdir + "./02-Dotplot_markers_anno.png")
之前教程01：http://www.reibang.com/p/3302c664e330 中遇到不支持多個cluster是同一種細(xì)胞類型的格式，看來這里又學(xué)習(xí)到了一種新的注釋方法啰脚！
手動注釋結(jié)果：
散點(diǎn)圖注釋后的結(jié)果：
sc.pl.umap(pbmc, color='cell type', legend_loc='on data',
frameon=False, legend_fontsize=10, legend_fontoutline=2)
pl.savefig(outdir + "./02-Dotplot_markers_anno_UMAP.png")
04 violin plot
探索這些標(biāo)記的另一種方法是用小提琴繪圖殷蛇。這里我們可以看到CD79A在集群5和8中的表達(dá)，以及MS4A1在集群5中的表達(dá)拣播。與點(diǎn)圖相比晾咪，小提琴圖給我們提供了基因表達(dá)值在細(xì)胞中的分布。
with rc_context({'figure.figsize': (4.5, 3)}):
sc.pl.violin(pbmc, ['CD79A', 'MS4A1'], groupby='clusters' )
pl.savefig(outdir + "./03-Violin_markers.png")
小提琴圖：
注意:小提琴繪圖還可以用于繪制存儲在.obs中的任何數(shù)值贮配。例如谍倦，這里用小提琴圖來比較不同集群之間的基因數(shù)量和線粒體基因的百分比
# use stripplot=False to remove the internal dots,
# inner='box' adds a boxplot inside violins
with rc_context({'figure.figsize': (4.5, 3)}):
sc.pl.violin(pbmc, ['n_genes', 'percent_mito'], groupby='clusters', stripplot=False, inner='box')
pl.savefig(outdir + "./03-Violin_n_genes-percent_mito.png")
結(jié)果如下：
05 stacked-violin plot
為了同時查看所有標(biāo)記基因的小提琴圖，我們使用sc.pl.stacked_violin泪勒。與前面一樣昼蛀，將一個樹形圖添加到類似的集群中。
ax = sc.pl.stacked_violin(pbmc, marker_genes_dict, groupby='clusters', swap_axes=False, dendrogram=True)
pl.savefig(outdir + "./04-stacked-violin.png")
結(jié)果如下：
06 matrixplot
將基因表達(dá)可視化的一個簡單方法是用矩陣圖圆存。這是按類別分組的每個基因的平均表達(dá)值的熱圖叼旋。這種類型圖顯示的信息基本上與dotplot中的顏色相同。
這里沦辙，基因的表達(dá)量歸一化為從0到1夫植，1表示最大的均值表達(dá)量，0表示最小的均值表達(dá)量
sc.pl.matrixplot(pbmc, marker_genes_dict, 'clusters', dendrogram=True, cmap='Blues', standard_scale='var', colorbar_title='column scaled\nexpression')
pl.savefig(outdir + "./05-matrixplot.png")
結(jié)果如下：
其他有用的選擇是使用sc.pp.scale歸一化基因表達(dá)。這里详民，我們將這些信息存儲在scale下延欠。然后我們調(diào)整了繪圖的最小值和最大值，并使用一個不同的顏色映射(在這種情況下沈跨，RdBu_r由捎，其中_r表示反轉(zhuǎn))。
# scale and store results in layer
pbmc.layers['scaled'] = sc.pp.scale(pbmc, copy=True).X
sc.pl.matrixplot(pbmc, marker_genes_dict, 'clusters', dendrogram=True, colorbar_title='mean z-score', layer='scaled', vmin=-2, vmax=2, cmap='RdBu_r')
pl.savefig(outdir + "./05-matrixplot-scaled.png")
結(jié)果：
07 合并圖
使用axis給繪圖以組合多個輸出饿凛，如下面的示例所示：
import matplotlib.pyplot as pl
fig, (ax1, ax2, ax3) = pl.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')
pl.savefig(outdir + "./05-plot_combined.png")
結(jié)果如下：
08 Heatmaps圖
熱圖不像以前的圖那樣歸類細(xì)胞狞玛。相反，每個細(xì)胞顯示在一行中(如果swap_axes=True則顯示在列中)涧窒⌒姆荆可以添加groupby信息，并使用與sc.pl.umap或任何其他嵌入相同的顏色代碼顯示纠吴。
ax = sc.pl.heatmap(pbmc, marker_genes_dict, groupby='clusters', cmap='viridis', dendrogram=True)
pl.savefig(outdir + "./06-Heatmaps.png")
結(jié)果圖：
熱圖也可以使用scaled數(shù)據(jù)繪制蒙畴。在下一幅圖中，類似于之前的矩陣圖呜象，最小值和最大值已經(jīng)被調(diào)整，并使用了一個不同的顏色映射
ax = sc.pl.heatmap(pbmc, marker_genes_dict, groupby='clusters', layer='scaled', vmin=-2, vmax=2, cmap='RdBu_r', dendrogram=True, swap_axes=True, figsize=(11,4))
pl.savefig(outdir + "./06-Heatmaps_scaled.png")
結(jié)果圖：
09 Tracksplot
軌跡圖顯示了與熱圖相同的信息碑隆，但是恭陡，基因表達(dá)用高度代替了顏色值
ax = sc.pl.tracksplot(pbmc, marker_genes_dict, groupby='clusters', dendrogram=True)
pl.savefig(outdir + "./07-tracksplot.png")
結(jié)果圖：
10 差異表達(dá)基因可視化
我們不像以前那樣通過已知的基因標(biāo)記來確定集群的特征，而是可以識別在集群或組中有差異表達(dá)的基因上煤。
為了識別差異表達(dá)的基因休玩，我們運(yùn)行sc.tl.rank_genes_groups。這個功能將取每組細(xì)胞劫狠，并將每一個基因在組內(nèi)的分布與不在組內(nèi)的所有其他細(xì)胞的分布進(jìn)行比較拴疤。在這里，我們將使用10倍給出的原始細(xì)胞標(biāo)記來識別這些細(xì)胞類型的標(biāo)記基因独泞。
在每個cluster中都展示差異表達(dá)基因：
sc.tl.rank_genes_groups(pbmc, groupby='clusters', method='wilcoxon')
sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=4)
pl.savefig(outdir + "./08-rank_genes_groups_dotplot.png")
氣泡圖可視化差異表達(dá)基因：每個cluster FC前4個差異基因
為了得到一個更好的表示呐矾，我們可以繪制對數(shù)log FC而不是基因表達(dá)。同時懦砂，我們想要關(guān)注在細(xì)胞類型表達(dá)和其他細(xì)胞之間具有l(wèi)og fold變化>= 3的基因蜒犯。
設(shè)置：values_to_plot='logfoldchanges' and min_logfoldchange=3
sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=4, values_to_plot='logfoldchanges', min_logfoldchange=3, vmax=7, vmin=-7, cmap='bwr')
pl.savefig(outdir + "./08-rank_genes_groups_dotplot_FC.png")
FC值可視化top4：
只畫某些類比如cluster1與clsuter5：
sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=30, values_to_plot='logfoldchanges', min_logfoldchange=4, vmax=7, vmin=-7, cmap='bwr', groups=['1', '5'])
pl.savefig(outdir + "./08-rank_genes_groups_dotplot_FC1.png")
結(jié)果：
使用matrixplot可視化差異表達(dá)基因：
sc.pl.rank_genes_groups_matrixplot(pbmc, n_genes=3, use_raw=False, vmin=-3, vmax=3, cmap='bwr', layer='scaled')
pl.savefig(outdir + "./08-rank_genes_groups_matrixplot.png")
結(jié)果圖：
使用stacked violin plots可視化差異表達(dá)基因：
sc.pl.rank_genes_groups_stacked_violin(pbmc, n_genes=3, cmap='viridis_r')
pl.savefig(outdir + "./08-rank_genes_groups_stacked_violin.png")
結(jié)果圖：
heatmap可視化差異表達(dá)基因：
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);
pl.savefig(outdir + "./08-rank_genes_groups_heatmap.png")
結(jié)果圖：
每個類別顯示10個基因，關(guān)閉基因標(biāo)簽并交換軸荞膘。請注意罚随，當(dāng)圖像交換時，類別的顏色代碼將出現(xiàn)羽资，而不是“括號”淘菩。
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')
pl.savefig(outdir + "./08-rank_genes_groups_heatmap10.png")
結(jié)果圖：
tracksplot可視化差異表達(dá)基因：
sc.pl.rank_genes_groups_tracksplot(pbmc, n_genes=3)
pl.savefig(outdir + "./08-rank_genes_groups_tracksplot.png")
結(jié)果圖：
split violin plots可視化差異表達(dá)基因：
with rc_context({'figure.figsize': (9, 1.5)}):
sc.pl.rank_genes_groups_violin(pbmc, n_genes=20, jitter=False)
結(jié)果圖：其中一個cluster
11 不同cluster之間的聚類樹
大多數(shù)可視化可以使用樹狀圖來排列類別。然而屠升，樹狀圖也可以單獨(dú)繪制如下：
# 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')
pl.savefig(outdir + "./09-dendrogram.png")
結(jié)果圖：
12 繪制相關(guān)性
與樹狀圖一起潮改，可以繪制出類別的相關(guān)性(默認(rèn)為pearson)
ax = sc.pl.correlation_matrix(pbmc, 'bulk_labels', figsize=(5,3.5))
pl.savefig(outdir + "./10-correlation_matrix.png")
結(jié)果圖：
