CITEseqとは

シングルセルトランスクリプトームと表面タンパク質発現量*1を同一細胞から計測する技術*2。概要としては以下のような手法である。

オリゴDNAを付与した抗体を使って細胞を標識。このオリゴDNAにはポリA配列*3、標的タンパク毎にに異なるバーコード配列が含まれる
ドロップレットの中に磁気ビーズ*4と細胞を封入
磁気ビーズ上でオリゴDNAとmRNAを増幅
シーケンシング
配列がどの抗体・細胞に由来していたかはバーコード情報から判別可能。これをもとに表面タンパク発現の定量を行う

引用元：Matula et al. 2020. “ Single-Cell Analysis Using Droplet Microfluidics.” Adv Biosyst. 2020 Jan;4(1):e1900188. doi: 10.1002/adbi.201900188.. CC-BY 4.0

10x Genomicsよりキットがリリースされており*5、CITEseq用のオリゴDNAを付与した抗体もBiolegend社から販売されている*6ため、少しづつではあるが報告例が増えてきている。また今年はIdo Amitのグループが細胞内タンパク質を抗体標識後、FACSで任意の細胞を分取しscRNA-seq解析を行う手法INs-seqを提案している*7。本手法は同時標識できるタンパク質の数は限定されているが、今後CITEseqのようにDNAバーコード技術を利用することで、より多くの細胞内タンパク質をシングルセルレベルで同時計測する技術が現れるものと予想される。

CITEseqデータの正規化

scRNA-seqデータの正規化手法は複数提案されているが、多くの研究では、Smart-seq2のようなFull-lengthの手法であればTPM (Transcript per million)、ChromiumのようなUMIベースの手法であればCP10k (Count per 10 thousands) などのシンプルな正規化手法が使われていると思う。これらの正規化手法が使われる背景には「一つの細胞内におけるmRNAの総量はほぼ等しい」という仮定が存在する。

一方でCITEseqはプロテオミクスのような網羅的計測でなく、せいぜい百種類の表面タンパク質のみが対象となる。この場合は「計測した表面タンパク質の総量が細胞間でほぼ等しい」という仮定は明らかに現実と乖離していることが予想されるため、TPMやCP10kのような合計値でのシンプルな正規化は行うことができない。このような場合にはどのように正規化をするのだろうか？色々と調べているうちに疑問に思ったので以下にメモを残す。

CITEseq 原著論文*8の場合

StoeckiusらのCITEseq原著論文ではCLR (Centered log-ratio) transformationを適用している。細胞i・タンパク質jのリードカウント数をXijとすると、変換後の値Yijは以下のようになる。

$Y_{i.} = clr(X_{i.}) = (ln(\frac{X_{i1}}{g(X_{i.})}), ln(\frac{X_{i2}}{g(X_{i.})}), ... , ln(\frac{X_{in}}{g(X_{i.})}))$
$g(X_{i.})$ ： $X_{i.}$ の幾何平均

この変換は組成データ(e.g. パーセント) によく用いられる。組成データが持つ、和が一定の値になる制約 (定数和制約) による不都合*9を解消するために使われるらしい。このあたりの話は太田らの資料が大変参考になった。

マテメソを読むと"we treated this data type as compositional data"と書いてある。結局「計測した表面タンパク質の総量が細胞間でほぼ等しい」という仮定を置いているということだと思うが、そのロジックは理解できない。明らかに間違った値へと変換される場合もあると思う*10。

Seuratの場合

Seuratでは関数NormalizeDataにCITEseqデータ用の正規化手法も実装されている。Vignetteによると、ここでもCLRが推奨されているが、"This is a slightly improved procedure from the original publication"と若干の改良がなされたことも記載されている*11。どのような改良がされているのだろうか？見当たらないのでソースコードを見てみた。

github.com

2326行目〜

      'CLR' = CustomNormalize(
        data = object,
        custom_function = function(x) {
          return(log1p(x = x / (exp(x = sum(log1p(x = x[x > 0]), na.rm = TRUE) / length(x = x)))))
        },
        margin = margin, # default margin = 1
        verbose = verbose
        # across = across
      )

3173行目〜

CustomNormalize <- function(data, custom_function, margin, verbose = TRUE) {
  # margin <- switch(
  #   EXPR = across,
  #   'cells' = 2,
  #   'features' = 1,
  #   stop("'across' must be either 'cells' or 'features'")
  # )
  norm.data <- myapply(
    X = data,
    MARGIN = margin,
    FUN = custom_function)
  if (margin == 1) {
    norm.data = t(x = norm.data)
  }
}

該当する箇所を上に示しました。どうやらデフォルトではFeatureごとにlog1p(x = x / (exp(x = sum(log1p(x = x[x > 0]), na.rm = TRUE) / length(x = x))))という変換を加えているようだ。

原著論文では細胞ごとに全タンパクの幾何平均を求めてCLR変換をしていたわけだが、ここではタンパク質ごとに全細胞の幾何平均を求めてCLR変換をしていた。根拠は不明。そしてlog1p (log(x+1))を使っているがこれも謎。同じところを気になった人はいるようだがSeuratチームの回答はありませんでした*12。

DSB

新しい手法DSBが提案されたようなのでチェック。手法の論文はまだbioRxivでしか公開されていないが、すでにDSBを使ってCITEseqデータを解析した事例がNature Medicineに出ていた*13。

mattpm.github.io

本手法はバックグラウンド(空のドロップレット)の情報を使う。確かにこれであればCLRのような組成データの仮定がないので問題が少なそう。

$Y_{ij} = \frac{log(X_{ij} + P) - \mu}{\sigma}$
$P$ ：pseudocount to stabilize variance of small values (default 10)
$\mu$ ：バックグラウンドにおけるタンパクjの平均リードカウント数
$\sigma$ ：バックグラウンドにおけるタンパクjのカウント数のSD