學(xué)習(xí)目標：

• 構(gòu)建 QC 指標和相關(guān)圖便于直觀地查看數(shù)據(jù)質(zhì)量
• 評估 QC 指標并設(shè)置過濾條件去掉低質(zhì)量細胞

單細胞 RNA-seq：質(zhì)量控制

image

此工作流程的每一步都有自己的目的和挑戰(zhàn)钉蒲。對于我們原始數(shù)據(jù)的 QC谦纱，它們包括：

目標：

• 篩選數(shù)據(jù)篙梢，使其僅包含高質(zhì)量的真實細胞，這樣當我們對細胞進行聚類時买雾，就更容易識別不同的細胞類群
• 識別任何不合格的樣本，并嘗試挽救數(shù)據(jù)或?qū)⑵鋸姆治鲋袆h除伤极，此外谎懦，還要嘗試了解樣本失敗的原因

挑戰(zhàn)：

• 從不太復(fù)雜的細胞中劃定質(zhì)量差的細胞
• 選擇合適的過濾閾值，保留高質(zhì)量細胞签赃，而不去除生物學(xué)相關(guān)細胞類型

建議：

• 在執(zhí)行QC之前谷异，要求您對細胞類型的期望有一個很好的了解。例如锦聊，您是否希望在您的樣本中有低復(fù)雜性的細胞或具有較高水平的線粒體表達的細胞歹嘹？如果是這樣，那么在評估我們的數(shù)據(jù)質(zhì)量時孔庭，我們需要考慮到這種生物學(xué)因素尺上。

生成質(zhì)控標準

將數(shù)據(jù)加載到 Seurat 并創(chuàng)建初始對象時，會為矩陣中的每個細胞組裝一些基本metadata史飞。要仔細查看此metadata尖昏，讓我們查看存儲在seurat對象meta.data數(shù)據(jù)：

# Explore merged metadata
View(merged_seurat@meta.data)

image

共有三列信息：

• orig.ident：通常包含樣本標識（已知）。默認是project在加載數(shù)據(jù)時提供的值
• nCount_RNA：此列表示每個細胞的 UMI 數(shù)量
• nFeature_RNA：這一列代表每個細胞檢測到的基因數(shù)量

為了為質(zhì)量控制分析創(chuàng)建適當?shù)膱D构资，我們需要計算一些額外的指標抽诉。這些包括：

• 每個 UMI 檢測到的基因數(shù)量：這個指標讓我們了解數(shù)據(jù)集的復(fù)雜性（每個 UMI 檢測到的基因越多，我們的數(shù)據(jù)就越復(fù)雜）
• 線粒體比率：該指標將為我們提供細胞中線粒體基因的read的百分比

每個細胞的每個UMI的基因數(shù)量非常容易計算吐绵，我們將對結(jié)果進行l(wèi)og10轉(zhuǎn)換迹淌，以便更好地在樣本之間進行比較。

# Add number of genes per UMI for each cell to metadata
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)

線粒體比率

Seurat有一個方便的功能己单，可以讓我們計算映射到線粒體基因的轉(zhuǎn)錄本的比例唉窃。PercentageFeatureSet()將采用一定模型并搜索基因標識符。此功能可以輕松計算屬于每個細胞的可能功能子集的所有計數(shù)的百分比纹笼。這里的計算只是將屬于該集合的要素的計數(shù)槽中存在的矩陣的列和除以所有要素的列和纹份，然后乘以100。由于我們要繪制比率值，所以我們將反轉(zhuǎn)這一步蔓涧，然后除以100

對于我們的分析件已，我們更愿意使用比率值而不是使用百分比值。因此元暴，我們將通過取輸出值并除以 100 來反轉(zhuǎn)函數(shù)執(zhí)行的最后一步篷扩。

# Compute percent mito ratio
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100

注意：提供的模式（“^MT-”）適用于人類基因名稱。您可能需要根據(jù)您感興趣的生物體調(diào)整模式參數(shù)茉盏。此外鉴未，如果您沒有使用基因名稱作為基因 ID，那么該函數(shù)將無法使用鸠姨，我們有代碼可用于您自己計算此指標铜秆。

其它metadata列

我們現(xiàn)在擁有評估數(shù)據(jù)所需的質(zhì)量指標。但是享怀，我們希望在元數(shù)據(jù)中包含一些有用的附加信息羽峰，包括cell ID 和條件信息。

當我們在上面的元數(shù)據(jù)文件中添加信息列時添瓷，我們只是使用$運算符將其直接添加到 Seurat 對象中的元數(shù)據(jù)中梅屉。對于接下來的幾列數(shù)據(jù)，我們可以繼續(xù)這樣做鳞贷，但我們會將數(shù)據(jù)幀提取到一個單獨的變量中坯汤。通過這種方式，我們可以將元數(shù)據(jù)數(shù)據(jù)幀作為與 seurat 對象分開的實體來使用搀愧，而不會影響存儲在對象內(nèi)的任何其他數(shù)據(jù)惰聂。

讓我們首先通過從 Seurat 對象中metadata提取meta.data槽來創(chuàng)建數(shù)據(jù)幀：
我們現(xiàn)在擁有評估數(shù)據(jù)所需的質(zhì)量指標，同時還需要將其他信息添加到QC指標的元數(shù)據(jù)中咱筛，例如cell ID搓幌、條件信息和其它各種指標。雖然使用$操作符將信息直接添加到Seurat對象的元數(shù)據(jù)槽非常容易迅箩，但是我們選擇把數(shù)據(jù)框提取到一個單獨的變量中溉愁。通過這種方式，我們可以繼續(xù)插入QC分析所需的其他指標饲趋，而不會有影響merded_seurat對象的風險拐揭。
我們將通過從seurat對象提取meta.data槽來創(chuàng)建元數(shù)據(jù)數(shù)據(jù)框：

# Create metadata dataframe
metadata <- merged_seurat@meta.data

接下來，我們將為細胞標識符添加一個新列奕塑。來源于元數(shù)據(jù)數(shù)據(jù)框的行名稱中堂污。我們將按原樣保留行名并將其復(fù)制到名為cells的新列中：

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

可以看到每個細胞 ID 都有一個ctrl_或stim_前綴，正如我們在合并 Seurat 對象時所指定的那樣龄砰。我們可以使用此前綴創(chuàng)建一個新列盟猖，表示每個細胞在哪種條件下分類讨衣。我們將此列稱為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"

最后，我們將重命名元數(shù)據(jù)數(shù)據(jù)框中的一些現(xiàn)有列扒披，以使其更直觀：

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

現(xiàn)在值依，您已設(shè)置好評估數(shù)據(jù)質(zhì)量所需的指標！您最終的元數(shù)據(jù)表將包含與每個細胞相對應(yīng)的行碟案，以及包含有關(guān)這些細胞信息的列：

image

將更新的metadata保存到 Seurat 對象中

在我們評估我們的指標之前，我們要將工作保存到 Seurat 對象颇蜡。我們可以通過簡單地將數(shù)據(jù)幀分配到meta.data中：

# Add metadata back to Seurat object
merged_seurat@meta.data <- metadata

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

評估質(zhì)量指標

現(xiàn)在我們已經(jīng)生成了要評估的各種指標价说，我們可以通過可視化來探索它們。我們將評估各種指標风秤，然后決定哪些細胞質(zhì)量低鳖目，應(yīng)從分析中刪除：

• 細胞計數(shù)
• 每個細胞的 UMI 計數(shù)
• 每個細胞檢測到的基因
• UMIs 與檢測到的基因
• 線粒體計數(shù)比
• Novelty

雙聯(lián)體細胞呢？在單細胞 RNA 測序?qū)嶒炛戌拖遥p聯(lián)體由兩個細胞產(chǎn)生领迈。它們通常是由于細胞分選或捕獲中的錯誤引起的，尤其是在涉及數(shù)千個細胞的基于液滴的protocol中碍沐。當目標是在單細胞水平上表征種群時狸捅，雙峰顯然是不可取的。特別是累提，它們可能會錯誤地暗示存在實際上并不存在的中間種群或過渡狀態(tài)尘喝。因此，最好去除雙峰庫斋陪，這樣它們就不會影響對結(jié)果的解釋朽褪。

為什么我們不檢查雙峰？許多工作流程使用 UMI 或基因的最大閾值无虚，認為檢測到的讀數(shù)或基因數(shù)量多得多缔赠，表明存在多個細胞。雖然這個理由看起來很直觀友题，但并不準確嗤堰。此外，許多用于檢測雙峰的工具往往會去除具有中間或連續(xù)表型的細胞咆爽，盡管它們可能適用于具有非常離散的細胞類型的數(shù)據(jù)集梁棠。Scrublet是一種流行的雙峰檢測工具，但我們還沒有對它進行充分的基準測試斗埂。目前符糊，我們建議此時不包括任何閾值。當我們確定了每個簇的標記后呛凶，我們建議探索這些標記以確定這些標記是否適用于一種以上的細胞類型男娄。

細胞計數(shù)

細胞計數(shù)由檢測到的唯一細胞barcode的數(shù)量決定。對于該實驗，預(yù)計有 12,000 -13,000 個細胞模闲。

在理想情況下建瘫，您希望barcode數(shù)量與細胞數(shù)量一一對應(yīng)。然而尸折，情況并非如此啰脚，因為細胞的捕獲率只是加載內(nèi)容的一部分。例如实夹，與10X 相比橄浓，inDrops 細胞捕獲效率更高（70-80%），而10X的效率在50-60% 之間 亮航。

** 注意:*如果用于庫制備的細胞濃度不準確, 捕獲效率可能會低得多荸实。細胞濃度不應(yīng)由 FACS 機器或生物分析儀測定（這些工具對于濃度測定不準確），而應(yīng)使用血細胞計數(shù)儀或自動細胞計數(shù)器來計算細胞濃度缴淋。

細胞編號也可能因方法而異准给，產(chǎn)生的細胞編號比我們加載的要高得多。例如重抖，在inDrops方法中露氮，細胞條形碼存在于水凝膠中，并與單個細胞和裂解/反應(yīng)混合物封裝在液滴中仇哆。雖然每個水凝膠都應(yīng)該有一個與之相關(guān)的細胞條形碼沦辙，但有時一個水凝膠可以有多個細胞條形碼。類似地讹剔，使用10X協(xié)議油讯，有可能只獲得乳液液滴(GEM)中的條形碼珠子，而沒有實際的細胞延欠。這兩種情況陌兑，加上死亡細胞的存在，都可能導(dǎo)致比細胞更多的細胞條形碼由捎。

# 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")

image

我們看到每個樣本有超過 15,000 個細胞兔综，這比預(yù)期的 12-13,000 個多得多。很明顯狞玛，我們可能存在一些垃圾“細胞”软驰。

每個細胞的 UMI 計數(shù)（轉(zhuǎn)錄本）

每個細胞的 UMI 計數(shù)通常應(yīng)在 500 以上，這是我們預(yù)期的低端心肪。如果 UMI 計數(shù)在 500-1000 計數(shù)之間锭亏，它是可用的，但可能應(yīng)該對細胞進行更深的測序硬鞍。

# 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)

image

我們可以看到兩個樣本中的大多數(shù)細胞都具有 1000 個或更多的 UMI慧瘤，這很好戴已。

每個細胞檢測到的基因

我們對基因檢測的期望與 UMI 檢測相似，盡管它可能比 UMI 低一點锅减。對于高質(zhì)量數(shù)據(jù)糖儡，比例直方圖應(yīng)包含代表被封裝細胞的單個大峰。如果我們看到主峰左側(cè)的小肩峰（不存在于我們的數(shù)據(jù)中）怔匣，或細胞的雙峰分布握联，則可以表明一些情況〗俸荩可能有一組細胞由于某種原因失敗了拴疤。也可能存在生物學(xué)上不同類型的細胞（即靜止細胞群、目標復(fù)雜度較低的細胞）独泞，或一種類型比另一種小得多（即具有高計數(shù)的細胞可能是尺寸較大的細胞）。因此苔埋，應(yīng)使用我們在本課中描述的其他指標來評估此閾值懦砂。

# 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")

image

image

UMIs 與基因檢測

通常一起評估的兩個指標是 UMI 的數(shù)量和每個細胞檢測到的基因數(shù)量。在這里组橄，我們繪制了基因數(shù)量與由線粒體read分數(shù)著色的 UMI 數(shù)量的關(guān)系圖荞膘。線粒體read分數(shù)僅在檢測到的基因很少（深色數(shù)據(jù)點）的特別低計數(shù)細胞中很高。這可能表明受損/垂死的細胞的細胞質(zhì) mRNA 通過破損的膜泄漏玉工，因此羽资，只有位于線粒體中的 mRNA 仍然是保守的。這些細胞被我們的計數(shù)和基因數(shù)量閾值過濾掉遵班。聯(lián)合可視化計數(shù)和基因閾值顯示聯(lián)合過濾效果屠升。

質(zhì)量差的細胞可能是每個細胞的基因和UMI都很低，并且對應(yīng)于圖左下象限中的數(shù)據(jù)點狭郑。好的細胞通常會表現(xiàn)出每個細胞有更多的基因和更高數(shù)量的UMI腹暖。

通過這個圖，我們還評估了線的斜率翰萨，以及圖右下象限中數(shù)據(jù)點的任何散點脏答。這些細胞具有大量的 UMI，但只有少數(shù)基因亩鬼。這些可能是垂死的細胞殖告，但也可能代表低復(fù)雜性細胞類型（即紅細胞）的群體。

# 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)

image

線粒體計數(shù)比

該指標可以確定死亡或瀕死細胞是否存在大量的線粒體污染雳锋。我們將線粒體計數(shù)的低質(zhì)量樣本定義為超過 0.2 線粒體比率標記的細胞黄绩，除非您希望在您的樣本中出現(xiàn)這種情況。

# 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)

image

復(fù)雜性

我們可以通過使用稱為新穎性評分的衡量標準魄缚，根據(jù) RNA 種類的復(fù)雜程度來評估每個細胞宝与。新穎性分數(shù)是通過 nGenes 與 nUMI 的比率來計算的焚廊。如果有很多捕獲的轉(zhuǎn)錄本（高 nUMI）和在細胞中檢測到的基因數(shù)量很少虐急，這可能意味著您只捕獲了少量基因妖碉，并且只是一遍又一遍地對這些較少數(shù)量的基因的轉(zhuǎn)錄本進行測序。這些低復(fù)雜性（低新穎性）細胞可能代表特定的細胞類型（即缺乏典型轉(zhuǎn)錄組的紅細胞）寄月，或者可能是由于一些其他奇怪的人工制品或污染诽里。通常袒餐，我們預(yù)計優(yōu)質(zhì)細胞的新穎性得分高于 0.80。

# 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)

image

注意： 每個細胞的reads是另一個可用于探索的指標谤狡；但是灸眼，使用的工作流程需要保存這些信息以進行評估。通常墓懂，使用此指標焰宣，您希望看到在每個細胞 10,000 到 100,000 個讀數(shù)之間的相對相同位置具有峰值的所有樣本。

過濾

總之捕仔，孤立地考慮這些 QC 指標中的任何一個都可能導(dǎo)致對細胞信號的誤解匕积。例如，線粒體計數(shù)相對較高的細胞可能參與呼吸過程榜跌，并且可能是您想要保留的細胞闪唆。同樣，其他指標可以有其他生物學(xué)解釋钓葫。因此悄蕾，在設(shè)置閾值時，必須始終考慮這些指標的聯(lián)合影響础浮，并將它們設(shè)置為盡可能寬松帆调，以避免無意中過濾掉活細胞群。

細胞水平過濾

現(xiàn)在我們已經(jīng)可視化了各種指標霸旗，我們可以決定要應(yīng)用的閾值贷帮，這將刪除低質(zhì)量的細胞。前面提到的建議通常是一個粗略的指導(dǎo)方針诱告，具體的實驗需要告知所選擇的確切閾值撵枢。我們將使用以下閾值：

• nUMI > 500
• nGene> 250
• log10GenesPerUMI > 0.8
• mitoRatio <0.2

在 Seurat 對象中使用該subset()函數(shù)過濾：

# Filter out low quality cells 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))

基因水平過濾

在我們的數(shù)據(jù)中，我們將有許多零計數(shù)的基因精居。這些基因可以明顯的降低細胞的平均表達锄禽，因此我們將從我們的數(shù)據(jù)中刪除它們。我們將首先確定每個細胞中哪些基因的計數(shù)為零：

# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")

# Output a logical matrix specifying for each gene on whether or not there are more than zero counts per cell
nonzero <- counts > 0

現(xiàn)在靴姿，我們將按prevalence進行過濾沃但。如果一個基因只在少數(shù)細胞中表達，它并不是特別有意義佛吓，因為它仍然會降低所有其他細胞的平均值宵晚。對于我們的數(shù)據(jù)垂攘，我們選擇只保留在 10 個或更多表達的基因細胞。通過使用這個過濾器淤刃，所有細胞中計數(shù)為零的基因?qū)⒈挥行У厝コ?/p>

# 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, ]

最后晒他，獲取這些過濾后的計數(shù)并創(chuàng)建一個新的 Seurat 對象用于下游分析。

# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)

重新評估 QC 指標

執(zhí)行過濾后逸贾，建議查看相關(guān)指標以確保您的數(shù)據(jù)符合您的期望并且有利于下游分析陨仅。

練習(xí)

1. 使用下面提供的代碼從過濾后的 Seurat 對象中提取新的元數(shù)據(jù)：

```
 # Save filtered subset to new metadata
 metadata_clean <- filtered_seurat@meta.data

```

2. 使用過濾后的數(shù)據(jù)執(zhí)行所有相同的 QC 繪圖。

3. 報告每個樣本剩余的細胞數(shù)铝侵，并評論移除的細胞數(shù)是高還是低灼伤。你能解釋為什么這個數(shù)字仍然不是~12K（這是為實驗加載了多少個細胞）？

4. 每個細胞的 nGene 過濾后咪鲜，您仍應(yīng)觀察到主峰右側(cè)的小肩峰狐赡。這個肩膀代表什么？

5. 當根據(jù) nUMI 繪制 nGene 時疟丙，您是否觀察到圖右下象限中的任何數(shù)據(jù)點猾警？你對這些被移除的細胞有什么想法？

保存過濾的細胞

根據(jù)這些 QC 指標隆敢，我們將識別任何不合格的樣本，并繼續(xù)處理我們過濾的細胞崔慧。通常我們會使用不同的過濾標準來遍歷 QC 指標拂蝎；它不一定是一個線性過程。當滿足過濾條件時惶室，我們將保存過濾后的細胞對象以進行聚類和標記識別温自。

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

注意：我們使用的數(shù)據(jù)質(zhì)量非常好。如果您有興趣了解在執(zhí)行 QC 時“不良”數(shù)據(jù)可能是什么樣子皇钞，我們在此處鏈接了一些材料悼泌，我們在其中探索了劣質(zhì)樣品的類似 QC 指標。