4. 質(zhì)控Quality Control

4.1 前言

選擇適合數(shù)據(jù)且不會(huì)過度校正或消除生物效應(yīng)的預(yù)處理方法至關(guān)重要。
由于新的測(cè)序技術(shù)以及捕獲的細(xì)胞、測(cè)量的基因和識(shí)別的細(xì)胞群數(shù)量的不斷增加，用于分析單細(xì)胞RNA測(cè)序數(shù)據(jù)的工具集正在快速發(fā)展。其中許多工具專用于預(yù)處理，旨在解決以下分析步驟：雙聯(lián)體doublet檢測(cè)、質(zhì)量控制、歸一化、特征feature選擇和降維。在本部分中選擇的工具，可能嚴(yán)重影響數(shù)據(jù)的下游分析和解釋。例如，如果在質(zhì)量控制過程中過濾掉過多的細(xì)胞，您可能會(huì)丟失稀有的細(xì)胞亞群，并錯(cuò)過對(duì)有趣的細(xì)胞生物學(xué)的深入了解。然而，如果指定的標(biāo)準(zhǔn)過于寬松，并且沒有在預(yù)處理步驟中排除質(zhì)量較差的細(xì)胞，則可能會(huì)很難注釋細(xì)胞。因此，選擇能夠提供最佳并被證明在下游任務(wù)方面優(yōu)于其他工具的方法非常重要。
本教程的起點(diǎn)是單細(xì)胞數(shù)據(jù)，已按照前面所述的部分進(jìn)行處理。將數(shù)據(jù)比對(duì)以獲得分子計(jì)數(shù)矩陣，即所謂的計(jì)數(shù)矩陣或讀數(shù)矩陣。計(jì)數(shù)矩陣和讀數(shù)矩陣之間的差異取決于單細(xì)胞文庫(kù)構(gòu)建方案中是否包含唯一分子標(biāo)識(shí)符(UMI)。讀數(shù)和計(jì)數(shù)矩陣的維數(shù)為條形碼數(shù) x 轉(zhuǎn)錄本數(shù)。值得注意的是，這里使用術(shù)語(yǔ)“條形碼”而不是“細(xì)胞”，因?yàn)闂l形碼可能錯(cuò)誤地標(biāo)記了多個(gè)細(xì)胞（雙聯(lián)體）或可能沒有標(biāo)記任何細(xì)胞（空滴/孔）。我們將在“雙聯(lián)體檢測(cè)”部分對(duì)此進(jìn)行詳細(xì)說明。

4.2 環(huán)境配置和數(shù)據(jù)

我們使用為2021年NeurIPS會(huì)議上的單細(xì)胞數(shù)據(jù)集成生成10x Multiome數(shù)據(jù)集[Luecken et al., 2021]。該數(shù)據(jù)集捕獲了在四個(gè)不同地點(diǎn)測(cè)量的12名健康人類捐贈(zèng)者的骨髓單核細(xì)胞的單細(xì)胞多組學(xué)數(shù)據(jù)，以獲得嵌套批次效應(yīng)。在本教程中，我們將使用一批上述數(shù)據(jù)集（供體8的樣本4）來(lái)展示scRNA-seq數(shù)據(jù)預(yù)處理的實(shí)戰(zhàn)流程。
第一步，我們首先使用scanpy加載數(shù)據(jù)集。

import numpy as np
import scanpy as sc
import seaborn as sns
from scipy.stats import median_abs_deviation

sc.settings.verbosity = 0
sc.settings.set_figure_params(
    dpi=80,
    facecolor="white",
    frameon=False,
)

adata = sc.read_10x_h5(
    filename="filtered_feature_bc_matrix.h5",
    backup_url="https://figshare.com/ndownloader/files/39546196",
)
adata

輸出結(jié)果：

Variable names are not unique. To make them unique, call `.var_names_make_unique`.
Variable names are not unique. To make them unique, call `.var_names_make_unique`.

AnnData object with n_obs × n_vars = 16934 × 36601
    var: 'gene_ids', 'feature_types', 'genome'

讀取數(shù)據(jù)后，scanpy會(huì)顯示一條警告，指出并非所有變量名稱都是唯一的。這表明某些變量（=基因）出現(xiàn)多次，這可能會(huì)導(dǎo)致下游分析任務(wù)出現(xiàn)錯(cuò)誤或意外行為。我們執(zhí)行建議的函數(shù)var_names_make_unique()，該函數(shù)通過將數(shù)字字符串附加到每個(gè)重復(fù)的索引元素：“1”、“2”等來(lái)使變量名稱唯一。

adata.var_names_make_unique()
adata

輸出結(jié)果：

AnnData object with n_obs × n_vars = 16934 × 36601
    var: 'gene_ids', 'feature_types', 'genome'

數(shù)據(jù)集的形狀為n_obs 16,934 x n_vars 36,601。這轉(zhuǎn)化為barcodes x number of transcripts。我們還檢查.var中有關(guān)gene_ids (Ensembl Id)、feature_types和基因組的更多信息。
大多數(shù)后續(xù)分析任務(wù)假設(shè)數(shù)據(jù)集中的每個(gè)觀測(cè)值代表來(lái)自一個(gè)完整單細(xì)胞的測(cè)量值。在某些情況下，低質(zhì)量細(xì)胞、無(wú)細(xì)胞RNA或雙聯(lián)體的污染可能會(huì)違反這一假設(shè)。本教程將指導(dǎo)您如何糾正和消除這種違規(guī)行為并獲得高質(zhì)量的數(shù)據(jù)集。

4.3 過濾低質(zhì)量的讀數(shù)reads

質(zhì)量控制的第一步是從數(shù)據(jù)集中刪除低質(zhì)量的讀數(shù)。當(dāng)細(xì)胞檢測(cè)到的基因數(shù)量較少、計(jì)數(shù)深度較低且線粒體計(jì)數(shù)較高時(shí)，這表明細(xì)胞膜可能破裂，細(xì)胞正在死亡。由于這些細(xì)胞通常不是我們分析的主要目標(biāo)，并且可能會(huì)扭曲我們的下游分析，因此我們?cè)谫|(zhì)量控制過程中將其去除。為了識(shí)別它們，我們定義了細(xì)胞質(zhì)量控制（QC）閾值。細(xì)胞質(zhì)控通常對(duì)以下三個(gè)質(zhì)控協(xié)變量進(jìn)行：

• 1.每個(gè)條形碼的計(jì)數(shù)數(shù)量（計(jì)數(shù)深度）
• 2.每個(gè)條形碼的基因數(shù)量
• 3.每個(gè)條形碼的線粒體基因計(jì)數(shù)分?jǐn)?shù)

在細(xì)胞QC中，這些協(xié)變量通過閾值過濾，因?yàn)樗鼈兛赡軐?duì)應(yīng)于死亡細(xì)胞。正如所指出的，它們可能反映了膜破裂的細(xì)胞，其細(xì)胞質(zhì)mRNA已泄漏，因此只有線粒體中的mRNA仍然存在。這些細(xì)胞可能會(huì)顯示出較低的計(jì)數(shù)深度、很少檢測(cè)到的基因以及較高比例的線粒體讀數(shù)。然而，共同考慮三個(gè)QC協(xié)變量至關(guān)重要，否則可能會(huì)導(dǎo)致細(xì)胞信號(hào)的誤解。例如，線粒體計(jì)數(shù)相對(duì)較高的細(xì)胞可能參與呼吸過程，不應(yīng)被過濾掉。計(jì)數(shù)低或高的細(xì)胞可能對(duì)應(yīng)于靜止細(xì)胞群或尺寸較大的細(xì)胞。因此，優(yōu)先考慮多個(gè)協(xié)變量再做出閾值決策。一般來(lái)說，建議排除較少的細(xì)胞，并盡可能避免過濾掉活細(xì)胞群或小亞群。
僅對(duì)少數(shù)或小型數(shù)據(jù)集進(jìn)行QC通常以手動(dòng)方式執(zhí)行，方法是查看不同QC協(xié)變量的分布并識(shí)別隨后將被過濾的異常值。然而，隨著數(shù)據(jù)集規(guī)模的增長(zhǎng)，這項(xiàng)任務(wù)變得越來(lái)越耗時(shí)，可能值得考慮通過MAD（median absolute deviations中值絕對(duì)偏差）進(jìn)行自動(dòng)閾值處理。如果細(xì)胞相差5 MADs，我們將其標(biāo)記為異常值，這是一種相對(duì)寬松的過濾策略。我們想強(qiáng)調(diào)的是，在細(xì)胞注釋后重新評(píng)估過濾可能是合理的。
在QC中，第一步是計(jì)算QC協(xié)變量或指標(biāo)。我們使用scanpy函數(shù)sc.pp.calculate_qc_metrics來(lái)計(jì)算這些值，該函數(shù)還可以計(jì)算特定基因群體的計(jì)數(shù)比例。因此，我們定義了線粒體、核糖體和血紅蛋白基因。值得注意的是，線粒體計(jì)數(shù)用前綴“mt-”或“MT-”注釋，具體取決于數(shù)據(jù)集中考慮的物種。如前所述，本篇中使用的數(shù)據(jù)集是人骨髓，因此線粒體計(jì)數(shù)用前綴“MT-”注釋。對(duì)于鼠標(biāo)數(shù)據(jù)集，前綴通常是小寫，即“mt-”。

# mitochondrial genes
adata.var["mt"] = adata.var_names.str.startswith("MT-")
# ribosomal genes
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
# hemoglobin genes.
adata.var["hb"] = adata.var_names.str.contains(("^HB[^(P)]"))

現(xiàn)在我們可以用scanpy計(jì)算對(duì)應(yīng)的QC指標(biāo)。

sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt", "ribo", "hb"], inplace=True, percent_top=[20], log1p=True
)
adata

輸出結(jié)果：

AnnData object with n_obs × n_vars = 16934 × 36601
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'

正如我們所看到的，該函數(shù)向.var和.obs添加了幾個(gè)附加列。我們想在這里重點(diǎn)介紹其中的一些，有關(guān)不同指標(biāo)的更多信息可以在scanpy文檔中找到：

• .obs中的n_genes_by_counts是細(xì)胞中計(jì)數(shù)為正的基因數(shù)量。
• total_counts是細(xì)胞的計(jì)數(shù)總數(shù)，這也可能稱為庫(kù)大小。
• pct_counts_mt是線粒體細(xì)胞總數(shù)的比例。
  現(xiàn)在，我們繪制每個(gè)樣本的三個(gè)QC協(xié)變量n_genes_by_counts、total_counts和pct_counts_mt，以評(píng)估各個(gè)細(xì)胞的捕獲情況。

p1 = sns.displot(adata.obs["total_counts"], bins=100, kde=False)
# sc.pl.violin(adata, 'total_counts')
p2 = sc.pl.violin(adata, "pct_counts_mt")
p3 = sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

輸出結(jié)果：

... storing 'feature_types' as categorical
... storing 'genome' as categorical

這些圖表明，一些讀數(shù)具有相對(duì)較高百分比的線粒體計(jì)數(shù)，這通常與細(xì)胞裂解相關(guān)。但由于每個(gè)細(xì)胞的計(jì)數(shù)數(shù)量足夠高，并且大多數(shù)細(xì)胞的線粒體讀數(shù)百分比低于 20%，我們?nèi)匀豢梢岳^續(xù)下一步分析數(shù)據(jù)。 基于這些圖，現(xiàn)在還可以定義用于過濾細(xì)胞的手動(dòng)閾值。在這里，我們將展示基于MAD的自動(dòng)閾值和過濾的QC。
首先，我們定義一個(gè)metric函數(shù)，即.obs中的一列和過濾策略中仍然允許的MAD(nmad) 數(shù)量。

def is_outlier(adata, metric: str, nmads: int):
    M = adata.obs[metric]
    outlier = (M < np.median(M) - nmads * median_abs_deviation(M)) | (
        np.median(M) + nmads * median_abs_deviation(M) < M
    )
    return outlier

我們現(xiàn)在將此函數(shù)應(yīng)用于log1p_total_counts、log1p_n_genes_by_counts和pct_counts_in_top_20_genesQC協(xié)變量，每個(gè)協(xié)變量的閾值為5 MADs。

adata.obs["outlier"] = (
    is_outlier(adata, "log1p_total_counts", 5)
    | is_outlier(adata, "log1p_n_genes_by_counts", 5)
    | is_outlier(adata, "pct_counts_in_top_20_genes", 5)
)
adata.obs.outlier.value_counts()

輸出結(jié)果：

False    16065
True       869
Name: outlier, dtype: int64

pct_counts_Mt使用3 MADs進(jìn)行過濾。此外，線粒體計(jì)數(shù)百分比超過8%的細(xì)胞也會(huì)被濾除。

adata.obs["mt_outlier"] = is_outlier(adata, "pct_counts_mt", 3) | (
    adata.obs["pct_counts_mt"] > 8
)
adata.obs.mt_outlier.value_counts()

輸出結(jié)果：

False    15240
True      1694
Name: mt_outlier, dtype: int64

現(xiàn)在，我們根據(jù)這兩個(gè)附加列來(lái)過濾AnnData對(duì)象。

print(f"Total number of cells: {adata.n_obs}")
adata = adata[(~adata.obs.outlier) & (~adata.obs.mt_outlier)].copy()

print(f"Number of cells after filtering of low quality cells: {adata.n_obs}")

輸出結(jié)果：

Total number of cells: 16934
Number of cells after filtering of low quality cells: 14814

p1 = sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

4.4 環(huán)境RNA校正

對(duì)于基于液滴的單細(xì)胞RNA-seq實(shí)驗(yàn)，稀釋液中存在一定量的背景mRNA，這些mRNA會(huì)與細(xì)胞一起分布到液滴中并與其一起測(cè)序。其最終效果是產(chǎn)生背景污染，該背景污染代表的表達(dá)不是來(lái)自液滴內(nèi)包含的細(xì)胞，而是來(lái)自包含細(xì)胞的溶液。
基于液滴的scRNA-seq生成跨多個(gè)細(xì)胞的基因的唯一分子標(biāo)識(shí)符 (UMI) 計(jì)數(shù)，旨在識(shí)別每個(gè)基因和每個(gè)細(xì)胞的分子數(shù)量。它假設(shè)每個(gè)液滴都含有來(lái)自單個(gè)細(xì)胞的mRNA。 雙聯(lián)體、空滴和無(wú)細(xì)胞RNA可能違反這一假設(shè)。無(wú)細(xì)胞mRNA分子代表稀釋液中存在的背景mRNA。這些分子沿著液滴分布，并與它們一起測(cè)序。輸入溶液中無(wú)細(xì)胞mRNA的這種污染通常稱為“湯”，是由細(xì)胞裂解產(chǎn)生的。

無(wú)細(xì)胞mRNA分子（也稱為環(huán)境RNA）可能會(huì)混淆觀察到的計(jì)數(shù)數(shù)量，并可被視為背景污染。糾正基于液滴的scRNA-seq數(shù)據(jù)集以獲得無(wú)細(xì)胞mRNA非常重要，因?yàn)樗赡軙?huì)扭曲我們下游分析中數(shù)據(jù)的解釋。一般來(lái)說，每個(gè)輸入溶液的湯都不同，并且取決于數(shù)據(jù)集中各個(gè)細(xì)胞的表達(dá)模式。去除環(huán)境mRNA的方法（例如SoupX和 DecontX）旨在估計(jì)湯的成分，并根據(jù)湯的表達(dá)校正計(jì)數(shù)表。
第一步，SoupX計(jì)算湯的profile。它根據(jù)未過濾的Cellranger矩陣給出的空液滴估計(jì)環(huán)境mRNA表達(dá)譜。接下來(lái)，SoupX估計(jì)細(xì)胞特定的污染分?jǐn)?shù)。最后，它根據(jù)環(huán)境mRNA表達(dá)譜和估計(jì)的污染來(lái)校正表達(dá)矩陣。
SoupX的輸出是修改后的計(jì)數(shù)矩陣，可用于任何下游分析工具。
我們現(xiàn)在加載運(yùn)行SoupX所需的python和R包。

import anndata2ri
import logging

import rpy2.rinterface_lib.callbacks as rcb
import rpy2.robjects as ro

rcb.logger.setLevel(logging.ERROR)
ro.pandas2ri.activate()
anndata2ri.activate()

%load_ext rpy2.ipython

%%R
library(SoupX)

SoupX可以在沒有聚類信息的情況下運(yùn)行，但是如果提供基本聚類，結(jié)果會(huì)更好。SoupX可以與cellranger生成的默認(rèn)cluster一起使用，也可以通過手動(dòng)定義cluster來(lái)使用。我們將在本篇中展示后者，因?yàn)镾oupX的結(jié)果對(duì)所使用的聚類并不強(qiáng)烈敏感。
現(xiàn)在，我們創(chuàng)建AnnData對(duì)象的副本，對(duì)其進(jìn)行標(biāo)準(zhǔn)化，降低其維度，并在處理后的副本上計(jì)算默認(rèn)leiden簇。后續(xù)章節(jié)將更詳細(xì)地介紹聚類。現(xiàn)在我們只需要知道leiden聚類為我們提供了數(shù)據(jù)集中細(xì)胞的分區(qū)（社區(qū)）。我們將獲得的簇保存為soupx_groups并刪除AnnData對(duì)象的副本以節(jié)省內(nèi)存。
首先，我們生成AnnData對(duì)象的副本，對(duì)其進(jìn)行標(biāo)準(zhǔn)化和log1p轉(zhuǎn)換。此時(shí)我們使用簡(jiǎn)單的移位對(duì)數(shù)歸一化。

adata_pp = adata.copy()
sc.pp.normalize_per_cell(adata_pp)
sc.pp.log1p(adata_pp)

接下來(lái)，我們計(jì)算數(shù)據(jù)的主成分以獲得較低維的形式。然后使用該形式來(lái)生成數(shù)據(jù)的鄰域圖并在KNN 圖上運(yùn)行l(wèi)eiden聚類。我們將簇作為soupx_groups添加到.obs并將它們保存為向量。

sc.pp.pca(adata_pp)
sc.pp.neighbors(adata_pp)
sc.tl.leiden(adata_pp, key_added="soupx_groups")

# Preprocess variables for SoupX
soupx_groups = adata_pp.obs["soupx_groups"]

現(xiàn)在，我們可以刪除AnnData對(duì)象的副本，因?yàn)槲覀兩闪丝稍趕oupX中使用的簇向量。

del adata_pp

接下來(lái)，我們保存細(xì)胞名稱、基因名稱和過濾后的cellranger輸出的數(shù)據(jù)矩陣。SoupX需要features x barcodes的矩陣，因此我們必須轉(zhuǎn)置.X。

cells = adata.obs_names
genes = adata.var_names
data = adata.X.T

SoupX還需要細(xì)胞矩陣的原始基因，該矩陣在cellranger輸出中通常稱為raw_feature_bc_matrix.h5。我們像之前一樣使用scanpy加載filtered_feature_bc_matrix.h5并在對(duì)象上運(yùn)行.var_names_make_unique()并轉(zhuǎn)置相應(yīng)的.X。

adata_raw = sc.read_10x_h5(
    filename="raw_feature_bc_matrix.h5",
    backup_url="https://figshare.com/ndownloader/files/39546217",
)
adata_raw.var_names_make_unique()
data_tod = adata_raw.X.T

輸出結(jié)果：

Variable names are not unique. To make them unique, call `.var_names_make_unique`.
Variable names are not unique. To make them unique, call `.var_names_make_unique`.

del adata_raw

現(xiàn)在，我們已準(zhǔn)備好運(yùn)行SoupX的一切。輸入是經(jīng)過過濾的條形碼 x 細(xì)胞的cellranger矩陣、來(lái)自cellranger的barcodes x droplets液滴原始表、基因和細(xì)胞名稱以及通過簡(jiǎn)單的leiden聚類獲得的簇。輸出將是校正后的計(jì)數(shù)矩陣。
我們首先從液滴表和細(xì)胞表構(gòu)建一個(gè)SoupChannel。接下來(lái)，我們將元數(shù)據(jù)添加到SoupChannel對(duì)象中，該對(duì)象可以是data.frame形式的任何元數(shù)據(jù)。 我們?cè)谶@里添加:

%%R -i data -i data_tod -i genes -i cells -i soupx_groups -o out 

# specify row and column names of data
rownames(data) = genes
colnames(data) = cells
# ensure correct sparse format for table of counts and table of droplets
data <- as(data, "sparseMatrix")
data_tod <- as(data_tod, "sparseMatrix")

# Generate SoupChannel Object for SoupX 
sc = SoupChannel(data_tod, data, calcSoupProfile = FALSE)

# Add extra meta data to the SoupChannel object
soupProf = data.frame(row.names = rownames(data), est = rowSums(data)/sum(data), counts = rowSums(data))
sc = setSoupProfile(sc, soupProf)
# Set cluster information in SoupChannel
sc = setClusters(sc, soupx_groups)

# Estimate contamination fraction
sc  = autoEstCont(sc, doPlot=FALSE)
# Infer corrected table of counts and rount to integer
out = adjustCounts(sc, roundToInt = TRUE)

SoupX成功推斷出校正后的計(jì)數(shù)，我們現(xiàn)在可以將其存儲(chǔ)為附加層layer。在以下所有分析步驟中，我們希望使用SoupX校正計(jì)數(shù)矩陣，因此我們用soupX校正矩陣覆蓋.X。

adata.layers["counts"] = adata.X
adata.layers["soupX_counts"] = out.T
adata.X = adata.layers["soupX_counts"]

接下來(lái)，我們另外過濾掉少于20個(gè)細(xì)胞中未檢測(cè)到的基因，因?yàn)檫@些基因沒有提供信息。

print(f"Total number of genes: {adata.n_vars}")

# Min 20 cells - filters out 0 count genes
sc.pp.filter_genes(adata, min_cells=20)
print(f"Number of genes after cell filter: {adata.n_vars}")

輸出結(jié)果：

Total number of genes: 36601
Number of genes after cell filter: 20171

4.5 雙聯(lián)體檢測(cè)

雙聯(lián)體被定義為在相同細(xì)胞條形碼下測(cè)序的兩個(gè)細(xì)胞，它們是在同一液滴中捕獲的。這就是為什么我們到目前為止使用術(shù)語(yǔ)“條形碼”而不是“細(xì)胞”。如果雙聯(lián)體由相同細(xì)胞類型（但來(lái)自不同個(gè)體）形成，則稱為同型，否則稱為異型。同型雙聯(lián)體不一定可以從計(jì)數(shù)矩陣中識(shí)別出來(lái)，并且通常被認(rèn)為是無(wú)害的，因?yàn)樗鼈兛梢酝ㄟ^細(xì)胞散列或SNP來(lái)識(shí)別。因此，它們的識(shí)別不是雙聯(lián)體檢測(cè)方法的主要目標(biāo)。
由不同細(xì)胞類型或狀態(tài)形成的雙聯(lián)體稱為異型。它們的識(shí)別至關(guān)重要，因?yàn)樗鼈兒芸赡鼙诲e(cuò)誤分類，并可能導(dǎo)致下游分析步驟的錯(cuò)誤。因此，雙聯(lián)體檢測(cè)和去除通常是初始預(yù)處理步驟。雙聯(lián)體可以通過其大量讀取和檢測(cè)到的特征來(lái)識(shí)別，也可以通過創(chuàng)建人工雙聯(lián)體并將其與數(shù)據(jù)集中存在的細(xì)胞進(jìn)行比較的方法來(lái)識(shí)別。雙聯(lián)體檢測(cè)方法計(jì)算效率高，并且存在多個(gè)用于此任務(wù)的軟件包。
在本教程中，我們將展示scDblFinder R包。scDblFinder隨機(jī)選擇兩個(gè)液滴，并通過平均它們的基因表達(dá)譜來(lái)創(chuàng)建人工雙聯(lián)體。然后，雙聯(lián)體得分定義為主成分空間中每個(gè)液滴的k最近鄰圖中人工雙聯(lián)體的比例。

我們首先載入一些python和R包。

%%R
library(Seurat)
library(scater)
library(scDblFinder)
library(BiocParallel)

data_mat = adata.X.T

現(xiàn)在，我們可以使用data_mat作為SingleCellExperiment中scDblFinder的輸入來(lái)啟動(dòng)雙聯(lián)體檢測(cè)。scBblFinder向sce的colData添加幾列。其中三個(gè)主要的組成：

• sce$scDblFinder.score：最終的雙聯(lián)體得分（越高，細(xì)胞越有可能是雙聯(lián)體）
• sce$scDblFinder.ratio：細(xì)胞鄰域中人工雙聯(lián)體的比例
• sce$scDblFinder.class：分類（雙聯(lián)體或單聯(lián)體）

我們只會(huì)輸出參數(shù)并將其存儲(chǔ)在.obs的AnnData對(duì)象中。 其他參數(shù)可以類似地添加到AnnData對(duì)象。

%%R -i data_mat -o doublet_score -o doublet_class

set.seed(123)
sce = scDblFinder(
    SingleCellExperiment(
        list(counts=data_mat),
    ) 
)
doublet_score = sce$scDblFinder.score
doublet_class = sce$scDblFinder.class

scDblFinder輸出具有分類Singlet(1) 和Doublet (2) 的類。我們將其添加到.obs中的AnnData對(duì)象中。

adata.obs["scDblFinder_score"] = doublet_score
adata.obs["scDblFinder_class"] = doublet_class
adata.obs.scDblFinder_class.value_counts()

輸出結(jié)果：

singlet    11956
doublet     2858
Name: scDblFinder_class, dtype: int64

我們建議暫時(shí)將已識(shí)別的雙聯(lián)體留在數(shù)據(jù)集中，并在可視化過程中檢查雙聯(lián)體。
在下游聚類過程中，重新評(píng)估質(zhì)量控制和所選參數(shù)可能會(huì)很有用，以過濾掉更多或更少的細(xì)胞。我們現(xiàn)在可以保存數(shù)據(jù)集并繼續(xù)之后的標(biāo)準(zhǔn)化過程。

adata.write("s4d8_quality_control.h5ad")