- 数据预处理
- 降维与聚类
我将分析 10X Genomics 免费提供的外周血单核细胞 (PBMC3k) 数据集。 数据为 2700 个单细胞在 Illumina NextSeq 500 上测序；
数据集下载地址：pbmc3k_filtered_gene_bc_matrices
下载后需解压得到以下目录：
|- data |- filtered_gene_bc_matrices |- hg19 |- barcodes.tsv |- genes.tsv |- matrix.mtx
导入相关包：
import numpy as np import pandas as pd import scanpy as sc
进行一些设置：
sc.settings.verbosity=3 # verbosity: errors (0), warnings (1), info (2), hints(提示) (3) sc.logging.print_header() sc.settings.set_figure_params(dpi=80,facecolor='white') result_file='./write/pbmc3k.h5ad' # 用于缓存分析结果
读入PBMC3k数据集：
adata=sc.read_10x_mtx(
'./data/filtered_gene_bc_matrices/hg19/', # .mtx文件的所在目录
var_names='gene_symbols', # 使用gene_symbols作为AnnData的特征名称
cache=True) # 写入缓存便于快速读写
adata.var_names_make_unique() # 如果在 sc.read_10x_mtx 中使用 var_names='gene_ids'，这是不必要的
print(adata)
"""
AnnData object with n_obs × n_vars = 2700 × 32738
var: 'gene_ids'
"""
针对AnnData的特征（基因），统计每个基因在所有细胞中的分布，并绘制这些基因分布的箱型图：
sc.pl.highest_expr_genes(adata,n_top=20)
使用两个工具：sc.pp.filter_cells进行细胞的过滤，sc.pp.filter_genes进行基因的过滤。这两个工具的说明可回顾Scanpy(一)AnnData数据结构与一些API用法介绍；过滤过程如下：
sc.pp.filter_cells(adata,min_genes=200) sc.pp.filter_genes(adata,min_cells=3)
上面过程为基因过滤（Basic filtering），下面我需要收集关于线粒体基因（mitochondrial genes）的信息，这些信息对质量控制（Quality Control，QC）很重要。
我们需要过滤包含线粒体基因和基因过多的细胞，因为：
- 检测出高比例的线粒体基因，表明细胞的测量质量较差；
- 表达基因过多可能是由于一个测量油滴包裹了多个细胞，从而检测出比正常情况要多的基因数，因此要过滤这些细胞。
首先，查看数据的特征量（AnnData.var）：
print(adata.var)
"""
gene_ids n_cells
AL627309.1 ENSG00000237683 9
... ... ...
SRSF10-1 ENSG00000215699 69
"""
对于var这个Dataframe的第0列，也就是index（基因名称），也被称为var_names，现在我们要为var中的基因做标记，判断哪些基因是线粒体相关的：
# startswith()方法用于检查字符串是否是以指定子字符串开头, 如果是则返回 True, 否则返回 False
adata.var['mt']=adata.var_names.str.startswith('MT-') # var新增一列'mt', 这是用于标记基因是否是线粒体相关的一个mask
使用 pp.calculate_qc_metrics，我们可以计算许多QC指标：
sc.pp.calculate_qc_metrics(adata,qc_vars=['mt'],percent_top=None,log1p=False,inplace=True)
扩展部分：QC与sc.pp.calculate_qc_metrics
在单细胞分析中，我们通常要对某些基因计算数据指标，根据指标过滤数据，以便保证后续计算的质量，这就是质量控制；
sc.pp.calculate_qc_metrics是用于计算QC指标的函数：
sc.pp.calculate_qc_metrics(data, qc_vars=(), percent_top=(50, 100, 200, 500), inplace=False, log1p=True)
参数说明：
- data：Anndata；
- qc_vars：用于标识要控制的特征（基因），布尔型元素，用于作为mask使用；
- percent_top：计算与常出现基因的比，percent_top=[50] 计算与第 50 个最常出现基因的比例，None则不计算；
- inplace：决定是否将计算指标添加到var和obs中；
- log1p：设置为False可以跳过转换到log1p空间的过程；log1p即log(1+number)，用于压缩数据并确保结果是一个正数；
实例如下：
import seaborn as sns
import matplotlib.pyplot as plt
pbmc=sc.datasets.pbmc3k() # 自动下载到./data/pbmc3k_raw.h5ad
# 标记线粒体基因, 并保存为 mask
pbmc.var['mito']=pbmc.var_names.str.startswith('MT-')
print(pbmc.var)
"""
gene_ids mito
index
MIR1302-10 ENSG00000243485 False
... ... ...
AC002321.1 ENSG00000215611 False
"""
print(pbmc.var.columns)
"""
Index(['gene_ids', 'mito'], dtype='object')
"""
print(pbmc.obs.columns)
"""
Index([], dtype='object')
此时obs只有index, obs的index为细胞的名称
"""
# 计算线粒体基因的 QC 指标
sc.pp.calculate_qc_metrics(pbmc,qc_vars=['mito'],inplace=True)
print(pbmc.var.columns) # 关于基因轴上的QC指标
"""
Index(['gene_ids', 'mito', 'n_cells_by_counts', 'mean_counts',
'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts',
'log1p_total_counts'],
dtype='object')
"""
print(pbmc.obs.columns) # 关于细胞轴上的QC指标
"""
Index(['n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts',
'log1p_total_counts', 'pct_counts_in_top_50_genes',
'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes',
'pct_counts_in_top_500_genes', 'total_counts_mito',
'log1p_total_counts_mito', 'pct_counts_mito'],
dtype='object')
"""
# 对obs中的计算指标可视化
sns.jointplot(data=pbmc.obs,
x="log1p_total_counts",
y="log1p_total_counts_mito",
kind="hex")
plt.show()
sns.histplot(pbmc.obs["pct_counts_mito"])
plt.show()
total_counts记录了对每个细胞，所包含基因的总数；
total_counts_mito记录了对每个细胞，所包含控制基因（qc_vars中列出的基因）的总数；
pct_counts_mito记录了对每个细胞，含线粒体的基因总数；（其中 ‘mito’ 来自var中的index）
现在，我们用小提琴图可视化两个QC指标：
# jitter 表示抖动 sc.pl.violin(adata, ['pct_counts_mt'], jitter=0.4) sc.pl.violin(adata, ['n_genes_by_counts'], jitter=0.4)
从上面两个指标可以看出，大部分细胞的线粒体基因数在5以下，大部分细胞的基因总数在2500以下，所以我们要过滤掉那些异常的细胞：
# 去除线粒体基因过多或基因总数过多的细胞 adata = adata[adata.obs.n_genes_by_counts < 2500, :] adata = adata[adata.obs.pct_counts_mt < 5, :]
然后对数据标准化，并转换到 log1p 空间：
# 标准化 sc.pp.normalize_total(adata, target_sum=1e4) # 压缩到logp1空间 sc.pp.log1p(adata)
确定高变基因并过滤数据：
# 确定高变基因
sc.pp.highly_variable_genes(adata)
"""
extracting highly variable genes
finished (0:00:00)
--> added
'highly_variable', boolean vector (adata.var) # 布尔型mask
'means', float vector (adata.var)
'dispersions', float vector (adata.var)
'dispersions_norm', float vector (adata.var)
"""
# 根据高变基因过滤数据
adata=adata[:,adata.var.highly_variable]
对数据进行回归并校正数据：
""" sc.pp.regress_out(data, keys) keys:要回归的data.obs中的数据列 """ # 回归 adata.obs 中的列 (columns) sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
将数据adata.X缩放到单位方差和零均值，对于缩放后的数据，在值为10处截断：
sc.pp.scale(adata, max_value=10)降维与聚类
我们可以用PCA进行降维：
# PCA降维并可视化 sc.tl.pca(adata, svd_solver='arpack') sc.pl.pca(adata, color='CST3') # 写入缓存 adata.write(result_file)
sc.tl.pca
使用PCA对数据进行变换；
sc.pl.pca
sc.pl.pca(adata,color=None)
注意参数color，我们可以传递adata中的'ann1' or ['ann1', 'ann2']，这些ann可以是基因，可以是聚类结果，所以我们可以根据基因标记颜色，也可以根据聚类算法结果标记颜色；
tl用于做变换，pl用于可视化（只取前2个维度），调用pl前需要先找到tl计算得到的数据空间；
比如上图，PCA将高维数据adata.X聚类到2维空间后得到的只是一些平面下的散点，但我们可以根据每个散点（细胞）中包含基因CST3的数量为散点标记颜色。
我们也可以用UMAP进行降维并可视化，我们现在使用的数据是经过PCA变换后的数据（保留了所有成分），使用这个数据并用UMAP方法做变换：
# UMAP 降维得到 embedding 并可视化 sc.pp.neighbors(adata,n_neighbors=10,n_pcs=40) sc.tl.umap(adata) sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])
现在，利用经过UMAP变换后的数据，使用 leiden 方法进行聚类：
# 用 leiden 方法聚类
sc.tl.leiden(adata)
"""
running Leiden clustering
finished: found 8 clusters and added
'leiden', the cluster labels (adata.obs, categorical) (0:00:00)
"""
print(adata.obs)
"""
n_genes n_genes_by_counts ... pct_counts_mt leiden
AAACATACAACCAC-1 781 779 ... 3.017776 0
AAACATTGAGCTAC-1 1352 1352 ... 3.793596 2
... ... ... ... ... ...
TTTGCATGCCTCAC-1 724 723 ... 0.806452 0
"""
# 用 UMAP 可视化聚类结果
sc.pl.umap(adata, color='leiden')
# 写入结果
adata.write(result_file)
