单细胞转录组流程三：超详细！Scanpy打通单细胞常规流程

有人可能会说：单细胞分析使用Seurat，monocle等R包会更加方便。但是实际分析中,测试情况是当细胞量大于5万时。一般小型服务器内存很容易不足，这时候请不要过多尝试使用Seurat来进行分析，monocle2更是。而基于python的单细胞转录分析包scanpy，能很好的解决内存不足的问题，亲测整合80万细胞量时，整个预处理流程在6小时左右能够完成。
scanpy相关python 包安装（安装好python3之后，终端运行），一定要注意的是，这里的python最好不要是3.9版本往上的，否则涉及多单细胞数据整合加载bbknn包会报错，numba将无法与当前环境适配。

pip install -i https://pypi.doubanio.com/simple/ scanpy==1.6.1 
##注意要加上-i参数，给pip加上一个豆瓣或者清华(https://pypi.tuna.tsinghua.edu.cn/simple)的镜像，否则下载起来你即使加了````--default-timeout=1000```也无济于事
#我个人建议在conda中新建一个环境，执行：
conda create -n scrna_test python=3.7.0 numba=0.55.2 pandas=1.1.5 llvmlite=0.38.1 numpy=1.21.6 bbknn=1.5.1
#然后在scanpy官网上下载scanpy-1.9.1-py3-none-any.whl（https://pypi.org/project/scanpy/#files）
wget https://files.pythonhosted.org/packages/51/87/a55c7992cba9b189de70eae37e9f1e2abe6fdaf3f087d30356f28698948e/scanpy-1.9.1-py3-none-any.whl
#下载好后
pip install scanpy-1.9.1-py3-none-any.whl 
#下载scanpy非常的困难，因为他对numba和 llvmlite有版本要求
#如果你在安装scanpy没有安装好pandas, numpy,numba,  llvmlite,事情会变得非常复杂，报错一个接一个。
#conda env create -n env_name pakage1=version package2=version... 非常好用,最好是通过文献看一下单细胞的开发环境，然后把他们复制过来。
#然后也可以看一下自己的conda环境下pip list都是哪些版本，也可以进行移植。

1. 运行python3，导入相关包及设置一些必要路径：

import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sc.settings.verbosity = 3 # verbosity即冗余 。设置日志等级: errors (0), warnings (1), info (2), hints (3)，取值表示测试日志结果显示的详细程度，数字越大越详细。
sc.logging.print_versions() # 输出版本号
sc.settings.set_figure_params(dpi=80) # set_figure_params 设置图片的分辨率/大小以及其他样式
import os #在服务器运行，习惯性会设置一个输出路径，用于保存pdf图片
os.getcwd()  #查看当前路径
os.chdir('./filtered_gene_bc_matrices/scanpy') #修改路径
os.getcwd()
results_file = 'pbmc3k.h5ad' #声明什么用于储存分析结果

2. 读取并查看数据：

scanpy可以读取.h5 .h5ad 以及10X标准化(features.tsv,barcodes.tsv, matrix.mtx)格式数据

data = sc.read_10x_h5("./pbmc3K.h5",genome=None, gex_only=True)
print(data)
AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'
"""
#cell数700 基因数32738的矩阵
#也可以是data = sc.read_10x_matrix('path',var_names = 'gene_symbols', cache = True)
#使用gene_symbols作为AnnData的特征名称
#cache=True , cache为True表示写入缓存(cache) 便于快速读写
#如果是h5ad格式，可以直接read('/.h5ad')
#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]
"""

AnnData数据结构

AnnData各部分数据类型

	功能	数据类型
data.X	矩阵数据	numpy,scipy sparse,matrix
data.obs	观察值数据	pandas dataframe
data.var	特征和高变基因数据	pandas dataframe
data.uns	非结构化数据	dict

3. 索引去重：

data.var_names_make_unique()
# # 如果在 sc.read_10x_mtx 中使用了var_names='gene_ids'，这一步就是不必要的

4. 质量控制：

质量控制3个指标：测到的转录本总数(total_counts)、测到的基因总数(total_cells)、来源于线粒体基因的转录本所占比例。

4.1基本过滤

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
"""
AnnData object with n_obs × n_vars = 2700 × 13714
    obs: 'n_genes'
    var: 'gene_ids', 'n_cells'
"""

4.2 计算质量控制指标：

*选择阈值去除高表达量的细胞，阈值很大程度上取决于对自己项目的了解程度，因为不同器官组织提取的单细胞，线粒体基因平均水平不一样。

## startswith()方法用于检查字符串是否是以指定子字符串开头, 如果是则返回 True, 否则返回 False
# mitochondrial genes
data.var['mt'] = data.var_names.str.startswitch('MT-')  
# hemoglobin genes.  血红蛋白基因
data.var['hb'] = data.var_names.str.contains('^HB[^P]')
# ribosomal genes
data.var['ribo'] = data.var_names.str.startswitch('RPS','RPL')
"""
AL627309.1       False
                 ...  
SRSF10-1         False
Name: mt, Length: 13714, dtype: bool
"""

sc.pp.calculate_qc_metrics(data, qc_vars=['mt',‘hb’,'ribo'], percent_top=None, log1p=False, inplace=True)
print(data)
#data：Anndata；
#qc_vars：用于标识要控制的特征（基因），布尔型元素，用于作为mask使用；
#percent_top：计算与常出现基因的比，percent_top=[50] 计算与第 50 个最常出现基因的比例，None则不计算；
#inplace：决定是否将计算指标添加到var和obs中；
#log1p：设置为False可以跳过转换到log1p空间的过程；log1p即log(1+number)，用于压缩数据并确保结果是一个正数；
"""
注释中多了很多QC计算得到的信息
AnnData object with n_obs × n_vars = 2700 × 13714
    obs: 'n_genes', 'nFeature_RNA', nCount_RNA', 'total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
"""

# 绘制高表达的前20个基因
sc.pl.highest_expr_genes(data, n_top=20, save='_pbmc3k.png')

n_top=20

"""
使用violinplot度量以下质量：
n_genes_bt_counts：每个细胞中，有表达的基因的个数；
total_counts：每个细胞的基因总计数（总表达量,umi数）；
pct_counts_mt：每个细胞中，线粒体基因表达量占该细胞所有基因表达量的百分比
pct_counts_hb:每个细胞中，血红蛋白基因表达量占该细胞所有基因表达量的百分比
pct_counts_ribo:每个细胞中，核糖体RNA基因表达量占该细胞所有基因表达量的百分比
"""
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)

1664981487608.png

# filter for percent mito
data = data[data.obs['pct_counts_mt'] < 20, :]
# filter for percent ribo > 0.05
data = data[data.obs['pct_count_ribo'] > 5,:]

由于线粒体和 MALAT1 基因的表达水平被认为主要是技术性的，因此除了计算每个细胞中线粒体基因，血红蛋白基因和核糖体基因所占的比例外，明智的做法是还要将线粒体基因，血红蛋白基因,MALAT1基因从数据集中直接删除，然后再进行任何进一步的分析
#delete var_names.str.startswitch['MT-'] 
#var_names.str.startswitch('MALAT1') var_names.str.contains('^HB[^(P)]')
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进行切片(类似于Seurat中的subset)

data = data[data.obs['n_genes_by_counts'] < 2500, :]
data = data[data.obs['n_genes_by_counts'] > 500,:]
"""
View of AnnData object with n_obs × n_vars = 2638 × 13714
    obs: 'n_genes', 'nFeature_RNA', nCount_RNA', 'total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
"""

我们先把矩阵提取出来看一下表达量的值

mat = pd.DataFrame(data=test.X.todense(),index=test.obs_names,columns=test.var_names)
mat

注意原始值都是整数

image.png

5. 数据标准化：

Normalize each cell by total counts over all genes, so that every cell has the same total count after normalization.
总计数标准化 ，以便细胞之间的基因表达量具有可比性

sc.pp.normalize_total(data, target_sum=1e4)
## normalize with total UMI count per cell

函数normalize_total可以对每个细胞进行标准化，以便每个细胞在标准化后沿着基因方向求和具有相同的总数 ````
target_sum：

data.X
array([[ 3.,  3.,  3.,  6.,  6.],
       [ 1.,  1.,  1.,  2.,  2.],
       [ 1., 22.,  1.,  2.,  2.]], dtype=float32)
# 设置 target_sum=1 标准化后
X_norm
array([[0.14, 0.14, 0.14, 0.29, 0.29],
       [0.14, 0.14, 0.14, 0.29, 0.29],
       [0.04, 0.79, 0.04, 0.07, 0.07]], dtype=float32)

看一下做完标准化的结果

image.png

没取标准化之前：

View of AnnData object with n_obs × n_vars = 15072 × 21655
    obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_hb', 'pct_counts_hb', 'total_counts_ribo', 'pct_counts_ribo', 'doublet_score', 'predicted_doublet'
    var: 'n_cells', 'mt', 'hb', 'ribo', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
    uns: 'scrublet'

标准化之后：

AnnData object with n_obs × n_vars = 15072 × 21655 
obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_hb', 'pct_counts_hb', 'total_counts_ribo', 'pct_counts_ribo', 'doublet_score', 'predicted_doublet' 
var: 'n_cells', 'mt', 'hb', 'ribo', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts' 
uns: 'scrublet'

标准化前后只是改变了data.X 没有增加obs或者var的内容。

#将数据压缩到log1p的空间
sc.pp.log1p(data)

log1p 也就是log(X+1) 也就是In(X+1)

image.png
image.png

5. 高变基因选取：

高变异基因就是highly variable features（HVGs），就是在细胞与细胞间进行比较，选择表达量差别最大的基因，Seurat使用FindVariableFeatures函数鉴定高变基因，这些基因在不同细胞之间的表达量差异很大（在一些细胞中高表达，在另一些细胞中低表达）。默认情况下，会返回2,000个高变基因用于下游的分析。
利用FindVariableFeatures函数，会计算一个mean-variance结果，也就是给出表达量均值和方差的关系并且得到top variable features，这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因，用于后续鉴定细胞类型。
标记基因 (marker gene)，是一种已知功能或已知序列的基因，能够起着特异性标记的作用。
标记基因通常是高变基因中很小的子集；

#识别高度可变基因，并进行过滤：
sc.pp.highly_variable_genes(data,min_mean = 0.0125,max_mean=3,min_disp=0.5)
sc.pl.highly_variable_genes(data)
#保存原始数据后再把data变成只有高变基因
data.raw = data
#过滤
##注意切片data[obs:var]
data = data[:,data.var['highly_variable']]

image.png

6. 归一化(将数据缩放到单位方差)：

>>>将数据data.X缩放到单位方差和零均值，对于缩放后的数据，在值为10处截断：

""" sc.pp.regress_out(data, keys) keys:要回归的data.obs中的数据列 """
# 回归 adata.obs 中的列 (columns) 
# 回归每个细胞的总计数和表达的线粒体基因的百分比的影响。
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

如果你这时候想看矩阵情况，需要进行转换

s = scipy.sparse.csr_matrix(data.X)
mat = pd.DataFrame(s.todense(),index = data.obs_names,colums = data.var_names)

image.png

7. 主成分分析：

通过运行主成分分析(PCA)来降低数据的维数，该分析揭示了变化的主轴并对数据进行去噪。

#svd_solver：指定奇异值分解SVD的方法，有4个可以选择的值：{‘auto’, ‘full’, ‘arpack’, ‘randomized’，}
#如果设为’arpack’，则使用scipy.sparse.linalg.svds计算SVD分解。这种方法严格要求 0 < n_components < min(样本数，特征数)。
sc.tl.pca(data,svd_solver = 'arpack')
sc.pl.pca(data, color = 'CST3')

PCA将高维数据data.X聚类到2维空间后得到的只是一些平面下的散点，但我们可以根据每个散点（细胞）中包含基因CST3的数量为散点标记颜色。

计算单个pc对数据总方差的贡献,这可以提供给我们应该考虑多少个 PC 以计算细胞的邻域关系的信息(resolution选择多少)，例如用于后续的聚类函数 sc.tl.louvain()

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

8.降维（对neighborhood graph进行embedding）

sc.pp.neighbors(data,n_neighbors=10,n_pcs=25)
"""
computing neighbors
    using 'X_pca' with n_pcs = 40
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:05)
"""
sc.tl.umap(data)

9. 聚类

scanpy提供leiden和louvain两种图聚类算法，图聚类在开始之前要先找到邻居。
首先计算neighborhood graph，我们使用数据矩阵的PCA表示来计算细胞的邻域图。可以在这里简单地使用默认值。为了可视化聚类的结果我们做一下umap降维，看看分群结果。

sc.tl.leiden(data,resolution = 0.3)
#用umap可视化聚类结果
sc.pl.umap(data,color = 'leiden', frameon=False, title='',
           legend_loc='right margin', legend_fontsize=12,
           legend_fontweight='normal', legend_fontoutline=6,
           palette=None, cmap=None)
#platelet是指给clutesr指定一组颜色。例如color=['red','blue'....]
#cmap是指颜色选择
#legend_fontweight 是指字体粗细
#legend_loc='on data'
#legend_fontoutline:图例字体轮廓的线宽，单位为pt。使用withStroke的路径效果绘制白色轮廓。
#. 当frameon=True的时候, 图示会被绘制在一个patch实体上;
# 否则, 如果frameon=False, 则图示会被直接绘制在图片上。
#这里, 讨论是否将图示绘制在一个patch实体上的意义在于,
#当把它绘制在一个patch实体上时, 
#我们才可以使用facecolor, edgecolor, framealpha, fancybox等参数来设置图示的背景(不是图片的背景)的颜色, 边框颜色, 透明度, 以及形状, 而当frameon=False的时候这些参数就会失效

10. Marker基因查找

通过文献或cellmarker查找marker 基因

PBMC common marker genes

#先设置分辨率以及长宽
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)

1665051072562.png

11. 通过差异表达基因寻找Marker基因

让我们计算每个cluster中高度差异基因的排名，最简单和最快的方法是t-test，当然还有wilcoxon。

sc.settings.set_figure_params(dpi=80,dpi_save=600,figsize=(10,10))
sc.tl.rank_genes_groups(data,'leiden',method='t-test')
#绘图指定每个cluster前多少个基因，每个cluster之间是否共享
sc.pl.rank_genes_groups(data,n_genes=25,sharey=False,fontsize=15)

#sharely ： 控制是否应共享每个panel的y轴，sharey =False表示每个panel都有自己的y轴范围。

image.png