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

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

质量控制

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

• 由于测序深度不够，测序数据中会存在很多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")

在这个案例中，首先对双细胞进行过滤（表达量过大），然后进行质控：过滤比如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)

环境RNA校正

对于基于液滴的单细胞RNA测序实验，在分配到含有细胞的液滴中存在一定量的背景mRNA，并且随着液滴一起被测序。这样做的结果是产生了一种背景污染，代表了不是来自液滴中包含的细胞，而是来自含有细胞的溶液中的RNA表达。

无细胞的 mRNA 分子，也被称为环境 RNA，混淆了观察到的计数的数量。对于无细胞 mRNA，纠正基于液滴的 scRNA-seq 数据集非常重要，因为它可能会扭曲我们下游分析中数据的解释。具体校正方法参考：soupX