前回、Seuratのハンズオンを用いてSeuratObjectの作り方から、セルのクラスタリングまでを解説いたしました。
前編:https://labo-code.com/bioinformatics/seurat-1/
後編:https://labo-code.com/bioinformatics/seurat-2/
この記事では、Seuratの一歩踏み込んだ使い方として、各クラスターをバイオマーカーを元に細胞種をラベリングしてみたいと思います。
SingleRライブラリを使うことで、自動で細胞種を割り当てることができ、バイアスを排除できます。是非挑戦してみましょう。
macOS Monterey(12.4), クアッドコアIntel Core i7, メモリ32GB
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説!
RとSeuratで始めるSingle Cell RNA-seq解析!
今回の解析の目的
前回のハンズオンの最後にUMAPプロットを行い、下記の出力を得たかと思います。
上記では、8つのクラスターに分けることができています。しかしながら、
- 各クラスターは何細胞なのか?
- 各クラスターに特異的に発現するバイオマーカーはあるのか?
といった疑問が生まれると思います。これらの問題もSeuratでは解析することが可能です。
用いるデータの確認
では早速、解析を行っていきます。今回用いるデータセットは以下になります。
データセット:GSE149689
この論文では重症COVID-19の炎症進行要因でタイプIインターフェロン応答の役割についてSingle Cell RNA-seq解析を用いています。またこの論文で用いられたGSE149689は、下記の論文で再解析されています。
Early peripheral blood MCEMP1 and HLA-DRA expression predicts COVID-19 prognosis
この論文では重症COVID-19の予後バイオマーカー: CD14+細胞のMCEMP1およびHLA-DRA遺伝子発現を見るために、GSE149689を含む7つのデータセットをあわせて再解析しています。
このように、公共データベースに存在するscRNA-seqデータは再解析することで新たな論文として発表することが可能です。
今回は、Early peripheral blood MCEMP1 and HLA-DRA expression predicts COVID-19 prognosisにある、こちらのFig4bをGSE149689データセットをSeuratで解析して再現してみましょう。完全に同じにすることは難しいですが、似たFigは作ることができます。
Seuratの基本処理
GSE149689をダウンロードします。こちらからアクセスしてください。
ページ下部に行くと、Supplementary fileがあると思います。ftpかhttpをクリックして3つともすべてダウンロードしてください。
ダウンロードしたら、名前を、下記のようにbarcodes, features, matrixの名前だけ残して変更してください。
Rstudioを起動し、Working Directoryをセットします。
今回は、デスクトップに「GSE149689」フォルダを作成し、そちらに先程ダウンロードしたファイルを格納し、DesktopフォルダをWorking Directoryとしました。後にデータをロードする際にきちんとパス指定ができるなら問題ないです。
次に、SeuratObjectを作成していきます。エディターに以下をコーディングしてください。
install.packages("dplyr")
install.packages("Seurat")
install.packages("patchwork")
library(dplyr)
library(Seurat)
library(patchwork)
# Load the dataset
pbmc.data <- Read10X(data.dir = "./GSE149689/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
CreateSeuratObject
する際の各引数の説明は以下の通りです。ほぼデフォルト値を使用しています。
counts
: セルごとの遺伝子発現データを含む行列。この例では、**pbmc.data
**が入力されています。project
: 作成されるSeuratオブジェクトに付けられる名前。この例では、”PBMC”が指定されています。min.cells
: この値以上のセルで検出された遺伝子のみを含めるための閾値。この例では、最低3つのセルで検出された遺伝子が含まれます。min.features
: この値以上の遺伝子が検出されたセルのみを含めるための閾値。この例では、少なくとも200の遺伝子が検出されたセルが含まれます。
次にQuality Control(QC)と正規化を行っていきます。条件は論文に書かれているので見ていきましょう。マテリアルメソッドの以下の項目です。
Single-cell data processing and analysis
Integrative analysis of single-cell data was performed using the Seurat R package (Version 3), and single-cell visualisation was performed using Uniform Manifold Approximation and Projection (UMAP). During quality control, cells with a mitochondrial gene ratio of more than 15% were removed, which may be potentially dead cells. Only those cells with gene numbers in the range of 200–5000 or RNA numbers detected between 1000 and 30,000 cells were retained. After quality control, 48,583 PBMCs and 78,666 BALF single cells were included in the subsequent analysis. Data were normalised using the Seurat package and a principal component analysis was performed taking the top 2000 genes with the largest coefficient of variation.
QC・正規化周りで要点をまとめると以下になります。
- ミトコンドリア遺伝子比率が15%以上の細胞は除外
- 遺伝子数が200~5000の範囲、または検出されたRNA数が1000~30,000の細胞のみ保持
- Seuratパッケージを使用してデータを正規化
大体の論文でQC・正規化方法について記載されていますので論文を確認しつつ進めていきましょう。
まずはミトコンドリア遺伝子比率が15%以上の細胞は除外し、遺伝子数が200~5000の範囲、または検出されたRNA数が1000~30,000の細胞のみ保持し、最後に正規化を行うコードを書きます。
# Calculate percent mitochondrial genes
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# Quality control
pbmc <- subset(pbmc, nFeature_RNA >= 200 & nFeature_RNA <= 5000 & percent.mt < 15)
# Normalize data
pbmc <- NormalizeData(pbmc)
データセット内の変動遺伝子(高変動遺伝子)を検出します。これらの遺伝子は、細胞間の異なる生物学的状態を捉えるために使用されます。
# Find variable features
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
FindVariableFeatures関数では以下の引数が使われています。
selection.method
: 変動遺伝子を検出する方法。ここでは “vst”(variance stabilizing transformation)が選択されています。nfeatures
: 検出する変動遺伝子の数。ここでは上位2000個が選択されています。
データのスケーリング(正規化)を行います。これにより、各遺伝子の平均が0、標準偏差が1になります。スケーリングは、データの構造を維持しながら、異なる遺伝子やサンプル間での比較を容易にするために行われます。
# Scale data
pbmc <- ScaleData(pbmc)
RunPCA関数により主成分分析(PCA)を実行して、シングルセルデータの次元削減を行います。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
データの変動を最大限捉える新しい座標軸を見つけることで、高次元データを低次元に縮約します。これにより、データの可視化や解析が容易になります。引数は以下のとおりです。
features
: PCAを実行する際に使用する遺伝子のリスト。ここでは、FindVariableFeaturesで検出した変動遺伝子が選択されています。
主成分分析により次元削減を行った後の、各主成分の寄与率を見ていきましょう。JackStraw
とScoreJackStraw
で確認していきます。JackStrawは少し時間がかかります。
# Choose top 20 PCs
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
実際にElbowPlot
を使って確認していきます。
ElbowPlot(pbmc)
以下の図は、主成分分析の寄与を縦軸で表しています。PC20まで比較的寄与率が下がり続けているプロットが確認できるため、今回はPC20までを用いて後続の処理を行っていきましょう。
それでは、実際にシングルセルデータをクラスター化させていきます。FindNeighbors関数とFindClusters関数を用いて行います。
# Cluster cells
pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)
FindNeighbors関数では、先程主成分分析のときに確認したPC20までの情報をdims=1:20
で指定しています。FindClusters関数のresolution
は、クラスタリングの解像度を制御するパラメータで、0から1の範囲の値を取ります。解像度が高いほど、より多くの細かいクラスターが検出されますが、過剰なクラスタリングのリスクも高まりますので注意しましょう。
以上までの工程が終わったら、実際にUMAP処理を行っていきましょう。dims=1:20
を指定します。
# Run UMAP
pbmc <- RunUMAP(pbmc, dims = 1:20)
DimPlotで可視化します。
# Visualization
DimPlot(pbmc, reduction = "umap", label = T)
以下のようなUMAPプロットが得られたら成功です。
クラスター毎に自動で細胞種のラベルをつける
ここまで、実際のデータを公共データからダウンロードし、UMAPで表示するところまで行いました。しかし、いくつかクラスターが確認できたものの、それぞれのクラスターが何を表しているかを特定する必要があります。
今回は、celldexとSingleRとSingleCellExperimentライブラリを用いて各細胞に自動で細胞種のラベルをつける方法をご紹介します。各ライブラリの説明は下記になります。
celldex
:一連の単一細胞RNAシーケンシング(scRNA-seq)データセットのコレクション
SingleR
:細胞タイプのアノテーションが付与された参照データセットと、アノテーションが不明なサンプルデータを比較し、細胞タイプを推定
SingleCellExperiment
:Rプログラミング言語とBioconductorプロジェクト用に開発されたパッケージで、シングルセルRNAシーケンシング(scRNA-seq)データの格納、操作、および解析を容易にするためのデータ構造を提供
実際のコードはこちらになります。
# Check if BiocManager is installed, if not, install it
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# Install required packages using BiocManager
BiocManager::install("SingleR")
BiocManager::install("celldex")
BiocManager::install("SingleCellExperiment")
# Load required libraries
library(SingleR)
library(celldex)
library(SingleCellExperiment)
# Load reference data from celldex
ref <- celldex::HumanPrimaryCellAtlasData()
# Run SingleR to infer cell types of pbmc dataset using reference data
results <- SingleR(test = as.SingleCellExperiment(pbmc), ref = ref, labels = ref$label.main)
# Add inferred cell type labels to pbmc object
pbmc$singlr_labels <- results$labels
# Visualize cell types in a UMAP plot with labels
DimPlot(pbmc, reduction = 'umap', group.by = 'singlr_labels', label = TRUE)
コードを実行し、以下のようなデータが得られれば成功です。
コードを解説します。
# Check if BiocManager is installed, if not, install it
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# Install required packages using BiocManager
BiocManager::install("SingleR")
BiocManager::install("celldex")
BiocManager::install("SingleCellExperiment")
# Load required libraries
library(SingleR)
library(celldex)
library(SingleCellExperiment)
ライブラリのインストールを行っております。
# Load reference data from celldex
ref <- celldex::HumanPrimaryCellAtlasData()
HumanPrimaryCellAtlasData
関数を使って、celldexパッケージからHuman Primary Cell Atlasの参照データを取得し、ref
という変数に格納しています。
# Run SingleR to infer cell types of pbmc dataset using reference data
results <- SingleR(test = as.SingleCellExperiment(pbmc), ref = ref, labels = ref$label.main)
SingleR
関数を使って、pbmc
データセットの細胞タイプを推定しています。この推定には、1行目で取得したref
を参照データとして使用し、細胞タイプラベルとしてref$label.main
を指定しています。推定結果はresults
という変数に格納されます。
# Add inferred cell type labels to pbmc object
pbmc$singlr_labels <- results$labels
SingleR
で推定された細胞タイプラベルをpbmc
オブジェクトのsinglr_labels
というメタデータに追加しています。これにより、細胞タイプ情報が可視化や後続の解析で利用可能になります。
# Visualize cell types in a UMAP plot with labels
DimPlot(pbmc, reduction = 'umap', group.by = 'singlr_labels', label = TRUE)
DimPlot
関数を使って、UMAP(Uniform Manifold Approximation and Projection)次元削減技術を用いたプロットを作成しています。プロット上の細胞は、SingleR
で推定された細胞タイプラベル(singlr_labels
)に基づいてグループ分けされ、各グループにラベルが付けられます。これにより、データセット内の細胞タイプの分布を視覚的に確認できます。
FeaturePlotで実際のマーカー発現を確認する
先程、各クラスターを細胞種ごとに割り当てることができました。しかし、本当に正しく細胞種が割り当てられたか確認することが大切です。FeaturePlotは各遺伝子毎の発現をUMAPプロット上で確認できるので、こちらでマーカーの発現を確認してみましょう。
今回はmonocytesのマーカーとしてCD14、T-cellのマーカーとしてCD3Eを使っていきます。下記を実行してください。
FeaturePlot(pbmc, features = c("CD14", "CD3E"))
下記の結果が得られれば成功です。
monocytesとして割り当てられたクラスターにはCD14の発現が見られ、T-cellとして割り当てられたクラスターにはCD3Eの発現が見られたことが確認できると思います。このように各種マーカーの発現が見られるか実際にプロットして確かめてみることが大切です。
最後に
いかがだったでしょうか。論文記載のFigの完全再現とまではいかないですが、やりたい解析はできたと思います。ぜひとも挑戦してみてください!