Trajectory inferenceとRNA velocityの違い

We first split multi-allelic sites to represent separate bi-allelic sites. All calls that did not pass the following hard filters were then set to no-call in our analysis: - For homozygous reference calls: - Genotype Quality < 20; Genotype Depth < 10; Genotype Depth > 200 - For heterozygous calls: - (A1 Depth + A2 Depth)/Total Depth < 0.9; A2 Depth/Total Depth < 0.2; Genotype likelihood[ref/ref] < 20; Genotype Depth < 10; Genotype Depth > 200 - For homozygous alternative calls: - (A1 Depth + A2 Depth)/Total Depth < 0.9; A2 Depth/Total Depth < 0.9; Genotype likelihood[ref/ref] < 20; Genotype Depth < 10; Genotype Depth > 200 # Trajectory inferenceとRNA velocityの違い ## まとめ - [ジエンブル社シングルセルRNA-seq解析](https://genble.co.jp/single-cell-rna-seq/#kaisekirei8) - [和光純薬時報　Vol.89 No.1(2021年1月号)にて鈴木絢子先生執筆](https://labchem-wako.fujifilm.com/jp/journal/docs/jiho891.pdf) から - RNA velocity - 細胞の分化や状態変化などの**変化する方向**が推定できる - Pseudotime trajectory - 分化経路に沿った**遺伝子発現の変化**を読み取れる ## Trajectory analysis #### [NIH Bioinformatics Training and Education Program](https://btep.ccr.cancer.gov/question/single_cell_rna_seq/can-you-say-something-about-trajectory-analysis-vs-velocity-analysis/) より > Can you say something about trajectory analysis vs. velocity analysis? - Trajectory analysis (trajectory inference, pseudotime analysis) - variable genesからクラスター間の「経路」を見つけようとする。 - (服部追記: 時間軸となる曲線のことを指す。例として[Ernst et al., 2019](https://www.nature.com/articles/s41467-019-09182-1)では, [Principal Curve](https://salad-bowl-of-knowledge.github.io/hp/statistics/2019/10/06/princurve2.html)を用いていたりするのがイメージとしてわかりやすい。) - 分化の方向は推定しないので、意味のないtrajectoryを作成することもできる。 - Velocity analysis - 分化の方向と速度が描写される。 - そのため、生物学的により正確なpseudotimeの作成ができる。 #### [A guide to trajectory inference and RNA velocity](https://www.biorxiv.org/content/10.1101/2021.12.22.473434v1.full) より - Trajectory inference - trajectory inferenceは単一次元への次元削減とみなすことができる。つまり、生物学的な変化をできるだけ多くとらえる単一次元を見つけたい。 - この時、PCAよりもより専門的なアプローチの方が良い結果をもたらす。 - PC1はシーケンス深度や、各細胞のzeroの割合と相関していることが多いらしい - 方向性を推論できない - 追加データやドメイン知識なしで、スタートを決めることができない。 ## RNA velocity ### the steady-state model (→ velocyto) [ La Manno et al. (Nature, 2018)](https://doi.org/10.1038/s41586-018-0414-6) → the original model published in [Nature](https://www.nature.com/articles/s41586-018-0414-6) [良いgifを見つけた](https://singlecellanalysistutorial.readthedocs.io/en/latest/notebooks/scVelo.html) <img width="600" alt="image.png (905.1 kB)" src="https://user-images.githubusercontent.com/31883718/80227452-eb822480-864d-11ea-9399-56886c5e2785.gif"> - 左の図 - 時間軸(転写促進(誘導?) → 転写抑制)に対してRNA量の変化をシミュレーションしたもの - 「転写(+) → 抑制/停止(-)」というプロセスを表していて、unspliced mRNA (紫色)はspliced mRNA (緑)より先に立ち上がる - 右の図 - Induction (上の弧を通り、右上の定常状態になる) - 転写が始まる前、unspliced/spliced mRNA共に存在しないので、左下にプロットされる - 転写が始まると、初めunsplicedの方が多いため上側に向かっていくが、徐々にsplicingによってspliced mRNAも増加して右側方向に移動する - 定常状態(左図で3 h経った頃あたり)に達したのが、右上の細胞集団 - Repression (下の弧通って、左下の定常状態) 上図の`steady-state ratio`と、各細胞の$\frac{unspliced\ mRNA}{spliced\ mRNA}$比からのdeviationを速度として推定する。しかし、このモデルは2つの大きな仮定を置いている。(1)全ての転写産物は同じsplicing rateであること、(2)発現レベルが定常状態にある集団が存在すること。 ### the dynamical model (→ scVelo) [Bergen et al. (Nature Biotechnology, 2020)](https://doi.org/10.1038/s41587-020-0591-3) [Bergen2020.pdf (10.1 MB)](https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/21/62720/40260081-173f-41af-a1a3-36971b26b735.pdf) - the steady-state modelの以下の2つの仮定を克服するためのモデル - a common splicing rate - the observation of the full splicing dynamics with steady-state mRNA levels <details><summary>詳細</summary> - **各遺伝子について**のモデル。 - 各遺伝子のvelocityから、それぞれの細胞の将来の状態を推定する <img width="2128" alt="image.png (925.3 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/26/62720/9a53b306-1e92-4fbc-9bd2-c4e5c9c027b6.png"> the likelihood-based dynamical model > Hence, we infer the parameters by EM, iteratively estimating the reaction rates and latent variables via maximum likelihood. <img width="2128" alt="image.png (1.6 MB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/26/62720/83348aaa-9680-484e-b620-ade82f62f9cd.png"> - Extended kinetic models - for complex gene dynamics - Stochastic models - for cell-specific dynamics - Multivariate models - towards system dynamics - Multi-omics models </details> ### scVelo (RNA velocity) : どのモデルを選べば良いのか・どう違うのか [scVeloのdocsより](https://scvelo.readthedocs.io/about/) > RNA velocity estimation can currently be tackled with three existing approaches: 1. steady-state / deterministic model (using steady-state residuals) 1. stochastic model (using second-order moments), 1. dynamical model (using a likelihood-based framework). → [`scv.tl.velocity(adata, mode='stochastic', **params)`](https://scvelo.readthedocs.io/scvelo.tl.velocity/#scvelo.tl.velocity)でvelocityの推定ができる。 - `mode`で利用できるmodelが選べる (default : `'stochastic'`)。 - `mode='deterministic'` - 決定論的 - `mode='stochastic'` - 確率的 - [`mode='dynamical'`](https://scvelo.readthedocs.io/DynamicalModeling/#Dynamical-Model) - 完全な動的モデルの解はこのmodelで得られる (the full transcriptional dynamics of splicing kinetics) 良い質問があった ['dynamical' vs 'stochastic' mode #125](https://github.com/theislab/scvelo/issues/125) #### → **"dynamical" modelの方が良い**。 - 1分以内で計算が終了する"stochastic" modelに比べて、"dynamical" modelは、遺伝子数に依って計算時間がより多くなる(通常20~30分程度)。 - したがって、velocityにのみ興味があるなら"stochastic" modelで十分だからデフォルトになっている。 - 多くのpreprintで用いられ、original paperでも示されていたように、"dynamical" modelの方が優れている。さらにvelocity以外にも多くの事が捉えられる(例えばlatent time、top-likelihood genes (putative drivers)、kinetic ratesなど) ### 文献 (オススメの) - Original - [ La Manno et al. (Nature, 2018)](https://doi.org/10.1038/s41586-018-0414-6) (→ **velocyto**) - RNA velocityの始まり。 - the steady-state model。 - [Bergen et al. (Nature Biotechnology, 2020)](https://doi.org/10.1038/s41587-020-0591-3) (→ **scVelo**) - これまで(Nature 2018)のモデルが改良されたもの。 - the dynamical model。 - scVeloはどちらのモデルも使用可能。 - Review系 - [A guide to trajectory inference and RNA velocity](https://www.biorxiv.org/content/10.1101/2021.12.22.473434v1.full) - Trajectory inference methods とRNA velocity両アプローチについて。 - [Bergen et al., 2021, *MSB*](https://doi.org/10.15252/msb.202110282) - RNA velocity—current challenges and future perspectives - Review - Blog系 - [RNA Velocity: The Cell’s Internal Compass](https://towardsdatascience.com/rna-velocity-the-cells-internal-compass-cf8d75bb2f89) - とりあえずこれを読めばRNA velocityがわかった気になれる。 - 長すぎず短すぎない解説 + scVeloのチュートリアル付き - [RNA velocity analysis with scVelo](https://smorabit.github.io/tutorials/8_velocyto/) - Seuratで解析済みのデータから、RNA velocityの下流解析まで。 - YouTube - [MIA: Volker Bergen, RNA velocity generalized to transient cell states through dynamical modeling](https://www.youtube.com/watch?v=QWfd_kSNPgk) - scVeloの作成者 (Theis groupから)による解説 - 8:29 ~ : Concept of RNA velocity ## Notes <details> Pseudotime解析は、細胞集団を擬時間という1次元の"軸"に変換することで、疑似的に細胞の分化を追跡できる。 > uncover global structure and convert this structure into smooth lineages represented by one-dimensional variables, called “pseudotime.” > > [Slingshot: Trajectory Inference for Single-Cell Data](https://bioconductor.org/packages/devel/bioc/vignettes/slingshot/inst/doc/vignette.html) より > the trajectory inference problem can be seen as a dimensionality reduction to a single dimension: we want to find a single dimension that captures as much biologically relevant variation as possible. > > [A guide to trajectory inference and RNA velocity](https://www.biorxiv.org/content/10.1101/2021.12.22.473434v1.full) より RNA velocityは、スプライスされていない未成熟な mRNA とスプライスされた成熟 mRNA の比（RNA 速度）に着目することにより、個々の細胞の運命(状態変化)を予測する。 > 他にも、scRNA-seq データよりさまざまな情報を抽出するための解析手法が開発されている。Trajectory解析は、遺伝子発現パターンの類似性から細胞を擬似時間軸で並べるものであり、Monocle 8）等の解析ツールを用いて実施する。また、RNA velocity 9）という、スプライスされていない未成熟な mRNA とスプライスされた成熟 mRNA の比（RNA 速度）に着目することにより、個々の細胞の運命を予測する解析も実施されている。これらの解析を行うことで、本来はある瞬間の計測データであるはずのものから疑似的に細胞の分化や状態変化を追跡することができる。 > > [和光純薬時報　Vol.89 No.1(2021年1月号)にて鈴木絢子先生執筆](https://labchem-wako.fujifilm.com/jp/journal/docs/jiho891.pdf) より [NIH Bioinformatics Training and Education Program](https://btep.ccr.cancer.gov/question/single_cell_rna_seq/can-you-say-something-about-trajectory-analysis-vs-velocity-analysis/) より > Can you say something about trajectory analysis vs. velocity analysis? - Trajectory analysis (trajectory inference, pseudotime analysis) - variable genesからクラスター間の「経路」を見つけようとする。 - (服部追記: 時間軸となる曲線のことを指す。例として[Ernst et al., 2019](https://www.nature.com/articles/s41467-019-09182-1)では, [Principal Curve](https://salad-bowl-of-knowledge.github.io/hp/statistics/2019/10/06/princurve2.html)を用いていたりするのがイメージとしてわかりやすい。) - 分化の方向は推定しないので、意味のないtrajectoryを作成することもできる。 - Velocity analysis - 分化の方向と速度が描写される。 - そのため、生物学的により正確なpseudotimeの作成ができる。 #### おまけ Velocityは、イントロンの保持率から推定されるが、10xなどのfull-lengthではないデータでも上手くいく。ただしsingle-nucleiデータでは、intron-retaining transcript readsを多く含むため、あまり上手くいかないらしい。 #### その他 [RNA velocity—current challenges and future perspectives](https://www.embopress.org/doi/full/10.15252/msb.202110282) [How to do single-cell RNA velocity analysis | Complete introduction](https://www.youtube.com/watch?v=AUiYxtGJYtg) [RNA Velocity: The Cell’s Internal Compass](https://towardsdatascience.com/rna-velocity-the-cells-internal-compass-cf8d75bb2f89) → 長すぎず短すぎない解説 + scVeloのチュートリアル付き ## 理論的側面 - Title - [A guide to trajectory inference and RNA velocity](https://www.biorxiv.org/content/10.1101/2021.12.22.473434v1.full) - 先行研究 - シングルセルのデータが爆発的に増えている。近年では、scRNA-seqから細胞分化や細胞周期、cell (de)activationなどの動的なプロセスに注目が集まっている。 - この論文の目的 - Trajectory inference methods とRNA velocity両アプローチについて、 - 概念的・理論的側面 - 両アプローチの組み合わせ方 - 実データを用いた使用例 - Trajectory inference - trajectory inferenceは単一次元への次元削減とみなすことができる。つまり、生物学的な変化をできるだけ多くとらえる単一次元を見つけたい。 - この時、PCAよりもより専門的なアプローチの方が良い結果をもたらす。 - PC1はシーケンス深度や、各細胞のzeroの割合と相関していることが多いらしい - 方向性を推論できない - 追加データやドメイン知識なしで、スタートを決めることができない。 - RNA velocity - splicing dynamicsは以下の3つに分類でき、これらの情報を利用する。それぞれの速度をα, β, γと置く。 - 転写 - スプライシング - 分解 - <img width="400" alt="image.png (160.5 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/14/62720/246f8722-9ebb-40cf-b6f5-2bad65e4540c.png"> - precursor mRNA ( *u* )とmRNA ( *s* )の時間変化を微分方程式でモデル化 - $\frac{du}{dt} = α - βu$ - $\frac{ds}{dt} = βu - γs$ - 定常状態(steady-state)における解析解 - <img width="400" alt="image.png (41.2 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/14/62720/e08e7792-2b2d-4bdd-aef0-92e307daf4e2.png"> - 上記のODEを用いてシミュレーションした結果 (mRNA存在量の時間変化) <img width="400" alt="image.png (905.1 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/14/62720/7109cef7-b6b8-47c2-8b95-0bd5901f7b1b.png"> [良いgifを見つけた](https://singlecellanalysistutorial.readthedocs.io/en/latest/notebooks/scVelo.html) <img width="600" alt="image.png (905.1 kB)" src="https://user-images.githubusercontent.com/31883718/80227452-eb822480-864d-11ea-9399-56886c5e2785.gif"> - unspliced mRNAはspliced mRNAより先に立ち上がる - 右図(gifの上の図のこと)は、定常状態に達する前に転写が終了した場合を示す - 左の図と違って、長い時間その状態が続かないため、定常状態には達しない。 → 「early switching」 - 下の図(phase portraits)は、上の弧を通り、右上の定常状態になる。転写抑制が始まると下の弧通る。左下と右上の集団が平衡状態にある細胞。 - モデルのパラメータを推定する方法は2つある - モデルのパラメータは以下のこと - 各細胞に割り当てられるtime point - 転写活性化(induction)から抑制に切り替わるタイミング(the time when the system switches from induction to repression) - Steady-state model → 2018年Natureに掲載されたoriginalのモデル - splicing dynamicsは定常状態に達していると仮定。位相ポートレート(phase portraits)において、左下と右上の集団が平衡状態にある細胞。 - **定常状態であることと、splicing rate βが全ての転写産物に共通であること**が必要。 - $\frac{ds}{dt} = βu - γs = 0$ より、 $u = \frac{γ}{β}s$ - $\frac{ds}{dt} = 0$であるということは、splicingによって産生されるmRNAと分解されるmRNAが釣り合っている状態 (→ steady-state) <img width="600" alt="image.png (596.0 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/14/62720/6f45a9b8-9486-42e4-ab01-dd1f99cb17f6.png"> - $\hat{γ} = \frac{γ}{β}$ とすると、一般化線形回帰モデルとして、この線より上に存在する細胞($u > \hat{γ}s$)は正の速度、つまり転写活性(induction)が起きていると考えることができる。 - なぜならunspliced mRNAの方が多いから (u - s > 0)、splicingを受けるmRNAより新たに転写されるmRNAの方が多い - - Dynamical model - transcriptome-wide assumptionsを作らないで、各遺伝子の全てのtranscriptional kineticsを直接解く(?, solve) - 推定するパラメータは、各遺伝子についてのα (transcription rate)・β (splicing rate)・γ (degradation rate)。しかし、実験的にtime point *t*とstate *k* (induction, repression, or steady state)は推定できない。 - 最尤法では推定できなかった。代わりにEMアルゴリズムを用いた。 - E-step : (i.e., expectation step) - M-step : (i.e., maximization step) - 定常状態(steady-state)における解析解 - <img width="400" alt="image.png (41.2 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/14/62720/e08e7792-2b2d-4bdd-aef0-92e307daf4e2.png"> - これをスピードアップのために、$s(t)$を以下のように近似 - <img width="859" alt="image.png (10.9 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/26/62720/736f6898-0ac4-46ac-965b-023e4b5e95e0.png"> ### EMアルゴリズム [【徹底解説】EMアルゴリズムをはじめからていねいに](https://academ-aid.com/ml/em) - EMアルゴリズムは**確率モデルの潜在変数・パラメータに関する最尤推定を行うため**の手法 ## scVelo - [Trajectory Analysis using 10x Genomics Single Cell Gene Expression Data](https://www.10xgenomics.com/resources/analysis-guides/trajectory-analysis-using-10x-Genomics-single-cell-gene-expression-data) <img width="1408" alt="image.png (59.6 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/25/62720/189de482-cd65-45b0-b5f3-4fd7cd8b2bd1.png"> RNA velocityでintronic readsが必要。`cellranger count`の際に`--include-introns`を付けないとintronic readsは捨てられる。 [Release notes for 7.0](https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/release-notes)で、`--include-introns=true`はデフォルトになった。 ### Usage - [Getting Started](https://scvelo.readthedocs.io/getting_started/) - #### Input data - two count matrices of - pre-mature (unspliced) - mature (spliced) abundances - **velocytoで作成した`.loom`ファイルを読み込む** (もしくはloompy/kallistoで作れる) - `.loom`ファイルは`cellranger count --include-introns`(version 6以前の場合) - #### その他 - 可視化にはmatplotlibが使われている - matrixは`adata`として保持する - spliced/unspliced counts、velocitiesは`adata.layers`に保持される <img width=500 alt="anndata.svg (40.1 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2022/07/27/62720/a5ff2295-f65a-4ae6-9c23-df46c5cec321.svg"> - The probabilities of one cell transitioning into another cell are computed using cosine correlation (between the potential cell transition and the velocity vector) ### 複数のモデルの利用法の違い(どのモデルを選べば良いのか・どう違うのか) [`scv.tl.velocity(adata, mode='stochastic', **params)`](https://scvelo.readthedocs.io/scvelo.tl.velocity/#scvelo.tl.velocity)でvelocityの推定ができる。 - `mode`で利用できるmodelが選べる (default : `'stochastic'`)。 - `mode='deterministic'` - 決定論的 - `mode='stochastic'` - 確率的 - [`mode='dynamical'`](https://scvelo.readthedocs.io/DynamicalModeling/#Dynamical-Model) - 完全な動的モデルの解はこのmodelで得られる (the full transcriptional dynamics of splicing kinetics) 良い質問があった ['dynamical' vs 'stochastic' mode #125](https://github.com/theislab/scvelo/issues/125) 質問 > I am wondering why the default mode in scvelo.tl.velocity being set to mode='stochastic' when the majority of the biorxiv paper seems to focus on the "dynamical" model. Is the "dynamical" model still under development? > > Also, from the paper, it seems to me that the "dynamical" model is superior to the "stochastic" model. Thus, should I always try to run scVelo in "dynamical" model instead? 回答 > As you mentioned the dynamical model is superior and also yields many additional insights beyond velocities such as latent time, top-likelihood genes (putative drivers), kinetic rates etc. > ... > If you're only interested in velocities, the stochastic model does a fairly good job yielding results approximating those from the dynamical model and usually runs within a minute. - 1分以内で計算が終了する"stochastic" modelに比べて、"dynamical" modelは、遺伝子数に依って計算時間がより多くなる(通常20~30分程度)。 - したがって、velocityにのみ興味があるなら"stochastic" modelで十分だからデフォルトになっている。 - 多くのpreprintで用いられ、original paperでも示されていたように、"dynamical" modelの方が優れている。さらにvelocity以外にも多くの事が捉えられる(例えばlatent time、top-likelihood genes (putative drivers)、kinetic ratesなど) ### 疑問 - [ ] each velocity vectorの方向はどうやって計算？スカラーはsteady-state ratioとの残差 </details> --- ## 模索 #### searched by 「"CDR1" sequence clustering」 - :star: [T cell receptor sequence clustering and antigen specificity](https://www.sciencedirect.com/science/article/pii/S2001037020303305) - We review the main TCR sequence clustering methods and the different similarity measures they use, and discuss their performance and possible improvement. - Currently there are a number of approaches that aim to cluster TCRs by extrapolating information from their primary sequences to study their specificites. - <img width="667" alt="image.png (244.8 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/01/27/62720/6e8940ee-2c51-42d5-b200-86c0a4261383.png"> - TCRの一次配列がわかっている場合、配列比較により距離行列TCR distance matrixを作成することができる。このマトリックスを用いて、配列の類似性に基づいて個々のTCRをクラスタリングし、生物学的類似性、すなわちエピトープ反応によってクラスタリングすることができる - [x] なんでCDR3のみが多い？ - :star: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8635517/ - https://github.com/immunoengineer/gliph - CDR3のアミノ酸残基はペプチドと密接に相互作用するから、特異性に重要。 - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5794212/ - TCRにおいて、CDR3がペプチドと接触する確率が高いことを示している(上, a) - CDRの接触の仕方で5つのクラスターに分類(下, b) - TCRs are clustered into five general contact modes according to contact profiles of all six CDRs. - <img width="800" alt="image.png (278.3 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/02/02/62720/f6fc1db8-2f56-49b7-8baa-6f8e26ebae33.png"> - [ ] どうやってCDR1を決定している？ #### searched by 「CDR3 clustering」 - [Clustering based approach for population level identification of condition-associated T-cell receptor β-chain CDR3 sequences](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04087-7) - ### [CellaRepertorium](http://www.bioconductor.org/packages/devel/bioc/vignettes/CellaRepertorium/inst/doc/cdr3_clustering.html) - [Clustering and differential usage of repertoire CDR3 sequences](http://www.bioconductor.org/packages/devel/bioc/vignettes/CellaRepertorium/inst/doc/cdr3_clustering.html) - <img width="1492" alt="image.png (53.8 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/01/27/62720/240589ba-2469-4f7c-89f6-ad08c5228f3e.png"> - <img width="1492" alt="image.png (108.7 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/01/27/62720/b41704b9-69e1-4989-8a25-a6a195820ee3.png"> - The min_length can be set somewhat smaller, but there is a lower limit for the cdhit algorithm. (ここでは"min_length = 5"となっていた) - ### [clusTCR](https://svalkiers.github.io/clusTCR/) - [ClusTCR: a python interface for rapid clustering of large sets of CDR3 sequences with unknown antigen specificity](https://academic.oup.com/bioinformatics/article/37/24/4865/6300511) - A two-step clustering approach that combines the speed of the Faiss Clustering Library with the accuracy of Markov Clustering Algorithm - <img width="2081" alt="image.png (415.6 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/01/27/62720/967245b5-d331-47bd-9247-ccb8e04c269a.png"> - #### [immunarch — Fast and Seamless Exploration of Single-cell and Bulk T-cell/Antibody Immune Repertoires in R](https://immunarch.com/articles/web_only/clustering.html) - B細胞受容体（BCR）およびT細胞受容体（TCR）の配列のクラスタリングは、B細胞およびT細胞の解析において一般的なステップである。 - B細胞レパートリーでは、このステップはクローン系譜を作るために必要である。B細胞のクローン性系統は、共通の起源（同じVDJ再配列イベントから生じる）と共通の祖先を持つと推定されるB細胞の集合を表しています。クローン系譜の情報により、体細胞超変異レベルを推定し、結合親和性変更変異を見つけることができる。 - 同一のエピトープを認識できるT細胞は、類似性の高いTCRを持つことが示されました（Glanville et al.）クラスタリングにより、同一の抗原を認識する可能性の高いレパートリーから、類似性の高いTCRを見つけることができます。 - クラスタリング解析には2つのステップがあります。最初のステップでは、配列間の距離を計算する必要があります。第二段階は、配列間の距離の情報を用いてクラスタリングを行う必要があります。 - クラスタリングの応用により、クロノタイプの概念を拡張することができる。クローノタイプとは、共通の祖先や性質を持つとされる細胞の集合のことである。クロノタイプの定義にはいくつかの方法が考えられる。例えば、ある研究者は、同じV遺伝子と同じヌクレオチドCDR3配列を持つ配列の集合をクロノタイプと定義している。また、同じV遺伝子、同じJ遺伝子、同じアミノ酸CDR3配列を持つ細胞集合をクロノタイプと定義する研究者もいる。クラスタリングは、あなたのビジョンに基づいた新しい配列のグループを作成する方法を提供します。 - 距離の計算 - CDR3領域は、BCRやTCRの中で最も可変性の高い領域です。この領域は、一般的に配列間の距離を計算するために使用されます。しかし、他の領域もこの距離の計算に使用することができます - フルヌクレオチドリードのように。カラムの選択は研究目的のみに依存するため、距離を計算するために任意のユーザー定義の配列のカラムを使用することができます。例えば、'nt'配列を使用すると、より保守的なクラスタリングが行われます。 - 距離の計算が終わったら、今度は配列のクラスタリングを行い、CDR3配列やその他の類似した領域を持つTCRやBCRを探します。CDR3配列の類似性のみで解析を行うことも可能ですが、この配列の類似性に限定して解析を行わないことをお勧めします。 - 数学的に言えば、クロノタイプはグラフの頂点である。V遺伝子とJ遺伝子が同じで、CDR3.aaがしきい値以下しか違わない場合、その頂点の間に辺を描く。クラスタとは、エッジで結ばれたクロノタイプの集合であり、連結成分である。 - クラスタリングは、クロノタイプの概念を拡張するためによく使われる。例えば、同じアミノ酸のCDR3配列を持つ配列は同じ性質を持つ可能性が高いので、クロノタイプを配列のグループと考えることができる。 - searched by 「immunarch "CDR1"」 > [For B-cell connoisseurs and T-cell fans we included support for CDR1 and CDR2 to #immunarch data format. Enjoy!](https://twitter.com/immunomind/status/1458125147106332673) #### searched by 「"CD-HIT" "CDR1"」 - https://mhlw-grants.niph.go.jp/system/files/2013/133151/201324073A/201324073A0006.pdf - 樹形図解析 - CDR1、CDR2領域についてIMGT/V-QUESTにより解析した配列情報についてClustal W2のKimura Methodを用いて樹形図を作成 - 健常者では発散して上手くクラスターを形成していない部分もある - 全体として疾患群でクラスターを形成する傾向が見られた <img width="200" alt="image.png (764.0 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/02/02/62720/cf746f6c-cc65-4867-af70-0ed3d344fe45.png"> <img width="200" alt="image.png (627.2 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/02/02/62720/7e75e05b-759a-461f-902d-1db7b61bb563.png"> ## 関連ツール - IMGT/V-QUEST - https://amelieff.esa.io/posts/2918#IMGT/HighV-QUEST - https://bio.tools/imgt_highv-quest - <img width="1573" alt="image.png (125.9 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/02/02/62720/5d2127c0-0e13-40b6-966a-20e6d638ffc1.png"> https://academic.oup.com/bioinformatics/article/37/24/4865/6300511 ClusTCRより > clustering methods for TCRs or CDR3 sequences: GLIPH2 ([Huang et al., 2020](http://www.ncbi.nlm.nih.gov/pubmed/28636589)), iSMART ([Zhang et al., 2020](http://www.ncbi.nlm.nih.gov/pubmed/31831563)) and TCRDist ([Glanville et al., 2017](http://www.ncbi.nlm.nih.gov/pubmed/28636589); [Mayer-Blackwell et al., 2020](https://pubmed.ncbi.nlm.nih.gov/33398288/)). - :star: [**ClusTCR**](https://github.com/svalkiers/clusTCR) - [docs](https://svalkiers.github.io/clusTCR/) - シンプル - [GLIPH2](http://50.255.35.37:8080) - [inputとoutput](http://50.255.35.37:8080/demo) - The format of TCR input file is tab-delimited, expecting the following columns in this order: **CDR3b Vb Jb CDR3a Subject:condition Frequency** - [ ] 選ばなかった理由 - [iSMART](https://github.com/s175573/iSMART) - iSMART performs a specially parameterized pairwise local alignment on T cell receptor CDR3 sequences to group them into antigen-specific clusters. - (欠点) Python2で書かれている - :star: [TCRDist](https://github.com/kmayerb/tcrdist3) - [docs](https://tcrdist3.readthedocs.io/en/latest/) - CellaRepertorium - シングルセル - CD-HITを利用 ## テスト解析 ### ([CellaRepertorium](https://r.mcdavid.dev)) CDR1でやっている論文は無さそう - CellaRepertorium - CD-HITを利用 - ContigCellDBオブジェクト - Many tools from the Immcantation suite can work directly on ContigCellDB() objects. - IMGT/V-QUEST → Change-O(挟む必要があるのかはわからない) → CellaRepertorium https://r.mcdavid.dev ``` mkdir CellaRepertorium_test && cd $_ R ``` ```R BiocManager::install('CellaRepertorium') ``` https://r.mcdavid.dev/articles/cdr3_clustering.html ```R library(CellaRepertorium) library(dplyr) library(ggplot2) library(readr) library(tidyr) library(stringr) library(purrr) ``` ```R data(contigs_qc) MIN_CDR3_AA = 6 cdb = ContigCellDB_10XVDJ(contigs_qc, contig_pk = c('barcode', 'pop', 'sample', 'contig_id'), cell_pk = c('barcode', 'pop', 'sample')) > names(cdb@contig_tbl) [1] "anno_file" "pop" "sample" "barcode" [5] "is_cell" "contig_id" "high_confidence" "length" [9] "chain" "v_gene" "d_gene" "j_gene" [13] "c_gene" "full_length" "productive" "cdr3" [17] "cdr3_nt" "reads" "umis" "raw_clonotype_id" [21] "raw_consensus_id" "celltype" ``` <details><summary>`str(cdb)`</summary> ``` Formal class 'ContigCellDB' [package "CellaRepertorium"] with 8 slots ..@ contig_tbl : tibble [1,508 × 22] (S3: tbl_df/tbl/data.frame) .. ..$ anno_file : chr [1:1508] "/Users/amcdavid/Box Sync/research/scRNAseq/CellaRepertorium/inst/extdata/all_contig_annotations_b6_4.csv.xz" "/Users/amcdavid/Box Sync/research/scRNAseq/CellaRepertorium/inst/extdata/all_contig_annotations_b6_4.csv.xz" "/Users/amcdavid/Box Sync/research/scRNAseq/CellaRepertorium/inst/extdata/all_contig_annotations_b6_4.csv.xz" "/Users/amcdavid/Box Sync/research/scRNAseq/CellaRepertorium/inst/extdata/all_contig_annotations_b6_4.csv.xz" ... .. ..$ pop : chr [1:1508] "b6" "b6" "b6" "b6" ... .. ..$ sample : chr [1:1508] "4" "4" "4" "4" ... .. ..$ barcode : chr [1:1508] "AAAGTAGTCGCGCCAA-1" "AAAGTAGTCGCGCCAA-1" "AAAGTAGTCGCGCCAA-1" "AACCATGCATTTGCCC-1" ... .. ..$ is_cell : logi [1:1508] TRUE TRUE TRUE TRUE TRUE TRUE ... .. ..$ contig_id : chr [1:1508] "AAAGTAGTCGCGCCAA-1_contig_1" "AAAGTAGTCGCGCCAA-1_contig_2" "AAAGTAGTCGCGCCAA-1_contig_4" "AACCATGCATTTGCCC-1_contig_3" ... .. ..$ high_confidence : logi [1:1508] TRUE TRUE TRUE TRUE TRUE TRUE ... .. ..$ length : num [1:1508] 611 609 538 799 634 923 693 658 558 614 ... .. ..$ chain : chr [1:1508] "TRB" "TRB" "TRA" "TRA" ... .. ..$ v_gene : chr [1:1508] "TRBV26" "TRBV31" "TRAV13D-2" "TRAV9-1" ... .. ..$ d_gene : chr [1:1508] "TRBD1" "None" "None" "None" ... .. ..$ j_gene : chr [1:1508] "TRBJ2-7" "TRBJ1-5" "TRAJ37" "TRAJ15" ... .. ..$ c_gene : chr [1:1508] "TRBC1" "TRBC1" "TRAC" "TRAC" ... .. ..$ full_length : logi [1:1508] TRUE TRUE TRUE TRUE TRUE TRUE ... .. ..$ productive : chr [1:1508] "True" "True" "True" "True" ... .. ..$ cdr3 : chr [1:1508] "CASSPTDYEQYF" "CAWSPGDNQAPLF" "CAIEAGNTGKLIF" "CAVSAYQGGRALIF" ... .. ..$ cdr3_nt : chr [1:1508] "TGTGCCAGCAGTCCGACAGACTATGAACAGTACTTC" "TGTGCCTGGAGTCCCGGGGACAACCAGGCTCCGCTTTTT" "TGTGCTATAGAGGCAGGCAATACCGGAAAACTCATCTTT" "TGTGCTGTGAGCGCATACCAGGGAGGCAGAGCTCTGATATTT" ... .. ..$ reads : num [1:1508] 34803 37898 11520 3799 33548 ... .. ..$ umis : num [1:1508] 6 6 2 3 6 3 6 4 1 9 ... .. ..$ raw_clonotype_id: chr [1:1508] "clonotype23" "clonotype23" "clonotype23" "clonotype27" ... .. ..$ raw_consensus_id: chr [1:1508] "clonotype23_consensus_3" "clonotype23_consensus_1" "clonotype23_consensus_2" "clonotype27_consensus_1" ... .. ..$ celltype : chr [1:1508] "T_ab" "T_ab" "T_ab" "T_ab" ... ..@ contig_pk : chr [1:4] "barcode" "pop" "sample" "contig_id" ..@ cell_tbl : tibble [832 × 3] (S3: tbl_df/tbl/data.frame) .. ..$ barcode: chr [1:832] "AAAGTAGTCGCGCCAA-1" "AACCATGCATTTGCCC-1" "AACTGGTGTCTGATCA-1" "AAGCCGCAGTAAGTAC-1" ... .. ..$ pop : chr [1:832] "b6" "b6" "b6" "b6" ... .. ..$ sample : chr [1:832] "4" "4" "4" "4" ... ..@ cell_pk : chr [1:3] "barcode" "pop" "sample" ..@ cluster_pk : chr(0) ..@ cluster_tbl : tibble [0 × 0] (S3: tbl_df/tbl/data.frame) Named list() ..@ cluster_type: chr(0) ..@ equalized : logi TRUE ``` </details> この時点で`"cdr1"`・`"cdr3_nt"`があれば良い ```R cdb$contig_tbl = dplyr::filter(cdb$contig_tbl, full_length, productive == 'True', high_confidence, chain != 'Multi', str_length(cdr3) > MIN_CDR3_AA) %>% mutate( fancy_name = fancy_name_contigs(., str_c(pop, '_', sample))) ``` ```R aa80 = cdhit_ccdb(cdb, sequence_key = 'cdr3', type = 'AA', cluster_pk = 'aa80', identity = .8, min_length = 5, G = 1) aa80 = fine_clustering(aa80, sequence_key = 'cdr3', type = 'AA', keep_clustering_details = TRUE) ``` これシングルセル用のツールだった。。。 ### [ClusTCR](https://github.com/svalkiers/clusTCR) ```bash conda create -n clusTCR conda activate clusTCR conda install clustcr -c svalkiers -c bioconda -c pytorch -c conda-forge python ``` CDR3ベースのクラスタリングに加え、ClusTCRはV遺伝子情報をクラスタリングプロセスに含める機能を提供します。V遺伝子クラスタリングを有効にすると、TCR配列はまずV遺伝子ファミリーでソートされ、そのV遺伝子ファミリーに属する配列の各グループ内でクラスタリングが適用されます。このようにすることで、クラスタリングの精度は向上しますが、クラスタリングの保持率は低下します（クラスタに入る配列が少なくなる）。 .fit()メソッドでinclude_vgene = Trueを設定することで、V遺伝子情報をクラスタリング処理に含めることができます。 ```py from clustcr import Clustering clustering = Clustering() from clustcr import datasets cdr3 = datasets.test_cdr3() ``` ```py >>> cdr3.head() junction_aa v_call 0 CASSYLPGQGDHYSNQPQHF TRBV13*01 1 CASSFEAGQGFFSNQPQHF TRBV13*01 2 CASSFEPGQGFYSNQPQHF TRBV13*01 3 CASSYEPGQVSHYSNQPQHF TRBV13*01 4 CASSFGVEDEQYF TRBV7-2*01 output = clustering.fit(cdr3["junction_aa"]) ``` https://github.com/s175573/iSMART `v_call`というのは`Variable gene name in Imgt format: TRBVXX-XX*XX`らしい(別のツールでたまたま見つけた) ```py >>> output.clusters_df.head() junction_aa cluster 0 CSATNRDRGLEQYF 0 1 CSATSRDRGLEQYF 0 2 CASSEAGGTEAFF 1 3 CASSETGGTEAFF 1 4 CASSARTSGGSDTQYF 2 ``` 各アミノ酸配列(CDR3)に対してクラスター番号が割り当てられた。 ```py >>> output.summary().head() size motif 0 2 CSAT[SN]RDRGLEQYF 1 2 CASSE[AT]GGTEAFF 2 9 CASSrRTSGG[TA]DTQYF 3 6 CA.[SG]TGDSNQPQHF 4 8 CASSSA.YGYTF ``` 各クラスターのモチーフ配列を抽出することもできる。 ## メモいろいろ参考になる。↓ https://amelieff.esa.io/posts/1522 - レパトア解析の勉強 - https://www.jstage.jst.go.jp/article/naika/108/11/108_2347/_pdf > BCRでは，遺伝子再構成後の体細胞超突然変異(somatic hypermutation)により，CDR3 に加え，CDR1，CDR2 にも多様性が獲得される. - <img width="1022" alt="image.png (414.4 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/01/27/62720/a577a498-2cfb-4a14-943d-b1f34623e4c2.png"> - こういう風に得られたとして、CDR1のアミノ酸配列の類似度でクラスタリングがしたい？ - [ ] UpSet plotは配列の重複で使おうとしている？mutation入ったら何を持って重複？ - https://genble.co.jp/single-cell-repertoire - /https://www.repertoire.co.jp/research/technology/repertoire/ - [ ] 何がしたいのか - [ ] どういう手法？ - [ ] 何を見たい？ - 上記のCDR3の配列の長さを見て、それぞれのアミノ酸配列の長さが気になった。(「cdr1 cdr2 cdr3 length」で画像検索) - https://www.researchgate.net/figure/The-characteristics-of-CDRs-in-the-IgH-repertoire-of-the-neonates-and-the-adults-a_fig3_322935750 より([元論文](https://www.frontiersin.org/articles/10.3389/fimmu.2018.00128/full)) - <img width="300" alt="image.png (222.9 kB)" src="https://amelieff-esa-bucket.s3-ap-northeast-1.amazonaws.com/uploads/production/attachments/2943/2023/01/27/62720/bacbdaa4-e81b-48cb-8a22-42933b197af4.png">

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.