1.基于python的单细胞数据预处理-质量控制

质量控制

原始的单细胞测序数据具有以下特点：

由于测序深度不够，测序数据中会存在很多0值，这被称为drop-out事件；这是具有迷惑性的，因为对于一个细胞的具体基因表达量为0，我们不能判断到底是测序深度不够导致的，还是本身这个基因就不表达。
由于细胞状态不同，我们必然会测量到一些即将死亡的细胞，这些衰老细胞的基因表达情况相对正常细胞的数据也是一个误差。
测序技术本身不是百分百精确，比如液滴技术有时候会一个液滴中包含两个细胞，这样的基因表达量会非常高。一般正常细胞的基因表达量在3000-4000左右。

质量控制方法被用于初步过滤以上异常情况。随着大量的工程积累，下面演示为目前大家认可的最佳质量控制步骤。

使用的数据：这是一个基于10x-Multiome技术的数据集。该数据集测量了来自12名健康人类受试者在4个不同site（骨髓抽取的位置：手臂，胸骨等）的骨髓单核细胞的单细胞多组学数据（包含了嵌套的批次效应）。这里使用的数据是受试者8的样本4。

下载地址为：https://figshare.com/ndownloader/files/39546196。

读取数据：

import omicverse as ov
import scanpy as sc

ov.utils.ov_plot_set()

adata = sc.read_10x_h5("./data/filtered_feature_bc_matrix.h5")
print(adata)

"""
AnnData object with n_obs × n_vars = 16934 × 36601
    var: 'gene_ids', 'feature_types', 'genome'
"""

由于原始数据中有的obs_names或者var_names会重名，我们执行unique，在重名names字符串后自动添加1,2等字符使变量名唯一：

adata.var_names_make_unique()
adata.obs_names_make_unique()

adata.X形状为16934 × 36601，分别对应着barcodes和number of transcripts，var里有三个元素，gene_ids为来自Ensembl的基因id，genome为基因组的名字：

Ensembl ID是生物学领域中用于唯一标识基因和其他生物学序列的一种标识符。
如果已知genome（基因组的名字），我们可以通过额外的注释数据获得基因组所在的坐标，然后就能为RNA特征和ATAC特征建立关联。

print(adata.var.head())
"""
                    gene_ids    feature_types  genome
MIR1302-2HG  ENSG00000243485  Gene Expression  GRCh38
FAM138A      ENSG00000237613  Gene Expression  GRCh38
OR4F5        ENSG00000186092  Gene Expression  GRCh38
AL627309.1   ENSG00000238009  Gene Expression  GRCh38
AL627309.3   ENSG00000239945  Gene Expression  GRCh38
"""

过滤低质量细胞的指南

这是质量控制的第一步，主要针对三个质控协变量：

1.每个barcode的计数数量（计数深度），barcode即一个细胞样本；
2.每个barcode的基因数量；
3.每个barcode的线粒体基因计数比例；

当检测到基因数量较少，计数深度较低，线粒体计数较高时，细胞膜可能会破裂，这说明细胞正在死亡，这种样本一般都不是我们分析的主要目标，从而误导下游分析，应该去除。

说明：如果一个细胞正在死亡，那么其mRNA被释放到内环境，导致线粒体基因的比例较高，所以可以通过线粒体基因的比例来过滤掉低质量的单细胞测序数据。但是如果仅考虑一个变量可能会造成生物学误差，共同考虑三个 QC 协变量至关重要。例如，线粒体计数相对较高的细胞可能参与呼吸过程，不应被过滤掉。然而，计数低或高的细胞可能对应于静止细胞群或尺寸较大的细胞。故在过滤低质量细胞的时候，要同时考虑不同的QC协变量之间的关系。

对于简单的数据，可以观察数据分布来确定协变量过滤的阈值。随着数据规模的增长，手动观察阈值不可行，我们可以通过计算MAD（median absolute deviations）设置阈值，如果计数大于5倍的MAD，我们可以认为这是一个异常值。

这里使用的数据集是人类骨髓，因此线粒体计数是以“ MT-”前缀注释的。对于鼠数据集，前缀通常是小写“ mt-”。

除了统计线粒体基因比例，scanpy还能统计指定基因的计数比例，比如核糖体，血红蛋白基因：

# 线粒体基因
adata.var["mt"] = adata.var_names.str.startswith("MT-")
# 核糖体基因
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
# 血红蛋白基因
adata.var["hb"] = adata.var_names.str.contains(("^HB[^(P)]"))

sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt", "ribo", "hb"], inplace=True, percent_top=[20], log1p=True
)

print(adata)
"""
AnnData object with n_obs × n_vars = 16934 × 36601
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'
"""

在obs中有三个变量：

n_genes_by_counts：一个细胞中发现的有效基因数量（即表达量不为0）
total_counts：一个细胞中发现的分子数量（UMI），通常也可以被认为是这个细胞的文库大小（adata.obs['total_counts'] = adata.X.toarray().sum(axis=1) ）
pct_counts_mt：一个细胞中线粒体基因的表达计数占比（adata.obs['pct_counts_mt'] = adata[:, adata.var["mt"]].X.toarray().sum(axis=1) / adata.obs['total_counts'].values）

绘制协变量分布图，可以直观看到三个协变量：

mito_filter = 15
n_counts_filter = 4300
fig, axs = plt.subplots(ncols = 2, figsize = (8,4))
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt',ax = axs[0], show=False)
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts',ax = axs[1], show = False)
#draw horizontal red lines indicating thresholds.
axs[0].hlines(y = mito_filter, xmin = 0, xmax = max(adata.obs['total_counts']), color = 'red', ls = 'dashed')
axs[1].hlines(y = n_counts_filter, xmin = 0, xmax = max(adata.obs['total_counts']), color = 'red', ls = 'dashed')
fig.tight_layout()
plt.savefig("./result/2-1.png")

fig1
在这个案例中，首先对双细胞进行过滤（表达量过大），然后进行质控：过滤比如total_counts小于500的细胞，比如n_genes_by_counts小于250的细胞，线粒体基因的计数比例不超过15%。请不要忽视质量控制的每一步，这对后续分析尤为重要。

双细胞过滤

双细胞是可能存在的，两个小尺寸细胞容易被同一个液滴捕获，每个液滴被barcode标记，这也是用barcode而不是cell的原因（没有完美单细胞测量，只能说分辨率几乎接近单细胞）。双细胞存在下面两种：

同型：同型通常被认为是不影响下游分析的，因为其是由一类相同的细胞中的两个所构成，所以这部分细胞不是我们所需要过滤的对象
异型：异型通常是由来自两类不同的细胞所构成的，异型的存在会使得我们后续的细胞分类出现错误

使用sc中的scrublet去除异型双细胞：

# Original QC
n0 = adata.shape[0]
print(f'Original cell number: {n0}')

print('Begin of post doublets removal and QC plot')
sc.external.pp.scrublet(adata, random_state=112)
adata = adata[adata.obs['predicted_doublet']==False, :].copy()
n1 = adata.shape[0]
print(f'Cells retained after scrublet: {n1}, {n0-n1} removed.')
print(f'End of post doublets removal and QC plots.')

"""
Cells retained after scrublet: 16931, 3 removed.
"""

在双细胞过滤后，按照指南，过滤低质量读数细胞，我们分别演示手动和自动过滤，所以，复制两份adata：

adata_manual=adata.copy()
adata_auto=adata.copy()

手动过滤低质量读数细胞

首先定义一个过滤字典，然后按照字典进行过滤：

import numpy as np
tresh={'mito_perc': 15, 'nUMIs': 500, 'detected_genes': 250}

adata_manual.obs['passing_mt'] = adata_manual.obs['pct_counts_mt'] < tresh['mito_perc']
adata_manual.obs['passing_nUMIs'] = adata_manual.obs['total_counts'] > tresh['nUMIs']
adata_manual.obs['passing_ngenes'] = adata_manual.obs['n_genes_by_counts'] > tresh['detected_genes']

print(f'Lower treshold, nUMIs: {tresh["nUMIs"]}; filtered-out-cells: {n1-np.sum(adata_manual.obs["passing_nUMIs"])}')
print(f'Lower treshold, n genes: {tresh["detected_genes"]}; filtered-out-cells: {n1-np.sum(adata_manual.obs["passing_ngenes"])}')
print(f'Lower treshold, mito %: {tresh["mito_perc"]}; filtered-out-cells: {n1-np.sum(adata_manual.obs["passing_mt"])}')

保留剩余的细胞的交集：

QC_test = (adata_manual.obs['passing_mt']) & (adata_manual.obs['passing_nUMIs']) & (adata_manual.obs['passing_ngenes'])
removed = QC_test.loc[lambda x : x == False]
print(f'Total cell filtered out with this last  QC (and its chosen options): {n1-np.sum(QC_test)}')
adata_manual = adata_manual[QC_test, :].copy()
n2 = adata_manual.shape[0]
   
print(f'Cells retained after scrublet and filtering: {n2}, {n0-n2} removed.')

最后，再添加一步，直接过滤掉一些从基因和细胞层面低计数的细胞/基因：

sc.pp.filter_cells(adata_manual, min_genes=200)
sc.pp.filter_genes(adata_manual, min_cells=3)
print(adata_manual)

自动过滤低质量读数细胞

自动过滤也需要设置基础最低阈值，MAD计算可以获得最高阈值：

tresh={'pct_counts_mt': 15, 'total_counts': 500, 'n_genes_by_counts': 250}
adata_auto.obs['passing_mt'] = adata_auto.obs['pct_counts_mt'] < tresh['pct_counts_mt']
adata_auto.obs['passing_nUMIs'] = ov.pp._qc.mads_test(adata_auto.obs, 'total_counts', nmads=5, lt=tresh)
adata_auto.obs['passing_ngenes'] = ov.pp._qc.mads_test(adata_auto.obs, 'n_genes_by_counts', nmads=5, lt=tresh)  

nUMIs_t = ov.pp._qc.mads(adata_auto.obs, 'total_counts', nmads=5, lt=tresh)
n_genes_t = ov.pp._qc.mads(adata_auto.obs, 'n_genes_by_counts', nmads=5, lt=tresh)

print(f'Tresholds used, nUMIs: ({nUMIs_t[0]}, {nUMIs_t[1]}); filtered-out-cells: {n1-np.sum(adata_auto.obs["passing_nUMIs"])}')
print(f'Tresholds used, n genes: ({n_genes_t[0]}, {n_genes_t[1]}); filtered-out-cells: {n1-np.sum(adata_auto.obs["passing_ngenes"])}')
print(f'Lower treshold, mito %: {tresh["pct_counts_mt"]}; filtered-out-cells: {n1-np.sum(adata_auto.obs["passing_mt"])}')

剩下步骤与手动注释一样：

QC_test = (adata_auto.obs['passing_mt']) & (adata_auto.obs['passing_nUMIs']) & (adata_auto.obs['passing_ngenes'])
removed = QC_test.loc[lambda x : x == False]
print(f'Total cell filtered out with this last  QC (and its chosen options): {n1-np.sum(QC_test)}')
adata_auto = adata_auto[QC_test, :].copy()
n2 = adata_auto.shape[0]

print(f'Cells retained after scrublet and filtering: {n2}, {n0-n2} removed.')

sc.pp.filter_cells(adata_auto, min_genes=200)
sc.pp.filter_genes(adata_auto, min_cells=3)
print(adata_auto)