😺

【シングルセル解析】[Seurat] scRNAseq公開データ読み込み例 ~ 発現マトリクスファイル ~

2023/05/19に公開

前置き

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

この記事は発現マトリクスファイルを主に扱う。その他ファイル形式のものはそれぞれの記事に記録する。


Matrixファイル

行に遺伝子、列に細胞名の単純なmatrix/data.frameからでもSeuratオブジェクトへ変換できる。
ただ、PCメモリに載りきらないサイズの行列をRに読み込ませる場合はエラーがでる。(Read10X機能で読み込まれるdgCMatrixオブジェクトは、カウントが0の箇所を.にするなど大きなmatrixでもメモリを節約して読み込んでいる。)


例) GSE155622

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

行列読み込み

この例では行に遺伝子、列に細胞名が入ったmatrixがtxtファイルとしてdepositされていた。
txtファイルであればread.delim()で以下のように読み込める。(たまにtxt拡張子で半角スペース区切りの場合もあるので、その場合はsep=" "引数を付ける。)

read.delimでtxtファイルを読み込み
mat <- read.delim(file = "GSE155622_raw_UMI_counts_1.txt.gz", row.names=1)

しかし、巨大なmatrixをローカルから読み込むのは時間がかかる。上記では10分もかかっている。

data.tableパッケージfread()機能を使えば速度は改善する。1分強で読み込むことができた。

freadでtxtファイルを読み込み
mat2 <- data.table::fread(file = "GSE155622_raw_UMI_counts_1.txt.gz")
# freadには行名を指定して読み込むオプションが無いようなので、1列目が行名となるように変換。
mat2 <- data.frame(mat2, row.names = 1)

dgCMatrixへ変換

ちなみに以下のようにMatrixパッケージを使用してdgCMatrixに変換することも可能。メモリ上のサイズを節約できる。

dgCMatrixへ変換
mat3 <- Matrix::Matrix(as.matrix(mat), sparse = TRUE)


オブジェクトのサイズ


Seuratオブジェクトへ変換

読み込んだ行列をCreateSeuratObject()機能でSeuratオブジェクトへ変換すればよい。

発現マトリクスからSeuratオブジェクトを作成
sc <- CreateSeuratObject(mat)



アノテーション情報の追加

このdepositデータには細胞名とmeta.dataの対応表もdepositされていた。

depositされていたmeta情報を取得
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オブジェクトにアノテーション情報を登録するには

  1. Seuratオブジェクト作成時にmeta.data引数にアノテーション情報を入れる
meta.data引数を使用してSeuratオブジェクトを作成
sc <- CreateSeuratObject(counts = mat, meta.data = meta)
  1. AddMetaData機能で作成済みのSeuratオブジェクトに追記する
AddMetaData機能でmeta.dataスロットに追記
sc <- CreateSeuratObject(counts = mat)
sc <- AddMetaData(sc, metadata = meta)



例) GSE139891

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=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;遺伝子名"になっていて見づらい。そこで遺伝子名だけになるように試みるも重複が生じるため行名に持ってこれない。



重複している遺伝子は意味の無さそうなものばかりだが、念のため重複遺伝子の中でも最もカウント値が高いものを残すようにしておく。具体的には遺伝子ごとの平均や合計を求めて、降順に並び替えたうえで重複除去を行えばよい。

発現マトリクスを並び替えてからduplicate遺伝子を除去
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]

Seuratオブジェクトへ変換
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するなら。

merge
sc_merge <- merge(x = sclist[[1]], y = sclist[2:length(sclist)])

Discussion