单细胞转录组分析---Scanpy（1）

我们前面也分享过很多篇关于Seurat的使用，尤其是在单细胞转录组分析上，今天我们来分享和学习另外一个和它齐名的工具：scanpy。也许有人会说单细胞分析使用Seurat，monocle等R包会更加方便。但是实际分析中,测试情况是当细胞量大于5万时。一般小型服务器内存很容易不足，这时候请不要过多尝试使用Seurat来进行分析，monocle2更是。而基于python的单细胞转录分析包scanpy，能很好的解决内存不足的问题，网上经验说在整合80万细胞量时，整个预处理流程在6小时左右能够完成。

scanpy相关python 包安装（安装好python3之后，终端运行）。参考官方文档：https://github.com/scverse/scanpy

======安装相关的库========

conda create -n scanpy    #我自己喜欢重新创建一个环境

source activate scanpy

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

======测试数据======

这儿我们还是使用了Seurat官网使用的pbmc3k数据进行测试（我是很早之前就下载过了）：

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，导入相关包及设置一些必要路径

import scanpy as sc

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

sc.settings.verbosity = 3 # verbosity即冗余 。设置日志等级

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'

读取数据：

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格式数据

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数2700 基因数32738的矩阵，而data的数据存储架构如下图所示：

data.X 存储count matrix，数据类型为稀疏矩阵scipy.sparse.csr.csr_matrix

data.obs 存储关于 obervations(cells) 的 metadata，数据类型为 pandas dataframe

data.var 存储关于 variables(genes) 的 metadata，数据类型为  pandas dataframe

#AnnData.uns 存储后续附加的其他非结构化信息

#data.obs_names 和 data.var_names是 index

#细胞名和基因名可分别通过 data.obs_names 和 data.var_names 查看。 AnnData 对象可以像 dataframe 一样进行切片操作，例如，data_subset = data[:, list_of_gene_names]

数据预处理：

读取数据之后，像Seurat一样，我们需要对数据进行预处理，比如测到的转录本总数(total_counts)、测到的基因总数(total_cells)、来源于线粒体基因的转录本所占比例等。

基础过滤

去除表达基因200以下的细胞；去除在3个细胞以下表达的基因。这也是通常Seurat通常用的默认过滤标准。

sc.pp.filter_cells：进行细胞的过滤，该函数保留至少有 min_genes 个基因（某个基因表达非0可判断存在该基因）的细胞，或者保留至多有 max_genes 个基因的细胞；

sc.pp.filter_genes：进行基因的过滤，该函数用于保留在至少 min_cells 个细胞中出现的基因，或者保留在至多 max_cells 个细胞中出现的基因；

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

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

data

质量过滤：

其实怎么过滤，取决于对自己数据的了解程度。

比如：可视化所有细胞中计数最多的基因。

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

plt.savefig("Highest_expr_genes.pdf") 

#我们可以看到前20个基因中有14个基因是属于核糖体基因，说明核糖体基因在这此数据集中表达丰度很高，有时候我们会去除核糖体占比过高的细胞，认为这些细胞质量低，这儿我们先不考虑核糖体基因，考虑线粒体相关基因。

例如，计算各种类型相关基因比例：线粒体相关基因，血红蛋白相关基因，核糖体相关基因，叶绿体相关基因等。

#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：计算与常出现基因的比，percent_top=[50] 计算与第 50 个最常出现基因的比例，None则不计算；

inplace：决定是否将计算指标添加到var和obs中；

log1p：设置为False可以跳过转换到log1p空间的过程；log1p即log(1+number)，用于压缩数据并确保结果是一个正数；

data

这样注释中就得到很多QC计算得到的信息。

使用violinplot度量查看各类质量，类似seurat：

n_genes_by_counts：每个细胞中，有表达的基因的个数；

total_counts：每个细胞的基因总计数（总表达量,umi数）；

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 基因的表达水平被认为主要是技术性的，因此除了计算每个细胞中线粒体基因，血红蛋白基因和核糖体基因所占的比例外，明智的做法是还要将线粒体基因，血红蛋白基因,MALAT1基因从数据集中直接删除，然后再进行任何进一步的分析。

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]

也可以根据基因总数进行进一步的过滤。

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

数据标准化：

数据预处理结束之后，就是和seurat一样进行标准化，挑选HVG基因了。

sc.pp.normalize_total(data, target_sum=1e4) ##标准化

sc.pp.log1p(data)  #将数据压缩到log1p的空间，就是取了个对数

下面就是类似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")

可以看到，执行完之后，就多出了uns:log1p和hvg的操作。

然后就是归一化，将数据放在单位方差。

sc.pp.regress_out(data, ['total_counts', 'pct_counts_mt']) #校正细胞基因计数和线粒体基因比例的影响。

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

data.write(results_file)

需要注意的是，在做完归一化后，data.X的数据格式从scipy.sparse.csr_matrix转换为numpy.ndarray。

下面，我们就开始做经典的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进行聚类，注意安装对应的python包，conda install -c conda-forge leidenalg ，当然也可使用其他的聚类算法，如louvain，sc.tl.louvain(data),sc.pl.umap(adata, color=['louvain'])，比较了一下，聚类结果差异不大

#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之间的表达差异，通过差异表达基因寻找Marker基因。

让我们计算每个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细胞相关的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")