【シングルセル解析】[Seurat] scRNAseq公開データ読み込み例 ~ 発現マトリクスファイル ~
前置き
scRNAseqデータはfastqファイルが巨大なためかfastqファイルでの公開は減ってきて、1次処理したファイル (cellrangerのoutputや単純な発現マトリクス、解析データも含めたannData、rdsファイルなど)がGEOやZenodo、各国管理のデポジット先に様々な形式で落ちている。
全くもって統一されていないので、まずは読み込むことに時間を喰う。本記事では自身が経験したデータ読み込み例をここに追記していく。
この記事は発現マトリクスファイルを主に扱う。その他ファイル形式のものはそれぞれの記事に記録する。
Matrixファイル
行に遺伝子、列に細胞名の単純なmatrix/data.frameからでもSeuratオブジェクトへ変換できる。
ただ、PCメモリに載りきらないサイズの行列をRに読み込ませる場合はエラーがでる。(Read10X
機能で読み込まれるdgCMatrixオブジェクトは、カウントが0の箇所を.にするなど大きなmatrixでもメモリを節約して読み込んでいる。)
例) GSE155622
行列読み込み
この例では行に遺伝子、列に細胞名が入ったmatrixがtxtファイルとしてdepositされていた。
txtファイルであればread.delim()
で以下のように読み込める。(たまにtxt拡張子で半角スペース区切りの場合もあるので、その場合はsep=" "
引数を付ける。)
mat <- read.delim(file = "GSE155622_raw_UMI_counts_1.txt.gz", row.names=1)
しかし、巨大なmatrixをローカルから読み込むのは時間がかかる。上記では10分もかかっている。
data.tableパッケージのfread()
機能を使えば速度は改善する。1分強で読み込むことができた。
mat2 <- data.table::fread(file = "GSE155622_raw_UMI_counts_1.txt.gz")
# freadには行名を指定して読み込むオプションが無いようなので、1列目が行名となるように変換。
mat2 <- data.frame(mat2, row.names = 1)
dgCMatrixへ変換
ちなみに以下のようにMatrixパッケージを使用してdgCMatrixに変換することも可能。メモリ上のサイズを節約できる。
mat3 <- Matrix::Matrix(as.matrix(mat), sparse = TRUE)
オブジェクトのサイズ
Seuratオブジェクトへ変換
読み込んだ行列をCreateSeuratObject()
機能でSeuratオブジェクトへ変換すればよい。
sc <- CreateSeuratObject(mat)
アノテーション情報の追加
このdepositデータには細胞名とmeta.dataの対応表もdepositされていた。
meta <- read.delim("GSE155622_raw_UMI_counts_1_metadata.txt.gz", row.names=1)
各細胞がどの細胞種かなどの情報もあり有用である。matrixの細胞数とアノテーション情報の細胞数も揃っているので、このオブジェクトをそのままSeuratオブジェクトに反映させると良さそう。
35894細胞
しかし、一部の細胞名のみが一致していて、大半は名前が合っていないようである。
9103細胞は名前が一致、26791細胞は名前が不一致
アノテーション情報の細胞名には半角スペースを持つものがあるが、matrixの細胞名は.
に変換されている。
データ処理を進める上で半角スペースは扱いづらいので、アノテーション情報の細胞名を変換しておく。これでSeuratオブジェクトでの細胞名とアノテーション情報の細胞名が一致する。
rownames(meta) <- gsub(pattern = " ", replacement = ".", x = rownames(meta))
Seuratオブジェクトにアノテーション情報を登録するには
-
Seuratオブジェクト作成時に
meta.data
引数にアノテーション情報を入れる
sc <- CreateSeuratObject(counts = mat, meta.data = meta)
-
AddMetaData
機能で作成済みのSeuratオブジェクトに追記する
sc <- CreateSeuratObject(counts = mat)
sc <- AddMetaData(sc, metadata = meta)
例) GSE139891
デポジットされているGSE139891_RAW.tarをダウンロード、解凍すると以下のようなファイルが含まれていた。txtファイルをgz圧縮している。
原著を見ると、6サンプルのscRNAseqとそのうち2サンプルはCITEseqを行っているとのこと。scRNAseqはumiという文字を含む6つのファイルが相当しそうだ。
以下はワーキングディレクトリをtxt.gzファイルがあるディレクトリに変更したうえで行う想定。
まずは1ファイル読み込んで形式を確認する。
mat <- data.table::fread(file = "GSM4148370_GC1_umi_grch38.txt.gz")
mat <- data.frame(mat, row.names = 1)
行に遺伝子、列に細胞がある形式は良いものの、遺伝子名が"Ensembl ID;遺伝子名"になっていて見づらい。そこで遺伝子名だけになるように試みるも重複が生じるため行名に持ってこれない。
重複している遺伝子は意味の無さそうなものばかりだが、念のため重複遺伝子の中でも最もカウント値が高いものを残すようにしておく。具体的には遺伝子ごとの平均や合計を求めて、降順に並び替えたうえで重複除去を行えばよい。
mat <- mat[ order(rowSums(mat), decreasing = TRUE), ]
mat <- mat[!duplicated(stringr::str_split(string = rownames(mat), pattern = ";", simplify = T)[,2]),]
rownames(mat) <- stringr::str_split(string = rownames(mat), pattern = ";", simplify = T)[,2]
これで遺伝子名を行名に持ってくることができた。
最終スクリプト
6サンプル分を読み込んだ最終版。
# umiという文字を含むファイル名一覧を取得
filepath <- list.files(pattern = "umi")
# txt.gzファイルを読み込み、行名処理したリストを取得
sclist <- lapply(filepath, function(x){
x <- data.table::fread(file = x)
x <- data.frame(x, row.names = 1)
x <- x[ order(rowSums(x), decreasing = TRUE), ]
x <- x[ !duplicated(stringr::str_split(string = rownames(x), pattern = ";", simplify = T)[,2]) , ]
rownames(x) <- stringr::str_split(string = rownames(x), pattern = ";", simplify = T)[,2]
return(x)
})
# リストにサンプル名をつける。
names(sclist) <- stringr::str_split(string = filepath, pattern = "_", simplify = T)[,3]
library(Seurat)
sclist <- mapply(function(x,y){
x <- CreateSeuratObject(counts = x, project = y, min.cells = 3, min.features = 200)
x <- RenameCells(x, add.cell.id=y)
return(x)
},
x = sclist, y= names(sclist), SIMPLIFY = FALSE
)
6検体をmergeするなら。
sc_merge <- merge(x = sclist[[1]], y = sclist[2:length(sclist)])
Discussion