由于网速的原因，昨天的测试数据没有下载成功，在我挂机下载了一整个晚上之后显示

image.png

我也只能作罢。。。
先看看第二步吧

Analyzing single-cell RNA-seq data containing read counts

1 Overview

• In this workflow, we use a relatively simple dataset (Lun et al. 2017) to introduce most of the concepts of scRNA-seq data analysis.

This dataset contains two plates of 416B cells (an immortalized mouse myeloid progenitor cell line), processed using the Smart-seq2 protocol (Picelli et al. 2014). A constant amount of spike-in RNA from the External RNA Controls Consortium (ERCC) was also added to each cell’s lysate prior to library preparation.

• Counts for all genes/transcripts in each cell are available from ArrayExpress using the accession number E-MTAB-5522. We download both the count tables (in the “processed files”) as well as the metadata file using the BiocFileCache package. This saves the files to a local cache (raw_data) and avoids re-downloading them if they are already present.

BiocManager::install("BiocFileCache")
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
lun.zip <- bfcrpath(bfc, 
                    file.path("https://www.ebi.ac.uk/arrayexpress/files",
                              "E-MTAB-5522/E-MTAB-5522.processed.1.zip"))
lun.sdrf <- bfcrpath(bfc, 
                     file.path("https://www.ebi.ac.uk/arrayexpress/files",
                               "E-MTAB-5522/E-MTAB-5522.sdrf.txt"))
unzip(lun.zip, exdir=tempdir())

2 Setting up the data

2.1 Loading in the count matrix

• 第一步 将计数矩阵加载到内存中。分别为研究中使用的每个细胞板生成了一个矩阵plate1、plate2。在每个矩阵中，每行代表基因或转录本，每列表示一个细胞（跟RNAseq表达矩阵类似，行名是基因名ensemble/entrizID，列名是样本名GSM）。随后，矩阵的每个条目中的计数表示映射到特定细胞中特定基因/转录的读取次数。

plate1 <- read.delim(file.path(tempdir(), "counts_Calero_20160113.tsv"), 
    header=TRUE, row.names=1, check.names=FALSE)
plate2 <- read.delim(file.path(tempdir(), "counts_Calero_20160325.tsv"), 
    header=TRUE, row.names=1, check.names=FALSE)

gene.lengths <- plate1$Length # First column is the gene length.
plate1 <- as.matrix(plate1[,-1]) # Discarding gene length (as it is not a cell).
plate2 <- as.matrix(plate2[,-1])
rbind(Plate1=dim(plate1), Plate2=dim(plate2))

• 第二步 将两个矩阵合并到一个对象中，以便进一步处理。ps(这是在验证两个矩阵之间的基因顺序相同后完成的。)

stopifnot(identical(rownames(plate1), rownames(plate2)))
all.counts <- cbind(plate1, plate2)

• 为方便起见，计数矩阵存储在SingleCellExperiment包中的SingleCellExperiment对象中。这允许将不同类型的行级和列级元数据与计数一起存储，以便在整个工作流中同步操作。(这一步也和DEGlist的获得方法类似)

library(SingleCellExperiment)
sce <- SingleCellExperiment(list(counts=all.counts))
rowData(sce)$GeneLength <- gene.lengths
sce

• 从行名称中标识与 ERCC spike-in 转录脚本对应的行。我们将此信息存储在 SingleCellExperiment对象中，以便将来使用。这是必需的，因为spike-in在下游步骤如规一化过程中需要特殊处理。

isSpike(sce, "ERCC") <- grepl("^ERCC", rownames(sce))
summary(isSpike(sce, "ERCC"))

• 此数据集稍有不寻常，因为它包含来自另一组spike-in转录本（Spike-In RNA Variants (SIRV) set)的信息。为简单起见，我们将仅在此分析中使用 ERCC spike-in。因此，在进一步分析之前，我们必须删除与 SIRV 转录本对应的行，只需通过子设置SingleCellExperiment对象即可完成。

is.sirv <- grepl("^SIRV", rownames(sce))
sce <- sce[!is.sirv,] 
summary(is.sirv)

Comments from Aaron:

• 某些要素计数工具将在计数矩阵中报告映射统计信息（例如，未对齐或未分配的读取数）。虽然这些值对质量控制很有用，但如果将其视为基因表达值，则它们会具有误导性。因此，在进一步分析之前，应将其移除（或至少移动到 colData）。
• 请注意，在计数矩阵的行名称为基因符号的人类数据中使用 ^ERCC 正则表达式。ERCC基因家族实际上存在于人类注释中，因此这将导致基因作为spike-in转录本的错误识别。通过发布具有标准标识符的计数矩阵（例如，Ensembl、Entrez），可以避免此问题。

2.2 Incorporating cell-based annotation

• 从 sdrf.txt 文件加载每个库/单元格的元数据。请务必检查元数据表的行与计数矩阵的列的顺序相同。否则，不正确的元数据将分配给每个单元格。

metadata <- read.delim(lun.sdrf, check.names=FALSE, header=TRUE)
m <- match(colnames(sce), metadata[["Source Name"]]) # Enforcing identical order.
stopifnot(all(!is.na(m))) # Checking that nothing's missing.
metadata <- metadata[m,]
head(colnames(metadata))

• 只保留相关的元数据字段，以避免在SingleCellExperiment对象的 colData 中存储不必要的信息。特别是，此处保留每个细胞的原产板 (i.e., block)和表型。第二个字段是有关联的，因为所有细胞都表达oncogene，CBFB-MYH11基因，但是这种oncogene基因只在细胞的一个亚型中诱导表达。

colData(sce)$Plate <- factor(metadata[["Factor Value[block]"]])
pheno <- metadata[["Factor Value[phenotype]"]]
levels(pheno) <- c("induced", "control")
colData(sce)$Oncogene <- pheno
table(colData(sce)$Oncogene, colData(sce)$Plate)

2.3 Incorporating gene-based annotation

• 特征计数工具(Feature-counting)通常根据 Ensembl 或 Entrez 的标准标识符报告基因。这些标识符的使用，因为它们是明确的，高度稳定。然而，与文献中较常用的基因符号( gene symbols )相比，它们很难理解。但是给定Ensembl标识符，方便我们使用注释包annotation packages（例如org.Mm.eg.db.）获得相应的基因符号。

library(org.Mm.eg.db)
symb <- mapIds(org.Mm.eg.db, keys=rownames(sce), keytype="ENSEMBL", column="SYMBOL")
rowData(sce)$ENSEMBL <- rownames(sce)
rowData(sce)$SYMBOL <- symb
head(rowData(sce))

• 下一步就是过滤去重
• 通常需要将 sce 的行名称重命名为基因符号，因为这些符号更易于解释。但是，这需要一些工作来注释缺省和重复的gene symbols。下面的代码将用 Ensembl ID替换NA的gene symbols，并将重复的gene symbols替换为 Ensembl ID。

library(scater)
rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ENSEMBL, rowData(sce)$SYMBOL)
head(rownames(sce))

• 我们还使用TxDb.Mmusculus.UCSC.mm10.ensGene包来确定每个基因的染色体位置。这在以后将非常有用，因为几个质量控制指标将从与线粒体基因对应的行中计算。

BiocManager::install("TxDb.Mmusculus.UCSC.mm10.ensGene")
library(TxDb.Mmusculus.UCSC.mm10.ensGene)
location <- mapIds(TxDb.Mmusculus.UCSC.mm10.ensGene, keys=rowData(sce)$ENSEMBL, 
                   column="CDSCHROM", keytype="GENEID")
rowData(sce)$CHR <- location
summary(location=="chrM")

或者，可以使用scater中的getBMFeatureAnnos函数将来自BioMart资源的注释直接添加到对象中。这可能比上面所示的方法更方便，但取决于 BioMart 数据库的可用性与互联网连接情况，这个网速就很玄妙了。

3 Quality control on the cells

3.1 Defining the quality control metrics

• 需要去除低质量细胞，以确保技术差异不会干扰下游分析结果。我们使用多个质量控制 （QC） 指标：

• 库的大小定义为所有要素（即基因和spike-in 转录本）的计数总和。小库的细胞质量低，一般是因为RNA在制备过程中没有被有效捕获（即未能转化为cDNA和扩增）。
• 每个细胞中的表征数定义为该细胞具有非零计数的表征数。任何表达基因很少的细胞都可能质量较差，因为不同的转录组群尚未被成功捕获。
• 映射到spike-in转录本的读取比例是相对于每个单元格的库大小计算的。高比例表明细胞质量低劣，其中内源RNA在加工过程中丢失（例如，由于细胞裂解或RNA降解）。每个细胞的spike-in RNA量相同，因此spike-in 计数中的富集是内源RNA损失的表现。
• 在没有spike-in转录本的情况下，也可以使用与线粒体基因组中基因映射的读取比例。高比例表明细胞质量低(Islam et al. 2014; Ilicic et al. 2016)，可能是因为穿孔细胞的细胞质RNA的丧失。其理由是线粒体比单个转录分子大，不太可能通过细胞膜中的撕裂而逃脱。

• 对于每个细胞，我们使用scater包（McCarthy 等人 2017 年）中的QCMetrics函数计算这些质量控制指标。这些质量控制指标存储在SingleCellExperiment的行和列元数据中，以供将来参考。
• 这些指标的分布如图1所示，按基因诱导状态和原板分层。目的是去除具有小库、低提取表征数量和高spike-in（或线粒体）比例的低质量细胞。这种细胞会干扰下游分析，例如，形成不同的簇，使结果解释复杂化。

sce$PlateOnco <- paste0(sce$Oncogene, ".", sce$Plate)
multiplot(
    plotColData(sce, y="total_counts", x="PlateOnco"),
    plotColData(sce, y="total_features_by_counts", x="PlateOnco"),
    plotColData(sce, y="pct_counts_ERCC", x="PlateOnco"),
    plotColData(sce, y="pct_counts_Mt", x="PlateOnco"),
    cols=2)

Figure 1Distributions of various QC metrics for all cells in the 416B dataset. This includes the library sizes, number of expressed genes, and proportion of reads mapped to spike-in transcripts or mitochondrial genes.

• 研究质量控制指标如何相互配合也是很有价值的（图2）。通常，它们将大致一致，即总计数低的细胞也将具有低表示特征的数量和高ERCC/线粒体比例。明显的差异可能对应于细胞批次之间的技术差异（见下文）或RNA含量的真正生物学差异。

par(mfrow=c(1,3))
plot(sce$total_features_by_counts, sce$total_counts/1e6, xlab="Number of expressed genes",
    ylab="Library size (millions)")
plot(sce$total_features_by_counts, sce$pct_counts_ERCC, xlab="Number of expressed genes",
    ylab="ERCC proportion (%)")
plot(sce$total_features_by_counts, sce$pct_counts_Mt, xlab="Number of expressed genes",
    ylab="Mitochondrial proportion (%)")

Figure 2: Behaviour of each QC metric compared to the total number of expressed features. Each point represents a cell in the 416B dataset.

3.2 Identifying outliers for each metric

• 为这些指标选择阈值并不简单，因为其绝对值取决于实验流程。例如，无论细胞的质量如何，到更深的测序深度都会导致更多的reads和更多expressed features。同样，在实验中使用更多的spike-in RNA将导致更高的spike-in比例。为了获得自适应阈值，我们假设大多数数据集由高质量细胞组成，并标识出各种 QC 指标异常的细胞。
• 异常值是根据所有细胞中每个指标的中值的绝对偏差 （MAD） 定义的。我们删除日志库大小低于日志库中位数 3 个以上的细胞。对数变换可提高较小值的分辨率，尤其是在原始值的 MAD 与中位数相当或大于中位数时。我们还去除了对数转换的表达基因数量低于中值3 MAD的细胞。

libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", 
    log=TRUE, batch=sce$PlateOnco)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", 
    log=TRUE, batch=sce$PlateOnco)

• batch=参数确保在指定板/基因的每个级别内异常值。这允许isoutlier() 适应不同板的 QC 指标的系统差异（图 1），这可能是由于处理技术差异（例如，测序深度的差异）而不是质量的任何变化而产生的。同样的推理也适用于肿瘤基因诱导状态，其中诱导细胞可能由于生物原因自然表达的基因较少。如果无法解释这些系统差异，就会使MAD估计值膨胀，并损害低质量细胞的去除。
• 我们以类似的方式识别基于比例的指标的异常值。在这里，不需要转换，因为我们正在识别大型异常值，在原始尺度上，区别应该相当清楚。我们不需要使用线粒体比例，因为我们已经具有此数据集的spike-in比例（用于类似目的）。这避免了由于细胞类型之间的线粒体含量的真正差异而产生的潜在问题，这些细胞类型可能会混淆异常值识别。

spike.drop <- isOutlier(sce$pct_counts_ERCC, nmads=3, type="higher",
    batch=sce$PlateOnco)

-按列进行设置将仅保留通过上述每个筛选器的高质量细胞。我们检查每个筛选器移除的细胞数以及保留的细胞总数。去除相当比例的细胞（> 10%）可能表示数据质量存在整体问题。

keep <- !(libsize.drop | feature.drop | spike.drop)
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop),
    BySpike=sum(spike.drop), Remaining=sum(keep))

结果是

  ByLibSize ByFeature BySpike Remaining
1         5         4       6       183

• 然后，我们对SingleCellExperiment对象进行子集化，以仅保留假定的高质量细胞。我们还将原始对象保存到文件中以供以后使用。

sce$PassQC <- keep
saveRDS(sce, file="416B_preQC.rds")
sce <- sce[,keep]
dim(sce)

Comments from Aaron:

• 有关基于异常值的低质量单元检测所依据的假设的更详细讨论，请参阅本节。
• isOutlier()还将返回每个指标的确切筛选器阈值（如果指定了 batch=，则在每个批次中）。这些对于检查自动选择的阈值是否合适可能很有用。

attr(libsize.drop, "thresholds")
attr(spike.drop, "thresholds")

4 Classification of cell cycle phase

• 我们使用Scialdone et al. (2015)描述的预测方法，根据基因表达数据将细胞分类为细胞周期阶段。使用训练数据集，计算了每对基因之间两个基因表达差异的符号。选择在细胞周期阶段与符号变化成对作为标记。然后，测试数据集中的细胞可以分类到适当的阶段，具体取决于每个标记对的观测符号是否与一个相位或另一个的阶段一致。
• 此方法用scran包中cyclone函数实现。该包包含一组预先训练的鼠的标准数据，我们可以在readRDS函数中加载这些数据。我们使用数据集中每个基因的 Ensembl 标识符来匹配标准数据的基因对集中的基因。

set.seed(100)
BiocManager::install("scran")
library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", 
    package="scran"))
assignments <- cyclone(sce, mm.pairs, gene.names=rowData(sce)$ENSEMBL)

• HSC 数据集中每个细胞的cyclone结果如图 3 所示。为每个细胞分配每个阶段的分数，分数越高，与细胞处于该阶段的概率较高对应。我们专注于 G1 和 G2/M 指标，因为这些是最具信息性的分类。

plot(assignments$score$G1, assignments$score$G2M, 
    xlab="G1 score", ylab="G2/M score", pch=16)

Figure 3: Cell cycle phase scores from applying the pair-based classifier on the 416B dataset. Each point represents a cell, plotted according to its scores for G1 and G2/M phases.

• 如果 G1 分数高于 0.5 且大于 G2/M 分数，则细胞被归类为 G1 阶段;G2/M 阶段如果 G2/M 分数高于 0.5 且大于 G1 分数;在S阶段，如果两个分数都高于0.5。在这里，绝大多数细胞被归类为G1阶段。我们将这些归类保存到 SingleCellExperiment对象中以供以后使用。

sce$phases <- assignments$phases
table(sce$phases)

结果如下

> table(sce$phases)

 G1 G2M   S 
 98  62  23 

• scran包中提供预先训练的数据集（Pre-trained classifiers），包括人类和小鼠数据。虽然这里使用的小鼠预先训练的数据集是根据胚胎干细胞的数据训练得来的，但对于其他细胞类型来说仍然准确（Scialdone et al. 2015）。这可能是由于与细胞周期相关的转录程序（Bertoli, Skotheim, and Bruin 2013; Conboy et al. 2007）。基于配对的方法也是一个非参数过程，对数据集之间的大多数技术差异是可靠的。
  Comments from Aaron:

• 为了消除由于细胞周期阶段造成的混淆效应，我们可以筛选细胞，只保留特定阶段（通常为 G1）的细胞，以便进行下游分析。或者，如果不可忽略的细胞数位于其他阶段，我们可以使用分配的阶段作为 blocking factor。这可以防止细胞循环效应，而不会丢弃信息，稍后将详细讨论。
• classifiers对于与训练集中使用的数据有很大不同的数据可能不准确，例如，由于使用不同的protocol。在这种情况下，用户可以使用sandbag函数从自己的训练数据构造自定义classifiers。对于没有预先训练的classifiers的其他模型生物体来说，这也是必要的。
• 在应用cyclone之前，不要过滤掉低丰度基因。即使基因未在任何细胞中表达，如果基因是相位特异性的，它仍然可用于分类。相对于其他基因，它的表达水平虽低但仍然会产生信息，而过滤掉它们会减少权重（reduce power，姑且先这么理解）。

5 Examining gene-level expression metrics

5.1 Inspecting the most highly expressed genes

• 富集分析高表达的基因（图4）。这通常应该由构成表达的转录本，如核糖体或线粒体蛋白的转录本主导。如果其他类别的功能与预期的生物学不一致，则它们的存在可能会引起关注。例如，包含许多spike-in转录本的数据集表明，在库制备过程中添加了过多的spike-in RNA，而没有核糖体蛋白和/或其假基因的存在则表明其排列不理想。

fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
plotHighestExprs(sce, n=50) + fontsize

图 4：分配给 416B 数据集中前 50 个高表达features的总计数百分比。对于每个features，每个条形表示分配给单个细胞的基因的百分比，而圆表示所有细胞的平均值。条形由每个单元格中的features总数着色，而圆形则根据features是否标记为control feature进行着色。

5.2 Filtering out low-abundance genes

• 低丰度基因存在问题，因为零或接近零计数不包含太多用于可靠统计推理的信息。在涉及假设检验的应用中，这些基因通常不能提供足够的证据来拒绝零假设，但它们仍然增加了多重检验校正的严重程度。计数的离散性也可能干扰统计过程，例如，通过影响连续近似的准确性。因此，在应用下游方法之前，许多RNA-seq分析策略中通常会去除低丰度基因。
• 筛选策略的"最优"选择取决于下游应用程序。与移除动力不足的测试所需的分立性相比，通常需要更激进的滤波器来去除离散性（e.g., for normalization）。对于假设检验，筛选器统计量还应独立于零假设下的检验统计量（Bourgon, Gentleman, and Huber 2010）。因此，我们（或相关函数）将根据需要在每个步骤进行筛选，而不是对整个分析应用单个筛选器（single filter）。
• 几个指标可以用来定义低丰度基因。最明显的是每个基因的平均计数，该数在数据集中的所有单元格中计算。我们使用calcAverage()函数计算这一点，该函数还对细胞之间的库大小差异执行一些调整。我们通常观察到中度表达基因的峰值，遵循低表达基因的plateau （从图上看像是一个平原）（图5）。

ave.counts <- calcAverage(sce, use_size_factors=FALSE)
hist(log10(ave.counts), breaks=100, main="", col="grey80", 
    xlab=expression(Log[10]~"average count"))

Figure 5: Histogram of log-average counts for all genes in the 416B dataset.

• 可以对此值应用最小阈值来筛选出低表示的基因。下面的示例演示了我们如何删除平均计数小于 1 的基因。demo.keep 中的 TRUE 值数对应于筛选后保留的行/基因数。

demo.keep <- ave.counts >= 1
filtered.sce <- sce[demo.keep,]
summary(demo.keep)

• 我们还检查表达每个基因的细胞数量。这与大多数基因的平均计数密切相关，因为基因在许多细胞中的表达将导致更高的平均值（图6）。在极少数细胞中表达的基因往往会不被关注，因为它们是由扩增伪影驱动的（尽管它们也可能来自稀有种群）。然后，我们可以去除在少于n个细胞中表达的基因。

num.cells <- nexprs(sce, byrow=TRUE)
smoothScatter(log10(ave.counts), num.cells, ylab="Number of cells", 
    xlab=expression(Log[10]~"average count"))

Figure 6: The number of cells expressing each gene in the 416B dataset, plotted against the log-average count. Intensity of colour corresponds to the number of genes at any given location.

• 如上所述，我们将在每个步骤中应用这些筛选器，而不是通过子设置 sce 全局应用它们。这可确保在每个应用程序中使用最合适的筛选器。尽管如此，我们删除任何细胞中不表达的基因，以减少下游步骤的计算工作。这些基因不提供信息，任何过滤策略都会被移除。

to.keep <- num.cells > 0
sce <- sce[to.keep,]
summary(to.keep)

6 Normalization of cell-specific biases

6.1 Using the deconvolution method to deal with zero counts

• 读取计数取决于细胞之间的捕获效率和测序深度的差异(Stegle, Teichmann, and Marioni 2015)。在下游定量分析之前，需要标准化来消除这些细胞特异性偏差。这通常是通过假设大多数基因在细胞之间没有差异表达（DE）来实现的。假定两个细胞之间非DE大多数基因之间的计数大小的任何系统差异都代表偏差，并通过缩放去除。更具体地说，"大小因子"的计算表示每个库中应缩放计数的程度。
• 大小因子可以用几种不同的方法计算，例如，使用DESeq2包中的estimateSizeFactorsFromMatrix函数 (Anders and Huber 2010; Love, Huber, and Anders 2014)，或用edgeR包的calcNormFactors函数。但是，单细胞数据对于这些基于批量数据的方法来说可能存在问题，因为低计数和零计数占主导地位。为了克服这一点，我们汇集了来自多个单细胞的计数，以增加计数的大小，以便精确估计大小因子(Lun, Bach, and Marioni 2016)。然后，基于池的大小因子被“deconvolved”到基于细胞的因子中，以规范化每个细胞的表达式配置文件。

sce <- computeSumFactors(sce)
summary(sizeFactors(sce))

• 大小因子与所有细胞的库大小密切相关（图 7）。这表明，细胞之间的大多数系统差异是由捕获效率或测序深度的差异引起的。细胞之间的任何 DE 都会在总计数和大小因子之间产生非线性趋势，并且/或增加围绕趋势的散射。我们在基因诱导后观察到一些证据，其中诱导后的大小因子有计划地降低。这与诱导后基因向上调节引入的成分偏差(Robinson and Oshlack 2010)是一致的。

plot(sce$total_counts/1e6, sizeFactors(sce), log="xy",
    xlab="Library size (millions)", ylab="Size factor",
    col=c("red", "black")[sce$Oncogene], pch=16)
legend("bottomright", col=c("red", "black"), pch=16, cex=1.2,
    legend=levels(sce$Oncogene))

到这一步的时候开始出错了，提示legend 为0，我单独运行了等号后面的是有东西的，不知道是不是前面哪里出了问题，我从新从头运行了一次然后出现前面3个图合在了一起，估计是这里的设置问题

test

但是仍然出错，

image.png

Error in legend("bottomright", col = c("red", "black"), pch = 16, cex = 1.2,  : 
  'legend' is of length 0

回头去查一下看看是哪里的问题

Comments from Aaron:

• 虽然deconvolution approach对scRNA-seq数据中的高频为零是可靠的，但如果计数过多为零，它最终将失败。这表现为manifests as negative size factors that are obviously nonsensical。为了避免这种情况，computeSumFactors函数在计算size factors 之前将自动删除低丰度基因。library size-adjusted的平均计数低于指定阈值（最小值）的基因将被忽略。对于读取计数数据，默认值 1 通常令人满意。
• 基于细胞的QC应始终在normalization前执行，以去除表达基因数量极低的细胞。否则，computeSumFactors()可能会为低质量细胞生成size factors。这是因为这些细胞中存在太多的零，从而降低了池以消除零的有效性。有关函数的更多详细信息，请参阅?computeSumFactors。
• sizes参数可用于指定用于计算sizes因子的池sizes数。sizes越多，估计量就越精确，但需要花费一些计算时间和内存。通常，所有sizes应大于 20 个细胞，以确保每个池中有足够的非零表达式值。用于有效池化的细胞总数也应至少为 100。
• 对于高度异构的数据集，建议对细胞执行粗略聚类以削弱非 DE 假设。这可以通过quickCluster()函数完成，并且通过 clusters参数传递给computeSumFactors()函数的结果。我们在下一个工作流中使用更大的数据集演示此方法。

6.2 Computing separate size factors for spike-in transcripts

• 从内源基因计数计算的Size factors通常不适合normalizing spike-in转录本的计数。考虑一个without library quantification，即，在池化和multiplexed sequencing之前，每个库中的 cDNA 量未相等。在这里，含有更多RNA的细胞对内源性基因的计数较大，因此， larger size factors to scale down those counts。然而，在库制备期间，每个细胞中都会添加相同数量的spike-in RNA。这意味着spike-in转录本的计数不受RNA含量的影响。尝试用基于基因的size factors使spike-in计数normalization将导致表达的过度normalization和不正确的量化。在进行库量化的情况下，也适用类似的推理。对于恒定的总量cDNA，内源RNA含量的任何增加都会抑制spike-in转录本的覆盖率。因此，spike-in计数中的偏差将与基于基因的size factors捕获的偏差相反。
• 为了确保正确执行normalization，我们为spike-in set 计算一组单独的size factors。对于每个细胞，特定于峰值的size factors定义为spike-in set中所有转录本的总数。这假定没有一个spike-in转录本本身差异地表示，这是合理的，因为每个细胞中应添加相同量和成分的spike-in RNA(Lun et al. 2017)。（See below for a more detailed discussion on spike-in normalization。）通过设置computeSpikeFactors中的 general.use=FALSE，这些size factors存储在SingleCellExperiment对象的单独字段中。这确保了它们将只与spike-in转录本一起使用，而不是与内源基因一起使用。

6.3 Applying the size factors to normalize gene expression

• count data用于计算normalized的log-expression值，用于下游分析。每个值定义为每个计数与相应细胞的size factors的，在计算前+1以避免零计数时值为空。按size factors划分可确保消除任何特定于细胞的偏差。如果sce中存在spike-in-specific size factors，则它们将自动应用于将spike-in转录本与内源基因分开normalize。

sce <- normalize(sce)

• log-transformation很有用，因为它意味着值中的任何差异都直接表示细胞之间表达差异的值的changes。这通常比绝对的覆盖范围差异更相关，需要根据总体丰度来解释这些差异。对数变换还提供一些方差稳定性度量 (Law et al. 2014)，因此具有较大方差的高丰度基因不会主导下游分析。计算值还存储为"logcounts"矩阵还可以成为用于其他分析的元素。

7 Modelling the technical noise in gene expression

• 基因间观察到的表达值的差异可能由真正的生物学差异或无意义的技术性噪声。为了区分这两个可能性，我们需要对每个基因的表达值方差的技术成分进行建模。我们使用一组spike-in转录本来执行此操作，这些spike-in转录本以相同的数量添加到每个细胞样品中。因此，spike-in转录本应没有生物变异性，即其计数的任何差异都应是技术来源。
• 我们使用trendVar()函数将均值相关趋势与spike-in转录本的对数表达式值的方差拟合。我们设置block=，以确保板之间的技术差异不会膨胀方差。鉴于一个基因的平均丰度，该趋势的拟合值将用作该基因技术成分的估计值。方差的生物分量最终通过decomposeVar函数的每个基因的总方差中减去技术分量来计算。

var.fit <- trendVar(sce, parametric=TRUE, block=sce$Plate,
    loess.args=list(span=0.3))
var.out <- decomposeVar(sce, var.fit)
head(var.out)

• 我们通过可视化检查趋势，确认它对应于峰值方差（图 8）。波浪状形状是log-expression值的平均方差趋势的典型。当平均值从零增加时，方差的线性增加可以被观察到，因为当计数增加时，大方差是可能的。在非常高的丰度下，采样噪声的影响会由于大量数定律而减小，从而导致方差减小。

plot(var.out$mean, var.out$total, pch=16, cex=0.6, xlab="Mean log-expression", 
    ylab="Variance of log-expression")
curve(var.fit$trend(x), col="dodgerblue", lwd=2, add=TRUE)
cur.spike <- isSpike(sce)
points(var.out$mean[cur.spike], var.out$total[cur.spike], col="red", pch=16)

Figure 8: Variance of normalized log-expression values for each gene in the 416B dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the spike-in transcripts (red).

• 我们检查具有最大生物成分的基因的表达值分布。这可确保方差估计不由一两个异常值的细胞驱动（图 9）。(这个小提琴图还是挺漂亮的像几个穿旗袍的美女，哈哈)

chosen.genes <- order(var.out$bio, decreasing=TRUE)[1:10]
plotExpression(sce, features=rownames(var.out)[chosen.genes]) + fontsize

Figure 9: Violin plots of normalized log-expression values for the top 10 genes with the largest biological components in the 416B dataset. Each point represents the log-expression value in a single cell.

• 在实践中，趋势拟合由于记录量少和丰度分布不均而变得复杂。有关如何优化拟合的更多详细信息，请参阅此处。
• 在没有spike-in的情况下，用户可以设置use.spikes=FALSE，以适合内源基因方差的趋势（参见此处）。或者，我们可以基于泊森技术噪声的假设创建趋势，如此处所述。
• Negative biological components通常来自decomposeVar。这些在直觉上是毫无意义的，因为基因不可能在技术噪声下具有总变异。尽管如此，这些值的发生是由于对总方差的估算不精确，尤其是对于数量少的细胞。
• decomposeVar还生成 p 值，可用于在错误发现率 （FDR） 的特定阈值下定义 HVG。稍后我们将更详细地讨论这一点，因为在数据探索期间，特征选择不需要 HVG 的正式检测。

8 Removing the batch effect

• 如前所述，数据是在两个板块上收集的。板间处理的微小无法控制的差异可能导致批次效应，即不同板上细胞之间的表达系统差异。这种差异并不有趣，可以通过从 limma 包中应用removeBatchEffect()函数(Ritchie et al. 2015)来删除。这消除了原板的影响，同时考虑了对子基因诱导的（感兴趣的）影响。

library(limma)
assay(sce, "corrected") <- removeBatchEffect(logcounts(sce), 
    design=model.matrix(~sce$Oncogene), batch=sce$Plate)
assayNames(sce)

• 对于非基于模型的下游过程（例如聚类和大多数尺寸减少形式），需要手动批处理校正。但是，如果分析方法可以接受设计矩阵，则最好在设计矩阵中阻止滋扰因子，而使用removeBatchEffect()。这是因为后者没有考虑到残余自由度的丧失，也没有考虑到估计blocking factor terms的不确定性。

Comments from Aaron:

• removeBatchEffect()执行线性回归，并将与blocking factors对应的系数设置为零。如果每个批次中的人口组成已知（设置design=）或跨批次相同，则此功能有效。此数据集的这种假设是合理的，涉及两个板上的均质细胞系总体。然而，在大多数scRNA-seq应用中，不同批次的变化因素并不相同，而且事先不知道。这促使使用更复杂的批处理校正方法，如mnnCorrect()。

9 Denoising expression values using PCA

• 一旦对技术噪声进行建模，我们就可以使用主要元件分析 （PCA） 来消除随机技术噪声。假设每个单元格表示高维表达式空间中的点，其中点的分布表示总方差。PCA 标识此空间中的坐标轴，这些轴可捕获尽可能多的此方差。每个轴都是一个主组件 （PC），其中任何早期的 PC 都会比后来的 PC 解释更多的差异。
• 我们假设，涉及共同调控的基因组的生物过程将占数据差异最大的部分。如果是这种情况，此过程应由一个或多个较早的 PC 表示。相反，随机技术噪声影响每个基因独立，将由以后的PC表示。denoisePCA()函数将删除后来的PC，直到丢弃的总方差等于PCA中使用的所有基因的技术成分之和。
• 该函数返回一个SingleCellExperiment对象，其中包含reducedDims插槽中每个单元格的 PC 分数。其目的是消除技术噪音，并丰富保留 PC 中的生物信号。这提高了下游过程（如聚类）中基础生物学的分辨率。

Comments from Aaron:

• denoisePCA()将只使用具有正生物成分的基因，即方差大于拟合趋势。这可确保要丢弃的总技术差异不会大于数据中的总方差。
• 对于technical=参数，该函数还将直接接受趋势函数（即var.fit$trend）或每个基因的技术组件的矢量。在这里，我们提供来自decomposeVar()的DataFrame，以允许函数在批处理校正后调整剩余自由度的损失。具体来说，批处理矩阵中的方差略低，需要对技术组件进行一些重新缩放以进行补偿。
• 此处不对丰度执行滤波，这可确保仍可检测到与罕见亚种群对应的PC。离散性问题较少，因为低丰度基因的方差也较低，从而降低了它们对 PCA 的贡献。
• 还可以获取原始表达式矩阵的低等级近似值，捕获与保留的 PC 等效的方差。这对于在需要基因表达值的下游过程之前进行去角非常有用。

10 Visualizing data in low-dimensional space

10.1 With PCA

• 我们通过构建前三个主成分的 PCA 图来可视化细胞之间的关系（图 10）。具有相似表达式配置文件的细胞应位于图中，而不同的细胞应相距很远。在这种情况下，我们观察到基于基因诱导状态的细胞的清晰分离，与转录组的预期效果一致。

plotReducedDim(sce, use_dimred="PCA", ncomponents=3, 
    colour_by="Oncogene") + fontsize

Figure 10: Pairwise PCA plots of the first three PCs in the 416B dataset, constructed from normalized log-expression values of genes with positive biological components. Each point represents a cell, coloured according to oncogene induction status.

-相比之下，我们观察到细胞没有按批次明确分离（图11）。这表明使用removeBatchEffect()的批处理更正步骤已成功。

plotReducedDim(sce, use_dimred="PCA", ncomponents=3, 
    colour_by="Plate") + fontsize

Figure 11: Pairwise PCA plots of the first three PCs in the 416B dataset, constructed from normalized log-expression values of genes with positive biological components. Each point represents a cell, coloured according to the plate of origin.

• 请注意，plotReducedDim()将使用已通过denoisePCA()存储在sce中的PCA结果。这使我们能够快速生成具有不同美学的新绘图，而无需重复整个 PCA 计算。同样，如果plotPCA()在SingleCellExperiment中可用，将使用现有结果，否则将重新计算它们。用户应设置rerun=TRUE，以在现有结果存在的情况下强制重新计算 PCs。

Comments from Aaron:

• 对于每种可视化方法，可以将其他特定于单元格的信息合并到每个点的大小或形状中。这是使用大多数绘图函数中的size_by=和shape_by=参数来完成的。
• 可以显示更多分量，但这些组件的信息量通常较少，因为它们解读的方差较少。它们通常也更难解读，因为它们被定义为与早期 PCs 正交（因此取决于在这些 PCs 中检测到的内容）。

10.2With t-SNE

• 另一种广泛使用的降维方法是t-stochastic neighbour embedding (t-SNE)方法(Van der Maaten and Hinton 2008)。t-SNE往往比PCA更好地在更多样化的种群分离细胞。这是因为前者可以直接捕获高维空间中的非线性关系，而后者必须在线性轴上表示非线性关系。但是，这种改进的代价是付出更多的计算工作量，并且要求用户考虑random seed and perplexity等参数（请参阅comments）。
• 我们使用plotTSNE()函数演示了图12中t-SNE图的生成。我们设置use_dimred="PCA"来对数据的低级近似执行t-SNE，允许算法利用前面的降噪步骤。

set.seed(100)
out5 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=5),
    colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 5")

set.seed(100)
out10 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=10),
    colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 10")

set.seed(100)
out20 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=20),
    colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 20")

multiplot(out5, out10, out20, cols=3)

又出现跟上次一样点图没有电的情况的figure12

• t-SNE 是一种随机方法，因此用户应多次运行该算法，以确保结果具有代表性。脚本应设置种子（通过set.seed()命令），以确保所选结果可重现。还建议测试"perplexity"参数的不同设置，因为这将影响低维空间中点的分布。
• 在这里，我们调用 20 的perplexity的runTSNE()， 将 t-SNE 结果存储在我们的SingleCellExperiment对象中。这样可以避免在我们想要使用plotTSNE()创建新绘图时重复计算，因为将改用存储的结果。同样，用户可以将rerun=TRUE设置为在存在存储的结果时强制重新计算。
• 这里没有考虑许多其他降维策略，但也可以使用，例如多维缩放、扩散图。它们各有优缺点，例如，扩散图（参见plotDiffusionMap）将细胞沿连续轨迹放置，并适合可视化分化等分级过程(Angerer et al. 2016)。

Comments from Aaron:

• 有关如何解释 t-SNE 绘图的良指南可在http://distill.pub/2016/misread-tsne/中找到。这演示了二维嵌入中的聚类之间的距离如何没有多大意义，并且群集的明显"size"（即spread）也没有什么意义。

11 Clustering cells into putative subpopulations

11.1Defining cell clusters from expression data

• 去噪声的log-expression值用于将细胞集中到假定的亚群中。具体而言，我们使用 Ward’s criterion来最小化每个群集内的总方差，对单元之间的欧氏距离执行分层聚类。这会产生一个树状图，将具有相似表达模式的细胞组合在一起，在选定的基因上。

pcs <- reducedDim(sce, "PCA")
my.dist <- dist(pcs)
my.tree <- hclust(my.dist, method="ward.D2")

• 聚类是通过将dynamic tree cut (Langfelder, Zhang, and Horvath 2008) 来显式定义的。这利用树突图中的分支形状来优化群集定义，并且比cutree更适合复杂树形图。通过在cutreeDynamic中手动指定cutHeight，可以更好地控制经验聚类。我们还将minClusterSize设置为低于默认值 20 的值，以避免远距离小群集的虚假聚合。

library(dynamicTreeCut)
my.clusters <- unname(cutreeDynamic(my.tree, distM=as.matrix(my.dist), 
    minClusterSize=10, verbose=0))

• 我们根据已知因素检查每个簇中的细胞分布。每个群集由两个批处理的单元组成，表示聚类不是由批处理效应驱动的。观察到每个簇的组成与Oncogene有关的差异，与基因诱导的生物效应一致。

table(my.clusters, sce$Plate)
table(my.clusters, sce$Oncogene)

• 我们可视化图 13 中 t-SNE 图上所有细胞的聚类分配。相邻细胞通常分配给同一群集，表示群集过程已正确应用。

sce$cluster <- factor(my.clusters)
plotTSNE(sce, colour_by="cluster") + fontsize

Figure 13: t-SNE plot of the denoised PCs of the 416B dataset. Each point represents a cell and is coloured according to the cluster identity to which it was assigned.

• 我们使用轮廓宽度检查聚类的分离度（图 14）。具有较大的正轮廓宽度的单元格比不同群集中的单元格更接近同一群集中的其他单元格。相反，宽度为负的单元格比分配给它的群集中的其他单元格更接近其他群集。理想情况下，每个群集将包含许多具有较大正宽度的单元格，这表明它与其他群集很好地分离。

  
  
  Figure 14: Barplot of silhouette widths for cells in each cluster. Each cluster is assigned a colour and cells with positive widths are coloured according to the colour of its assigned cluster. Any cell with a negative width is coloured according to the colour of the cluster that it is closest to. The average width for all cells in each cluster is shown, along with the average width for all cells in the dataset.
• 轮廓宽度可用于确定参数值，以最大化群集之间的分隔。例如，我们可以在 cutreeDynamic中更改剪切高度或分割深度，以最大化所有单元格的平均轮廓宽度。这通常为进一步检查提供令人满意的初始聚类。但是，请记住，聚类的粒度与显微镜上的放大倍率非常类似。不同的数据视图可以使用不同的粒度获得，其中一些在分离测量上可能不够理想。如果群集不能为特定生物问题提供所需的粒度，则用户不应将聚类与最大分离性集中在一起。
• 图 14 中大多数单元的轮廓正宽度相对较小，表明簇之间的分离较弱。这可能是过度聚类的症状，其中在子基因归纳状态上明确定义的聚类进一步被分割成分离得不太良好的子集。尽管如此，我们将在my.cluster中继续当前聚类方案，因为它为进一步描述异质性提供了合理的分区。

Comments from Aaron:

• 另一种聚类策略是使用从相关关系派生的距离矩阵（例如，如quickCluster）。这对噪声和规范化错误更加可靠，但对表达式配置文件中的细微变化也不太敏感。
• 沃德的标准和完整的连杆产生球形，紧凑的集群。特别是，完全联系有利于形成直径相同的簇。在某些情况下，这可能是可取的，但当亚种群的方差不同时，这不太合适。因此，我们通常使用沃德的标准进行初始聚类。当然，尝试其他方法（建议）很简单（建议），前提是执行一些评估，例如，使用轮廓宽度。

11.2Detecting marker genes between clusters

• 一旦通过聚类识别假定的亚群体，我们就可以使用findMarkers函数识别每个聚类的标记基因。这将在每个基因的对数表达值和每对簇之间执行Welch t-tests(Soneson and Robinson 2018)。目的是测试每个群集中的 DE 与其他群集相比，同时阻止无意义因素，如原板差异（详情请参阅此处）。top DE基因可能是很好的候选标记，因为它们可以有效地区分不同簇中的细胞。

markers <- findMarkers(sce, my.clusters, block=sce$Plate)

• 对于每个群集，相关比较的 DE 结果合并到单个输出表中。这允许通过从每个集群之间的成对比较中获取Top DE 基因来轻松定义一组标记基因。例如，要从每个比较的前 10 个基因中构造群集 1 的标记集，将筛选marker.set 保留 Top小于或等于 10 的行。还报告了每个基因的其他统计信息，包括调整后的 p 值（见下文）和相对于其他每个群集的对数折叠变化。

marker.set <- markers[["1"]]
head(marker.set, 10)

• 我们保存候选标记基因列表以供进一步检查。

write.table(marker.set, file="416B_marker_1.tsv", sep="\t", 
    quote=FALSE, col.names=NA)

• 我们可视化top候选项的表达式配置文件，以验证 DE 签名是否健壮（图 15）。与部分或所有其他群集相比，大多数top 标记在cluster1 的单元中具有强且一致的向上或向下调节。粗略地检查热图表明，cluster1包含的肿瘤基因诱导细胞，DNA复制和细胞周期基因有很强的低调节。这与衰老作为抗肿瘤反应的潜在诱导是一致的(Wajapeyee et al. 2010)。通过基因集浓缩分析，例如，使用来自limma的kegga或goana，可以对这些标记的功能进行更全面的调查。

top.markers <- rownames(marker.set)[marker.set$Top <= 10]
plotHeatmap(sce, features=top.markers, columns=order(sce$cluster), 
            colour_columns_by=c("cluster", "Plate", "Oncogene"),
            cluster_cols=FALSE, center=TRUE, symmetric=TRUE, zlim=c(-5, 5), show_colnames = F) 

Figure 15: Heatmap of mean-centred and normalized log-expression values for the top set of markers for cluster 1 in the 416B dataset. Column colours represent the cluster to which each cell is assigned, the plate of origin or the oncogene induction status of each cell, as indicated by the legend.

• 图 15 中的许多标记在所选群集中不是唯一向上或向下调节的。对unique DE的测试往往过于严格，因为它忽略了在两个或多个簇中表达的重要基因。例如，在仅CD4+、仅CD8+、双阳性和双阴性T细胞的混合种群中，Cd4或Cd8都不会被检测为亚群特异性标记，因为每个基因在两个亚种群中表达。通过我们的方法，这两个基因将被拾取为候选标记，因为它们将在至少一对亚种群之间被认为是DE。然后可以选择标记的组合来描述亚群，这比尝试找到unique的DE基因更灵活。
• 我们强烈建议选择一些标记用于具有独立复制细胞群的验证研究。其目的是确定表达上调节标记但不表示向下调节标记的细胞的相应子集。理想情况下，在验证过程中也会使用不同的表达量化技术，例如荧光原位杂交或定量PCR。这证实了亚群确实存在，并且不仅仅scRNA-seq或计算分析的结果。

Comments from Aaron:

• 通过设置direction="up", findMarkers将只返回每个cluster中与其他集群相比被调控的基因。这在高度异质的人群中非常方便，可以专注于能够立即识别每个cluster的基因。虽然未表达基因可能也有信息，但它对positive identification用处不大。
• findMarker也可以定向查找所选聚类和所有其他群集之间的 DE 基因。这应该通过设置pval.type="all"来完成，该值将每个基因的 p 值定义为涉及所选群集的所有对对比较中的最大值。结合direction="up"，这可用于标识每个群集的唯一标记。但是，这对过度聚类很敏感，因为如果将聚类拆分为两个较小的子群集，则唯一标记基因将不再存在。
• 必须强调，此处计算的（adjusted）p 值不能正确解释为重要度量。这是因为群集是从数据中根据经验识别的，请参阅此处。

12 Concluding remarks

• 基本分析完成后，将SingleCellExperiment对象保存到本地使用saveRDS 函数存档通常很有用。然后，可以使用readRDS函数轻松地将该对象还原到新的 R 任务中。这允许进行进一步的工作，而无需重复上述所有处理步骤。

saveRDS(file="416B_data.rds", sce)

• 有多种方法可用于对已处理的表达式数据执行更复杂的分析。例如，细胞可以在伪时间（例如，沿着分化途径的进度）与monocle (Trapnell et al. 2014)或TSCAN(Ji and Ji 2016)排序;细胞状态层次结构可以用sincell package (Julia, Telenti, and Rausell 2015)来表示;和振荡行为可以用Oscope (Leng et al. 2015)。HVGs 可用于基因集扩充分析，以识别具有异构活性的生物通路和过程，使用专为 topGO 等批量数据设计的包，或使用专用的 scde (Fan et al. 2016)。这些分析的完整描述超出了此工作流的范围，因此建议感兴趣的用户查阅相关文档。

All software packages used in this workflow are publicly available from the Comprehensive R Archive Network (https://cran.r-project.org) or the Bioconductor project (http://bioconductor.org). The specific version numbers of the packages used are shown below, along with the version of the R installation.

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cluster_2.1.0                          dynamicTreeCut_1.63-1                 
 [3] limma_3.40.6                           scran_1.12.1                          
 [5] TxDb.Mmusculus.UCSC.mm10.ensGene_3.4.0 GenomicFeatures_1.36.4                
 [7] scater_1.12.2                          ggplot2_3.2.1                         
 [9] org.Mm.eg.db_3.6.0                     AnnotationDbi_1.46.0                  
[11] SingleCellExperiment_1.6.0             SummarizedExperiment_1.14.0           
[13] DelayedArray_0.10.0                    BiocParallel_1.18.0                   
[15] matrixStats_0.54.0                     Biobase_2.44.0                        
[17] GenomicRanges_1.36.0                   GenomeInfoDb_1.20.0                   
[19] IRanges_2.18.1                         S4Vectors_0.22.0                      
[21] BiocGenerics_0.30.0                    BiocFileCache_1.8.0                   
[23] dbplyr_1.4.2                           Rsubread_1.34.7                       

loaded via a namespace (and not attached):
 [1] bitops_1.0-6             bit64_0.9-7              RColorBrewer_1.1-2       progress_1.2.2          
 [5] httr_1.4.1               tools_3.6.1              backports_1.1.4          R6_2.4.0                
 [9] irlba_2.3.3              KernSmooth_2.23-15       vipor_0.4.5              DBI_1.0.0               
[13] lazyeval_0.2.2           colorspace_1.4-1         withr_2.1.2              tidyselect_0.2.5        
[17] gridExtra_2.3            prettyunits_1.0.2        bit_1.1-14               curl_4.0                
[21] compiler_3.6.1           BiocNeighbors_1.2.0      rtracklayer_1.44.4       labeling_0.3            
[25] scales_1.0.0             rappdirs_0.3.1           stringr_1.4.0            digest_0.6.20           
[29] Rsamtools_2.0.0          XVector_0.24.0           pkgconfig_2.0.2          rlang_0.4.0             
[33] rstudioapi_0.10          RSQLite_2.1.2            DelayedMatrixStats_1.6.0 dplyr_0.8.3             
[37] RCurl_1.95-4.12          magrittr_1.5             BiocSingular_1.0.0       GenomeInfoDbData_1.2.1  
[41] Matrix_1.2-17            Rcpp_1.0.2               ggbeeswarm_0.6.0         munsell_0.5.0           
[45] viridis_0.5.1            stringi_1.4.3            edgeR_3.26.8             zlibbioc_1.30.0         
[49] Rtsne_0.15               plyr_1.8.4               grid_3.6.1               blob_1.2.0              
[53] dqrng_0.2.1              crayon_1.3.4             lattice_0.20-38          Biostrings_2.52.0       
[57] cowplot_1.0.0            hms_0.5.1                locfit_1.5-9.1           zeallot_0.1.0           
[61] pillar_1.4.2             igraph_1.2.4.1           reshape2_1.4.3           biomaRt_2.40.4          
[65] XML_3.98-1.20            glue_1.3.1               packrat_0.5.0            BiocManager_1.30.4      
[69] vctrs_0.2.0              gtable_0.3.0             purrr_0.3.2              assertthat_0.2.1        
[73] rsvd_1.0.2               viridisLite_0.3.0        pheatmap_1.0.12          tibble_2.1.3            
[77] GenomicAlignments_1.20.1 beeswarm_0.2.3           memoise_1.1.0            statmod_1.4.32