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

最新推荐文章于 2025-03-28 17:55:49 发布

tzc_fly

最新推荐文章于 2025-03-28 17:55:49 发布

阅读量2k

点赞数 29

分类专栏：生物计算工具文章标签： python 开发语言

本文链接：https://blog.csdn.net/qq_40943760/article/details/138585941

版权

质量控制

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

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