我们前面也分享过很多篇关于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")
网友评论