【scRNA-seq】Seuratを用いてscRNA-seq解析を始める方法(前編)【Seurat】

【scRNA-seq】Seuratを用いてscRNA-seq解析を始める方法(前編)【Seurat】

single cell RNA-seq解析(scRNA-seq)とはRNA-seq解析を個々の細胞集団ごとに行う手法として注目を集めています。

今回はscRNA-seqのRパッケージであるSeuratのハンズオンに取り組むことで、scRNA-seqを始める前準備を行おうと思います。
本記事はSeuratのハンズオンに則って行いますので使い方の概要を理解するが可能です。ぜひ取り組んでみましょう。

動作検証済み環境
macOS Monterey(12.4), クアッドコアIntel Core i7, メモリ32GB

Seuratとは? 

SeuratはscRNA-seqデータ解析を行うRパッケージのことで、scRNA-seqデータの正規化やビジュアライズなどを行う事ができます。

Seuratには製作者たちによってハンズオンが用意されています。こちらに取り組むことでSeuratの理解を深めることができますが、Rに詳しくない方が取り組むには難しいところがあるため、こちらの記事では補足説明を加えつつ解説していきたいと思います。

R環境の用意

Seuratのハンズオンを行うにはRを実行する環境が必要となります。私的には、Rを実行する環境としてRstudioを使うことをおすすめします。

各OS環境を元にこちらのサイトよりインストーラーをダウンロードしてください。
https://cran.r-project.org/

インストールするバージョンは最新版を選んでください。一番トップにあるものが最新版です。

インストール画面では「続ける」を押して完了させてください。

Rstudioでは左上のエリアにコードを書いていきます。Runを押すと行実行ができます。
Sourceを押すとスクリプト全体が実行できますので覚えておいてください。

Seuratハンズオンをやってみる(前編)

ハンズオンの準備

それではSeuratのハンズオンをやっていきたいと思います。まずは必要なライブラリをインポートしてください。

install.packages("dplyr")
install.packages("Seurat")
install.packages("patchwork")

インストールが終わったら、以下の記述を記載してライブラリが使える状態にします。

library(dplyr)
library(Seurat)
library(patchwork)

生データをハンズオンサイトからダウンロードして来ます。
こちらにアクセスし、図のhereからダウンロードしてください。
https://satijalab.org/seurat/archive/v3.2/pbmc3k_tutorial.html#:~:text=data%20can%20be%20found%20here

ダウンロードしてきたファイルをデスクトップ上で解凍してください。その後Rstudioの右下画面でデスクトップまで進み、

  1. ダウンロードしてきたファイルにチェック
  2. Moreをクリック
  3. Set As Working Directory

を順に選びます。この作業を間違うと後のコードが一切通らなくなるのでお気をつけください。

以下のコマンドを叩いて最後にDesktopで終わるようなパスが出てきたらきちんと設定がなされています。

> getwd()
[1] "/Users/[ユーザー名]/Desktop"

以下のコマンドを実行し、データをRに読み込んでいきます。これでハンズオンの準備は完了です。

pbmc.data <- Read10X(data.dir = "./filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

クオリティチェック(QC)

scRNA-seqではバイオリンプロットを用いてクオリティチェックを行い、外れ値に存在する細胞がいるかを確認していきます。

細胞が死んでいるなど低いクオリイティのセルは遺伝子発現量が極端に低かったりミトコンドリア由来のコンタミが懸念されます。また遺伝子数カウントする際に細胞がダブる(シングルセルではない)と異常な遺伝子数でカウントされます。

# MT-で始まる全ての遺伝子の集合をミトコンドリア遺伝子の集合として使用する
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")


VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

こちらを打ち込むと以下のような図がRstudioの右下に現れると思います。

バイオリンプロットの膨らみが確認できる場所の細胞群を用いると全体としてもっともらしい遺伝子数カウントがされると判断します。今回のサンプルでのQCの参考値はいかになります。

  • nFeature_RNAより、ユニークフィーチャー数が2,500以上、または200以下の細胞をフィルタリングします。
  • percent.mtよりミトコンドリア数が5%以上の細胞をフィルタリングします。

それぞれの相関関係を見るために散布図を利用することができます。下記のコードを実行してみてください。

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

percent.mtnCount_RNAでは相関が見られません(-0.13)が、nFeature_RNAnCount_RNAとは相関関係がある(0.95)ことがわかるかと思います。

QC後に決まった値を再度サンプルpbmcに入れておきます。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

データの正規化

QCによってデータセットから不要な細胞を除去した後は、データの正規化を行います。

デフォルトでは、各細胞の特徴的な発現量を総発現量で正規化し、これにスケールファクター(デフォルトでは10,000)を掛け、結果を対数変換するグローバルスケーリング正規化手法 “LogNormalize” が採用されています。
データの正規化後に決まった値を再度サンプルpbmcに入れておきます。

# LogNormalize
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

# pbmcに再度入れる
pbmc <- NormalizeData(pbmc)

ここまでできれば、scRNA-seqデータの前処理は完了になります。後の処理は正規化されたデータの可視化方法の解説になります。

変動量の大きい遺伝子を抽出

ある細胞では高発現、別の細胞では低発現しているデータセット中の細胞間変動が大きいサブセットを計算します。デフォルトでは、1データセットあたり2,000のVariable countをを返すようになっています。

# FindVariableFeatures関数を使ってデータを用意
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# top10の遺伝子名を記述のため、headを取得
top10 <- head(VariableFeatures(pbmc), 10)

# プロット
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

もしコマンド通りに叩いてもエラーが出る場合は…

コマンド通りに叩いてもエラーが出る場合は、pbmc に入る値が間違っています。エラーが出ても慌てずに、クオリティチェック→正規化…とコマンドを最初から叩き直してください。順番通りきちんと入力すればエラーなく通ります。

本記事では、一旦ここまでにしようと思います。

最後に

いかがだったでしょうか。ハンズオンは折返しになります。次の記事も取り組んで一通りSeuratを使えるようにしてみましょう。

コメントを残す

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