シングルセルシーケンス解析ツールSeuratの使い方
1. Seuratのインストール
まず、RにSeuratをインストールしないといけません。Seuratの他にHdf5ファイルを読み込むためのhdf5rもインストールしましょう。インストールのコマンドは以下です。
Biocondaでインストール（おすすめ）
何かをインストールするとき、Biocondaでインストールできればそれが一番簡単です。ターミナルの画面で、R環境に 入ってない 状態で、以下のコマンドを入力します。
R環境下でインストール
Biocondaじゃなくても、Rのinstall.packagesコマンドでSeuratをインストール可能です。ターミナルで”R”と入力すると、Rの環境に入るはずです。R環境下で以下のコマンドを入力します。
上記のいずれかで、Seuratがちゃんとインストールされたかを確認するにはR環境下で、以下のコマンドを入力します。
このコマンドを入力して何もエラーが出なければちゃんとSeuratがインストールされ、呼び出されています。
いよいよ本番の解析です。現在、singleセルシーケンスの多くは10x Genomics社のChromiumプラットフォームを用いて行われています。なので、今回は10x Genomics社の解析ツールCellrangerで解析された出力ファイルを用いて解析を行います。
練習用データをダウンロード
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から抜けてください。その後に、
これで5k_pbmc_v3_filtered_feature_bc_matrix.h5というファイルがダウンロードされました。再びR環境下で、そのファイルを読み込みます。R環境下に入ったら
hdf5ファイルを読み込む
Seuratではhdf5形式（〇〇.h5）のファイルを読み込むための Read10X_h5という 関数が用意されています。 それを使用します。Xは大文字です。
ちゃんと読み込まれているかを確認するためには、dataを確認します。以下のようなデータが出力されれば大丈夫です。
Seurat objectに変換
前回まではhdf5ファイルを読み込むだけでしたが、この形式からseurat objectと呼ばれる形式に変換します。この形式が何なのかは今のところ知らなくて大丈夫です。Seurat objectに変換するには、CreateSeuratObjectで行います。
Normalization
Seuratのチュートリアルでは、正規化の前にTrimmingなどがありますが、一応cellrangerの時点で、不要なデータは取り除かれていることと、本質の解析ではないのでここではTrimmingは省いてNormalization（正規化）に進みます。正規化するには以下のコマンドで行います。
変化の激しい遺伝子を同定して、今後の解析に使用します。FindVariableFeaturesで行います。nfeatures(特異遺伝子)を2000に設定していますが、この数は変更してもらって構いません。ほとんどの場合2000で大丈夫だと思います。
遺伝子の発現量は遺伝子によってまちまちなのでスケール化をします。以下のコマンドで行います。まず、遺伝子名を抜き出して、それを特徴量として使います。
all.genes <- rownames(seurat_object) seurat_object <- ScaleData(seurat_object, features = all.genes)
スケール化したデータを用いてPCA解析を行います。
1. Dimentional Plot
1.1 基本的なDimentional Plotの表示
DimPlot(seurat_clusters)
1.2 Legendを消したいとき
DimPlot(seurat_clusters) + theme(legend.position="none")
#もしくは
DimPlot(seurat_clusters) + NoLegend()
1.3 Labelを表示
DimPlot(seurat_clusters, label=TRUE)
#labelの大きさを設定する
DimPlot(seurat_clusters, label=TRUE, label.size=10)
1.4 クラスターのグループを変える
まず、どのグループが使えるかは、metadataを表示すればわかります。
seurat_clusters@meta.data
DimPlot(seurat_clusters, group.by="RNA_snn_res.0.5")
DimPlot(seurat_clusters, group.by="RNA_snn_res.0.5")