1. Seuratのインストール
まず、RにSeuratをインストールしないといけません。Seuratの他にHdf5ファイルを読み込むためのhdf5rもインストールしましょう。インストールのコマンドは以下です。
1.1 Biocondaでインストール（おすすめ）
何かをインストールするとき、Biocondaでインストールできればそれが一番簡単です。ターミナルの画面で、R環境に入ってない状態で、以下のコマンドを入力します。
conda install -c conda-forge r-hdf5r -y
conda install -c bioconda r-seurat
1.2 R環境下でインストール
Biocondaじゃなくても、Rのinstall.packagesコマンドでSeuratをインストール可能です。ターミナルで”R”と入力すると、Rの環境に入るはずです。R環境下で以下のコマンドを入力します。
install.packages('Seurat')
install.packages("hdf5r")
上記のいずれかで、Seuratがちゃんとインストールされたかを確認するにはR環境下で、以下のコマンドを入力します。
library(Seurat)
このコマンドを入力して何もエラーが出なければちゃんとSeuratがインストールされ、呼び出されています。
2. hdf5ファイルをSeurat objectに変換
いよいよ本番の解析です。現在、singleセルシーケンスの多くは10x Genomics社のChromiumプラットフォームを用いて行われています。なので、今回は10x Genomics社の解析ツールCellrangerで解析された出力ファイルを用いて解析を行います。
2.1 練習用データをダウンロード
10x Genomics社のウェブサイトにデータセットがあるので、今回は5k Peripheral Blood Mononuclear Cells (PBMCs) from a Healthy Donor (v3 chemistry)のデータを使いたいと思います。必要なデータは、Feature / cell matrix HDF5 (filtered)です。hdf5形式と呼ばれるファイルです。まずこのファイルをwgerもしくはcurlでダウンロードします。こちらはR環境下ではうまくいきませんので、q()でRから抜けてください。その後に、
#wget
wget https://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_v3/5k_pbmc_v3_filtered_feature_bc_matrix.h5
#もしくはcurl
curl -O https://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_v3/5k_pbmc_v3_filtered_feature_bc_matrix.h5
これで5k_pbmc_v3_filtered_feature_bc_matrix.h5というファイルがダウンロードされました。再びR環境下で、そのファイルを読み込みます。R環境下に入ったら
# もう一度Seuratを呼び込む
library(Seurat)
# ファイルのパスを指定
file <- "5k_pbmc_v3_filtered_feature_bc_matrix.h5"
2.2 hdf5ファイルを読み込む
Seuratではhdf5形式（〇〇.h5）のファイルを読み込むためのRead10X_h5という関数が用意されています。 それを使用します。Xは大文字です。
# fileは上記のファイルパス
data <- Read10X_h5(file)
ちゃんと読み込まれているかを確認するためには、dataを確認します。以下のようなデータが出力されれば大丈夫です。
head(data)
# 出力されるデータ
MIR1302-2HG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FAM138A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
OR4F5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
AL627309.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
AL627309.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
AL627309.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Seurat objectに変換
前回まではhdf5ファイルを読み込むだけでしたが、この形式からseurat objectと呼ばれる形式に変換します。この形式が何なのかは今のところ知らなくて大丈夫です。Seurat objectに変換するには、CreateSeuratObjectで行います。
# seurat objectに変換
seurat_object <- CreateSeuratObject(data)
# seurat_objectと入力して以下のようなデータが出力されればOK
seurat_object
# 出力データ
An object of class Seurat
33538 features across 5025 samples within 1 assay
Active assay: RNA (33538 features, 0 variable features)
３. クラスタリングの計算
Seurat objectに変換しただけでは、ほぼ生データなので、これから正規化、クラスタリングを行います。
3.1 Normalization
Seuratのチュートリアルでは、正規化の前にTrimmingなどがありますが、一応cellrangerの時点で、不要なデータは取り除かれていることと、本質の解析ではないのでここではTrimmingは省いてNormalization（正規化）に進みます。正規化するには以下のコマンドで行います。
seurat_object <- NormalizeData(seurat_object)
# 以下のようなデータが出力される
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
3.2 特異遺伝子を検出
変化の激しい遺伝子を同定して、今後の解析に使用します。FindVariableFeaturesで行います。nfeatures(特異遺伝子)を2000に設定していますが、この数は変更してもらって構いません。ほとんどの場合2000で大丈夫だと思います。
seurat_object <- FindVariableFeatures(seurat_object, selection.method = "vst", nfeatures = 2000)
# 以下のようなデータが出力される
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
3.3 データのスケール化
遺伝子の発現量は遺伝子によってまちまちなのでスケール化をします。以下のコマンドで行います。まず、遺伝子名を抜き出して、それを特徴量として使います。
all.genes <- rownames(seurat_object)
seurat_object <- ScaleData(seurat_object, features = all.genes)
# 出力
Centering and scaling data matrix
|======================================================================| 100%
3.4 PCA
スケール化したデータを用いてPCA解析を行います。
seurat_object <- RunPCA(seurat_object, features = VariableFeatures(object = seurat_object))
# 以下のような出力
PC_ 1
Positive: LTB, TRBC2, TRAC, IL7R, IL32, CD69, CD7, CD2, TRBC1, CD247
CCR7, TRABD2A, LEF1, SYNE2, ITM2A, RORA, MAL, TRAT1, HIST1H4C, GZMM
AQP3, CD40LG, MT-CO3, CXCR4, MYC, FHIT, LRRN3, DDIT4, ID3, ZNF331
Negative: FCN1, MNDA, LYZ, CST3, SERPINA1, NCF2, CSTA, FGL2, S100A9, CTSS
PSAP, GRN, CD68, CPVL, SPI1, CLEC12A, MPEG1, LST1, TNFAIP2, VCAN
MS4A6A, S100A8, TMEM176B, CYBB, IGSF6, CD14, AIF1, TYROBP, CFP, FCER1G
seurat_object <- FindNeighbors(seurat_object, dims = 1:10)
# 出力
Computing nearest neighbor graph
Computing SNN
seurat_object <- FindClusters(seurat_object, resolution = 0.5)
# 出力
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5025
Number of edges: 165355
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9073
Number of communities: 13
Elapsed time: 0 seconds
3.5 UMAP
最後にUMAPを計算します。
seurat_object <- RunUMAP(seurat_object, dims = 1:10)
# 出力
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
16:39:35 UMAP embedding parameters a = 0.9922 b = 1.112
16:39:35 Read 5025 rows and found 10 numeric columns
16:39:35 Using Annoy for neighbor search, n_neighbors = 30
16:39:35 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:39:36 Writing NN index file to temp file /tmp/Rtmph2s6ZQ/file505fe23b8461a
16:39:36 Searching Annoy index using 1 thread, search_k = 3000
16:39:38 Annoy recall = 100%
16:39:39 Commencing smooth kNN distance calibration using 1 thread
16:39:39 Initializing from normalized Laplacian + noise
16:39:40 Commencing optimization for 500 epochs, with 201484 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:39:47 Optimization finished
ここまでで基本的な計算は終了しました。出来上がったseurat_objectを用いてデータを可視化していきます。
それはまた次回。