哈佛大学单细胞课程|笔记汇总 （四）

（四）Single-cell RNA-seq: Quality control

在QC步骤中，我们的目标和挑战主要包括：

目标：

• 筛选数据以仅包括高质量的真实细胞，以便在对细胞进行聚类时，更容易识别不同的细胞类型。

• 识别任何失败的样本并尝试挽救数据或从分析中删除，此外还试图了解样本失败的原因。

挑战：

• 从不太复杂的细胞（不太复杂只细胞中转录本种类构成）中划定质量较差的细胞。

• 选择合适的阈值进行过滤，保留高质量的细胞而不去除生物学相关的细胞类型。

建议：

在执行质量控制之前，需要对细胞类型进行一个预估。例如，您是否需要样品中的复杂性较低的细胞或线粒体表达水平较高的细胞？如果是，那么在评估我们的数据质量时，我们需要考虑这种生物学特征。

生成质控指标：

除了我们上面提到过的orig.ident、nCount_RNA和nFeature_RNA外，我们还可以计算其他QC指标如number of genes detected per UMI（能反映数据的复杂度，UMI数越大数据复杂度越高）和mitochondrial ratio（可以得知来源于线粒体基因的reads）。

#将每个细胞的UMI数量进行log10转换并加入到metadata
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA)/log10(merged_seurat$nCount_RNA)

#计算线粒体相关基因比例，注意！("^MT-")只限于人哦。。。，其它物种自己寻找线粒体基因或其模式
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
merged_seurat$mitoRatio <- [email protected]$mitoRatio / 100

当然，我们也可以使用 $ 添加其他质控指标到metadata中（R变量索引 - 什么时候使用 @或$），这样不会影响merged_seurat对象。我们可以直接从Seurat对象中提取meta.dataslot 创建metadata：

# Create metadata dataframe
metadata <- [email protected]

你应该看到每个细胞ID都有一个ctrl_或stim_前缀，这是我们合并Seurat对象时指定的，它们应与orig.ident中列出的样本匹配。我们可以先添加带有细胞ID的列，并更改当前列名称以使其更加直观：

# Add cell IDs to metadata
metadata$cells <- rownames(metadata)

# Rename columns
metadata <- metadata %>%
        dplyr::rename(seq_folder = orig.ident,
                      nUMI = nCount_RNA,
                      nGene = nFeature_RNA)

基于细胞前缀得到每个细胞的样本名。给sample这一列进行命名：

# Create sample column
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^ctrl_"))] <- "ctrl"
metadata$sample[which(str_detect(metadata$cells, "^stim_"))] <- "stim"

最终的metadata，包含了用来质控的指标:

# 将新的metadata重新添加回seurat
[email protected] <- metadata

#保存
save(merged_seurat, file="data/merged_filtered_seurat.RData")

评估质量指标

我们将对多个质量指标进行联合筛选并去除不合格的细胞，这些指标包括：

• 细胞计数（Cell counts）

• 每个细胞的UMI（UMI counts per cell）

• 每个细胞观察到的基因（Genes detected per cell）

• UMIs vs. genes detected

• 线粒体counts比例（Mitochondrial counts ratio）

• Novelty

什么是doublets？简单的说就是两个细胞混在一起，可能发生在细胞捕获过程中，并且可能会误导认为是两种细胞类型的过渡态（transitory states），所以应该去除（单细胞预测Doublets软件包汇总-过渡态细胞是真的吗？）。

我们为什么不检查doublets呢？许多的工作流程都是通过设置UMI或genes的最大阈值进行的，其原理为大量的reads或基因表明存在着多个细胞，尽管此理由似乎很直观，但并不准确。Scrublet (https://github.com/AllonKleinLab/scrublet)是检测doublets的重要工具，但是我们还没有对它进行基准测试。我们不建议使用阈值进行筛选，*当我们确定了每个cluster的marker后，我们建议探索这些标记，以确定这些marker是否适用于一种以上的细胞类型*。

细胞计数

细胞计数由检测到的barcode的数量确定。对于此实验，预计将有12,000-13,000个细胞。inDrops的细胞捕获效率略高（70-80％），而10X则为50-60％。注意：如果用于文库制备的细胞浓度不准确，捕获效率可能会低得多。我们不应使用FACS或Bioanalyzer确定细胞浓度（这些工具无法准确确定浓度），而应使用血细胞计数器或自动细胞计数器来计算细胞浓度。

同时在inDrops和10X中，也会发生得到的细胞数目（即细胞barcodes数）大于实际细胞数的情况。

# Visualize the number of cell counts per sample
metadata %>%
      ggplot(aes(x=sample, fill=sample)) +
      geom_bar() +
      theme_classic() +
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
      theme(plot.title = element_text(hjust=0.5, face="bold")) +
      ggtitle("NCells")

从上图可以看出每个sample的细胞数均在15000以上。

每个细胞的UMI

每个细胞的UMI计数通常应高于500，500是预期的下限。如果UMI在500-1000计数之间，虽然可以使用，但可能应该对细胞进行更深的测序。

# Visualize the number UMIs/transcripts per cell
metadata %>%
      ggplot(aes(color=sample, x=nUMI, fill= sample)) +
      geom_density(alpha = 0.2) +
      scale_x_log10() +
      theme_classic() +
      ylab("Cell density") +
      geom_vline(xintercept = 500)

我们可以看到大多数细胞的UMI数均在1000以上。

每个细胞观察到的基因

正常情况应该和UMI相同，也是有一个大大的峰，如果我们在主峰的左侧也看到了一个小峰，出现的原因也是多种多样的，可能是由于捕获单个细胞失败，或者具有某些生物学意义，需要仔细评估。

对于高质量数据，基因数目分布直方图应包含一个大峰，该峰代表完整的细胞。如果我们在主要峰的左侧看到一个小峰（我们使用的数据中不存在这种情况），或者细胞的双峰分布，则可能表明有两件事。一个可能是某些原因导致一组细胞失败。一个可能是存在生物学上不同类型的细胞（例如，静止细胞群体，不太复杂的目标细胞），和/或一种类型比另一种类型小得多（即，具有高计数的细胞可能是尺寸较大的细胞） ）的情况。

因此应该使用本课中介绍的其他指标来评估此阈值。

# Visualize the distribution of genes detected per cell via histogram
metadata %>%
      ggplot(aes(color=sample, x=nGene, fill= sample)) +
      geom_density(alpha = 0.2) +
      theme_classic() +
      scale_x_log10() +
      geom_vline(xintercept = 300)

# Visualize the distribution of genes detected per cell via boxplot
metadata %>%
      ggplot(aes(x=sample, y=log10(nGene), fill=sample)) +
      geom_boxplot() +
      theme_classic() +
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
      theme(plot.title = element_text(hjust=0.5, face="bold")) +
      ggtitle("NCells vs NGenes")

UMIs vs. genes

通常我们会一起评估UMI的数量和每个细胞检测到的基因数量。下图整合了线粒体reads、基因数和UMI数之间的关系，显示出线粒体reads比例在低基因数和低总reads数细胞中很高。

为什么我们需要去除线粒体基因占比较高的细胞呢，主要是由于细胞破损的胞质中的mRNA会由于细胞破损游离出来，而线粒体中的被保留了下来，导致线粒体基因比例高。这样的细胞为低质量细胞，需要通过count数、基因数、线粒体比例联合去除。

质量较差的细胞的基因和UMI可能很低，对应于该图左下象限中的数据点。好的细胞通常会同时显示更多的基因和更高数量的UMI。

通过此图，我们还可以评估线的斜率 （斜率在1左右最符合预期，斜率过高说明细胞中基因数多需要加大测序深度，斜率低说明细胞的复杂性低或捕获的基因少）以及该图右下象限中数据点的分布。具有大量的UMI，但只有少数基因的细胞可能是垂死的细胞，但也可能是低复杂度的细胞类型（即红细胞）。

# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
metadata %>%
      ggplot(aes(x=nUMI, y=nGene, color=mitoRatio)) +
      geom_point() +
    scale_colour_gradient(low = "gray90", high = "black") +
      stat_smooth(method=lm) +
      scale_x_log10() +
      scale_y_log10() +
      theme_classic() +
      geom_vline(xintercept = 500) +
      geom_hline(yintercept = 250) +
      facet_wrap(~sample)

线粒体counts比例

通常死细胞或垂死细胞存在大量线粒体污染。除非必要 （如癌细胞或心肌细胞，线粒体比例自然高），我们将线粒体比率超过0.2的细胞定义为线粒体计数质量差的样品 (对一篇单细胞RNA综述的评述：细胞和基因质控参数的选择)。

# Visualize the distribution of mitochondrial gene expression detected per cell
metadata %>%
      ggplot(aes(color=sample, x=mitoRatio, fill=sample)) +
      geom_density(alpha = 0.2) +
      scale_x_log10() +
      theme_classic() +
      geom_vline(xintercept = 0.2)

转录本复杂度

虽然每个细胞都是低深度测序，但其转录本组成具有高度的复杂性。样品中的异常细胞也可能是其RNA复杂度较低，我们可以通过该指标检测筛选低复杂度细胞（如红细胞）的污染。一般情况下，我们把评分（novelty score 每个UMI对应的平均基因数）设置为0.8。

# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI
metadata %>%
      ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
      geom_density(alpha = 0.2) +
      theme_classic() +
      geom_vline(xintercept = 0.8)

NOTE: Reads per cell is another metric that can be useful to explore; however, the workflow used would need to save this information to assess. Generally, with this metric you hope to see all of the samples with peaks in relatively the same location between 10,000 and 100,000 reads per cell.

过滤

线粒体数相对高的细胞也可能参与了呼吸过程，是某研究中正需要的。因此作者认为在设置条件时应综合考虑，尽量宽松，防止有意义的生物学变化被删除。

细胞层面过滤

• nUMI > 500

• nGene > 250

• log10GenesPerUMI > 0.8

• mitoRatio < 0.2

用subset()函数对Seurat对象进行设置：

# Filter out low quality reads using selected thresholds - these will change with experiment
filtered_seurat <- subset(x = merged_seurat,
                         subset= (nUMI >= 500) &
                           (nGene >= 250) &
                           (log10GenesPerUMI > 0.80) &
                           (mitoRatio < 0.20))

基因层面过滤

删除0表达值的基因和在少于10个细胞中表达的基因。

# Output a logical vector for every gene on whether the more than zero counts per cell
# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")

# Output a logical vector for every gene on whether the more than zero counts per cell
nonzero <- counts > 0

# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10

# Only keeping those genes expressed in more than 10 cells
filtered_counts <- counts[keep_genes, ]

# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = [email protected])

Re-assess QC metrics

过滤之后可以回顾一下过滤前后的效果，以评判过滤指标是否合适。

Saving filtered cells

# Create .RData object to load at any time
save(filtered_seurat, file="data/seurat_filtered.RData")