Zenn
🙆

[R][Seurat] scRNAseq公開データ読み込み例 ~ h5ファイル ~

2023/05/19に公開・約4,700字
R
#
seurat
#
bioinfomatics
#
scrnaseq
#
公共データ
tech

前置き

scRNAseqデータはfastqファイルが巨大なためかfastqファイルでの公開は減ってきて、1次処理したファイル (cellrangerのoutputや単純な発現マトリクス、解析データも含めたannData、rdsファイルなど)がGEOやZenodo、各国管理のデポジット先に様々な形式で落ちている。
全くもって統一されていないので、まずは読み込むことに時間を喰う。本記事では自身が経験したデータ読み込み例をここに追記していく。

この記事はh5ファイルを主に扱う。その他ファイル形式のものはそれぞれの記事に記録する。


h5ファイル

cellrangerの出力に.h5という拡張子のファイルがあり、こちらを使ってデータを読み込むことも可能。自身でcellranger countを行うとraw_feature_bc_matrix.h5とfiltered_feature_bc_matrix.h5が確認できる。


例) GSE203552

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE203552

GEOからGSE203552_RAW.tarをダウンロードし、解凍してみると以下のようにh5形式のシングルセルデータとその他BCRレパトアのcsvなどが格納されている。ここではscRNAseqのデータのみを読み込む流れを紹介する。


ファイルパスの取得

visiumはvisium.h5、scRNAseqはrna.h5とsuffixを付けてデポジットされているようなので、rna.h5をkeyにしてファイルパス一覧を取得する。

ファイルパスのベクトル作成
(paths <- dir(path = ".", pattern = "rna.h5"))


h5ファイルの読み込み

Read10X_h5機能にファイルパスを指定するだけで読み込める。

1ファイルの読み込み例
mat <- Read10X_h5(paths[1])


Read10X機能を使った読み込みと同様のdgCMatrixが得られるので、以後のSeuratオブジェクトへ変換していく流れは同じである。


データの読み込みとSeuratオブジェクトの作成を一括処理

scRNAseqデータが計9つあり、1つ1つ読み込むのは面倒である。ファイルパスのベクトルを使って、lapply文で繰り返し処理を行う。
今回はこれら9つのデータをmergeしたりintegrateすることを想定して、細胞バーコードにサンプル名を付ける処理を行っている。(細胞バーコードは1サンプル内では重複は無いが、複数サンプルを混ぜて解析する際には重複が生じる。後にややこしいことになる前に、細胞バーコードにサンプル名を付与しているだけである。)

ファイルパスベクトルを使って全ファイル読み込み
sclist <- lapply(paths, function(x){
  # 細胞バーコードに付与するサンプル名の用意。ファイルパスの中からサンプル名の箇所を取得している。
  fname <- stringr::str_split(string = x, pattern = "matrix_", simplify = T)[,2]
  fname <- gsub(pattern = "_rna.h5", replacement = "", x = fname)
  
  # h5ファイルの読み込み
  x <- Read10X_h5(x)
  
  # Seuratオブジェクトの作成
  x <- CreateSeuratObject(counts = x, project = fname)
  # 細胞バーコードにサンプル名を付与
  x <- RenameCells(x, fname)
  return(x)
})


データのmerge

sclistには9つのSeuratオブジェクトが保存されている。これをmerge機能を使って1つのSeuratオブジェクトへまとめる。
mergeは2つのオブジェクトしか引数に取れないので、sclistから1つ取り出して、mergeしていくというのを繰り返している。sc_merge <- merge(sclist[[1]], sclist[2:length(sclist)])と書いてもmergeできる。

データのmerge
sc_merge <- sclist[[1]]
for(i in 2:length(sclist)){
  sc_merge <- merge(x=sc_merge, y=sclist[[i]])
}

合計97399細胞の1つのSeuratオブジェクトが作成された。


meta.dataの付与

GSE203552には論文中で使用したサンプル名や細胞名の対応表も同時にdepositされている。
これをRで読み込んでsc_mergeのmeta.dataスロットへ付与したり、論文で使われた細胞バーコードのみ抜き出して解析に使用する。

depositされていたmeta情報を読み込み
anno_sample <- read.csv("GSE203552_annotation_sample.csv.gz", row.names = 1)
anno_clusterbc <- read.csv("GSE203552_annotation_cluster_bc.csv.gz", row.names = 1)
anno_cluster <- read.csv("GSE203552_annotation_cluster.csv.gz", row.names = 1)

行名に細胞バーコードが入っていた。これらが値も並びも全く同じかをidentical機能を使って確認。

3つのアノテーションファイルは行名が全く同じであることが確認できたので、1つのdata.frameとしてまとめる。

meta情報を1つのdata.frameにまとめる。
anno <- cbind(anno_sample, anno_cluster, anno_clusterbc)
colnames(anno) <- c("anno_sample", "anno_cluster", "anno_clusterbc")

sc_mergeは9万細胞以上あったが、実際に論文に使われたものは73896細胞とのこと。

この情報をSeuratオブジェクトに反映させるには、Seuratオブジェクトに保存されている細胞バーコード名と揃える必要があるが、アノテーションファイルの細胞バーコードはこのままでは使えない。

このようにアノテーションファイルの細胞バーコード名を編集する。

細胞バーコード名を編集
# 文字列中の"rna:"を削除
rownames(anno) <- gsub(x = rownames(anno), pattern = "rna:", replacement = "")
# 文字列末尾の"x"を"-1"に変換
rownames(anno) <- gsub(x = rownames(anno), pattern = "x$", replacement = "-1")

アノテーションファイルの細胞バーコード全てがsc_merge中に存在することが確認できた。

Seuratオブジェクトはindex操作で遺伝子名、細胞バーコード名をkeyにデータを抜き出すことができる。
アノテーションファイルにある細胞バーコードをkeyにして以下のようにSeuratオブジェクトのsubsetが可能である。

sc_merge <- sc_merge[,rownames(anno)]

AddMetaData機能を使って、meta.dataスロットにアノテーション情報を付加

Seuratオブジェクトにmeta情報を追加
sc_merge <- AddMetaData(object = sc_merge, metadata = anno)

このようにmeta.dataスロットに新規列ができている。


Discussion

ログインするとコメントできます