其實(shí)這一部分在前面就已經(jīng)涉及到一些，不過官網(wǎng)既然把這部分拿出來單獨(dú)作為一大塊講解，可能也是因?yàn)檫@一部分可供選擇的可視化方法有很多。對(duì)于圖片的優(yōu)化上也有比較詳細(xì)的介紹。
官網(wǎng)這部分講解的地址：https://scanpy-tutorials.readthedocs.io/en/latest/visualizing-marker-genes.html
開始之前需要注意一點(diǎn)的是，這部分官網(wǎng)分別使用了marker 基因和高變基因，來進(jìn)行可視化。所以看上去好像來來回回都是那么幾張圖，但實(shí)際上前半部分展示的是特定細(xì)胞群的marker基因的可視化。而后半部分展示的是每一個(gè)cluster里的高變基因的可視化。
#調(diào)用軟件，設(shè)置參數(shù)
>>> 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
這一部分練習(xí)的數(shù)據(jù)來自10x PBMC 68k數(shù)據(jù)庫(kù)。這個(gè)dataset已經(jīng)被過濾了，包含700個(gè)細(xì)胞和765個(gè)高變基因。另外，這個(gè)dataset也已經(jīng)經(jīng)過PCA和UMAP計(jì)算。louvain方法聚類和細(xì)胞周期檢測(cè)的結(jié)果都存放在了pbmc.obs里:
#你不需要下載，直接導(dǎo)入數(shù)據(jù)
>>> 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'
可以看到這個(gè)AnnData對(duì)象已經(jīng)是一個(gè)預(yù)處理相當(dāng)全面的dataset了，這一章也是主要介紹經(jīng)過前期處理后，有哪些方法可以讓你的data很漂亮的展示出來。
可以看一下UMAP的圖：
#用rcParams修改圖片的默認(rèn)大小
>>> rcParams['figure.figsize'] = 4, 4
>>> sc.pl.umap(pbmc, color=['bulk_labels'], s=50)
給定一個(gè)marker基因的list:
>>> marker_genes = ['CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1','FCGR3A', 'FCER1A', 'CST3']
小提琴圖
畫在每一個(gè)cluster里這些marker基因的表達(dá)情況：
>>> 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軸交換，這樣就是每一個(gè)marker基因在所有cluster里的表達(dá)情況：
>>> 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)
點(diǎn)圖
只有當(dāng)你的count矩陣?yán)镉校暗臅r(shí)候才有意義，點(diǎn)圖不適用于那些已經(jīng)被矯正過不存在０的矩陣。
Marker基因可以是列表，也可以是字典。如果存為字典，你畫出的圖將會(huì)被分組并標(biāo)記:
>>> 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')
為了讓上面的圖更加能清楚的表達(dá)我們data的內(nèi)容，可以加點(diǎn)東西進(jìn)去：
（１）dendrogram=True：在圖里加入樹狀圖，也就是分支圖。
（２）dot_max=0.5 最大的點(diǎn)設(shè)置為0.5，基因表達(dá)量大于50%的都用這個(gè)型號(hào)的點(diǎn)表示
（３）dot_min=0.3 最小的點(diǎn)設(shè)置為0.3，基因表達(dá)量小于30%的都用這個(gè)型號(hào)
（４）standard_scale=’var’ 對(duì)矩陣進(jìn)行scale矯正（０－１）
>>> ax = sc.pl.dotplot(pbmc, marker_genes, groupby='bulk_labels', dendrogram=True, dot_max=0.5, dot_min=0.3, standard_scale='var')
上圖還可以換顏色，改變圖片大?。?/p>
>>> 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標(biāo)簽：
>>> 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中基因的平均表達(dá)。與點(diǎn)圖相比，矩陣圖可以用于校正的矩陣。默認(rèn)情況下使用原始count。
>>> gs = sc.pl.matrixplot(pbmc, marker_genes_dict, groupby='bulk_labels')
接下來我們給圖里加入系統(tǒng)樹，并用standar_scale=var把矩陣進(jìn)行scale到０－１：
>>> gs = sc.pl.matrixplot(pbmc, marker_genes_dict, groupby='bulk_labels', dendrogram=True, standard_scale='var')
如果我們用use_raw=False畫圖呢？使用pbmc里存儲(chǔ)的已經(jīng)經(jīng)過sc.pp.scale處理后的矩陣來畫圖，并將坐標(biāo)軸進(jìn)行交換：
# 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))
熱圖
熱圖是將每一個(gè)細(xì)胞作為一行或者一列，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和熱圖是一樣的，只不過基因的表達(dá)不是用標(biāo)尺來表示，而是用高度來表示：
# Track plot 用來表示非log的count值會(huì)更好
>>> 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工具，我們可以根據(jù)一些條件來選擇性的可視化marker基因，比如說，在一個(gè)cluster里，選擇那些變化倍數(shù)（fold change）至少是相對(duì)于其他cluster里的２倍以上，并且這個(gè)cluster里的marker基因變化２倍以上的細(xì)胞數(shù)至少要50%以上。
建立一個(gè)新的category（名為broad_type），把所有的T細(xì)胞放進(jìn)去：
>>> 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)
#這里官網(wǎng)寫錯(cuò)了?。?！他們把20個(gè)基因?qū)懗闪?00個(gè)！還是那句話：不要輕易相信任何人的代碼，除非自己跑一遍?。。?！
>>> 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)
可以設(shè)置更為嚴(yán)格的marker基因過濾標(biāo)準(zhǔn)：
>>> sc.tl.filter_rank_genes_groups(pbmc,
min_in_group_fraction=0.5,
max_out_group_fraction=0.3,
min_fold_change=1.5)
然后比較一下過濾前和過濾后的點(diǎn)圖：
#過濾前
>>> 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可視化高變基因
其實(shí)這個(gè)方法在這一系列筆記的第一篇里就提到了。不知道為什么這里又講了一遍。我認(rèn)為區(qū)別在于上面的圖，我們都是根據(jù)給定的基因列表（marker gene）來做的，下面的圖是根據(jù)每一個(gè)cluster里的高變基因做的。也就是說每一個(gè)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)
把上面的圖轉(zhuǎn)換成點(diǎn)圖，每個(gè)cluster里顯示４個(gè)高變基因:
>>> sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=4)
上面是顯示了全部的cluster，現(xiàn)在可以讓它只顯示其中的兩個(gè)，每個(gè)cluster顯示15個(gè)基因：
>>> axs = sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=15, groups=['Dendritic', 'CD19+ B'])
用矩陣點(diǎn)圖可視化高變基因
>>> 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')
用小提琴圖展示高變基因
# 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)
把所有的“小提琴”設(shè)置成同樣的顏色：
>>> sc.pl.rank_genes_groups_stacked_violin(ad, n_genes=3, row_palette='slateblue')
坐標(biāo)軸交換：
>>> sc.pl.rank_genes_groups_stacked_violin(ad, n_genes=3, swap_axes=True, figsize=(6, 10), width=0.4)
用熱圖可視化高變基因
#展示前３個(gè)高變基因
>>> sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=3, standard_scale='var')
同樣的，你可以交換坐標(biāo)軸，更改顏色，圖不放了，放代碼：
>>> #坐標(biāo)軸交換
sc.pl.rank_genes_groups_heatmap(pbmc, n_genes=3, use_raw=False, swap_axes=True, vmin=-3, vmax=3, cmap='bwr')
下面可以顯示每個(gè)cluster里的top10高變基因，去掉基因標(biāo)簽，將圖旋轉(zhuǎn)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')
你也許注意到了這個(gè)熱圖有點(diǎn)奇怪，感覺兩邊缺兩條，實(shí)際不是缺，我感覺可能是其他原因。而且下面的cluster也沒有和熱圖對(duì)其。但是并不妨礙結(jié)果的展示，因?yàn)槿绻氵\(yùn)行出來的圖也有這樣的問題，不用擔(dān)心，你可以將圖片保存為svg格式的圖片，然后放到illustrator里進(jìn)行編輯（illustrator是什么就不用說了），把下面的bar的長(zhǎng)度調(diào)整到正確長(zhǎng)度，兩邊的黑線往里面拉一拉就好了。
tracksplot可視化高變基因
>>> sc.pl.rank_genes_groups_tracksplot(ad, n_genes=3)
Plot相關(guān)性
>>> sc.tl.dendrogram(pbmc, 'bulk_labels', n_pcs=30)
>>> ax = sc.pl.correlation_matrix(pbmc, 'bulk_labels')
