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を用いてデータを可視化していきます。

それはまた次回。