5.1 数据集构建和质控
5.1.1 简介
基因表达定量后形成一个表达矩阵，其中每一行对应一个基因（或转录本），每一列对应一个细胞。下一步我们将检查矩阵以去除低质量的细胞。在此阶段如果未能去除低质量细胞，可能会增加技术噪音，从而掩盖下游分析中感兴趣的生物信号。
由于目前没有进行scRNA-seq的标准方法，下文介绍的各种QC手段的期望值在不同实验方法之间可能会有很大差异。我们将寻找相对于数据集其余部分异常的细胞，而不是与独立的QC标准进行比较。因此，在比较使用不同实验方法测序的数据集的质控指标时应小心谨慎。
5.1.2 Tung数据集
为了说明细胞质控，我们使用芝加哥大学Yoav Gilad实验室的三个不同个体产生的诱导多能干细胞数据集（Tung等，2017）。实验在Fluidigm C1平台上进行，为了便于定量，使用UMI和外源RNA对照（ERCC）。由于基于液滴的方法使用的快速增长，ERCC不再被广泛使用；但是它们可以作为低通量方法的参照。数据文件位于工作目录中的tung文件夹中（下载地址：https://singlecellcourse.cog.sanger.ac.uk/index.html?shared=data/）。
使用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")
> rowData(umi)$SYMBOL <- gene_names
> table(is.na(gene_names))
FALSE TRUE
18079 859
删除所有没有找到名称的基因：
> umi <- umi[! is.na(rowData(umi)$SYMBOL),]
检查是否可以在新注释的名称中找到线粒体蛋白质。
> grep("^MT-",rowData(umi)$SYMBOL,value = T)
奇怪的是，该命令什么也没返回。而查找核糖体蛋白（以RPL或RPS开头）的类似命令按预期工作：
> grep("^RP[LS]",rowData(umi)$SYMBOL,value = T)
ENSG00000116251 ENSG00000142676 ENSG00000117676 ENSG00000142937 ENSG00000122406
"RPL22" "RPL11" "RPS6KA1" "RPS8" "RPL5"
ENSG00000177954 ENSG00000136643 ENSG00000138326 ENSG00000177600 ENSG00000166441
"RPS27" "RPS6KC1" "RPS24" "RPLP2" "RPL27A"
快速搜索线粒体蛋白ATP8（MT-ATP8），显示名称中不包含“MT-”。但是，我们的注释中存在该基因（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_id
> is_mito <- rownames(umi) %in% MT_names
> table(is_mito)
is_mito
FALSE TRUE
18066 13
5.1.3 基本QC
scater函数允许我们添加对数据集评估有用的指标。对于细胞最常用的指标是总counts数（UMI）、检测到的基因总数、线粒体counts总数、线粒体counts百分比等。
> 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 61706 8242 4883 13
NA19098.r1.A02 62298 8113 3732 13
NA19098.r1.A03 42211 7188 3089 13
NA19098.r1.A04 52323 7862 3606 13
NA19098.r1.A05 69193 8495 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.91333 1187 31
NA19098.r1.A02 5.99056 1277 31
NA19098.r1.A03 7.31800 1121 28
NA19098.r1.A04 6.89181 1240 30
NA19098.r1.A05 6.33157 1262 33
NA19098.r1.A06 4.87632 1308 30
altexps_ERCC_percent total
<numeric> <numeric>
NA19098.r1.A01 1.88733 62893
NA19098.r1.A02 2.00865 63575
NA19098.r1.A03 2.58700 43332
NA19098.r1.A04 2.31503 53563
NA19098.r1.A05 1.79121 70455
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")
当很难确定一个明显的过滤阈值时，自适应阈值可以帮助我们识别距离中位数超过3个中位数绝对偏差（median absolute deviations, MAD）的点，但需要注意确定数据偏离正常范围的正确方向。具体来说，如果一个细胞检测到的基因数量较少，但线粒体基因（MT基因）的比例较高，这是低质量细胞的标志：
> qc.lib2 <- isOutlier(umi_cell$sum, log=TRUE, type="lower")
> attr(qc.lib2, "thresholds")
lower higher
23594.12 Inf
> qc.nexprs2 <- isOutlier(umi_cell$detected, log=TRUE, type="lower")
> attr(qc.nexprs2, "thresholds")
lower higher
6253.217 Inf
> qc.spike2 <- isOutlier(umi_cell$altexps_ERCC_percent, type="higher")
> attr(qc.spike2, "thresholds")
lower higher
-Inf 3.6197
> qc.mito2 <- isOutlier(umi_cell$subsets_Mito_percent, type="higher")
> attr(qc.mito2, "thresholds")
lower higher
-Inf 9.29524
> discard2 <- qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2
> DataFrame(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
绘制多个coldata（细胞元数据）结果之间的相互关系，能够阐明它们之间的依赖关系。例如，线粒体含量高的细胞通常被认为是死亡的或正在死亡；这些细胞通常也具有较低的总UMI数和检测到的基因数量。
> plotColData(umi, x="sum", y="subsets_Mito_percent", colour_by="discard")
> plotColData(umi, x="sum", y="detected", colour_by="discard")
> plotColData(umi, x="altexps_ERCC_percent", y="subsets_Mito_percent",colour_by="discard")
还可以绘制分批次的coldata图表，查看是否存在显著的批次差异：
> library(scales)
> plotColData(umi, x="sum", y="detected", colour_by="discard", other_fields = "individual") +
facet_wrap(~individual) + scale_x_continuous(labels = unit_format(unit = "k", scale = 1e-3))
> plotColData(umi, x="sum", y="detected", colour_by="discard", other_fields = "replicate") +
facet_wrap(~replicate) + scale_x_continuous(labels = unit_format(unit = "k", scale = 1e-3))
5.1.4 高表达基因
让我们看一下整个数据集中表达最高的基因，使用上面获得的基因名称表示。大多数基因都是线粒体或核糖体蛋白，这对于大多数scRNA-seq数据集来说很典型。
> plotHighestExprs(umi, exprs_values = "counts",
feature_names_to_plot = "SYMBOL", colour_cells_by="detected")
保留在2个或更多细胞中检测到的基因（表达值>1），这将过滤掉大约4,000个表达较低的基因。
> keep_feature <- nexprs(umi,byrow = TRUE,detection_limit = 1) >= 2
> rowData(umi)$discard <- ! keep_feature
> table(rowData(umi)$discard)
FALSE TRUE
13876 4203
接下来创建一个新的assay，logcounts_raw，它将包含log2转换的counts并添加伪计数1。
> assay(umi, "logcounts_raw") <- log2(counts(umi) + 1)
最后，保存SingleCellExperiment对象，其中包含我们添加到细胞元数据的所有字段以及新assay（logcounts_raw）：
> saveRDS(umi, file = "data/tung/umi.rds")
往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——4. 使用Bioconductor进行scRNA-seq分析（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——4. 使用Bioconductor进行scRNA-seq分析（2）
