我們前面也分享過很多篇關(guān)于Seurat的使用，尤其是在單細胞轉(zhuǎn)錄組分析上，今天我們來分享和學(xué)習(xí)另外一個和它齊名的工具：scanpy。也許有人會說單細胞分析使用Seurat，monocle等R包會更加方便。但是實際分析中,測試情況是當細胞量大于5萬時。一般小型服務(wù)器內(nèi)存很容易不足，這時候請不要過多嘗試使用Seurat來進行分析，monocle2更是。而基于python的單細胞轉(zhuǎn)錄分析包scanpy，能很好的解決內(nèi)存不足的問題，網(wǎng)上經(jīng)驗說在整合80萬細胞量時，整個預(yù)處理流程在6小時左右能夠完成。

scanpy相關(guān)python 包安裝（安裝好python3之后，終端運行）。參考官方文檔：https://github.com/scverse/scanpy

======安裝相關(guān)的庫========

conda create -n scanpy? ? #我自己喜歡重新創(chuàng)建一個環(huán)境

source activate scanpy

conda install -c conda-forge scanpy python-igraph leidenalg

======測試數(shù)據(jù)======

這兒我們還是使用了Seurat官網(wǎng)使用的pbmc3k數(shù)據(jù)進行測試（我是很早之前就下載過了）：

wget http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -O data/pbmc3k_filtered_gene_bc_matrices.tar.gz

tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz

=======程序測試======

運行python3，導(dǎo)入相關(guān)包及設(shè)置一些必要路徑

import scanpy as sc

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

sc.settings.verbosity = 3?# verbosity即冗余 。設(shè)置日志等級

sc.logging.print_versions() # 輸出版本號

sc.settings.set_figure_params(dpi=80)

import os

os.getcwd()

os.chdir("./filtered_gene_bc_matrices/scanpy")

os.getcwd()

results_file = 'pbmc3k.h5ad'

讀取數(shù)據(jù)：

data=sc.read_10x_mtx('./filtered_gene_bc_matrices/hg19/',var_names='gene_symbols', cache=True)

data.var_names_make_unique() # 索引去重，若上一步中使用 `var_names='gene_ids'` 則這一步非必須進行

注意：當然scanpy可以直接讀取10Xgenomics的.h5格式數(shù)據(jù)

data=sc.read_10x_h5("./pbmc3K.h5",genome=None,gex_only=True)

data.var_names_make_unique()

>>> print(data)

AnnData object with n_obs × n_vars = 2700 × 32738

? ? var: 'gene_ids'

#可以看出data對象是，cell數(shù)2700 基因數(shù)32738的矩陣，而data的數(shù)據(jù)存儲架構(gòu)如下圖所示：

data.X 存儲count matrix，數(shù)據(jù)類型為稀疏矩陣scipy.sparse.csr.csr_matrix

data.obs 存儲關(guān)于 obervations(cells) 的 metadata，數(shù)據(jù)類型為 pandas dataframe

data.var 存儲關(guān)于 variables(genes) 的 metadata，數(shù)據(jù)類型為? pandas dataframe

#AnnData.uns 存儲后續(xù)附加的其他非結(jié)構(gòu)化信息

#data.obs_names 和 data.var_names是 index

#細胞名和基因名可分別通過 data.obs_names 和 data.var_names 查看。 AnnData 對象可以像 dataframe 一樣進行切片操作，例如，data_subset = data[:, list_of_gene_names]

數(shù)據(jù)預(yù)處理：

讀取數(shù)據(jù)之后，像Seurat一樣，我們需要對數(shù)據(jù)進行預(yù)處理，比如測到的轉(zhuǎn)錄本總數(shù)(total_counts)、測到的基因總數(shù)(total_cells)、來源于線粒體基因的轉(zhuǎn)錄本所占比例等。

基礎(chǔ)過濾

去除表達基因200以下的細胞；去除在3個細胞以下表達的基因。這也是通常Seurat通常用的默認過濾標準。

sc.pp.filter_cells：進行細胞的過濾，該函數(shù)保留至少有 min_genes 個基因（某個基因表達非0可判斷存在該基因）的細胞，或者保留至多有 max_genes 個基因的細胞；

sc.pp.filter_genes：進行基因的過濾，該函數(shù)用于保留在至少 min_cells 個細胞中出現(xiàn)的基因，或者保留在至多 max_cells 個細胞中出現(xiàn)的基因；

sc.pp.filter_cells(data,min_genes=200)

sc.pp.filter_genes(data,min_cells=3)

data

質(zhì)量過濾：

其實怎么過濾，取決于對自己數(shù)據(jù)的了解程度。

比如：可視化所有細胞中計數(shù)最多的基因。

sc.pl.highest_expr_genes(data,n_top=20)

plt.savefig("Highest_expr_genes.pdf")?

#我們可以看到前20個基因中有14個基因是屬于核糖體基因，說明核糖體基因在這此數(shù)據(jù)集中表達豐度很高，有時候我們會去除核糖體占比過高的細胞，認為這些細胞質(zhì)量低，這兒我們先不考慮核糖體基因，考慮線粒體相關(guān)基因。

例如，計算各種類型相關(guān)基因比例：線粒體相關(guān)基因，血紅蛋白相關(guān)基因，核糖體相關(guān)基因，葉綠體相關(guān)基因等。

#mitochondrial genes

data.var['mt'] = data.var_names.str.startswith('MT-')?

#hemoglobin genes.? 血紅蛋白基因

data.var['hb'] = data.var_names.str.contains('^HB[^P]')

#ribosomal genes

data.var['ribo'] = data.var_names.str.startswith('RPS','RPL')

data

sc.pp.calculate_qc_metrics(data, qc_vars=['mt','hb','ribo'], percent_top=None, log1p=False, inplace=True)

其中：

qc_vars：用于標識要控制的特征（基因），布爾型元素，用于作為mask使用；

percent_top：計算與常出現(xiàn)基因的比，percent_top=[50] 計算與第 50 個最常出現(xiàn)基因的比例，None則不計算；

inplace：決定是否將計算指標添加到var和obs中；

log1p：設(shè)置為False可以跳過轉(zhuǎn)換到log1p空間的過程；log1p即log(1+number)，用于壓縮數(shù)據(jù)并確保結(jié)果是一個正數(shù)；

data

這樣注釋中就得到很多QC計算得到的信息。

使用violinplot度量查看各類質(zhì)量，類似seurat：

n_genes_by_counts：每個細胞中，有表達的基因的個數(shù)；

total_counts：每個細胞的基因總計數(shù)（總表達量,umi數(shù)）；

pct_counts_mt：每個細胞中，線粒體基因表達量占該細胞所有基因表達量的百分比

pct_counts_hb:每個細胞中，血紅蛋白基因表達量占該細胞所有基因表達量的百分比

pct_counts_ribo:每個細胞中，核糖體RNA基因表達量占該細胞所有基因表達量的百分比

植物的話，還需要查看pct_counts_ct：每個細胞中，葉綠體基因表達量占該細胞所有基因表達量的百分比

sc.pl.violin(data, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt','pct_counts_ribo','pct_counts_hb'],jitter=0.4, multi_panel=True,save="QC_violin.pdf")

由于線粒體和 MALAT1 基因的表達水平被認為主要是技術(shù)性的，因此除了計算每個細胞中線粒體基因，血紅蛋白基因和核糖體基因所占的比例外，明智的做法是還要將線粒體基因，血紅蛋白基因,MALAT1基因從數(shù)據(jù)集中直接刪除，然后再進行任何進一步的分析。

mito_genes = data.var_names.str.startswith('MT-')

malat1 = data.var_names.str.startswith('MALAT1')

hb_genes = data.var_names.str.contains('^HB[^(P)]')

remove = np.add(mito_genes, malat1)

remove = np.add(remove, hb_genes)

keep = np.invert(remove)

data = data[:,keep]

也可以根據(jù)基因總數(shù)進行進一步的過濾。

data = data[data.obs['n_genes_by_counts'] < 2500, :]

數(shù)據(jù)標準化：

數(shù)據(jù)預(yù)處理結(jié)束之后，就是和seurat一樣進行標準化，挑選HVG基因了。

sc.pp.normalize_total(data, target_sum=1e4) ##標準化

sc.pp.log1p(data)? #將數(shù)據(jù)壓縮到log1p的空間，就是取了個對數(shù)

下面就是類似seurat的挑選高變異的基因HVG。

data.raw = data

sc.pp.highly_variable_genes(data, min_mean=0.0125, max_mean=3, min_disp=0.5)

sc.pl.highly_variable_genes(data)

plt.savefig("highly_variable_genes.pdf")

可以看到，執(zhí)行完之后，就多出了uns:log1p和hvg的操作。

然后就是歸一化，將數(shù)據(jù)放在單位方差。

sc.pp.regress_out(data, ['total_counts', 'pct_counts_mt']) #校正細胞基因計數(shù)和線粒體基因比例的影響。

sc.pp.scale(data, max_value=10)

data.write(results_file)

需要注意的是，在做完歸一化后，data.X的數(shù)據(jù)格式從scipy.sparse.csr_matrix轉(zhuǎn)換為numpy.ndarray。

下面，我們就開始做經(jīng)典的PCA聚類，降維和UMAP分析等。

sc.tl.pca(data, svd_solver='arpack')# svd_solver 指定奇異值分解 SVD 的方法

#sc.tl.pca(data,color="CT3")

sc.pl.pca_variance_ratio(data, log=True)

plt.savefig("PCA.pdf")

data.write(results_file)

降維（對neighborhood graph進行embedding）

sc.pp.neighbors(data, n_neighbors=10, n_pcs=25)

聚類：

sc.tl.leiden(data)? ?#使用leiden進行聚類，注意安裝對應(yīng)的python包，conda install -c conda-forge leidenalg ，當然也可使用其他的聚類算法，如louvain，sc.tl.louvain(data),sc.pl.umap(adata, color=['louvain'])，比較了一下，聚類結(jié)果差異不大

#sc.tl.leiden(data,resolution=0.3)

sc.tl.umap(data)

sc.pl.umap(data, color=['leiden'])

plt.savefig("Umap.pdf")

我們先查看一下常見的marker的表達情況。

sc.settings.set_figure_params(dpi=50, dpi_save=600, figsize=(5,5))

marker_names = ['IL7R','CCR7',

? ? ? ? ? ? ? ? 'CD14','LYZ',

? ? ? ? ? ? ? 'IL7R','S100A4',

? ? ? ? ? ? ? 'MS4A1',

? ? ? ? ? ? ? 'CD8A',

? ? ? ? ? ? ? 'FCGR3A','MS4A7',

? ? ? ? ? ? ? 'GNLY','NKG7',

? ? ? ? ? ? ? 'FCER1A','CST3',

? ? ? ? ? ? ? 'PPBP']

sc.pl.umap(data,color = marker_names, ncols=2)

plt.savefig("marker.pdf")

也可以通過對比不同cluster之間的表達差異，通過差異表達基因?qū)ふ襇arker基因。

讓我們計算每個cluster中高度差異基因的排名，最簡單和最快的方法是t-test，當然還有wilcoxon。

sc.tl.rank_genes_groups(data, 'leiden', method='wilcoxon')

sc.pl.rank_genes_groups(data, n_genes=30, sharey=False)

plt.savefig("dif_gene.pdf")

pd.DataFrame(data.uns['rank_genes_groups']['names']).head(5)

data.uns['rank_genes_groups'].keys()

#dict_keys(['logfoldchanges', 'names', 'params', 'pvals', 'pvals_adj', 'scores'])

result = data.uns['rank_genes_groups']

groups = result['names'].dtype.names

pd.DataFrame(? ? {group + '_' + key[:1]: result[key][group]? ? for group in groups for key in ['names', 'pvals']}).head(5)

res = pd.DataFrame( {group + '_' + key: result[key][group] for group in groups for key in ['names', 'pvals','logfoldchanges','pvals_adj','scores']})

res.to_csv("dif.csv") #基因差異情況輸出到本地保存

當然任意兩個簇之間也可以比較差異:

sc.tl.rank_genes_groups(data, 'leiden', groups=['0'], reference='1', method='wilcoxon') #cluster0 與 cluster1之間進行比較

sc.pl.rank_genes_groups(data, groups=['0'], n_genes=20)

plt.savefig("gene_exp_C0_C1_score.pdf")

sc.pl.rank_genes_groups_violin(data, groups='0', n_genes=8)

plt.savefig("gene_exp_C0_C1_vio.pdf")

sc.pl.rank_genes_groups_violin(data, groups='0', n_genes=8)

plt.savefig("gene_exp_C0_all_vio.pdf")

也可以像seurat一樣繪制任意基因的小提琴圖和UMAP圖。

sc.pl.violin(data, ['MS4A1', 'NKG7', 'PPBP'], groupby='leiden')

plt.savefig("marker_exp_Vio.pdf")

sc.pl.umap(data, color=['MS4A1', 'NKG7', 'PPBP'])

plt.savefig("marker_exp_umap.pdf")

從前面marker基因的分布來看，我們可以清晰的看到B細胞相關(guān)的MS4A1在簇4中表達，因此我們可以定于這簇為B細胞，同樣的CD8A在簇5表達最高，簇5為CD8 T；CD14在cluster0中表達，所以0位CD14+Mono；FCER1A在cluster 8中表達，cluster 8位DC；MS4A7則在cluster6中表達，所以cluster6為FCGR3A Monocytes；GNLY和NKG7在cluster7中表達，所以cluster 7位NK。

new_cluster_names = ['CD14 Monocytes','Memory CD41','Memory CD42','Naive CD4 T','B','CD8 T','FCGR3A Monocytes','NK','DC']

data.rename_categories('leiden', new_cluster_names)

sc.pl.umap(data, color='leiden', legend_loc='on data', title='', frameon=True)

plt.savefig("umap_celltype.pdf")

ax = sc.pl.dotplot(data, marker_names, groupby='leiden')

plt.savefig("gene_Dot_celltype.V2.pdf")

ax = sc.pl.stacked_violin(data, marker_names, groupby='leiden', rotation=90)

plt.savefig("gene_Vio_celltype.V2.pdf")