Aurora blog

バイオインフォ・バイオテクノロジーにまつわる情報

シングルセル解析⑦：RNA velocity

Bioinfo Single-cell

はじめに

Trajectory inferenceはscRNA-seqデータから細胞状態の遷移パターンを推定する解析手法である。一般的にはこれらの手法は、細胞の遺伝子発現データの分布から、どの細胞種（細胞サブセット）が連続的な関係にあるかを推定することができるが、「どの細胞がどの細胞へ向かっているのか」といった状態遷移の方向性については推定できず、解析者の生物学的な知識で解釈する *1 必要があった。

引用元：Wikipedia: Trajectory inference. https://en.wikipedia.org/wiki/Trajectory_inference CC-BY 4.0

その中で2018年にシングルセルデータから細胞状態遷移の方向性を推定する手法RNA velocityがLa Mannoらによって提案された。後述の通り、scRNA-seqデータ以外のデータを必要とせず、転写物に含まれるイントロンに着目するだけで推定できる画期的なアイデアであったが、しばしば誤った結果が得られることもあり使いづらい印象だった。RNA velocityの現状と課題についてのレビューが最近Theis labから公開されていたので改めてまとめてみる。

La Manno et al. 2018. Nature

RNA velocityのオリジナル論文。Gieole La Manno、Sten Linnarsson、Peter Kharchenkoのグループによる報告。RNA velocityでは「イントロンが取り除かれていない転写物は転写してからの時間が間もない」→「イントロンを含む転写産物（=新しい転写産物）が多い遺伝子ほど転写が活発に行われている」→「その遺伝子の発現量は次の時点で上昇する」という仮定を置く（下図B）ことで、イントロンを含有する転写産物の量から「近い未来の遺伝子発現の変化（細胞状態遷移）の方向性」を推定する。イントロンの情報を利用する本手法は多様なscRNA-seqデータに利用できる。論文では、Smart-seq2などの転写産物の全長を解析する手法だけでなく、Chromiumなどの3'端（または5'端）のみを解析する手法においてもイントロンが検出できることが示されている（下図A）。

RNA velocityではイントロンを含む転写産物 (Unspliced mRNA) と含まない転写産物 (Spliced mRNA) について遺伝子ごとに以下のようなモデルを置く（下図b）。

  • u(t): 時間tにおけるUnspliced mRNAの量
  • s(t): 時間tにおけるSpliced mRNAの量
  • \alpha(t): 時間tにおける転写スピード
  • \beta(t): 時間tにおけるスプライシングスピード
  • \gamma(t): 時間tにおけるmRNA分解スピード
  • \frac{du}{dt} = \alpha(t) - \beta(t)u(t): Unspliced mRNA量の変化
  • \frac{ds}{dt} = \beta(t)u(t) - \gamma(t)s(t): Spliced mRNA量の変化

引用元：La Manno et al. “RNA velocity in single cells” bioRxiv. 2018. https://doi.org/10.1101/206052 CC-BY 4.0

RNA velocityでは以下の仮定を置くことで、さらにシンプルなモデルを考える。

  • 転写スピード \alphaとmRNA分解スピード \gammaは時間によらず一定である
  • スプライシングスピード \betaは遺伝子・時間によらず一定である
  • スプライシングスピード \betaを \beta = 1とする*2
  • \frac{du}{dt} = \alpha - u(t): Unspliced mRNA量の変化
  • \frac{ds}{dt} = u(t) - \gamma s(t): Spliced mRNA量の変化

ある遺伝子の発現量が安定した状態 (定常状態: \frac{du}{dt} = 0, \frac{ds}{dt} = 0) の場合は以下の関係が成り立つ。従って定常状態にあるデータを用いればSpliced/Unspliced mRNAの量から転写スピード \alphaと分解スピード \gammaを推定することができる。

  • \alpha = u(t)
  • u(t) = \gamma s(t)

RNA velocityでは定常状態にあるSpliced/Unspliced mRNAのデータ s, uを使い最小二乗法でパラメータ \gammaを求める。遺伝子によってはエクソン・イントロンの情報が完全ではないのでOffset項 oを設けてその影響を考慮する。RNA velocityでは uと sは以下のように求められる。安定したデータを出すために近傍k個の細胞を考慮する。

  • u \sim \gamma * s + o
  • u = U / N_{u}, s = S / N_{s}
    • U, S: Unspliced/Spliced read count of the gene in k-nearest neighbor cells
    • N_{u}, N_{s}: Total counts of unspliced/spliced reads in k-nearest neighbor cells

ここで定常状態にあるデータをどのように見つけるか、が問題になるが、 s, uが最大・最小に近い細胞において (上図Dの右上・左下) に位置する細胞が定常状態 ( \frac{du}{dt} = 0, \frac{ds}{dt} = 0)) に近いと仮定する。RNA velocityでは s, uが上位 x%に位置する細胞を用いて上記の回帰を行う。この手法 (Extreme quantile regression) は全データを使った回帰と比較して安定した結果が出ることを示している *3。求めたパラメータ \gammaと \alphaを使うことで定常状態にない細胞における各遺伝子の転写スピード \frac{du}{dt}, \frac{ds}{dt}を計算できる。このように求めた転写スピード \frac{ds}{dt}を使って各細胞がどの細胞状態へと近づいているかを推定する。

www.nature.com

Bergen et al. 2020. Nat Biotechnol

Fabian Theisグループから発表されたRNA velocityの派生手法。Extreme quantile regressionにおける定常状態がデータに含まれている、という仮定に従わない*4場合にも適用できる手法を考案。ここでは定常状態のデータを使った線形回帰ではなく、 u(t), s(t)のカーブにフィッティングする方法をとる。

  • \frac{du}{dt} = \alpha(t) - \beta(t)u(t)
  • \frac{ds}{dt} = \beta(t)u(t) - \gamma(t)s(t)

オリジナル論文と同じ上式を変形し、本論文は以下のモデルを考える。オリジナルの手法は転写スピード \alphaが時間によらず一定であると仮定していたが、このモデルでは \alphaに4つのState (Induction ( k=1), Repression ( k=0), Steady state at maximum ( k=ss_1), Steady state at minimum ( k=ss_0) があり、各細胞が潜在的にいずれかのState k・時間 tに該当すると考える。*5

  • u(t) = u_{0}e^{-\beta\tau} + \frac{\alpha^{(k)}}{\beta} (1 - e^{-\beta\tau})
  • s(t) = s_{0}e^{-\gamma\tau} + \frac{\alpha^{(k)}}{\gamma} (1 - e^{-\gamma\tau}) + \frac{\alpha^{(k)} - \beta u_{0}}{\gamma - \beta} (e^{-\gamma\tau} - e^{-\beta\tau})
  • \tau = t - t_{0}^{(k)}

この手法では各細胞のパラメータ ( k, t) と遺伝子発現のパラメータ ( \alpha^{(k)}, \beta, \gamma) を推定する問題を解く。現実のデータ u^{obs}と s^{obs}と近い値を返すモデルのパラメータを探す。ここでは実測値との残差 eを x=(u, s)の符号付きユークリッド距離として、これができるだけ小さくなるような調整を行う。

  • e_i = sign(s_i^{obs} - \hat{s}_{t_i}) \cdot |x_{i}^{obs} - \hat{x}_{t_i}(\theta)|: 細胞iの残差

ここでは残差が正規分布に従うと考え、以下の尤度 L(\theta) が最大 (=negative log-likelihood l(\theta) が最小) になるパラメータを探索する。

  • \theta = ( \alpha^{(k)}, \beta, \gamma )
  • L(\theta) = \frac{1}{\sqrt{2\pi}\sigma} exp(\frac{1}{2n} \sum_{i}^{n} \frac{|x_{i}^{obs} - \hat{x}_{t_i}(\theta)|^{2}}{\sigma^{2}})
  • l(\theta) = \frac{1}{2} log(2\pi\sigma^{2}) + \frac{1}{2 \sigma^{2}} \frac{1}{n} \sum_{i}^{n} |x_{i}^{obs} - \hat{x}_{t_i}(\theta)|^{2}

各細胞の状態パラメータ ( k, t) と遺伝子発現のパラメータ ( \alpha^{(k)}, \beta, \gamma) を解析的に推定することは困難なので、本手法では以下のようなEMアルゴリズムを使う。

  • E step: \thetaを固定して tを最適化する
  • M step: tを固定して t_{0}^{k} (パラメータkの最適化と同義)と \thetaを最適化する
    • Nelder-Mead法で最適化をするらしい *6

引用元：Bergen et al. “Generalizing RNA velocity to transient cell states through dynamical modeling” bioRxiv. 2019. https://doi.org/10.1101/820936 CC-BY 4.0

www.nature.com

Bergen et al. 2021. Mol Sys Biol

同グループから最近発表されたレビュー論文。Fig. 2で紹介されているRNA velocityがうまく機能しない場合の事例が興味深い。以下のようなケース・遺伝子はうまく上記のモデルが適用できないらしい。

  1. 細胞種ごとに遺伝子発現のKineticsが異なる場合
    • 一般性が高そうな印象。全ての細胞を使ってKineticsを推定するのではなく、ある程度状態遷移の流れを見たい細胞集団に絞って解析をする方が良いと思われる。
  2. Transcriptional boost (ある状態時に発現が急上昇する) の場合
    • 時間によらず一定と仮定されているKinetics *7 が一定でない場合はうまくいかない
  3. 安定した状態の細胞しかいない場合
    • PBMCのような分化済みの細胞が中心の場合は発現変動のKineticsを推定しようとしてもランダムな結果が出やすい

引用元：Bergen et al. “RNA velocity—current challenges and future perspectives” Mol Syst Biol. 2021. doi: 10.15252/msb.202110282 CC-BY 4.0

pubmed.ncbi.nlm.nih.gov

使い方

ツールvelocyto.pyを使って、本解析のインプットであるUnspliced mRNAとSpliced mRNAのリード数をカウントする。ChromiumのデータであればCellrangerのマッピング結果(bam file)をそのまま転用することができる。リピート領域のデータはUCSC genome browserでダウンロード可能。GTFファイルはCellrangerのリファレンスゲノムに付随しているものを使うことを推奨。velocyto run10xというオプションも用意されているが、"This is an older version of cellranger"というエラーが出てきたので、開発者が推奨*8するようにvelocyto runを使用した。

velocyto run \
   --mask repeat_mask.gtf \  # リファレンスゲノムのリピート領域の情報
   --samtools-threads 4 \      # samtools sortで使うスレッド数
   --samtools-memory 16000 \    # samtools sort 1スレッドあたり使うメモリ数 (MB)
   path/to/cellranger_output \ # 入力データ Cellrangerの出力 `outs`などを含むフォルダを指定する
   path/to/gtf_file                    # gtfファイル

計算が終わるとSpliced/Unspliced mRNAのカウントデータが格納されたloom形式のファイルが得られる。これを解析するためのツールもvelocyto.pyに用意されているが、ツールの使いやすさと機能の充実度 *9 を考えるとscVeloを使う方が便利な印象。ツールの使い方については本家のTutorialが充実しているのでそちらを参照。

*1:細胞Aは未分化の細胞であるから状態遷移の起点にあるべき、というような解釈

*2:スプライシングスピードあたりの転写スピード、mRNA分解スピードを求める形になる

*3:Supplementary Note 2, Fig. 2

*4:対象の細胞種・細胞分化で発現が上昇・下降し続ける遺伝子など

*5:遺伝子発現の変動は転写調節による影響が強いはずで、オリジナル論文の仮定には明らかに無理がある、という彼らの主張だと思われる。

*6:https://www.kimoton.com/entry/20200624/1592944325が分かりやすそう

*7:Bergenらのモデルでは4状態を仮定

*8:https://github.com/velocyto-team/velocyto.py/issues/121

*9:上述のモデルをすべて使うことが可能

シングルセル解析⑥： Ambient RNA contamination除去

Bioinfo Single-cell

はじめに

10X ChromiumをはじめとするDroplet-basedのシングルセル解析を行う場合にAmbient RNA contaminationが問題になることがある。Ambient RNAは細胞に由来しないRNAであり、処理中に壊れた細胞から溶出したRNAや、コンタミネーションなどが由来となる。Ambient RNAはDropletに混入して、あたかも細胞に由来するRNAのように検出されるため、遺伝子発現データを歪める恐れがある。ここではAmbient RNAの影響をデータから除外するための方法についてまとめる。

SoupX

SoupXは2020年にGigaScienceで発表された手法。以下のステップでAmbient RNAの影響を取り除く。

  1. Empty dropletの情報を使ってAmbient RNAのプロファイルを推定する
    • 総UMI数が閾値よりも少ないDroplet*1をEmpty dropletと仮定*2
    • Empty droplet由来のUMI数からAmbient RNAの遺伝子プロファイル b_.を求める
      • b_g = \dfrac{\sum_{c \in D} n_{gc}}{\sum_{g \in G} \sum_{c \in D} n_{gc}}
      • D: Empty droplets
      • G: 全遺伝子
      • n_{gc}: Droplet cにおける遺伝子gのUMI数
  2. 細胞種特異的なExclusive marker geneset \Omegaを同定
    • 「細胞種Aでは発現していない遺伝子B」のような遺伝子セット*3を決める
    • 可能であれば既知情報に基づいてユーザーがマニュアルで設定する
    • マニュアルで決めるのが困難な場合は以下の方法で決める
      • 特定のクラスタに特異性の高い(tf-idfが1以上)の遺伝子を選び、次の基準で \Omegaを求める
      • P_{c} = \sum_{k=n_{gc}}^{\infty} \dfrac{\lambda^{k}e^{-\lambda}}{k!}
      • \lambda = \rho_{max}b_{g}n_{.c}: Ambient RNAが \rho_{max}混入している場合の遺伝子gの平均UMI数
      • \rho_{max}: Ambient RNAの最大混入量 (Default: 1.0)
      • ある細胞クラスタに含まれる全細胞の P_{c}が0.05以上となる遺伝子セットをExclusive marker \Omegaとする
  3. Ambient RNAの混入率 \rhoを細胞ごとに推定する
    • Exclusive marker geneの発現情報をもとに推定する
      • \rho = \dfrac{\sum_{g,c \in \Omega} n_{gc}}{\sum_{g,c \in \Omega} n_{.c} b_{g}}
  4. 各細胞の発現プロファイルを補正する

引用元：Young et al. “ SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data” GigaScience. 2020. https://doi.org/10.1093/gigascience/giaa151. CC-BY 4.0

https://academic.oup.com/gigascience/article/9/12/giaa151/6049831academic.oup.com

DecontX

DecontXは2020年にGenome Biologyで発表された手法。トピックモデリングのような以下の階層ベイズモデルを使ってAmbient RNAを分離する。入力データはUMIカウントデータで、細胞のクラスタリングを事前情報として与える必要がある。

  • p(X,Z,Y,\theta|\phi,a_{1},a_{2}) = \Pi_{j=1}^M p(\theta_{j} | a_{1},a_{2}) \Pi_{t=1}^{N_{j}} ( (p(y_{jt}=1 | \theta_j) p(x_{jt}=g | \phi_{z_j}))^{I(y_{jt}=1)} ( (p(y_{jt}=0 | \theta_j) p(x_{jt}=g | \phi_{-z_j}))^{I(y_{jt}=0)} )

  • k: 細胞種 (細胞クラスタ) の数

  • M: 細胞数
  • N_j: 細胞jのUMI数
  • p(y_{jt}=1 | \theta_j) = Bernoulli(\theta_j): 細胞jのUMI tが細胞に由来する確率
  • p(x_{jt}=g | \phi_{z_j}) = Categorical(\phi_{z_j}): 細胞jのUMI tが細胞種 z_jに由来するとき、遺伝子gに由来する確率
  • p(y_{jt}=0 | \theta_j) = Bernoulli(\theta_j): 細胞jのUMI tがAmbient RNAに由来する確率
  • p(x_{jt}=g | \phi_{-z_j}) = Categorical(\eta_{z_j}): 細胞jのUMI tが細胞種 z_j以外のAmbient RNAに由来するとき、遺伝子gに由来する確率
  • \eta_{z_j}) = \sum_{k \neq z_{j}} w_k \phi_{k}: Ambient RNAに含まれる各遺伝子の割合
  • \phi \sim Dir(\beta) : パラメータ \phiの事前分布
  • \theta \sim Beta(\alpha_1, \alpha_2): パラメータ \thetaの事前分布

引用元：Yang et al. “Decontamination of ambient RNA in single-cell RNA-seq with DecontX” Genome Biol. 2020. https://doi.org/10.1186/s13059-020-1950-6. CC-BY 4.0

genomebiology.biomedcentral.com

CellBender remove-background

2019年にbioRxivに公開された方法。以下のようにUMIカウントデータをモデリングしてバックグラウンドの推定を行う。この手法ではバックグラウンドがAmbient RNAに由来するだけでなく、Barcode swapping (PCRでキメラ配列が形成されることで誤ったバーコードが付与される) についても考慮している点が特徴である。

  • \chi_{ng} = NN_{\chi}(z_{n}): Droplet nにおけるノイズ除去後の遺伝子gの割合 (NN=Neural network)
  • z_{n} \sim Normal(0, 1): Droplet nを表す潜在ベクトル (N(0, 1)の正規分布を事前分布に設定)
  • d_{n}^{drop} \sim LogNormal(d_{\mu}^{drop}, d_{\sigma}^{drop}): Size factor (=バックグラウンド由来のUMI量)
  • d_{n}^{cell} \sim LogNormal(d_{\mu}^{cell}, d_{\sigma}^{cell})]: Size factor (=細胞由来のUMI量)
    • d_{\mu}^{drop}, d_{\sigma}^{drop}, d_{\mu}^{cell}, d_{\sigma}^{cell}: ハイパーパラメータ *5
  • y_n \sim Bernoulli(p): Droplet nに細胞が含まれていたか否か
    • p: 予想解析細胞数をもとに決める
  • \rho_n \sim Beta(\rho_{\alpha}, \rho_{\beta}): Droplet nにおける他の細胞由来 (上述のBarcode swappingによる) UMIの割合
  • \phi \sim Gamma(\phi_{\alpha}, \phi_{\beta})
  • c_{ng} \sim NB( (1-\rho_n)(y_{n}d_{n}^{cell} \chi_{ng} + d_{n}^{drop} \chi_{g}^{a}) + \rho_n (y_{n}d_{n}^{cell} + d_{n}^{drop}) \bar{\chi_{g}}, \phi): 細胞nにおける遺伝子gのUMI数
    • \bar{\chi_{g}}: 全細胞での遺伝子gの平均割合 (Barcode swappingでは他の全細胞からランダムに転写物が混入するので)
    • \chi_{g}^{a}: Ambient RNAにおける遺伝子gの割合 (これは推定する)

引用元：Fleming et al. “CellBender remove-background: a deep generative model for unsupervised removal of background noise from scRNA-seq datasets” bioRxiv. 2019. https://doi.org/10.1101/791699 . CC-BY 4.0

www.biorxiv.org

*1:論文では10 UMI以下

*2:SoupXのデフォルトは100 UMIが閾値に設定されている

*3:例示されていたのはIg関連遺伝子やヘモグロビン関連遺伝子

*4:m_gcが負の場合は0とする

*5:データから自動決定するらしい

シングルセル解析⑤：ディープラーニングはシングルセル解析に有用なのか

Single-cell Bioinfo

はじめに

シングルセル解析にディープラーニングを利用して、データのノイズ除去(≒インピューテーション)、細胞種の推定、バッチ補正などを行う事例が出てきている。本当に有用であるか些か疑問に感じていたので*1、まとめたいと思う。

研究事例

scVI

2018年にNir Yosefのグループが発表。VAE (Variational autoencoder)をシングルセルデータへ適用した最初期の手法。シングルセルデータのDenoising (Dropoutの除去)、バッチ補正、次元圧縮を実施できる。

www.nature.com

モデル

入力は遺伝子カウントデータ x_{ng} とBatch ID (どの実験・サンプルに由来するか) s_n *2。 l Encoder部分でSize factor l_n (各細胞のSequence depth) と潜在変数 z_n の平均(mean)と標準偏差(std)をNNで計算。

Decoder部分では各遺伝子の頻度 w_{ng} とドロップアウト率 h_{ng} をNNで計算。 w_{ng}・ h_{ng}・ l_nをパラメータとするZINB*3モデルで入力の遺伝子カウントデータをモデル化。

  • z_n \sim Normal(0, I)
  • l_n \sim logNormal(l_{\mu}, l_{\sigma})
    • l_{\mu} : 1細胞あたり総カウント数の平均
    • l_{\sigma} : 1細胞あたり総カウント数の標準偏差
  • \rho_{n} = f_w(z_n, s_n))
    • s_n: Batch ID of cell n (該当するBatchのindexのみ1、その他は0のOne-hot vector)
    • f_w: Neural network
  • w_{ng} \sim Gamma(\rho_{ng}, \theta_{g})
  • y_{ng} \sim Poisson(l_n w_{ng})
  • h_{ng} \sim Bernoulli(f_h^g(z_n, s_n))
    • f_h: Neural network
  • x_{ng} = y_{ng} (h_{ng} = 0)
  • x_{ng} = 0 (h_{ng} = 1)

引用元：Lopez et al. “Bayesian Inference for a Generative Model of Transcriptome Profiles from Single-cell RNA Sequencing” bioRxiv. 2018 https://doi.org/10.1101/292037. CC-BY 4.0

示したデータ

  • 7種類のシングルセルデータを使った精度評価
  • Imputation/Denoisingの精度評価
    • 一部の値を0に置換したシングルセルデータを使いImputaionの評価を実施
    • ZINB distributionでモデル化した手法が良い成績を残したことを紹介 (ディープラーニングベースであるかに関わらず) (Fig S3)
    • その一方でゼロのほとんどはNegative binomial distributionで説明でき、Zero inflationがほぼ認められないことも報告している
  • 次元圧縮手法としての有用性
    • 人力による細胞種アノテーションが潜在空間に反映されているかをSilhouette widthやARIで評価
    • 他の次元圧縮手法と比較しても遜色のない結果であることを示している (Fig S9)
    • SIMLRと違いscVIで得られた潜在空間は細胞分化のTrajectoryを反映していることを一部のデータで示している (Fig. 2)
  • バッチ補正手法としての有用性
    • MNNベースの方法よりも良好 (Fig. 2)

メモ

  • 一つ一つのタスク (Denoising・次元圧縮・バッチ補正) は既存の手法を超えるパフォーマンスではない (同等程度)
  • 単一のモデルで上記のことを一度にある程度の精度で実現できることが優れた点なのだろう

scGEN

2019年にFabian Theisグループが発表。VAEで得られる潜在空間がある摂動に対する細胞の反応の予想に使えることを主張した論文。イメージとしては、潜在空間で「細胞A(摂動あり) - 細胞A(摂動なし) + 細胞B(摂動なし) 」を計算すると「細胞B(摂動あり)」に近い細胞状態が得られる、という内容。

www.nature.com

引用元：Lotfollahi et al. “Generative modeling and latent space arithmetics predict single-cell perturbation response across cell types, studies and species” bioRxiv. 2018 https://doi.org/10.1101/478503. CC-BY 4.0

モデル

示したデータ

  • 潜在空間上で「摂動あり」と「摂動なし」の細胞が線形分離できることを2つ (IFNb刺激、細菌感染) のデータで示した (Fig S1, Note S1)
  • "IFNb刺激あり"の細胞種Aだけを学習から除外し、本モデルで細胞種AにおけるIFNbの影響 (発現変動) を推定できるか評価 (Fig. 2)
    • 細菌感染腸管上皮細胞のデータでも同様の結果 (Fig. 3)
    • 細胞種特異的なIFNbへの反応も予測可能 (Fig. 2f)
    • CVAE (Conditional VAE)やStyle transfer GANとの比較を実施 (Fig. S2-4)
  • Data Aの"刺激あり/刺激なし"のデータを使って、"刺激なし"のデータのみしか存在しないData Bから"刺激あり"の状態を予測 (Fig. 4)
    • Organism AのデータからOrganism Bのデータを予測する問題にも挑戦 (Fig. 5)

メモ

  • 面白い試みだが、どの程度一般化できる内容であるかは不明
  • Fig. 4とFig. 5で完全に実現可能性・有用性が評価できているのか疑問に感じた。

SAVER-X

Nancy R Zhangグループが2019年に発表した手法。シングルセルデータで学習したVAEモデルが新規のシングルセルデータのDenoisingに使えることを示した。

引用元：Wang et al. “Data Denoising with transfer learning in single-cell transcriptomics” bioRxiv. 2019 https://doi.org/10.1101/457879. CC-BY 4.0

www.nature.com

モデル

入力は観測された遺伝子カウントデータ y_{gc} のみ。以下の階層モデルで細胞cにおける遺伝子gの割合 x_{gc} を推定する。scVIやscGenとは違いZero inflationを考慮しないNegative binomial分布でモデル化している。

  • \Lambda_{.c} = f(y_{.c})
    • f: Autoencoder
  • x_{gc} \sim Gamma(\Lambda_{gc}, \theta_{g} \Lambda_{gc}^2)
  • y_{gc} \sim Poisson(l_c x_{gc})
    • Size factor l_c = \sum_{g} y_{gc}

示したデータ

  • HCAと10X GenomicsのPBMCデータで事前学習したモデルを使ってPBMCデータのDenoisingを実施 (Fig. 2)
    • Denoising後の遺伝子発現量をCITEseqで得られるタンパク発現量と比較、他のDenoising手法よりも良い結果 (Fig. 2c)
  • マウス脳シングルセルデータで学習したモデルでヒト脳シングルセルデータのDenoising、種を超えた転移学習の可能性を示唆 (Fig. 3)

メモ

  • 従来手法 (e.g. MAGIC, scImpute, DCA) は入力データのみを利用するDenoisingがほとんど。NNを利用した転移学習により過去のデータを利用して良い結果が出せる点は重要。
  • テストケースが少ない (2例) なので結果の一般性については疑問

totalVI

2021年にNir Yosefグループが発表。表面タンパクと遺伝子発現を同時計測するCITEseqデータを対象とするVAEモデルを提案。複数の異なる種類のデータを同時に用いて次元圧縮・クラスタリング・Denoisingすることが可能。

www.nature.com

モデル

同ラボが開発したscVIをベースとしたモデル。入力は遺伝子カウントデータ x_{ng}とタンパクカウントデータ y_{nt}とBatch ID s_n *4。

Encoder部分で遺伝子とタンパクのSize factor l_n・ \beta_nと潜在変数 z_n の平均(mean)と標準偏差(std)をNNで計算。

Decoder部分は各遺伝子の頻度 \rho_{ng} をNNで計算。 w_{ng}・ l_nをパラメータとするNegative binomialモデルで入力の遺伝子カウントデータ x_{ng}をモデル化*5。 タンパクデータに関してはバックグラウンドとフォアグラウンドのシグナルが混在していると仮定したモデル。学習後にフォアグラウンドのシグナルを推定することが可能。

  • z_n \sim LogisticNormal(0, I)

RNA-seq

  • l_n \sim logNormal(l_{\mu}^T s_n, l_{\sigma}^T s_n) *6
    • l_{\mu}^T s_n : Batch s_nにおける1細胞あたり総カウント数の平均
    • l_{\sigma}^T s_n : Batch s_nにおける1細胞あたり総カウント数の標準偏差
    • s_n: Batch ID of cell n (該当するBatchのindexのみ1、その他は0のOne-hot vector)
  • \rho_{n.} = f_w(z_n, s_n)) #
    • f_w: Neural network
  • w_{ng} \sim Gamma(\theta_{g}, l_n \rho_{ng})
  • x_{ng} \sim Poisson(w_{ng}): UMI counts

Protein

  • \beta_n \sim logNormal(c_{\mu}^T s_n, d_{\sigma}^T s_n): Background signalの平均
  • v_{nt} \sim Bernoulli(\pi_{nt})
  • \pi_{nt} = h_{\pi}(z_n, s_n)): Background signalの割合
    • h_{\pi}: Neural network
  • v_{nt} \sim Bernoulli(\pi_{nt})
  • r_{nt} \sim Gamma(\phi_{t}, v_{nt} \beta_{nt} + (1 - v_{nt}) \beta_{nt} \alpha_{nt})
    • \beta_{nt} \alpha_{nt}: Foreground signalの平均
  • \alpha_{nt} = g_{\alpha}(z_n, s_n))
    • g_{\alpha}: Neural network (出力が1以上になるような条件付き)
  • y_{nt} \sim Poisson(r_{nt}): UMI counts for protein t

引用元：Gyaoso et al. “Joint probabilistic modeling of paired transcriptome and proteome measurements in single cells” bioRxiv. 2020. https://doi.org/10.1101/2020.05.08.083337. CC-BY 4.0

示したデータ

  • 表面タンパクデータのDenoising
    • Foreground probabilityが陽性細胞を分離する上で有用 (Fig. 2a-f)
  • 異なる抗体パネルのデータも統合可能 (Fig. 3)
    • 通常のscRNA-seqデータをCITEseqデータと統合し、タンパク発現量をImputationできることを示している (Fig. 3g)
  • Archetypal analysis

メモ

  • タンパク発現量の前処理でよく使われているCLR normalizationよりもいくつかのメリットがある
    • BackgroundとForeground signalの分離が可能
    • 異なる抗体パネルのデータを比較可能 (CLRはCompositional dataと仮定するので計測項目が違うと比較できない)
    • 一方で、CLRと違ってSequence depthが考慮されない気がするが、この辺りはどうなのだろうか？Denoising後の数値は細胞間で単純に比較できないのでは？
  • あるオミクスデータから別のオミクスデータを生成できる可能性を示しているのは興味深い

scArches

2021年にFabian Theisグループが発表。シングルセルデータ用VAEモデルに対する効率的な転移学習の方法を提案。 scArchesは先行研究で提案されているモデルを転移学習に利用する方法を提案している。Nir Yosefグループとの共同開発のようで、ツールscArchesを使えば、TheisグループとYosefグループがこれまでに開発したモデル (scVI, totalVI, scgen, trVAE)を転移学習へ応用できる。Yosefグループが開発するscvi-toolsでもすでにscArchesの機能は実装されており、利用することができる。*7

www.nature.com

手法

scArchesではArchitecture surgeryという転移学習の方法を提案している。

データ1~N (Reference data)を使って学習したモデルを転移学習へ利用する場合を考える。このとき上記のモデル(e.g. totalVI)ではBatch ID tex: s_i (i=1~N)をEncoderとDecoderの入力に用いる。Architecture surgeryでは、学習済みモデルへ新たなデータ(N+1)~M(Query data)を適用・再学習する場合に、このOne-hot vectorに新たな項目 s_i (i=(N+1)~M)を追加する。このときパラメータの更新は、Encoder/Decoderの第一層の s_i (i=(N+1)~M)に紐づいているパラメータのみ行い、そのほかのパラメータについてはReference dataで学習した時点から変更はしない。これによってある程度の性能が示せることを論文(Fig. 2)で報告している。

引用元：Lotfollahi et al. “Query to reference single-cell integration with transfer learning” bioRxiv. 2020. https://doi.org/10.1101/2020.07.16.205997 . CC-BY 4.0

メモ

  • 同じ転移学習をテーマにしているSAVER-Xとは違いバッチ間差異を考慮して、学習済みモデルを新たなデータへ適用できることが特徴
  • totalVIで学習した他モダリティ(Protein)のデータを生成した事例(Fig.4)は有用だと思うが、実用性についてはデータが不十分な印象を受けた

*1:なんでもかんでもディープラーニング使っとけばよいという風潮を筆者が嫌いなだけ

*2:nは細胞、gは遺伝子のindex

*3:Zero-inflated negative binomial、ドロップアウトを含むシングルセルデータのモデリングによく用いられる

*4:nは細胞、gは遺伝子、tはタンパクのindex

*5:ZINBを選んでいない

*6:scVIと違いバッチごとの平均・標準偏差を使うようにしている

*7:こちらの方がマニュアル等充実しており使いやすそうな印象も受ける

マイクロバイオーム解析②： 16Sデータの前処理について

Bioinfo Microbiome Omics

はじめに

16Sアンプリコンシーケンスデータの前処理・一次解析のプロセスを自動化しようとしたら、色々とわかりにくいポイントがあったので以下にまとめる。

シーケンシングまでの過程

データ前処理のプロセスを最適化するには、アンプリコンシーケンスの実験過程を理解する必要がある。ここではIlluminaのシーケンサーで解析をする場合を考える。この場合、以下のようなライブラリを作製する。

f:id:auroratummy:20210802230652p:plain:w400

P5とP7はフローセルに結合するために必要な領域である。P5とP7がフローセルに結合することでブリッジPCRによるクラスタ形成 *1 が行われる。

R1SPとR2SPはシーケンスプライマー配列である。ブリッジPCRによる増幅の後に、R1SP・R2SPに対応するプライマー配列を入れ、Sequencing By Synthesis法を行うことで、Insert DNAの両端を解読することができる。

f:id:auroratummy:20210802232512p:plain:w400

In1はインデックス (バーコード) 配列である。In1にサンプルごとに異なる配列を入れ、それを解読することで、シーケンスがどのサンプルに由来しているかをデータから判別することができる。上図はR2SP側にのみバーコードが入っているが、R1SP側・R2SP側の両方にバーコードを入れるDual barcodingもよく行われる。

f:id:auroratummy:20210802232812p:plain:w400

アンプリコンシーケンスのライブラリ作製

アンプリコンシーケンスの場合は上記のようなライブラリを作製するためにtailed PCRを行う。Illuminaの資料によると、いくつかの方法があるようだが、そのうちtailed PCRを2ステップ挟む方法は以下のような流れで行われる。

f:id:auroratummy:20210803000149p:plain:w400

上図の場合はR1SPとR2SPに挟まれた領域にプライマー領域が含まれるので、プライマー領域を先頭に含んだ状態で配列が解読される。

したがって、この場合、解読される配列の先頭部分の多様性が低くなる。配列の先頭部分の多様性が低いと、フローセル上でクラスタを同定する精度が悪くなり、結果的に最終的な配列情報の質が低くなる*2。この問題を改善するために以下のような方法が採られることがある。

  1. PhiXのようなSpike-inを添加して配列の多様性を上げる方法
  2. 多様性を担保した配列 (Heterogeneity spacer) をPCRで挿入する方法
  3. Sequencing primerでプライマー領域をカバーする方法

2番目の方法はHeterogeneity spacerと呼ばれるランダムな配列を挟んだプライマーを1st PCRで使うことで、解読される配列の先頭部分の多様性を高めるものである。

f:id:auroratummy:20210803003144p:plain:w400

3番目はEarth microbiome project (EMP) などのプロトコルで採用されている方法で、以下のようにプライマー領域までSequencing primerでカバーしてしまうことで、配列の先頭部分をターゲット配列にするものである。

f:id:auroratummy:20210803052751p:plain:w400

16Sデータ前処理の自動化

DADA2によるDenoisingを行う場合は以下のポイントに留意する必要がある。

  1. Demultiplexing
  2. アダプタ配列 (e.g. プライマー) の除去
  3. 配列のクオリティを確認 & 低クオリティ領域の除去

Demultiplexing

上記のように、In1に入るバーコード領域をサンプルごとに変えることで、複数のサンプルを同時に解析できる。この場合はデータをサンプルごとに分割するステップ (Demultiplexing) が必要となる。このステップに関しては、バーコードの情報さえ用意すれば、QIIME2のドキュメント通りに実施することができる。

しかし、Forward/Reverse両方のリードにバーコード配列を付与するDual barcodingにはQIIME2は現段階で非対応なので、より汎用的なパイプラインを作るのであればIlluminaのbcl2fastqでDemultiplexingを行った方が無難なように感じる。

アダプタ配列の除去

DADA2でオプション--p-trim-left-f Nと--p-trim-left-r Nを指定するとForward readとReverse readの先頭N文字を除くことができる。除外すべき配列のサイズが分かっている場合はそれを指定すれば良い。

しかし、上述のようにデータによってはHeterogeneity spacerが入っていたり、EMPのようにプライマーを読まない手法が使われていたりするので、自動化する上では不便な設計になっている。

DADA2の公式ドキュメントには以下のような記載がある。これを踏まえるとcutadaptやtrimmomaticでアダプタ配列の除去を行う方が一般化しやすいように感じる。

If primers are at the start of your reads and are a constant length (a common scenario) you can use the trimLeft = c(FWD_PRIMER_LEN, REV_PRIMER_LEN) argument of dada2’s filtering functions to remove the primers. For more complex situations (e.g. heterogenity spacers, variable length amplicons such as the ITS region) we recommend the cutadapt tool or the trimmomatic tool.

cutadaptはQIIME2で使うことができる。cutadaptでプライマー配列を除去する場合は以下のようなコマンドで実行できる。以下のコマンドでcutadaptは${FORWARD_PRIMER} ${REVERSE_PRIMER}以前の5'領域を削ることができる。

qiime cutadapt trim-paired \
  --i-demultiplexed-sequences input.qza \ # input
  --o-trimmed-sequences trimmed.qza \ # output
  --p-front-f ${FORWARD_PRIMER} \ # forward primer
  --p-front-r ${REVERSE_PRIMER} \ # forward primer
  --p-error-rate 0.1 \ # maximum allowed error rate
  --p-cores 1

クオリティの確認と低クオリティ領域の除去

DADA2では3'端の低クオリティ領域をトリムすることを推奨している*3。DADA2でオプション--p-trunc-len-f Nと--p-trunc-len-r Nを指定すると長さNになるように全リードの3'端をTruncationすることができる。

以下のコマンドで配列のクオリティを確認するとFASTQCで出力されるような配列のクオリティ情報が出力される。この情報をもとに上記のパラメータNを決める。

# クオリティ情報の計算
qiime demux summarize \
  --i-data input.qza \
  --o-visualization quality.qzv

# 可視化
qiime tools view quality.qzv

f:id:auroratummy:20210803074954p:plain:w400

ちなみにQIIME固有の可視化ファイル(.qzv)は上のようにqiime tools viewを実行するとブラウザ上で閲覧することができる仕様になっているが、unzipするとhtmlファイルが出力される*4ので、そちらから結果にアクセスする方が楽。

# 解凍
unzip quality.qzv

QIIME 2のドキュメントやコミュニティページによるとこのプロセスは特に標準化されておらず、解析者の主観に委ねられている印象を受ける。例えば以下のスレッドでは、パラメータ--p-trunc-len-についてQIIME 2開発者がコメントしているが、それぞれバラバラなことを言っており、統一見解は存在しない模様。

DADA2では配列のTruncationをした後にForward/Reverse readsをマージするステップがあるが、Truncationによりリードが短くなりすぎるとマージできなくなる (マージのステップで大部分のリードが失われる) ので注意する必要がある。リードのマージには20bp以上の重複が推奨される*5とのこと(その一方でDADA2のマニュアルには12bpと書かれている)*6。あまりにシーケンスのクオリティが悪くてどうしようもない場合は、Reverse readの情報を捨てて、Forward readだけで解析することが推奨されるようだ。

*1:参考: https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2013_illumina_techsupport_session19.pdf

*2:参考: https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/webinar/2019_ilmn_webinar_Sequencing%20for%20low%20diversity%20sample_20191120.pdf

*3:https://benjjneb.github.io/dada2/tutorial.html

*4:https://docs.qiime2.org/2021.4/tutorials/qiime2-for-experienced-microbiome-researchers/#pro-tips-for-power-users

*5:https://forum.qiime2.org/t/dada2-truncation-lengths-and-features-number/1940/8

*6:https://docs.qiime2.org/2021.4/plugins/available/dada2/denoise-paired/

シングルセル解析④: Compass: Flux balance analysisを使った代謝活性予測

Single-cell Bioinfo

シングルセルデータを使った代謝活性予測

先週のCellに以下の興味深い論文が出ていたので手法を読み解いていく。この論文では、Flux balance analysisをscRNA-seqデータへ応用することで各細胞の代謝状態を推定するCompassというツールを開発し、Th17のサブセット間で代謝状態が如何に異なるかを推定、評価している。

doi.org

github.com

Flux balance analysis (FBA)

FBAは化学量論 (Stoichiometry) と代謝流量 (Metabolic flux) のみに着目して代謝状態をモデリングする方法。

Stoichiometryとは「Acetyl-CoA + L-Glutamate <=> CoA + N-Acetyl-L-glutamate」のような反応式で表される、代謝反応の入力・出力の量的関係にまつわる情報を指す。

FBAでは対象生物種の全ての代謝反応を同時にモデリングする。上述の論文ではRecon 2というヒト代謝マップのリファレンスデータを使っている。

www.nature.com

定式化

FBAでは全代謝反応のStoichiometryとMetabolic fluxを下図Bのように行列Sとベクトルvで表す。行列Sの行は代謝物、列は代謝反応に対応しており、各行列成分には反応によって産生される物質量が格納される。例えば代謝反応 r_{i}が「A + D → B」という反応であれば、産生される代謝物Bの行に1、消費される代謝物AとDの行に-1が格納される。

FBAでは、対象となる代謝マップを持つ生物において、以下の条件・仮定をおいた場合の「任意の代謝活性の最大値 (または最小値)」を求める。例えば、下図の代謝物Eの産生量 (細胞外への放出量) の最大値を知りたい場合は、下記の条件・仮定下での反応 r_7の流速 v_7の最大値を求める問題を解く。

仮定・条件
  1. 全ての代謝物の産生量と消費量が一致する状態 (定常状態) を仮定： Sv = \vec{0}
  2. 反応特性など代謝の事前情報に基づく条件：e.g. 反応 r_3が不可逆である場合 v_3 \geqq 0
  3. 実験条件： e.g. 反応 r_3を担う遺伝子Xを欠失した場合 v_3 = 0

引用元：Cuevas et al. “From DNA to FBA: How to build your own genome-scale metabolic model” Front Microbiol. 2016 Jun 17;7:907. doi: 10.3389/fmicb.2016.00907. CC-BY 4.0

上記の条件を満たすvは一意には決まらない。これはすなわち対象とする生物種の代謝状態に複数の定常状態があることを意味するが、FBAでは、この複数の定常状態の中から、目的関数が最大・最低になる場合の状態を求める。

線形計画法

これらの目的関数・制約条件は全て1次式なので、FBAでは上記の最適化問題を線形計画法 (Linear programming) を用いて解く。例えば次のような2変数の最適化問題を考える。

  • 目的関数
    • f(x,y) = x + y
  • 制約条件
    • 2x + y \leqq 6
    • x + 2y \leqq 6
    • x \geqq 0
    • y \geqq 0

このとき下図のように問題を捉えることができる。点線で囲まれた斜線部は制約条件を満たす領域である。オレンジの直線はf(x,y)の等高線である。このような領域の中で、ある線形な目的関数 f(x,y) = x + yが最大・最小となる点はいずれかの頂点であることがわかる。

f:id:auroratummy:20210713011111p:plain:w300

線形計画法ではこのような目的関数の最大・最小値に該当する頂点を効率よく探索する方法が提案されているのだが、詳しいアルゴリズムに関しては以下の参考文献がわかりやすかった。CompassではCplexというIBMの最適化問題用のPythonライブラリを用いてFBAを実施している。*1

はじめての最適化 | 関口 良行 | 数学 | Kindleストア | Amazon

Compass

ここで肝心の論文に話を戻す。この論文では、上述したFBAと、scRNA-seqでの代謝関連遺伝子の発現情報を利用して、各細胞の代謝状態を推定する方法を提案している。

最大流速 v_r^{opt} を求める

最初のステップでは単純なFBAにより各代謝反応 v_rの最大流速 v_r^{opt}を求める。条件(i)は定常状態、(ii)は各代謝反応の可逆性*2・制限*3を表す。ここでは遺伝子発現データは考慮していないため、全ての細胞で v_r^{opt}は等しい。

引用元：Wagner et al. “In Silico Modeling of Metabolic State in Single Th17 Cells Reveals Novel Regulators of Inflammation and Autoimmunity” bioRxiv 2020.01.23.912717; doi: https://doi.org/10.1101/2020.01.23.912717 CC-BY 4.0

遺伝子発現を考慮した最適化

遺伝子発現データを元に以下の代謝活性行列 Rとペナルティ行列 Pを求める。ペナルティ行列の行は代謝反応、列は細胞に対応する。各成分 P_{ij}を以下の通り求める。

  • R_{ij} = log(1 + X_{ij})

  • P_{ij} = 1/(1 + R_{ij})

X_{ij}は代謝iに対応する遺伝子の発現量を表す。*4

ペナルティ P_{ij}は代謝に対応する遺伝子が発現が弱いほど1、強いほど0に近づく(論文ではResistance (抵抗)と表現されている)。ここでは各代謝 r、各細胞 cごとに y_{rc} = P_{c.} \cdot vを最小化する以下の最適化問題を解く。(iv)は対象とする代謝rの流速が最大流速 v_r^{opt}の \omega倍 *5以上という制約条件なので、代謝r周辺のペナルティ(抵抗)が高いほど y_{rc}は高い値となる。

引用元：Wagner et al. “In Silico Modeling of Metabolic State in Single Th17 Cells Reveals Novel Regulators of Inflammation and Autoimmunity” bioRxiv 2020.01.23.912717; doi: https://doi.org/10.1101/2020.01.23.912717 CC-BY 4.0

最後に C^{raw}を以下の要領で変換したものをReaction scoreとして解析に利用する。

引用元：Wagner et al. “In Silico Modeling of Metabolic State in Single Th17 Cells Reveals Novel Regulators of Inflammation and Autoimmunity” bioRxiv 2020.01.23.912717; doi: https://doi.org/10.1101/2020.01.23.912717 CC-BY 4.0

*1:商用利用には有償ライセンスが必要らしい。PuLPなどオープンソースのものも存在しているのでそれを使おう

*2:a>0であれば不可逆反応

*3:Recon 2における事前知識に基づく

*4:複数遺伝子が一つの代謝に対応する場合は平均値を採用

*5:デフォルトだと0.95倍

データベース調査③： シングルセルデータベース

Bioinfo Omics Single-cell Public data

HCA Data Portal

Human Cell Atlasで開発されたデータベース。2021年6月時点では55種類のヒト組織に由来する92プロジェクトのシングルセルデータが登録されている。

data.humancellatlas.org

HCA Data Portalには基本的にユーザーにより登録されたデータが公開されている。Contributeページによると、対象となるシングルセルデータにはいくつかの基準があるようだ。

  • Consent: 誰でも自由にアクセス可能なデータが優先
  • Sample type: 組織サンプル由来のデータが優先(オルガノイドやCell lineよりも)
  • Health status: 健康人由来データが優先
  • Organism: ヒト由来データが優先
  • Data processing pipeline support: HCAが提供するデータ解析パイプラインに対応しているデータが優先

Exploreページからデータを検索できるが、実際に健康なヒトの組織に由来するデータが多く含まれている印象。メタデータは独自の辞書で標準化しており*1、疾患名や実験条件(e.g. シーケンサー、ライブラリ作製手法、サンプルの保存状態)などで比較的フレキシブルにデータを抽出可能。

ダウンロードできるデータの種類も豊富である。遺伝子発現データ(loom形式)、メタデータ(csv形式)に加えて、生シーケンスデータ(fastq形式)や、HCAのパイプラインで処理されたデータであればマッピングの結果(bam形式)も提供されている。APIでのダウンロードも可能であるらしい。*2

HCAではSmart-seq等の転写産物の全長を読む手法のパイプラインと、Microfluidics/Dropletベースで転写産物の3'端を読む手法(10X Chromium, inDrop, Drop-seq)用のパイプライン(Optimus pipeline)が開発・利用されている。パイプラインはBroad Institute発のワークフロー言語であるWDLで記述されたものが公開されているので、Cromwellなどのワークフローエンジンを使えば誰でもすぐに利用できる。

github.com

Single Cell Expression Atlas

EMBL-EBIが開発したデータベース。2021年6月時点では18生物種に由来する217プロジェクトのデータが登録されている。*3

www.ebi.ac.uk

ArrayExpressやGEOに登録されたシングルセルデータを対象としているた。幅広い生物種のデータが登録されている点が特徴。サンプルや細胞のメタデータ(e.g. 生物種名・疾患名、細胞種アノテーション)はEFO(Experimental Factor Ontology)やCell ontologyのターミノロジーを使うことで標準化されている。一方でデータの検索機能はHCA Data Portalの方が充実しており、生物種名やライブラリ作製法で検索できる一方で、疾患や組織名では検索できない。

簡単なデータ可視化が可能であるところは特徴。例えば、"Browse Experiments"から各プロジェクトのページに飛ぶと、どのような細胞サブセットが含まれているかを簡単に調べることができる。

引用元：Papatheodorou et al. “Expression Atlas update: from tissues to single cells” Nucleic Acids Research, Volume 48, Issue D1, 08 January 2020, Pages D77–D83, https://doi.org/10.1093/nar/gkz947. CC-BY 4.0

データは独自のパイプラインで処理されている。HCAと同様にFull-length/3'用にそれぞれ実装されている。パイプラインはNextflowで実装されており、以下のリンクからダウンロード可能。

github.com

PanglaoDB

カロリンスカ研究所のグループが開発したデータベース。ヒト・マウス由来のシングルセルデータが登録されている。

panglaodb.se

PanglaoDBはSRAから取得したシングルセルデータを対象としている。SRAよりダウンロードした生のシーケンスデータに対して独自のパイプラインでデータを処理。他のパイプラインとの違いは細胞種アノテーションを独自開発したツールalonaで自動的に実施するという点。alonaの細胞種アノテーションのアルゴリズムについては以下の論文に記載されている。*4

academic.oup.com

データ横断的な情報の検索が可能である点が特徴である。Searchページでは、遺伝子名を入力すると、その遺伝子を発現している*5クラスタを全サンプルから探索し、どの細胞種とアノテーションされているクラスタで発現しているかが表示される。

引用元：Oscar et al. "PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data", Database, Volume 2019, 2019, baz046, https://doi.org/10.1093/database/baz046. CC-BY 4.0

Cell type markerページへ飛び、興味のある細胞種の名前を指定すると、彼らが文献調査で収集したマーカー遺伝子と、シングルセルデータを元に計算したマーカー遺伝子の情報にアクセスすることができる。

Ubiquitous index(UI)という指標が計算されている。UIは、PanglaoDBの全クラスタで発現している*6Ubiquitousnessの高い遺伝子が1、ごくわずかなクラスタでしか発現していない遺伝子が0となるような指標である。絶対発現量が高い遺伝子の方がOverestimateされそうではあるが、興味深い。

Single Cell Portal

Broad Instituteで開発されたデータベース。2021年6月現在で、ヒトに限らない様々な生物種に由来する343プロジェクトのデータが公開されている。

HCA Data Portalと目的が被っているような気がするが、もともとはBRAIN (Brain Research through Advancing Innovative Neurotechnologies) initiativeで開発・利用されていたものの、クオリティが高かったからなのか、別のデータにも利用されるようになった、という経緯らしい。

singlecell.broadinstitute.org

フレキシビリティの高いデータ閲覧が可能であり、tSNE/UMAPなどのEmbeddingに加え、遺伝子名を入力すればViolin plot、Dot plot、Heatmapなどが作成可能である。tSNE/UMAP plot上で特定の細胞を選んでアノテーションをすることも可能。

引用元：https://github.com/broadinstitute/single_cell_portal/wiki/Annotations

github.com

REST APIも提供されている。

singlecell.broadinstitute.org

*1:https://data.humancellatlas.org/metadata

*2:https://data.humancellatlas.org/apis マニュアルのリンクが死んでいる

*3:うちヒト由来データは99プロジェクト

*4:そのうち細胞種アノテーションの状況についてもまとめたい

*5:クラスタにおける発現量のMedianが>0であれば発現、という定義らしい

*6:UIの詳細について記述が見つからなかったが、おそらく前述のMedian > 0を「発現あり」とする基準だと思われる

データベース調査②： メタボロームデータ (HMDB, Metabolomics Workbench, MetaboLights etc)

Omics Public data Bioinfo

代謝産物のデータベース

GC-MSやLC-MSなどでヒト資料から取得したメタボロームデータのUntargeted解析*1を実施すると、全ピークのうち既知の物質へ帰属できるものは僅か1.8%しかないという報告がある*2。(素人目に見ても)ある意味宝の山とも言えるメタボロームのデータはどのように整備されているのだろうか。

HMDB

Human Metabolome Database (HMDB)は2007年に公開された代謝物データベースである。名前の通りヒト由来の代謝物に特化したデータベースであり以下の情報を含む。

例: Butyric acid in HMDB

academic.oup.com

MetaboLights, Metabolomics Workbench

MetaboLightsはEMBL-EBIによって開発されたメタボロームデータベースである。Metabolomics WorkbenchはUCSDとNIHが開発したメタボロームデータ。メタボロームデータ(1次データ、2次データ)のレポジトリとして多くの研究者に利用されている。代謝物のアノテーションデータも整理されているが、情報量はHMDBの方が充実している印象を受ける。Metabolomics Workbenchでは公開されているデータを簡単に解析することも可能*4。

academic.oup.com

academic.oup.com

参考資料

理研・津川 裕司先生の講演が統合TVにアップされていて非常に参考になりました。

togotv.dbcls.jp

*1:解析者が対象となる代謝物の情報を与えずにスペクトルを帰属、代謝物を同定する方法

*2:https://www.pnas.org/content/112/41/12549.long

*3:HMDB独自のオントロジーが使われている

*4:触ってみたが、使いこなすのが難しい印象