scRNA-seqの解析に用いられるRパッケージのSeuratについて、ホームページにあるチュートリアルに沿って解説（和訳）していきます。ちゃんと書いたら長くなってしまいました。
あくまで自分の理解のためのものです。足らない所はご了承ください。
参考、というか大元→Sajita Lab
参考→やってみよう！ single cell RNA-seq 解析ことはじめ
##0. インストール
まずR環境を用意します。Rstudioもあると便利です。
インストール方法、使い方はこちら→統計用言語Rの使い方
RもしくはRstudioで、Seuratをインストールしてみましょう。
# Enter commands in R (or R studio, if installed)
> install.packages('Seurat')
> library(Seurat)
warning messageが出てきてもyと打てば大丈夫だそうです。
##1. セットアップ(データの読み込み)
今回のチュートリアルでは、PBMCが提供する10X Genomics のデータを使います。上のリンク先のページからダウンロードして決めた所においておきましょう。
10XのデータはRead10X関数を用いて読み込むことができます。
##先に必要なライブラリを読み込みます
##dplyrはデータフレームの操作、patchworkはプロット用のライブラリです
> library(dplyr)
> library(Seurat)
> library(patchwork)
# Load the PBMC dataset
## データがあるディレクトリを指定
> pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
これによって読み込んだデータpbmc.dataは行列になっていて、行は遺伝子、列は各細胞に対応しています。ちなみに中をみるとこんな感じです。
> pbmc.data[1:3, 1:3]
3 x 3 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1
MIR1302-10 . . .
FAM138A . . .
OR4F5 . . .
HDF5ファイル~.h5に対応したRead10X_h5というのもあるのですが、HDF5ファイルが読み込めるようにしておく必要があります。(エラーが出なかったら大丈夫)。
（私はダメだったので、(なぜか？)Homebrewでインストールしました。）
# in terminal
> brew install hdf5
# in R/R studio
> install.packages("hdf5r")
> pbmc.counts <- Read10X_h5("path to the file.h5")
データの読み込みは10X以外のものでも、組み込み関数を使って読み込むことが可能です。(やってみよう！ single cell RNA-seq 解析ことはじめより)
> df <- read.table(file ='PATH/TO/YOUR/GENE/EXPRESSION/MATRIX/FILE', #ファイルの場所を指定
sep = ',', #csvファイルなので、カンマ区切りを指定
#タブ区切りの場合は、”¥t”
header = TRUE, #先頭の行が列名の時は、TRUEに
row.names = 1, #1列目が行名であることを指定
stringsAsFactors = FALSE) #文字列が誤変換されないように
次に、読み込んだデータをSeurat用のオブジェクトに変換する必要があります。
# Initialize the Seurat object with the raw (non-normalized data)
> pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
これにより、データを保持しながらSeuratで解析が可能となります。
オプションの説明は以下の通りです。
counts →元データ
project →データの名前を指定
min.cells →細胞による遺伝子フィルター。=3であれば、3つ以上の細胞で発現が確認された遺伝子のみを解析に用いる。
min.features →遺伝子による細胞フィルター。=200であれば、200種類以上の遺伝子の発現がみられた細胞のみを解析に用いる。
その他、様々な調整オプションがありますが、詳しく知りたい場合は、RstudioでPackageタブ(右下の画面にある)からSeuratを探し、その中にあるCreateSeuratObject()をクリックすると、ドキュメンテーションが確認できます。
Seuratオブジェクトがどうなっているかは以下のようにして知ることができます。
> pbmc
An object of class Seurat
13714 features across 2700 samples within 1 assay
Active assay: RNA (13714 features, 0 variable features)
# Cell-level meta data is stored as a data frame
# Standard data frame functions work on the meta data data frame
##metadataを表示
##meta dataはデータフレームの形で格納されています。
##nCount -> RNAカウント数
##nFeature -> 検出された遺伝子数
> colnames(x = pbmc[[]])
[1] "orig.ident" "nCount_RNA" "nFeature_RNA"
> head(x = pbmc[[]])
orig.ident nCount_RNA nFeature_RNA
AAACATACAACCAC-1 pmbc3k 2419 779
AAACATTGAGCTAC-1 pmbc3k 4903 1352
AAACATTGATCAGC-1 pmbc3k 3147 1129
AAACCGTGCTTCCG-1 pmbc3k 2639 960
AAACCGTGTATGCG-1 pmbc3k 980 521
AAACGCACTGGTAC-1 pmbc3k 2163 781
##一部分を見ることも可能です
> head(x = pbmc[[c('orig.ident', 'nCount_RNA')]])
orig.ident nCount_RNA
AAACATACAACCAC-1 pmbc3k 2419
AAACATTGAGCTAC-1 pmbc3k 4903
AAACATTGATCAGC-1 pmbc3k 3147
AAACCGTGCTTCCG-1 pmbc3k 2639
AAACCGTGTATGCG-1 pmbc3k 980
AAACGCACTGGTAC-1 pmbc3k 2163
##count matrixを見る
##普通にやると一気に表示されてしまうので今回は5×5にしています。
> pbmc[["RNA"]]@counts[1:5, 1:5]
5 x 5 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
AL627309.1 . . . . .
AP006222.2 . . . . .
RP11-206L10.2 . . . . .
RP11-206L10.9 . . . . .
LINC00115 . . . . .
##以下GitHubにのっているコマンドです。
> dim(x = pbmc)
[1] 13714 2700
# In addtion to rownames and colnames, one can use dimnames
# which provides a two-length list with both rownames and colnames
##各行各列について、先頭5つの名前を表示
> head(x = rownames(x = pbmc))
[1] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" "LINC00115"
[6] "NOC2L"
> head(x = colnames(x = pbmc))
[1] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1"
[5] "AAACCGTGTATGCG-1" "AAACGCACTGGTAC-1"
# A vector of names of associated objects can be had with the names function
# These can be passed to the double [[ extract operator to pull them from the Seurat object
##データや解析結果は以下のような形で格納されています。
> names(x = pbmc)
[1] "RNA"
> pbmc[["RNA"]]
Assay data with 13714 features for 2700 cells
First 10 features:
AL627309.1, AP006222.2, RP11-206L10.2, RP11-206L10.9, LINC00115, NOC2L, KLHL17,
PLEKHN1, RP11-54O7.17, HES4
#FindVariableFeatures()後の結果が出ています、、
現段階で見られるのはこの辺りでしょうか。詳しくは(コード中にも書いていますが) Git Hubを見てください。
##2. 前処理
データを読み込んだら、解析の前に前処理を行います。
###2-1. Quality Check
Seuratでは様々なQC用の関数が用意されています。
single cellシークエンスではドロップレット中に1つの細胞が入るという前提のもとで配列を読んでいきますが、うまくいかないものもあります。
・ドロップレット中に細胞がなく、遺伝子が検出できない
・逆に1つのドロップレットに複数の細胞が入ってしまい、発現数が異常にあがる
・死んだ細胞が取り込まれてしまった
これらのデータをmeta dataの情報から取り除きます。
先に、ミトコンドリアゲノムの割合の情報をmetadataに追加しましょう。これは上でいう3番目の死んだ細胞を検出するのに用いられます。死んだ細胞ではミトコンドリアのゲノムが比較的多く検出されるそうです。
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
##PercentageFeatureSet()は特定の遺伝子の発現数の割合を計算する関数です。
##今回はミトコンドリアゲノムの割合を調べるので遺伝子名がMT-から始まるものを選びます。
> pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
percent.mtデータができたか確認してみましょう。
> head(pbmc@meta.data, 5)
orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1 pmbc3k 2419 779 3.0177759
AAACATTGAGCTAC-1 pmbc3k 4903 1352 3.7935958
AAACATTGATCAGC-1 pmbc3k 3147 1129 0.8897363
AAACCGTGCTTCCG-1 pmbc3k 2639 960 1.7430845
AAACCGTGTATGCG-1 pmbc3k 980 521 1.2244898
これらは全てのmetadetaです。このデータから不要なものをとりのぞいていきましょう！
まず、metadataの概要を掴むため、Vlnplot で可視化してみましょう。Vlnplotは特徴量の分布を見ることができます。
# Visualize QC metrics as a violin plot
> VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
このコマンドを打つと以下のような図が右側に表示されます(Rstudio)。
それぞれの特徴量間の相関もみてみましょう。
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
> plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
> plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
> plot1 + plot2
この可視化を行い、使う特徴量や境界を決めます。
とりあえずチュートリアルと同じようにnFeature_RNAが(200, 2500)のものと、percent.mtが5未満のものを選びます。
> pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
###2-2. 正規化
次に正規化(normalize)を行います。使用する関数はNormalizeData()です。デフォルトの手法(Method)はLogNormilizeを使用します。(ドキュメント)これは各遺伝子の発現量をその細胞の総発現量でわり、scale factor(デフォルトは10000)を掛けたあと、対数に直すというものです。Methodは他にもCLR(有心対数変換)やRC(対数変換しない)があります。
> pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
##単純にNormalizeData（）でも同じ
正規化されたデータはpbmc[["RNA"]]@dataに格納されます。
###2-3. Feature selection
細胞ごとに発現量が変動しやすい遺伝子はPCAやクラスタリング、細胞の決定などで重要な役割を持ちます。この遺伝子を見つけてくれるのがFindVariableFeatures（）になります(ドキュメント)。デフォルトは手法がvstでnfeatures = 2000の遺伝子を抽出します。
ドキュメントによると、vstは
vst: First, fits a line to the relationship of log(variance) and log(mean) using local polynomial regression (loess). Then standardizes the feature values using the observed mean and expected variance (given by the fitted line). Feature variance is then calculated on the standardized values after clipping to a maximum (see clip.max parameter).
とあるのですが、平均と分散についての回帰を行い、これをもとに分散を予測して大きいものを選ぶということでしょうか、、
とりあえず、実行とその結果をプロットしたものが以下になります。
> pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
このプロットでは特に分散の大きい10遺伝子の名前も表示されています。
###2-4. Scaling
次にスケーリングを行います。これは後々の次元削減の作業で必要になります。
使う関数としてはScaleData()で、
・発現量の平均を0、分散を1にする。
・計算後のデータをpbmc[["RNA"]]@scale.dataに格納
を行います。
> all.genes <- rownames(pbmc)
> pbmc <- ScaleData(pbmc, features = all.genes)
featuresのデフォルトはfeatures=2000になっていて、2-3で抽出した2000個の遺伝子について行います。今回は全ての遺伝子についてスケーリングしました。
結果は以下のコマンドで確認できます。
> GetAssayData(object = pbmc, slot = "scale.data")
##3. PCA
ここから実際のデータ解析に入ります。まずPCAを行うRunPCA()をやってみます。2-3で抽出した遺伝子でPCAをします。
> pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
# Examine and visualize PCA results a few different ways
##PCAの結果を確認します。
##上での実行結果はpbmc[["pca"]]に格納されています。
##各主成分での特徴量を表示
> print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
PC_ 1
Positive: CST3, TYROBP, LST1, AIF1, FTL
Negative: MALAT1, LTB, IL32, IL7R, CD2
PC_ 2
Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
Negative: NKG7, PRF1, CST7, GZMB, GZMA
PC_ 3
Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
Negative: PPBP, PF4, SDPR, SPARC, GNG11
PC_ 4
Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
Negative: VIM, IL7R, S100A6, IL32, S100A8
PC_ 5
Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
Negative: LTB, IL7R, CKB, VIM, MS4A7
結果を可視化する場合は以下の関数を利用します。
##第1、2主成分での特徴量とその係数を表示
> VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
##各細胞を第1主成分、第2主成分でプロット
> DimPlot(pbmc, reduction = "pca")
各主成分での特徴量(遺伝子)と各細胞での発現量を見るには、HeatMapを作成します。
> DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
dimsは表示する主成分の指定です。= 1:15とすれば一気に15主成分見ることができます。
cells = 500というのは発現のコントラストが大きいものから500個という指定になります。
このHeatmapから関連のある遺伝子(群)を見つけることもできそうです。
balancedは
Plot an equal number of genes with both + and - scores.
ということらしいです。(ドキュメントより)
主成分分析はできましたが、クラスタリングをする際には第何主成分までとるかというのが重要になりますよね。
これを調べるには、以下のコマンドで可能です。
# NOTE: This process can take a long time for big datasets, comment out for expediency. More
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce computation time
> pbmc <- JackStraw(pbmc, num.replicate = 100)
> pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
##可視化
> JackStrawPlot(pbmc, dims = 1:15)
これはJackstraw法(これについてはJackStraw Plot (QQ-plot) を理解するを見るとよくわかります)という手法を用いているのですが、コードの注意に書かれている通り、時間がかかります。
ので、データが多い時は(いつもの？) ElbowPlot()が楽だよ、ということみたいです(?)。
基本的に使うデータセットを引数(object=)にすればOKです。
> ElbowPlot(pbmc)
これを見ると9-10あたりまで使うといいねという感じ(だとチュートリアルには書いている)です。
次元数の決定は難しく、チュートリアルにのっているアドバイスとして、
・PC12,13になって特徴的になる細胞があったりする
・いろんな次元数で解析を行う(10,20,50など)
・なるべく多めの次元数で解析を行う
とありました。
##4. クラスタリング
Seuratではグラフを用いたクラスタリングを行っています。
まずはじめにPCA空間でユークリッド距離に基づくKNNグラフを作成し、近隣にあたる2細胞のエッジの重みを更新します(FindNeighbors())。これには前で決めた次元数のデータを用います。
クラスタリングには Louvain algorithm (デフォルト) やSLMといった手法を用いて行われます。使う関数のFindClusters()は resolution パラメータでクラスターの数を決めることができます。 3000個の細胞データをクラスタリングするときは 0.4-1.2ぐらいがいいそうです。
> pbmc <- FindNeighbors(pbmc, dims = 1:10)
Computing nearest neighbor graph
Computing SNN
> pbmc <- FindClusters(pbmc, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 2638
Number of edges: 96033
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8720
Number of communities: 9
Elapsed time: 0 seconds
上の操作でできたクラスターはIdents()で見ることができます。
# Look at cluster IDs of the first 5 cells
> head(Idents(pbmc), 5)
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
1 3 1 2 6
Levels: 0 1 2 3 4 5 6 7 8
##5. 非線形次元圧縮(UMAP, tSNE)
精度がいいとされているUMAPやtSNEにも対応しています。inputには3.で抽出した主成分情報を用います。
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages ='umap-learn')
> pbmc <- RunUMAP(pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters
DimPlot(pbmc, reduction = "umap")
##6. クラスターのバイオマーカーを見つける
クラスタリングで得られたクラスターで特徴的な遺伝子を見つけます。これをすることで、そのクラスタがどの細胞にあたるかを予測できます。
# find all markers of cluster 1
> cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
> head(cluster1.markers, n = 5)
p_val avg_logFC pct.1 pct.2 p_val_adj
IL32 1.894810e-92 0.8373872 0.948 0.464 2.598542e-88
LTB 7.953303e-89 0.8921170 0.981 0.642 1.090716e-84
CD3D 1.655937e-70 0.6436286 0.919 0.431 2.270951e-66
IL7R 3.688893e-68 0.8147082 0.747 0.325 5.058947e-64
LDHB 2.292819e-67 0.6253110 0.950 0.613 3.144372e-63
ident.1でクラスターの指定、min.pctで使う遺伝子の基準を設定(2つのグループのどちらかで、割合がmin.pct以上の発現が確認できた遺伝子のみ用いる)しています。似たようなパラメータは他にも、thresh.test max.cells.per.identがあります。これらは計算高速化に関わります。
ident.2で比較対象のクラスターを指定します。
# find all markers distinguishing cluster 5 from clusters 0 and 3
> cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
> head(cluster5.markers, n = 5)
p_val avg_logFC pct.1 pct.2 p_val_adj
FCGR3A 7.583625e-209 2.963144 0.975 0.037 1.040018e-204
IFITM3 2.500844e-199 2.698187 0.975 0.046 3.429657e-195
CFD 1.763722e-195 2.362381 0.938 0.037 2.418768e-191
CD68 4.612171e-192 2.087366 0.926 0.036 6.325132e-188
RP11-290F20.3 1.846215e-188 1.886288 0.840 0.016 2.531900e-184
全部のクラスターを一気に調べたい時は、FindAllMarkers()を使います。
logfc.threshold()も高速化用と思ってください。
# find markers for every cluster compared to all remaining cells, report only the positive ones
> pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
> pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
# A tibble: 18 x 7
# Groups: cluster [9]
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
1 1.96e-107 0.730 0.901 0.594 2.69e-103 0 LDHB
2 1.61e- 82 0.922 0.436 0.11 2.20e- 78 0 CCR7
3 7.95e- 89 0.892 0.981 0.642 1.09e- 84 1 LTB
4 1.85e- 60 0.859 0.422 0.11 2.54e- 56 1 AQP3
5 0. 3.86 0.996 0.215 0. 2 S100A9
6 0. 3.80 0.975 0.121 0. 2 S100A8
7 0. 2.99 0.936 0.041 0. 3 CD79A
8 9.48e-271 2.49 0.622 0.022 1.30e-266 3 TCL1A
9 2.96e-189 2.12 0.985 0.24 4.06e-185 4 CCL5
10 2.57e-158 2.05 0.587 0.059 3.52e-154 4 GZMK
11 3.51e-184 2.30 0.975 0.134 4.82e-180 5 FCGR3A
12 2.03e-125 2.14 1 0.315 2.78e-121 5 LST1
13 7.95e-269 3.35 0.961 0.068 1.09e-264 6 GZMB
14 3.13e-191 3.69 0.961 0.131 4.30e-187 6 GNLY
15 1.48e-220 2.68 0.812 0.011 2.03e-216 7 FCER1A
16 1.67e- 21 1.99 1 0.513 2.28e- 17 7 HLA-DPB1
17 7.73e-200 5.02 1 0.01 1.06e-195 8 PF4
18 3.68e-110 5.94 1 0.024 5.05e-106 8 PPBP
これらのdifferential expression解析の精度を評価するためのパラメータとしてtest.use()があります。これを使うと
> cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
> head(cluster0.markers, n = 5)
myAUC avg_diff power pct.1 pct.2
RPS12 0.822 0.5029988 0.644 1.000 0.991
RPS27 0.822 0.5020359 0.644 0.999 0.992
RPS6 0.820 0.4673635 0.640 1.000 0.995
RPL32 0.815 0.4242773 0.630 0.999 0.995
RPS14 0.807 0.4283480 0.614 1.000 0.994
という風にpower(0~1, 高い方がいい)で評価してくれます。
differential expressionについてはここに詳しくまとめられています。
可視化も様々な関数が用意されています。
> VlnPlot(pbmc, features = c("CD79A", "TCL1A"))
> FeaturePlot(pbmc, features = c("LDHB", "LTB", "S100A9", "CD79A", "CCL5", "FCGR3A", "GZMB", "FCER1A", "PF4"))
選んだ遺伝子は先ほどの、differential expressionで選ばれた、各クラスターで特徴的な遺伝子です。
Heatmapもクラスターごとで可能です。
> top5 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_logFC)
> DoHeatmap(pbmc, features = top5$gene) + NoLegend()
チュートリアルは10でやってたけど、多すぎてプロットが密すぎたので5でやったらこんな感じになりました、、
##7.細胞種の同定
最後に、differential expressionで選ばれた、各クラスターで特徴的な遺伝子から、そのクラスターがどの細胞種なのか推定します。
ただし、手動です。(Seuratにはデータベースを参照するとかいうのはない)
その遺伝子についてどのような機能に関わっているかを調べましょう。
今回はチュートリアルのをそのままいただいて、各クラスターに名前を設定し、可視化します。
> new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
> names(new.cluster.ids) <- levels(pbmc)
> pbmc <- RenameIdents(pbmc, new.cluster.ids)
> DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
完成です。
※上記のコマンドをまとめたもの、データ取得コマンドなど→コマンドリスト