单细胞最新入门教程系列（三）：单细胞数据分析中的质量控制

01、单细胞数据集构建和质量控制 

一旦基因表达被量化，它被总结为一个表达矩阵，每一行对应一个基因(或转录本)，每一列对应一个细胞。下一步，应检查矩阵以去除质量差的细胞。如果不能在这一阶段去除低质量的细胞，可能会增加技术噪声，从而有可能模糊下游分析中感兴趣的生物信号。

由于目前没有执行scRNA-seq的标准方法，因此本文将介绍的各种质量控制QC的期望值可能因实验而异。因此，为了执行QC，我们将寻找相对于数据集其余部分的异常值，而不是与独立的质量标准进行比较。因此，在比较使用不同测序协议的数据集的质量时，应该分情况讨论。      

02、Tung 数据集

为了说明细胞质量控制，我们考虑了芝加哥大学Yoav Gilad实验室中由三个不同个体（Tung等人，2017）产生的诱导多能干细胞数据集 （http://jdblischak.github.io/singleCellSeq/analysis/） 。实验在Fluidigm C1平台上进行，为了便于定量，使用了唯一分子标识符（UMI）和ERCC加标。数据文件位于工作目录中的 tung 文件夹中。这些文件是在 15/03/16 上制作的原始文件的副本。我们将出于可重复性目的使用这些 副本。

我们将使用 scater 包，以及 AnnotationDbi org.Hs.eg.db 将ENSEMBL ID转换为基因名称（符号）。

library(scater)library(SingleCellExperiment)library(AnnotationDbi)library(org.Hs.eg.db)library(EnsDb.Hsapiens.v86)

接下来，我们将读取矩阵和每个细胞的注释。后者转换为因子。

molecules<-read.delim("data/tung/molecules.txt",row.names=1)annotation<- read.delim("data/tung/annotation.txt",stringsAsFactors = T)

快速浏览一下数据集：

head(molecules[,1:3]) ## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03## ENSG00000237683 0 0 0## ENSG00000187634 0 0 0## ENSG00000188976 3 6 1## ENSG00000187961 0 0 0## ENSG00000187583 0 0 0## ENSG00000187642 0 0 0 head(annotation) ## individual replicate well batch sample_id## 1 NA19098 r1 A01 NA19098.r1 NA19098.r1.A01## 2 NA19098 r1 A02 NA19098.r1 NA19098.r1.A02## 3 NA19098 r1 A03 NA19098.r1 NA19098.r1.A03## 4 NA19098 r1 A04 NA19098.r1 NA19098.r1.A04## 5 NA19098 r1 A05 NA19098.r1 NA19098.r1.A05## 6 NA19098 r1 A06 NA19098.r1 NA19098.r1.A06

在这里，我们设置 altExp 包含 ERCC，从主对象中删除 ERCC 特征：

umi <- SingleCellExperiment(assays = list(counts = as.matrix(molecules)), colData = annotation)altExp(umi,"ERCC") <- umi[grep("^ERCC-",rownames(umi)), ]umi <- umi[grep("^ERCC-",rownames(umi),invert = T), ]

现在，让我们将ENSEMBL ID映射到基因符号。从命令中 table ，我们可以看到大多数基因都被注释了;但是，846返回了“NA”。默认情况下， mapIds 每个 ID 还原一个符号;可以使用参数 multiVals 更改此行为。

gene_names<-mapIds(org.Hs.eg.db, keys=rownames(umi), keytype="ENSEMBL", columns="SYMBOL",column="SYMBOL") ## 'select()' returned 1:many mapping between keys and columns rowData(umi)$SYMBOL<-gene_namestable(is.na(gene_names)) ## ## FALSE TRUE ## 18078 860

删除所有没有找到符号的基因：

umi<-umi[!is.na(rowData(umi)$SYMBOL),]

检查一下是否可以在新注释的符号中找到线粒体蛋白。

grep("^MT-",rowData(umi)$SYMBOL,value = T) ## named character(0)

奇怪的是，这什么也没返回。查找核糖体蛋白（以RPL或RPS开头）的类似命令按预期工作：

grep("^RP[LS]",rowData(umi)$SYMBOL,value = T)

快速搜索线粒体蛋白ATP8（也称为MT-ATP8）显示该名称不包含“MT-”。但是，正确的特征（ENSEMBL ID ENSG00000228253）存在于我们的注释中。

grep("ATP8",rowData(umi)$SYMBOL,value = T) ## ENSG00000143515 ENSG00000132932 ENSG00000104043 ENSG00000081923 ENSG00000130270 ## "ATP8B2" "ATP8A2" "ATP8B4" "ATP8B1" "ATP8B3" ## ENSG00000124406 ENSG00000228253 ## "ATP8A1" "ATP8"

大多数现代注释，例如 使用的 Cell Ranger 注释，将具有以 MT- 开头的线粒体基因名称。出于某种原因，我们发现的那个没有。注释问题通常很常见，应始终仔细考虑。在我们的例子中，我们也找不到基因的位置，因为染色体不受支持 org.Hs.eg.db ——这个数据库中没有基因组位置列： 

columns(org.Hs.eg.db) ## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" ## [11] "GENETYPE" "GO" "GOALL" "IPI" "MAP" ## [16] "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" ## [21] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG" ## [26] "UNIPROT"

让我们尝试一个不同的、更详细的数据库 - EnsDb.Hsapiens.v86 。利用这一资源，我们可以找到位于线粒体中的13个蛋白质编码基因：

ensdb_genes<-genes(EnsDb.Hsapiens.v86)MT_names<-ensdb_genes[seqnames(ensdb_genes) == "MT"]$gene_idis_mito<-rownames(umi) %in% MT_namestable(is_mito) ## is_mito## FALSE TRUE ## 18065 13

03、基本的质量控制

以下  scater  函数允许我们添加对数据集评估有用的每个细胞和每个基因的指标。每个细胞最受欢迎的指标是总数 （UMI）、检测到的基因总数、线粒体计数总数、线粒体计数百分比等。 umi_cell<-perCellQCMetrics(umi,subsets=list(Mito=is_mito))umi_feature<-perFeatureQCMetrics(umi)head(umi_cell) ## DataFrame with 6 rows and 9 columns## sum detected subsets_Mito_sum subsets_Mito_detected## <numeric> <numeric> <numeric> <numeric>## NA19098.r1.A01 61707 8242 4883 13## NA19098.r1.A02 62300 8115 3732 13## NA19098.r1.A03 42212 7189 3089 13## NA19098.r1.A04 52324 7863 3606 13## NA19098.r1.A05 69192 8494 4381 13## NA19098.r1.A06 66341 8535 3235 13## subsets_Mito_percent altexps_ERCC_sum altexps_ERCC_detected## <numeric> <numeric> <numeric>## NA19098.r1.A01 7.91320 1187 31## NA19098.r1.A02 5.99037 1277 31## NA19098.r1.A03 7.31782 1121 28## NA19098.r1.A04 6.89167 1240 30## NA19098.r1.A05 6.33166 1262 33## NA19098.r1.A06 4.87632 1308 30## altexps_ERCC_percent total## <numeric> <numeric>## NA19098.r1.A01 1.88730 62894## NA19098.r1.A02 2.00859 63577## NA19098.r1.A03 2.58694 43333## NA19098.r1.A04 2.31499 53564## NA19098.r1.A05 1.79124 70454## NA19098.r1.A06 1.93351 67649 head(umi_feature) ## DataFrame with 6 rows and 2 columns## mean detected## <numeric> <numeric>## ENSG00000187634 0.0300926 2.77778## ENSG00000188976 2.6388889 84.25926## ENSG00000187961 0.2384259 20.60185## ENSG00000187583 0.0115741 1.15741## ENSG00000187642 0.0127315 1.27315## ENSG00000188290 0.0243056 2.31481

现在，我们可以使用将上面计算的指标添加到每个细胞和每个基因元数据的函数：

umi<-addPerCellQC(umi, subsets=list(Mito=is_mito))umi<-addPerFeatureQC(umi)

手动过滤可以使用我们选择任何截止值。为了找到一个好的值，最好看一下分布：

hist( umi$total, breaks = 100)abline(v = 25000, col = "red")

hist( umi_cell$detected, breaks = 100)abline(v = 7000, col = "red")

有时很难想出一个明显的过滤截止点。在这种情况下，自适应阈值可以帮助我们识别在我们用于QC的任何变量中与中位数绝对偏差（MAD）相差超过3个中位数的点。注意指定偏差的正确方向：事实上，检测到的基因数量少，但MT基因百分比高，是低质量细胞的标志：

qc.lib2<-isOutlier(umi_cell$sum, log=TRUE, type="lower")attr(qc.lib2, "thresholds") ## lower higher## 23588.23 Inf qc.nexprs2<-isOutlier(umi_cell$detected, log=TRUE, type="lower")attr(qc.nexprs2, "thresholds") ## lower higher ## 6252.451 Inf qc.spike2<-isOutlier(umi_cell$altexps_ERCC_percent, type="higher")attr(qc.spike2, "thresholds") ## lower higher ## -Inf 3.619558 qc.mito2<-isOutlier(umi_cell$subsets_Mito_percent, type="higher")attr(qc.mito2, "thresholds") ## lower higher ## -Inf 9.294928 discard2<-qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2DataFrame(LibSize=sum(qc.lib2), NExprs=sum(qc.nexprs2), SpikeProp=sum(qc.spike2), MitoProp=sum(qc.mito2), Total=sum(discard2)) ## DataFrame with 1 row and 5 columns## LibSize NExprs SpikeProp MitoProp Total## <integer> <integer> <integer> <integer> <integer>## 1 47 65 137 75 194 上面执行的所有操作都可以在一个 scater 命令中完成： quickPerCellQC reasons<-quickPerCellQC(umi_cell, sub.fields=c("subsets_Mito_percent", "altexps_ERCC_percent"))colSums(as.matrix(reasons)) ## low_lib_size low_n_features high_subsets_Mito_percent ## 47 65 75 ## high_altexps_ERCC_percent discard ## 137 194

让我们添加另一个元数据列，该列将保留有关单元格是否被丢弃的信息： umi$discard<-reasons$discard

打赏