【シングルセル】Harmonyを使ったSingle Cell RNA-seqデータのバッチエフェクトの補正【Single Cell RNA-seq】

【シングルセル】Harmonyを使ったSingle Cell RNA-seqデータのバッチエフェクトの補正【Single Cell RNA-seq】

シングルセルのバッチエフェクトの補正のやり方はご存知でしょうか?

この記事では、シングルセルRNA-seqデータの異なる条件からの統合についてHarmonyという手法を使用します。以前に紹介したバッチ補正の手法では、Seuratの正準相関解析法を使用してデータを統合し、バッチエフェクトを補正しました。しかし、今回はHarmonyを使用してみます。

是非挑戦してみましょう。

Tested Environment

macOS Monterey (12.4), Quad-Core Intel Core i7, Memory 32GB seurat 4.3.0 harmony 1.2.1

バッチ補正とは?


バッチエフェクトとは、非生物学的な要因、つまり実験のセットアップや手順、機器の違いなどによって生じるデータのバリエーションのことを指します。これは、異なるバッチで生成されたデータを比較するときに特に問題となります。バッチエフェクトは、結果の解釈を歪め、偽の相関を生じさせる可能性があります。したがって、scRNA-seqデータを解析する際には、バッチエフェクトを適切に補正することが非常に重要です。

Harmonyとは?


Harmonyは、各イテレーションで各クラスタのセルが可能な限り多くのバッチから来るようにする反復クラスタリングメソッドを使用し、各セルに適用する補正係数を計算します。このプロセスを繰り返してイテレーションごとに洗練し、収束するまで繰り返します。

https://portals.broadinstitute.org/harmony/articles/quickstart.html

Quick start to Harmony Korsunsky et al.: Fast, sensitive, and accurate integration of single cell data with Harmony

Harmonyで使用するメリットは主に以下です。

  1. バッチエフェクトの補正: Harmonyは、多標本シングルセルデータのバッチエフェクトを補正し、データの一貫性を改善します。これは、生物学的な変動を維持しながら、非生物学的なバッチエフェクトを除去するために非常に重要です。
  2. データ統合: Harmonyは、異なるバッチや実験条件から得られた複数のデータセットを統合する能力があります。これにより、大規模なデータセットの解析が可能になり、細胞の種類や状態を横断的に比較することができます。

データセットの紹介


今日使用するデータセットはループスデータセットです。ループスは自己免疫疾患であり、免疫系の機能に影響を与え、免疫系が体の臓器や組織を攻撃することがあります。

今日使用するデータセットは、8人のループス患者からの末梢血単核細胞です。これらのデータセットは刺激されたグループと対照グループに分割されており、刺激されたグループはインターフェロンベータで処理されています。

今回の分析の目的は、条件ごとにデータを統合し、両方の条件で類似した細胞を重ね合わせることです。今回は以下の流れでHarmonyを解説します。

  • データを統合するためにHarmonyを使用する。
  • SeuratObjectまたはシングルセル実験オブジェクトを直接使用する。
  • データセットへのアクセスにはSeuratDataパッケージを使用する。
  • Seurat、Tidyverse、ggplot2などの通常のパッケージを使用する。

Harmonyによるバッチエフェクト補正のやり方


それではHarmonyによるバッチエフェクト補正のやり方を解説していきます。

まずは必要なライブラリをインストール&呼び出しをしてください。

install.packages("harmony")
install.packages("Seurat")
install.packages("SeuratData")
install.packages("tidyverse")
install.packages("ggplot2")
install.packages("devtools")
devtools::install_github("satijalab/seurat-data")
# set seed for reproducibility
set.seed(1234)

library(harmony)
library(Seurat)
library(SeuratData)
library(tidyverse)
library(ggplot2)

SeuratData利用可能なデータは以下です。今回はifnb を使っていきます。

> # get data -------------------------
> AvailableData()
                                   Dataset Version
adiposeref.SeuratData           adiposeref   1.0.0
bmcite.SeuratData                   bmcite   0.3.0
bonemarrowref.SeuratData     bonemarrowref   1.0.0
cbmc.SeuratData                       cbmc   3.1.4
celegans.embryo.SeuratData celegans.embryo   0.1.0
fetusref.SeuratData               fetusref   1.0.0
hcabm40k.SeuratData               hcabm40k   3.0.0
heartref.SeuratData               heartref   1.0.0
humancortexref.SeuratData   humancortexref   1.0.0
ifnb.SeuratData                       ifnb   3.1.0
kidneyref.SeuratData             kidneyref   1.0.1
lungref.SeuratData                 lungref   2.0.0
mousecortexref.SeuratData   mousecortexref   1.0.0
panc8.SeuratData                     panc8   3.0.2
pancreasref.SeuratData         pancreasref   1.0.0
pbmc3k.SeuratData                   pbmc3k   3.1.4

データロードを行います。

InstallData("ifnb")
ifnb <- LoadData("ifnb")
ifnb <- UpdateSeuratObject(ifnb)
str(ifnb)

ここからは解析コードを一気に流していってください。

# QC and filtering
ifnb$mito.percent <- PercentageFeatureSet(ifnb, pattern = '^MT-')
View(ifnb@meta.data)
# explore QC

# filter
ifnb
ifnb.filtered <- subset(ifnb, subset = nCount_RNA > 800 &
                          nFeature_RNA > 200 & 
                          mito.percent < 5)

# standard workflow steps
ifnb.filtered <- NormalizeData(ifnb.filtered)
ifnb.filtered <- FindVariableFeatures(ifnb.filtered)
ifnb.filtered <- ScaleData(ifnb.filtered)
ifnb.filtered <- RunPCA(ifnb.filtered)
ElbowPlot(ifnb.filtered)
ifnb.filtered <- RunUMAP(ifnb.filtered, dims = 1:20, reduction = 'pca')

before <- DimPlot(ifnb.filtered, reduction = 'umap', group.by = 'stim')

# run Harmony -----------
ifnb.harmony <- ifnb.filtered %>%
  RunHarmony(group.by.vars = 'stim', plot_convergence = FALSE)

ifnb.harmony@reductions

ifnb.harmony.embed <- Embeddings(ifnb.harmony, "harmony")
ifnb.harmony.embed[1:10,1:10]

# Do UMAP and clustering using ** Harmony embeddings instead of PCA **
ifnb.harmony <- ifnb.harmony %>%
  RunUMAP(reduction = 'harmony', dims = 1:20) %>%
  FindNeighbors(reduction = "harmony", dims = 1:20) %>%
  FindClusters(resolution = 0.5)

# visualize 
after <- DimPlot(ifnb.harmony, reduction = 'umap', group.by = 'stim')

before|after

以下の出力となれば成功です!

コードの解説


1. QC(品質管理)とフィルタリング

ifnb$mito.percent <- PercentageFeatureSet(ifnb, pattern = '^MT-')
View(ifnb@meta.data)
  • PercentageFeatureSet(): ミトコンドリア遺伝子(MT-で始まる)を持つ遺伝子の割合を計算します。mito.percentという新しいメタデータ列に格納されます。これは、細胞の品質管理において重要な指標です。ミトコンドリア遺伝子の割合が高いと、低品質な細胞(細胞死やストレスの影響を受けたもの)の可能性が高くなります。
  • View(): メタデータを視覚化するために使用します。@meta.dataスロットは、各細胞に関連するメタデータ(例: 細胞バーコード、UMIカウント、ミトコンドリア遺伝子の割合など)を含んでいます。

2. フィルタリング

ifnb.filtered <- subset(ifnb, subset = nCount_RNA > 800 &
                          nFeature_RNA > 200 &
                          mito.percent < 5)
  • subset(): 細胞のフィルタリングを行います。この例では、次の基準に基づいて細胞をフィルタリングしています。
    • nCount_RNA > 800: UMIカウント(細胞ごとに検出されたmRNAの数)が800より多い細胞のみを保持します。
    • nFeature_RNA > 200: 検出された遺伝子の数が200より多い細胞のみを保持します。
    • mito.percent < 5: ミトコンドリア遺伝子の割合が5%未満の細胞のみを保持します。

このステップで、低品質の細胞やストレスを受けた細胞を除去します。

3. 標準的な解析ワークフロー

ifnb.filtered <- NormalizeData(ifnb.filtered)
ifnb.filtered <- FindVariableFeatures(ifnb.filtered)
ifnb.filtered <- ScaleData(ifnb.filtered)
ifnb.filtered <- RunPCA(ifnb.filtered)
ElbowPlot(ifnb.filtered)
  • NormalizeData(): データを正規化します。これは、細胞ごとの読み取り数の差を補正するために行われます。
  • FindVariableFeatures(): 解析のために重要な変動遺伝子を検出します。このステップでは、バイオロジカルな違いを捉えるために、変動の大きい遺伝子を特定します。
  • ScaleData(): データをスケーリング(標準化)します。これは、遺伝子ごとの平均値をゼロに、標準偏差を1に揃える処理です。
  • RunPCA(): 主成分分析(PCA)を実行します。データの次元を削減し、主要な情報を保持します。
  • ElbowPlot(): PCAの次元数を選ぶためのプロットを表示します。このプロットで「ひじ」部分が見える箇所が、適切な次元数の選択基準となります。

4. Harmonyによるバッチ効果補正

ifnb.harmony <- ifnb.filtered %>%
  RunHarmony(group.by.vars = 'stim', plot_convergence = FALSE)
  • RunHarmony(): バッチ効果補正のためにHarmonyを使用します。この例では、stim(刺激条件)によってグループ分けされている細胞群のバッチ効果を補正します。plot_convergence = FALSEで、収束プロットを表示しないようにしています。
ifnb.harmony.embed <- Embeddings(ifnb.harmony, "harmony")
ifnb.harmony.embed[1:10,1:10]
  • Embeddings(): Harmonyによって生成された低次元埋め込み(embedding)を取得します。PCAのように、データの次元が圧縮されていますが、バッチ効果が取り除かれています。

5. UMAPとクラスタリング(Harmonyの埋め込みを使用)

ifnb.harmony <- ifnb.harmony %>%
  RunUMAP(reduction = 'harmony', dims = 1:20) %>%
  FindNeighbors(reduction = "harmony", dims = 1:20) %>%
  FindClusters(resolution = 0.5)
  • RunUMAP(): Harmonyの埋め込みを使用して、UMAP(次元削減手法)を実行し、細胞を低次元空間に可視化します。
  • FindNeighbors(): 細胞の近傍(隣接関係)を計算し、クラスタリングの準備をします。
  • FindClusters(): 細胞をクラスタリングします。resolution = 0.5はクラスタリングの粒度を指定します。値を大きくすると細かく、小さくすると大きなクラスタになります。

6. 可視化

before <- DimPlot(ifnb.filtered, reduction = 'umap', group.by = 'stim')
after <- DimPlot(ifnb.harmony, reduction = 'umap', group.by = 'stim')
  • DimPlot(): UMAPによる次元削減の結果をプロットします。group.by = 'stim'は、stim(刺激条件)ごとに色分けして可視化します。
  • beforeafter の比較: beforeにはバッチ効果が含まれたデータ、afterにはHarmonyでバッチ効果が補正されたデータが表示されます。この比較でバッチ効果がどれだけ除去されたかを視覚的に確認できます。

7. プロットの比較

before | after
  • before | after: 2つのプロットを並べて表示します。Harmonyを適用する前後のUMAPプロットを比較することで、クラスタの変化やバッチ効果の影響を評価します

最後に


いかがだったでしょうか?特にHarmonyを用いたバッチ効果補正は、異なる条件や実験バッチ間での系統的な違いを取り除くために重要なステップです。しっかりマスターしていきましょう。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です