とりあえずSeuratというRのパッケージで解析すると良いらしい。
ちなみに10xのCell Rangerでマッピングまではやっている前提で。
こちらの情報を参考にやってみることにする。ちょっと古いけど大丈夫かな？
Seuratを駆使する会 ① - ばいばいバイオ
Rはほんとうに敬遠していて何もわからないのだけど。
まずはパッケージをインストールから
わからないなりに環境だけいっちょ前にLinux上に構築してRStudioサーバを建ててあり、ネット越しにブラウザでアクセスして解析ができるようにしていたりする。なにもわからないのに。
Rのバージョンは4.4.0である。
> install.packages('Seuart')とするといつ終わるかもわからないくらいずらずらとメッセージが流れる。
幸い、致命的エラーは出ずにインストールできたらしい。
> library(Seuart)
ライブラリインポートも問題なくできたっぽい。
さて、チュートリアルを一通りやっておこう。
ConsoleタブからTerminalタブに切り替えて、
$ mkdir scrnaseqtest (ngs) kuro@guilty6-x2110:~$ cd scrnaseqtest/ (ngs) kuro@guilty6-x2110:~/scrnaseqtest$ mkdir immune_alignment_matrix (ngs) kuro@guilty6-x2110:~/scrnaseqtest$ cd immune_alignment_matrix/ (ngs) kuro@guilty6-x2110:~/scrnaseqtest/immune_alignment_matrix$ wget https://www.dropbox.com/s/79q6dttg8yl20zg/immune_alignment_expression_matrices.zip `immune_alignment_expression_matrices.zip' へ保存完了 [21329741/21329741] (ngs) kuro@guilty6-x2110:~/scrnaseqtest/immune_alignment_matrix$ unzip immune_alignment_expression_matrices.zip Archive: immune_alignment_expression_matrices.zip inflating: immune_control_expression_matrix.txt.gz creating: __MACOSX/ inflating: __MACOSX/._immune_control_expression_matrix.txt.gz inflating: immune_stimulated_expression_matrix.txt.gz inflating: __MACOSX/._immune_stimulated_expression_matrix.txt.gz (ngs) kuro@guilty6-x2110:~/scrnaseqtest/immune_alignment_matrix$ rm immune_alignment_expression_matrices.zip (ngs) kuro@guilty6-x2110:~/scrnaseqtest/immune_alignment_matrix$ gunzip immune_control_expression_matrix.txt.gz (ngs) kuro@guilty6-x2110:~/scrnaseqtest/immune_alignment_matrix$ gunzip immune_stimulated_expression_matrix.txt.gz
と、これで下準備OK
さてConsoleに戻って
> ctrl.data <- read.table("immune_control_expression_matrix.txt", sep = '\t', row.names = 1)
file(file, "rt") でエラー: コネクションを開くことができません
追加情報: 警告メッセージ:
file(file, "rt") で:
ファイル 'immune_control_expression_matrix.txt' を開くことができません: そのようなファイルやディレクトリはありませんああ、そうかディレクトリを移動しないといけないのか。どうやって？（このレベルです）
> setwd("~/scrnaseqtest/immune_alignment_matrix/")
> ctrl.data <- read.table("immune_control_expression_matrix.txt", sep = '\t', row.names = 1)
> stim.data <- read.table("immune_stimulated_expression_matrix.txt", sep = "\t", row.names = 1)
> ctrl <- CreateSeuratObject(raw.data = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)
CreateSeuratObject(raw.data = ctrl.data, project = "IMMUNE_CTRL", でエラー:
引数 "counts" がありませんし、省略時既定値もありませんまたひっかかる。どうもバージョン3のSeuratでなにか変わっているらしい。なのでチャットGPTの出番。
エラーの原因は、CreateSeuratObject 関数の引数の名前が変更されたためです。Seuratのバージョン3以降では、raw.data 引数が counts 引数に変更されました。
> ctrl <- CreateSeuratObject(counts = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)
警告: Feature names cannot have underscores ('_'), replacing with dashes ('-')
警告: Data is of class data.frame. Coercing to dgCMatrix.なんか警告は出たけどいいらしい。
続いてmetadataの付与
> ctrl@meta.data$stim <- "CTRL" > head(ctrl@meta.data$stim) [1] "CTRL" "CTRL" "CTRL" "CTRL" "CTRL" "CTRL"
フィルタリング
> ctrl <- FilterCells(ctrl, subset.names = "nGene", low.thresholds = 500, high.thresholds = Inf) FilterCells(ctrl, subset.names = "nGene", low.thresholds = 500, でエラー: 関数 "FilterCells" を見つけることができませんでした
chatGPT
FilterCells 関数は、Seuratのバージョン3以降では廃止されました。代わりに、subset 関数を使用します。
らしい。
> ctrl <- subset(ctrl, subset = nFeature_RNA > 500)
正規化
> ctrl <- NormalizeData(ctrl, scale.factor = 10000) Normalizing layer: counts Performing log-normalization 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************|
z-scoreの算出
> ctrl <- ScaleData(ctrl) Centering and scaling data matrix |=================================================================================================================| 100%
同様にstimulatedのほうも
> stim <- CreateSeuratObject(counts = stim.data, project = "IMMUNE_STIM", min.cells = 5)
警告: Feature names cannot have underscores ('_'), replacing with dashes ('-')
警告: Data is of class data.frame. Coercing to dgCMatrix.
> stim@meta.data$stim <- "STIM"
> stim <- subset(stim, subset = nFeature_RNA > 500)
> stim <- NormalizeData(stim)
Normalizing layer: counts
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> stim <- ScaleData(stim)
Centering and scaling data matrix
|=================================================================================================================| 100%分散の大きい（発現変動の大きい）遺伝子を抽出
> ctrl <- FindVariableGenes(ctrl) FindVariableGenes(ctrl) でエラー: 関数 "FindVariableGenes" を見つけることができませんでした > ctrl <- FindVariableFeatures(ctrl) Finding variable features for layer counts 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% [----|----|----|----|----|----|----|----|----|----| **************************************************|
FindVariableGenes 関数も、Seuratのバージョン3以降では廃止されました。代わりに FindVariableFeatures 関数を使用します。
> g.1 <- head(rownames(ctrl@hvg.info), 1000)
dimnames(x)[[1L]] %||% if (do.NULL) NULL else { でエラー:
名前 "hvg.info" というスロットが、クラス "Seurat" のこのオブジェクトには存在しません
> # 上位1000の変動遺伝子を取得
> g.1 <- head(VariableFeatures(ctrl), 1000)
> g.2 <- head(VariableFeatures(stim), 1000)hvg.info スロットも Seurat のバージョン3以降では変更されています。変動遺伝子の情報は VariableFeatures スロットに格納されています。
g1、g2で共通して発現変動の大きい遺伝子をunique関数で抽出
> genes.use <- unique(c(g.1, g.2))
> genes.use <- intersect(genes.use, rownames(ctrl@scale.data))
dimnames(x)[[1L]] %||% if (do.NULL) NULL else { でエラー:
名前 "scale.data" というスロットが、クラス "Seurat" のこのオブジェクトには存在しませんscale.data スロットはSeuratのバージョン3以降では存在しません。代わりに ScaleData 関数を使用してデータをスケーリングし、ctrl"RNA"@scale.data からスケールされたデータにアクセスできます。
だそうで
> genes.use <- intersect(genes.use, rownames(ctrl[["RNA"]]@scale.data))
dimnames(x)[[1L]] %||% if (do.NULL) NULL else { でエラー:
名前 "scale.data" というスロットが、クラス "Assay5" のこのオブジェクトには存在しませんまだだめ。
scale.data スロットが存在しないというエラーは、Seuratのバージョンによる変更の結果です。現在のSeuratバージョンでは、スケールされたデータは scale.data スロットではなく、ScaleData 関数が返すオブジェクト内に格納されます。
> scaled_genes <- rownames(ctrl[["RNA"]]@data)
dimnames(x)[[1L]] %||% if (do.NULL) NULL else { でエラー:
名前 "data" というスロットが、クラス "Assay5" のこのオブジェクトには存在しませんSeuratの最新バージョンでは、スケールされたデータは ScaleData 関数の出力として直接取得できないため、正しいスロットにアクセスする必要があります。スケーリングされたデータは通常、ctrl"RNA"@scale.data ではなく、ctrl"RNA"@scale.data に格納されます。
しかし、最新のバージョンではスケールデータが標準的なアクセス方法ではなく、GetAssayData 関数を使う必要があります。
> scaled_genes <- rownames(GetAssayData(ctrl, slot = "scale.data")) 警告メッセージ: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0. ℹ Please use the `layer` argument instead. This warning is displayed once every 8 hours. Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated. > genes.use <- intersect(genes.use, scaled_genes)
なんかあやしくなってきたぞ？
一旦停止して整理する。
