ダブレットとは

シングルセル解析で現在主流になりつつあるのは10x Genomics Chromiumのようなドロップレットベースの方法である。この方法はドロップレットの中で磁気ビーズと細胞を封入し、mRNAの逆転写・増幅を行う (下図b)。ビーズにはポリT領域と、ビーズごとにユニークなバーコード配列がついたオリゴDNA (図c) が生えており、この上で細胞由来のmRNAが逆転写・増幅される仕組みである。この時にバーコードがついた状態でmRNAは逆転写・増幅されるので、配列解析時に、同じバーコードを持つリードは同一細胞に由来すると判別することができる。

引用元：Zheng et al. “Massively parallel digital transcriptional profiling of single cells” Nat Commun. 2017 Jan 16;8:14049. doi: 10.1038/ncomms14049. CC-BY 4.0

一つのドロップレットの中にどのように一つの細胞を閉じ込めているのか？これは機械的に制御されているわけではない。ドロップレットの中に封入される細胞の数は確率的に決まり、ポアソン分布に従う。例えば、一つのドロップレットあたり平均0.5個の細胞が封入されるような濃度で実験を行えば、確率的に約60%, 30%, 8%のドロップレットには0個, 1個, 2個の細胞が封入される *1。すなわち細胞の濃度をある程度下げることで、ほとんどの細胞は単独でドロップレットの中に封入されることとなり、問題なくシングルセルレベルでの計測をすることができる、という算段である。

しかし実際には数%のドロップレットには複数個 (多くは2個) の細胞が入ってしまう (下図A)。この場合は、同じビーズ上でmRNAが処理され、同一のバーコードが付けられてしまう。したがって、配列解析時に異なる細胞由来のmRNAであることが区別できないため、複数個の細胞の遺伝子発現が、一つの細胞として、混ざった状態で検出されることになる。これがいわゆるダブレットと呼ばれるアーティファクトである。

ダブレットは解析結果を歪めかねない。一見するとユニークな細胞集団に見えて、実は単なるアーティフィシャルであることも多い (下図B)。

引用元：Wolock et al. “Scrublet: computational identification of cell doublets in single-cell transcriptomic data” bioRxiv. 2018 https://doi.org/10.1101/357368 CC-BY 4.0

Chromiumのマニュアル*2を見ると、１サンプルあたり10000細胞を解析する条件では、約7.6%がダブレットになるようだ。これは結構大きな影響がありそうなものだが、多くの論文ではあまりケアされていないように感じる。ダブレットを除くために、検出遺伝子数 (もしくはUMI数) が極端に多い細胞を解析から除外する前処理をよく目にするが、これは理想的とは言えない。検出遺伝子数は細胞種間で不均一であり、活性化した細胞や癌細胞は多くの遺伝子を発現しているからだ。

そんな中で2019年ごろから、よりエレガントにダブレットを解析から除外するための方法論が提案されてきた。ここではそれらの方法についてレビューする。

Scrublet, DoubletFinder

最初に報告された2つの手法 Scrublet*3とDoubletFinder*4。2019年4月の同時期にCell Systemsに公開されたこれらの手法はいずれも似たアプローチを提案している。基本的には以下の流れ。

入力: UMIカウント行列 (細胞 x 遺伝子)
1からランダムに細胞ペアをN組選出
2のペアのカウントデータをミックス仮想ダブレットデータを作成
3を1と混ぜる
Highly-variable geneの抽出 & PCAで次元圧縮 *5
主成分得点を使ってkNN graph *6 を描画
6の情報をもとにスコアを計算

Scrubletは以下の式でスコアが求まる*7。ダブレットはこのスコアが高くなる。

$score_{i}=\frac{q_{i} \hat{\rho} / r}{ 1 - \hat{\rho} - q_{i}(1 - \hat{\rho} - \hat{\rho} / r) }$
$q_{i} = \frac{k_{d}(i) + 1}{k_{adj} + 2}$
$\hat{\rho}$ : 予想されるダブレットの割合
$r = M_{d}/M_{obs}$
$k_{adj}$ : k
$k_{d}(i)$ : 細胞iのkNNに含まれる仮想ダブレットの数
$M_{obs}$ : 全細胞数
$M_{d}$ : 作成した仮想ダブレットの数

$\hat{\rho}$ はユーザーが解析の際に指定する。（このパラメータが結果に影響することがgithubに記載されていた）

シングルセルデータにはダブレットとシングレットが混ざっているのでスコアの分布は二峰性を示す。この分布をもとにユーザーがカットオフをマニュアルで設定してダブレットを特定することが推奨されている。

github.com

DoubletDecon

2019年に報告された手法。Scrublet/DoubletFinderが $\hat{\rho}$ などのパラメータをユーザーが指定しなければならないことを問題視 *8。仮想ダブレットを作成するプロセスは共通しているが、発現変動遺伝子の情報を使って、過度に細胞を除きすぎないようにするプロセス (Rescue step) を設けている点が特徴。

入力: 正規化済み発現行列 (細胞 x 遺伝子)
前処理 (低発現遺伝子除外 & 質の低い細胞を除外)
クラスタリング *9
3のクラスタごとに代表発現プロファイル (Centroid) を計算
Cluster merging: 類似したクラスタをマージする *10
仮想ダブレットを作成: 5でできたクラスタのペアごとにランダムに細胞を選んで仮想ダブレットを作成
4をリファレンスに各細胞/仮想ダブレットの遺伝子発現プロファイルをDeconvolute。下図Eのような組成を計算する
7で得られた組成が仮想ダブレットの組成と類似している細胞は除外
8で除外した細胞を対象に3と同様の方法でクラスタリング
3と9の全クラスタ間で発現変動遺伝子を抽出
Rescue step: 特異的遺伝子が検出されたクラスタは9からレスキューする

引用元：DePasquale et al. “DoubletDecon: Cell-State Aware Removal of Single-Cell RNA-Seq Doublets” bioRxiv. 2019 doi: https://doi.org/10.1101/364810 CC-BY 4.0

イントロではScrublet/DoubletFinderのパラメータ設定を問題視していたが、DoubletDeconも5, などのステップでad hocにユーザーがパラメータを決めなければならない様子。 Rescue stepは2つの細胞が混ざってできたダブレットであれば特異的に発現が動いている遺伝子はないという仮定に基づいている。この考え方は実装も楽だし、普遍的に使えそう。

github.com

SCDS

2020年に公開された手法。Co-expression based doublet scoring (cxds)とBinary classification based doublet scoring (bcds)という2つのスコアをもとにダブレットを推定する。

Co-expression based doublet scoring (cxds)

ある遺伝子 $i$ が検出される細胞の割合を $p_i$ とする。仮に2つの遺伝子 $i$ と $j$ の発現が独立している場合は、2つのうちどちらかが検出される細胞の数 $n_{ij}$ は次の二項分布に従う。これを利用すると2つの遺伝子の発現の排他性 (互いにExclusiveか) を検定することができる。

$n_{ij} \sim Bin(n, p_{i}(1 - p_{j}) + p_{j}(1 - p_{i}))$

cxdsは以下のように求められる。通常は排他的に発現しているはずの遺伝子ペアが同時に発現しているほどスコアが高くなる仕組み。

$cxds(i) = \sum_{k}\sum_{j} B_{ki}B_{ji}S_{ij}$
$B_{ki}$ : 細胞iにおいて遺伝子kが検出されていれば1, そうでなければ0のバイナリ行列
$S_{ij}$ : 上記の二項検定で得られたp-valueを-log10変換したもの (遺伝子iとjの排他性)