論文要約シリーズ第1弾！
論文じゃなかったレビューだった。

Single-cell RNA sequencing technologies and bioinformatics pipelines
Experimental & Molecular Medicinevolume 50, Article number: 96

概要

Single-cell RNA sequencing (scRNA-seq) は、がんゲノム分野は勿論のこと、マイクロバイオーム分野でもbulkのRNA-seqと比較して個々の細胞の多様性を観測することができる点で優れており、希少な細胞集団や、遺伝子間相互作用、分化中の細胞系譜等を明らかにすることができる。
同様のセルタイプであってもヘテロ性を持っていることが様々な研究から明らかになっている。このレビューでは、single cellの単離、ライブラリの調整、データ解析に使用できる計算解析ツール等の技術的な難点に焦点を当てる。分子、細胞生物学レベルでのさらなる技術向上及び利用可能なバイオインフォマティクスツールの増加により、シーケンシング技術を基礎科学や医療応用へ役立てることができるだろう。

1細胞単離技術

a ~ cまでの古典的な手法では、インプットが10,000細胞より少ないような場合、有用な情報を得ることが困難だったため、インプットの試料が大量に必要であった。また、FACSの場合、関心のあるタンパクに対するモノクローナル抗体が必須であった。

f:id:kimoppy126:20190209084545p:plain

限界希釈法
細胞濃度を薄めることで単一のクローンを取得する。収率は全プレートの1/3程度。効率・精度共に悪い。
顕微操作 (Micromanipulation)
日本語にすると、顕微操作とでも言うのだろうか。顕微鏡画像から単一細胞の初期杯を取得する。時間がかかる上に、低スループット。
flow-activated cell sorting (FACS)
現状、最もよく使われている手法。マーカー遺伝子の発現量が少ない場合等に特に有用。目的の細胞を蛍光モノクローナル抗体によりタギングすることで、populationの分離が可能となる。蛍光抗体で染色した細胞を流体に混合し、レーザー光などの光を照射することで、細胞を解析するとともに、特定の細胞を無菌的に分取することができる。あるいは、電荷を細胞に付与し、静電偏向を利用した振るい分けを行うことも可能。
レーザーマイクロダイセクション法
コンピュータにより制御されたレーザーを使用し、固体試料から細胞を単離する技術。詳細はレーザーマイクロダイセクション法（櫻田陽右ほか）
マイクロ流体力学を使用
サンプル量が少なく、低コストで済むことから近年利用者が増えている。ナノリットル規模の細胞しか使用しないため、コンタミの恐れが低いことが特徴。初期は生化学のアッセイに使用されていたが、近年はより拡張性の高いアレイも商用販売されており、広く利用されているFluidigm C1は、1つのチップで細胞溶解、RNAの単離、cDNA合成までを最大800細胞について並列に行うことができる。特徴として、擬陽性が低く、バイアスの少ないデータを取得することが可能。
抗体と接合した磁石を使用
がん患者の末梢血に存在する末梢血循環腫瘍細胞(CTC)は、転移の過程を知るのに重要。この細胞の単離には、抗体と接合した磁石を使用する。 Veridex社のCellSearchシステムはFDA承認を受けた唯一の検出系であり、CTCの解析に使用されている。
dropletを使用微小液滴高いキャプチャ効率とスループットにより商用では、10× GenomicsのChromiumシステムが有名。臨床サンプルは

CTCの単離技術に関して、CELLSearch社のかっこいい動画があったので、参考にどうぞ。

www.youtube.com

ライブラリの調整

一般的にライブラリを調整する際は、

細胞の溶解
↓
逆転写反応（cDNAの合成）
↓
cDNAの増幅

という過程を経る。

ここではポリdTプライマーを使用してポリdA配列のセレクションを行うが、転写産物が正常に逆転写される配列は10～20%のみである。
既存のscRNA-seqのプロトコルにおいて、この低いmRNAキャプチャ効率は重要な課題となっており、高い細胞溶解を行うことでこの効率を上げる必要が生じている。

一本鎖cDNAの準備には、低いリボヌクレアーゼ活性と体温性を持った人工的なMoloneyマウス白血病ウイルスの逆転写酵素が一般的に使用される。その後の二本差cDNAの合成には、ポリA鎖を使用した方法、若しくはSMART（Switching Mechanism At 5’End of RNA Template）法が取られる。後者は、均一なカバレッジを得ることができる点で前者より優れている。少量の合成cDNAはその後、PCR、若しくはin vitro transcription反応により増幅される。
in vitro transcription反応を使用すれば、線形にテンプレートを増幅させることができるが、さらに逆転写反応の工程を必要とするため、時間がかかってしまう。これにより3'末端のカバレッジにバイアスが生じることがある。

Smart-seq2では、完全長のcDNA配列を得ることができるため、選択的スプライシングや1塩基多型を使用したアリル特異的な発現を観測したい際に有用である。現在ではシーケンスにイルミナの製品が広く使用されている（Hiseq4000やNextSeq500）。 Miseqを使用すれば、3000万ペアエンドリードのシーケンスを1日で終えることができる。

depthの観点から、トランスクリープ解析には多くの細胞発現プロファイルを必要となる。シーケンスコストを抑えるため、以前の手法は5'若しくは3'のみの転写産物に焦点を当てていた。最近では、unique molecular identifiers (UMIs)と呼ばれる4~8bpのランダム配列を逆転写時に使用することで、PCRにあるバイアスを受けることなく各リードが由来する細胞を特定できるようになった。
しかし、UMIタグを使用した手法は5'若しくは3'のみの配列をシーケンスしているため、アリル特異的な発現やアイソフォームの解析には向いていない。

↓各プロトコルの特徴 f:id:kimoppy126:20190211192252p:plain

scRNA-seqデータ解析の課題

scRNA-seq解析に使用できる実験的な手法は徐々に増えているが、生データを扱えるパイプラインは限られている。
また、商用のものとしては10X genomicsやFluidigmからオープンソースの解析ツールが発表されてはいるが、未だgold-standard的な手法が確立されていない。

データ前処理

クオリティーコントロール(QC)

リードのクオリティーコントロールにはFastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)が使用される。
ここでは、 3′末端の低クオリティー塩基及びアダプター配列を除去する。

アライメント

BulkなRNA-seq解析で使用されていたアライメントツールとして、Burrows-Wheeler Aligner (BWA)や STARが挙げられる。
UMI (Unique Molexular Identifer)が使用されている場合、これらの配列はアライメント前に除いておく必要がある。

RNA-seQCを使用すると、アライメント後の結果のサマリを見ることができる。ここで、RNA劣化及び非効率な細胞溶解を検出する。

発現マトリクスの正規化

アライメント後、各リードはエクソン領域、イントロン領域、遺伝子間領域などのアノテーションが割り当てられる。発現マトリクス(N (cells)×m (genes))の生成には、エクソン領域のリードのみが使用される。 dropoutや偏りのある遺伝子発現のため、scRNA-seqのデータは欠損値の多いデータとなる。この特徴から、細胞特異的なバイアスを除去するために、正規化を行う必要がある。遺伝子特異的な発現量及び細胞特異的なスケール因子に比例する。キャプチャ効率及び逆転写効率を考慮したこれらの確率変数は推定が難しく、一般に固定係数を使用してモデリングが行われる。

正規化には

RPKM (Reads Per Kilobase of exon per Million mapped reads)
FPKM (fragments per kilobase of exon per million reads mapped)
TPM (transcripts per kilobase million)

等が使用される。

RPKM (Reads Per Kilobase of exon per Million mapped reads)

リファレンスにマッピングできた全リード数を $N$ 、そのうち転写産物 $i$ の領域にマッピングされたリード数を $Y_{i}$ 、転写産物 $i$ の長さを $L_{i}$ とすると、転写産物 $i$ の $RPKM_i$ は以下の式で与えられる。

$RPKM_i = Y_{i}\frac{10^{3}}{L_{i}}\frac{10^{6}}{N}$

FPKM (fragments per kilobase of exon per million reads mapped)

FPKMもRPKMと計算式自体は一緒。
ペアエンドシーケンスを行った場合に、片側のリードのみアライメントした場合をカウントする点がRPKMと異なる。
転写産物 $i$ の $FPKM_i$ は以下の式で与えられる。

$FPKM_{i} = Y_{i}\frac{10^{3}}{L_{i}}\frac{10^{6}}{N}$

TPM (Transcripts per million)

TPMは、異なるサンプル間でTPMの合計が均一になる点でRPKMと異なる。
転写産物 $i$ の領域にマッピングされたリードの数を $n_i$ 、長さを $l_i$ としたとき、

${RPK_i} = \frac{n_i}{l_i}$

これを用いて、

$TPM_{i} =\frac{RPK_i}{\sum_{j=0}^nRPK_j} ×10^{6}$

と表される。

ここまでの正規化手法はサンプル内での正規化に使用される。この手法でサンプル間の正規化を行うと、発現量が変動していないのに負に発現制御されたなどの誤解が生じてしまう。

サンプル間の正規化手法として、

TMM (trimmed mean of M-values)
DESeq

が挙げられる。これらの手法は発現変動の大きい遺伝子が発現量の多くを占めていることを前提としている。

TMM (trimmed mean of M-values)

まず、TMMではリファレンスのサンプルを選出する。各遺伝子についてのM-valueはリファレンスに対するサンプルの発現量の割合を対数変換した値となる。極端な値のM-valueの遺伝子を除いた後、各サンプルについてこれらのM-valueの重みづけ平均が設定される。
ライブラリkを構成している遺伝子g=1,2,...Gについて、リファレンスのライブラリをrとするとき、M-valueは以下のように定義される。

$M_{g,k}^{r} = \log_{2} \frac{r_{g,k}/N_k}{r_{g,r}/N_k}$

DESeq

TMMと同様に、DESeqは全サンプルにおける幾何平均に対する特定のサンプルにおける各遺伝子のリードカウントの割合の中央値として size factor を計算する。
m 個のライブラリーが存在するとき、ライブラリー j の size factor $s_{j}$ は以下のように定義される。ここで、kij はライブラリー j における遺伝子 i のカウントデータを表す。

$s_{j} = median \frac{k_{i,j}}{(\prod_{v=1}^{m}k_{i,v})^\frac{1}{m}}$

このとき、カウントデータ $q_{iρ}$ は以下で定義される。

$q_{iρ}=\frac{k_{i,j}}{s_j}$

いずれの手法も、欠損値の多いデータの場合うまく機能しない。発現値のプーリングをベースにした手法は、有意に発現している遺伝子に決まって存在する欠損値の影響を受けないように開発された。遺伝子発現変動の大きい遺伝子の選定は正規化の大きな影響を大きく受ける。この遺伝子の選定により最終的に決定される細胞のヘテロ性が決まってしまうため、サンプル内及びサンプル間の正規化手法の組み合わせに関する潜在的な影響はまだ大きくわかっておらず、今後も厳しい検証が行われていくだろう。

交絡因子の除去

正規化が終わると、次のステップは交絡因子を推定することになる。観測されるリードカウントは、生物学的な変動のみならず technical noise と呼ばれる生物学的には何の関係もないノイズを含んでいることがある。
この影響は、Ambion社のERCC Spike-In Mix 等のspike-inを使用して補正される。複数の環境において有意に発現している遺伝子を比較するbulk RNA-seqに比べ、scRNA-seqでは1つの環境における細胞群を比較するため、生物学的な発現変動と関係のないbatch effectが生じやすい。 Technical repliateとして、同じ実験を繰り返し行うことで、この影響を観測することが可能だが、労力とコストがかかってしまう。

この影響を取り除くため、scLVMという手法が開発された。これを用いると、潜在変数の影響を除くことができる。これを使用することで、cell cycle covariatesが存在する場合不可能であったTh2細胞の分化に必須の遺伝子を同定することが可能となった。

既知、若しくは未知の変数に対しては、ランダムノイズを盛り込んだ線形結合によっても対処することができる。

セルタイプの同定

次元削減

疾患進行中、分化中にわたって細胞は数多の状態を遷移する一方で、数少ないマーカー遺伝子が免疫細胞における分化抗原群（CD）等の良く確立されたマーカー遺伝子でさえ、隠された多様性が存在する。
各遺伝子による「次元の呪い」に対処するためには、正規化の後、次元削減を行う必要がある。広く利用されている線形変換による次元削減手法としてPrincipal component analysis (PCA) が存在する。次元削減を行い細胞を2次元空間に射影することで、解釈が容易になる。

非線形変換による手法としては、

tSNE (t-distributed stochastic neighbor embedding)
LLE (locally linear embedding)
Isomap (Isometric Mapping)

等が挙げられ、t-SNEはCell Ranger pipeline (10× Genomics) やSeurat (http://satijalab.org/seurat/)にも実装されている。 LLEやIsomapはマイクロアレイデータの解析では良い結果を示したため、scRNA-seqデータでももっと検証されるべきである。

一方で、これらの次元削減を行うと本当に見たい生物学的な情報を失うことになる可能性があることから、次元削減は慎重に行う必要がある。

クラスタリング

クラスタリングはミトコンドリアRNAが有意に発現しているクラスターを同定することで細胞のQCにも使用することができる。この手法は、ミトコンドリアRNAの割合が高く、細胞質RNAの割合が失われている細胞は細胞膜が破裂していることを示す研究に基づいている。クラスタリングが終わると、次のステップは異なるクラスタで有意に発現しているマーカー遺伝子を同定することになる。

マーカー遺伝子探索

最もシンプルなモデルとしてはポアソン分布が使用される。scRNA-seqデータの様々な因子によるノイズを考慮したモデルとしては、負の二項分布が挙げられる。負の二項分布では、過分散を考慮したモデルを構築することができる。代わりに、dropout等のtechnical noiseを考慮したエラーモデルを当てはめることもできる。 single cell発現解析に使用されるプラットフォームでは、2つの確率論的プロセスを組み合わせる。正しく増幅され、実際の量と相関のある発現量についてであり、もう1つは正常に増幅されなかった、若しくは検出されなかった発現量についてである。

モデリング手法としては、混合モデルは単峰性のモデルに比べ有利ではあるが、ヘテロな細胞分布は二峰性の分布を示すことが多い。

制御ネットワーク推論

制御ネットワーク (GRNs) の解明は、生細胞内の複雑な細胞プロセスを理解するのに役立つ。
また、これらの解明は遺伝子、タンパク間の相互作用の理解にも役立つことが多い。
scRNA-seqが出現するまでは、マイクロアレイデータを使用したRNA-seqがこの用途で使用されていたが、scRNA-seqの出現により、GRNsの推論がより容易になった。細胞のヘテロ性の高さと遺伝子間相互作用の多様さから、GRNの決定にはまだ課題が残っている。

GRN推定のアルゴリズムは、機械学習ベースのもの、共発現ベースのもの、統計モデルベースのもの、情報理論ベースのものなど多岐に渡る。
中でも共発現ベースのものは、推定上の関係を同定するのにおそらく最もシンプルな手法である。
しかし、これらの手法は細胞システムの正確なダイナミクスを表現することはできない。
Bayesian network等のモデルベースの手法では、多くのパラメータを使用し、時間のかかる作業となる。
加えて、グラフィックモデルは、多くの遺伝子について可能性として存在しうるパスウェイを全て検索する必要があるため、NP困難な問題である。情報理論ベースの手法は仮定を置かず、非線形な遺伝子間の相互作用を測定することができることから、人気が出ている。

上述したように、technical noiseは生物学的な変動情報との分離が難しい。また、残りの変動因子に関してはほとんどわかっていない。しかしながら、複数のセルタイプが非同期的に存在するsingle cellデータの性質からは、制御ネットワークの推定に必要な固有の統計的ばらつきがを観測することができる。これまでいくつかの著名な手法がsingle-cellデータからGRNを同定するために開発され、T細胞に適用されることで新たな知見が明らかになっている。

転写産物の変化は永続するわけではないので、制御ネットワークの検出は適切なタイムスケールで可能であるべきである。加えて、遺伝子間ネットワークの方向性についても検証されなければならない。

細胞階層の再構築

各細胞は、動的なプロセスを受け続けており、様々な環境刺激に対する応答を示している。この応答の速度は、一瞬で完了するものから何年もかかるものまで多岐に渡る。これをbulkのRNA-seqで行おうとすると、特別な手法で細胞を同期させる必要がある。ところが、single-cellの場合、細胞系譜に沿った異なる瞬間の細胞を観測することができる。さらにアルゴリズムを適用することで、これらの細胞を使用し、分化や細胞サイクルといった動的な細胞系譜を再構築することができる。

“pseudotime”という概念はMonocleにより取り入れられた。この“pseudotime”は、実際の時間とは異なり、

Monocleはまず、nodeが細胞を示し、edgeが各細胞のペアに対応するようなグラフを描画する。 edgeの重みは、independent component analysis (ICA)による次元削減結果から得られた細胞間距離に基づいて計算される。その後、 minimum spanning tree (MST) が最も長い系譜を検索するのに使用される。より新しいバージョンのMonocle2では、計算速度がより早く、正確な結果を出力する。また、データドリブンな教師なし学習を実装している。

時間に関する情報が使用できる場合は、教師あり学習がベースの手法を使用する方がより正確である。例えば、Single-cell clustering using bifurcation analysis (SCUBA)は分岐解析を実装しており、早期のマウス胚の系譜を複数のタイムポイントで取得した発現マトリクスから再構築する際等に使用される。

scRNA-seqは、ニューロン新生中の細胞系譜の再構築にも利用される。この技術の適用例として、Div-Seqは単離された核を直接シーケンスすることで、組織を分離する手間が必要ない。酵素を用いた組織の分離を行うと、RNAの構造が壊れてしまうため、複雑な組織から取得した細胞を研究する際にはこの技術は必要不可欠である。線形系譜の構築を推定することから始めるが、最近の研究成果では枝分かれの概念を取り入れている。これは細胞の動的なシステムを理解する際に欠かせない。

↓擬似系譜解析で使用されるツールとその概要 f:id:kimoppy126:20190213002049p:plain

将来の展望

scRNA-seq により、生物学の発展に革命的な影響を及ぼした。この技術を用いると、

がん細胞が薬剤耐性を得ていく過程を計測する
特定のパスウェイについて発現プロファイルを分析する
転写動態をモデリングすることで、クローンと系統発生の関係を再構築する

といったことが可能になる。

近年、血中CTCの解析はDNAを医療診断マーカーとして活用することに着目から、“liquid biopsy”の最盛期を迎えた。

系統追跡実験（linage tracing）は生物学において、長年の疑問であった。1細胞胚がいかにして、様々なセルタイプを持った組織に分化するのかカリフォルニア工科大学の研究者が最近、多くの世代にわたって、連続的なシーケンス結果から系譜を再構築する手法を開発した（https://www.ncbi.nlm.nih.gov/pubmed/27869821?dopt=Abstract）。もう1つ興味深いこととして、肝細胞の制御ネットワークに関与する遺伝子を同定することだ。幹細胞が機能的な細胞に分化するトリガーがどのように引かれているのか、といったことは理解され始めたばかりであり、ヒトの健康、疾患に潜在する生物学的なプロセスを理解するのに欠かせない情報である。

シーケンシングのコストは年々下がっており、次の5年間で100万細胞以上を日常的に解析することが可能になるといわれている。

Human Cell Atlasでは、ヒトを構成する35兆個の細胞の機能推定を行うことを目的としているが、 3千万から1兆の細胞をシーケンシングした結果の発現プロファイルを使用し、クラス分け及び新規のセルタイプの同定を行う。高度に多様化した免疫細胞のscRNA-seqは、免疫細胞に固有のヘテロ性をより深く理解するのに役立つ。脳細胞についての研究は、神経関連疾患に関与している新規のパスウェイを同定するのに役立ち、バイオマーカー探索における新規の治療ターゲットの提供に役立つ。生体医学の研究におけるscRNA-seqの応用によって様々な組織、器官における構造と機能の生理的な関係がわかると見越している。標準的なバイオインフォマティクスのパイプラインにより、新規の生物学的システムの解明及び治療開発機会の創出に繋がるだろう。