🐥

[R][Seurat] scRNAseq公開データ読み込み例

2023/02/22に公開・約18,800字

seurat

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

本記事では自身が経験したデータ読み込み例をここに追記していく。

cellranger countの3つの出力ファイル
h5ファイル
Matrixファイル
h5adファイル

cellranger countの3つの出力ファイル

cellranger countでfastqファイルの処理をすると、filtered_feature_bc_matrixフォルダに以下の3つのファイルが生成される。

barcodes.tsv.gz
features.tsv.gz
matrix.mtx.gz

SeuratのRead10X機能ではこの3つのファイルが入ったフォルダのパスを指定する。

GEOなどの公共データからダウンロードした場合、prefixにサンプルごとの識別名がつくことがあるが、面倒なことにRead10Xはファイル名がcellrangerのデフォルト出力と一字一句同じでなければ認識しない。

例) GSE181304

1サンプル分のcellranger outputがdepositされていた。3つのファイルのprefixにGSM番号が付与されている。

Read10Xでの読み込み

このままではファイルがmissingであるとエラーがでる。

簡単な対策としてはファイル名を編集してGSM番号の箇所を消せば良い。この例ではprefixが_で繋がっていたので、_で文字列を分割して、2番目の要素を使用した。

Rからファイル名を変更

fromname <- dir(path = ".", pattern = ".gz")
toname <- stringr::str_split(string = fromname, pattern = "_", simplify = T)[,2]
file.rename(from = fromname, toname)

mat <- Read10X(data.dir = getwd())

ReadMtx機能で1ファイルずつ指定して読み込み

この機能であればファイル名のルールに縛られなくてよい。

ReadMtx機能で3つのファイルを明示的に指定して読み込み

mat2 <- ReadMtx(mtx = "GSM5494566_matrix.mtx.gz",
                cells = "GSM5494566_barcodes.tsv.gz",
                features = "GSM5494566_features.tsv.gz")

ReadMtxで3つのファイルを指定する引数以外のオプションはデフォルトのままでRead10Xと同じ出力が得られた。

例) heiData, VRJUNV

https://doi.org/10.11588/data/VRJUNV

cellranger countの3つの出力ファイルが入ったzipフォルダがサンプル数 (12個)分depositされている。
Read10Xはzipファイルを読めなかったので、zipは解凍しておく。

depositされていたsample_sheetファイルを見ると、cellrangerのversionが2.1.1で処理されたものと3.0.2で処理されたものが混ざっている。

cellrangerはversionによってファイル名が変わっており、現行のfeatures.tsvファイルが古いversionではgenes.tsvファイルとなっている。

DLBCL1フォルダの中身

DLBCL3フォルダの中身

またcellranger countの3つの出力ファイルはgz圧縮された状態のはずだが、このdeposit dataでは解凍されている。

実際にRead10Xでサンプルフォルダを引数に読み込んでみると、何故かcellranger ver 2.1.1の出力はgz無しでもRead10Xで読めたのだが、ver 3.0.2の出力はgz無しだとRead10Xでは読めなかった。

ver2.1.1のDLBCL1の読み込みとver3.0.2のDLBCL3の読み込み画面

しかし、ReadMtx機能ではどちらのcellranger versionでも読むことができた。

最終コマンド

最終的に以下の流れで読み込み成功。(zipファイルから開始する想定)

data_dir <- "ダウンロードした12検体分のデータが入ったフォルダへのパス"

# 各サンプルのパスを取得
samplepaths <- dir(path = data_dir, 
                   pattern = "zip",
                   full.names = T)

# 各サンプルの識別名を作成
samplenames <- dir(path = data_dir, 
                    pattern = "zip")
(samplenames <- gsub(pattern = ".zip", replacement = "", x = samplenames))

# 各サンプルをfor loopで読み込んで保存する先の空リストを用意
sclist <- vector(mode = "list", length = length(samplenames))
names(sclist) <- samplenames

# meta情報もdepositされていたのでmeta.dataスロットへ反映させるために使用する。
meta <- read.csv(file = paste(data_dir, "sample_sheet.csv",sep = "/"))

# for loopで以下を全サンプルをSeuratオブジェクトに変換する。
for (i in 1:length(samplenames)) {
  # zipを解凍
  unzip(samplepaths[i], overwrite = T)
  
  # ReadMtxの引数を用意。(3つのファイルへのパス)
  mtx <- paste(samplenames[i], "matrix.mtx", sep = "/")
  cells <- paste(samplenames[i], "barcodes.tsv", sep = "/")
  features <- paste(samplenames[i], "genes.tsv", sep = "/")
  
  # cellranger versionの違いでfeatures.tsvかgenes.tsvか変わるので、if文で対応
  if (file.exists(features) == FALSE){
    features <- paste(samplenames[i], "features.tsv", sep = "/") 
  }
  
  # 発現マトリクスを読み込み
  x <- ReadMtx(mtx = mtx, cells = cells, features = features)
  
  # zip解凍したフォルダを削除
  unlink(samplenames[i], recursive = TRUE)
  
  # Seuratオブジェクトに変換
  x <- CreateSeuratObject(x, project = samplenames[i])
  # meta.dataスロットにメタ情報を記録
  x <- AddMetaData(object = x, metadata = rep(x = meta[meta$SampleName == samplenames[i],], times = ncol(x)))
  
  # リストのi番目の要素として保存
  sclist[[i]] <- x
}

h5ファイル

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

例) 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)
  # 細胞バーコードにサンプル名を付与
  x <- RenameCells(x, fname)
  return(x)
})

データのmerge

sclistには9つのSeuratオブジェクトが保存されている。これをmerge機能を使って1つのSeuratオブジェクトへまとめる。
mergeは2つのオブジェクトしか引数に取れないので、sclistから1つ取り出して、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スロットに新規列ができている。

Matrixファイル

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

例) GSE155622

この例では行に遺伝子、列に細胞名が入った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をローカルから読み込むのは時間がかかる。
data.tableパッケージのfread機能を使えば速度は改善する。

freadでtxtファイルを読み込み

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

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

dgCMatrixへ変換

mat3 <- Matrix::Matrix(as.matrix(mat), sparse = TRUE)

オブジェクトのサイズ

発現マトリクスから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オブジェクトに反映させると良さそう。

しかし、一部の細胞名は一致していて、大半は名前が合っていないようである。

アノテーション情報の細胞名には半角スペースを持つものがあるが、matrixの細胞名は.に変換されている。

データ処理を進める上で半角スペースは扱いづらいので、アノテーション情報の細胞名を変換しておく。

細胞名を編集

rownames(meta) <- gsub(pattern = " ", replacement = ".", x = rownames(meta))

Seuratオブジェクト作成時に、meta.data引数にアノテーション情報を入れることもできるし、AddMetaData機能で作成済みのSeuratオブジェクトに追記することもできる。

meta.data引数を使用してSeuratオブジェクトを作成

sc <- CreateSeuratObject(counts = mat, meta.data = meta)

AddMetaData機能でmeta.dataスロットに追記

sc <- CreateSeuratObject(counts = mat)
sc <- AddMetaData(sc, metadata = meta)

h5adファイル

Pyhton Scanpyで解析されたAnnDataクラスのデータをR Seuratで読み込む。
SeuratDiskパッケージが必要となる。まだ入れてなければ以下のようにインストールする。

SeuratDiskのインストール

if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
remotes::install_github("mojaveazure/seurat-disk")

例) PRJCA001063

この例ではZenodoにアノテーション情報つきの解析データがh5ad拡張子でdepositされている。

読み込むにはSeuratDiskパッケージのConver機能とLoadH5Seurat機能を使用する。

# h5adファイルからh5seuratファイルを作成。(ローカルのフォルダに作成される。)
SeuratDisk::Convert(source = "StdWf1_PRJCA001063_CRC_besca2.annotated.h5ad", 
		    dest="h5seurat",
		    overwrite=TRUE)
# h5seurat拡張子のものを読み込み		    
sc <- SeuratDisk::LoadH5Seurat(file = "StdWf1_PRJCA001063_CRC_besca2.annotated.h5seurat")

Seuratオブジェクトとして読み込まれる。meta.dataも再現。

さらに次元削減の情報なども残っていたので、そのまま論文通りのplotが再現できた。

DimPlot(sc, group.by = "Cell_type")

例) Human Cell Atlas, AdultHemOrgans

Human Cell Atlasからは複数のデータ形式でのダウンロードが選択可能であるが、今回は既に解析済みのデータが入っていることを期待してh5adファイルを使用してみた。

h5adのダウンロードボタンをクリックし、移動先のページで得られるcurlコマンドをコピペしてダウンロード。(私はwindowsでWSL ubuntu環境を構築しているが、curlを使える環境作りは各自OSで調べてほしい。)

このデータも以下のコマンドでSeuratオブジェクトとして読み込むことができた。

SeuratDisk::Convert(source = "20210204_S_SUBS4_4x_annotated_donors_UMI_counts_HCA.h5ad", 
		    dest="h5seurat",
		    overwrite=TRUE)
sc <- SeuratDisk::LoadH5Seurat(file = "20210204_S_SUBS4_4x_annotated_donors_UMI_counts_HCA.h5seurat")

患者番号や由来組織などのメタ情報が既に入った状態なのはありがたい。

Seuratオブジェクトのmeta.dataスロットを確認

str(sc@meta.data)

'data.frame': 66003 obs. of 4 variables:
$ batch : Factor w/ 4 levels "0","1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ library: Factor w/ 11 levels "SIGAA12","SIGAB10",..: 5 5 5 5 5 5 5 5 5 5 ...
$ donor : Factor w/ 4 levels "DOD1","DOD2",..: 1 1 1 1 1 1 1 1 1 1 ...
$ organ : Factor w/ 3 levels "BM","PB","SPL": 2 2 2 2 2 2 2 2 2 2 ...

ただ、慣例のごとくデータのフィルタリングを行おうとすると、nCount_RNAやnFeature_RNAが無いことに気づいた。

sc <- PercentageFeatureSet(sc, pattern = "^MT-", col.name = "Percent_mito")
VlnPlot(sc, features = c("nFeature_RNA", "nCount_RNA", "Percent_mito"), pt.size = 0)

Error: Cannot find 'nCount_RNA' in this Seurat object

そこで、以下のように自身でnCount_RNAとnFeature_RNAの情報を追加。

# 発現マトリクスを取得
mat <- GetAssayData(sc, slot = "counts", assay = "RNA")

# 細胞ごとの合計発現カウントをnCount_RNAとして保存
sc$nCount_RNA <- colSums(mat)

# 細胞ごとに発現が0より大きい遺伝子の数を集計してnFeature_RNAに保存
sc$nFeature_RNA <- colSums(mat>0)

無事にいつものデータ分布の確認ができた。

sc <- PercentageFeatureSet(sc, pattern = "^MT-", col.name = "Percent_mito")
VlnPlot(sc, features = c("nFeature_RNA", "nCount_RNA", "Percent_mito"), pt.size = 0)

例) Human Cell Atlas, HematopoieticImmuneCellAtlas

このdepositデータも幾つかのデータ形式 (fastq, bam, loom, h5, h5ad)が選べるが、h5ad以外はファイル数が100を超えており、これをダウンロードして処理すのは気が遠くなる。
色々と処理済みのデータが入っていることを願ってh5adをダウンロードした。

以下のように取得したh5adファイルを読み込もうとすると、levelsとvaluesの情報が無いというエラーが出た。

SeuratDisk::Convert(source = "BL_hashing.h5ad", 
	    dest="h5seurat",
	    overwrite=TRUE)
sc <- SeuratDisk::LoadH5Seurat(file = "BL_hashing.h5seurat")

Error: Missing required datasets 'levels' and 'values'

検索してみると、scanpyとanndataのversion更新により、新しいh5adファイルでは同様のエラーが出ているようだ。

対策として、LoadH5Seuratのmeta.data引数とmisc引数をFALSEにするとひとまずエラー無く読み込めるとのことである。

seuratobj <- LoadH5Seurat("BL_hashing.h5seurat", meta.data = F, misc = F)

確かにSeuratオブジェクトを作ることはできたが、meta.dataスロットは空である。一方reductionスロットは読み込めている。

anndataパッケージを使うと、h5ad拡張子をRで扱うことができる。

Rでanndataを使うには

おそらくanndataを使うにはpythonのanndataも必要となる。

python環境が作られていなければ、anndataの機能を使おうとした際にMinicondaを入れるか多分聞いてくるのでそれに従うとよい。
何もエラーが無ければ問題なし。

既にAnaconda3などを入れておりパスが通っていればその環境を使おうとする。

ただ既にpythonでの解析用に環境構築している場合、Rから同じ環境を使いたくないかもしれない。

私の場合はreticulateパッケージの紹介のようにr-reticulate環境を作成して、そこにanndataをインストールした。

reticulateパッケージで仮想環境を作成

library(reticulate)

# r-reticulate環境の構築
conda_create("r-reticulate")
conda_list() # 環境が作られたか確認

# r-reticulate環境に切り替え
use_condaenv("r-reticulate") # もしfailed to initialize requested version of Pythonとでれば、Rを立ち上げ直すとよい

# r-reticulate環境にanndataパッケージをインストール
conda_install(envname = "r-reticulate",packages = "anndata")

もし書き込み権限でインストールできないようであれば、R/Rstudioを管理者権限で立ち上げてからやってみてほしい。

read_h5ad機能を使って、depositされていたデータ構造を見てみる。

# install.packages("anndata")
library(anndata)

h5ad <- read_h5ad("BL_hashing.h5ad")
h5ad

h5ad
AnnData object with n_obs × n_vars = 9278 × 19651
obs: 'assignment', 'n_genes', 'n_counts', 'percent_mito', 'scale', 'leiden_labels', 'doublet_score', 'pred_dbl', 'anno'
var: 'featureid', 'n_cells', 'percent_cells', 'robust', 'highly_variable_features', 'mean', 'var', 'hvf_loess', 'hvf_rank'
uns: 'PCs', '_attr2type', 'genome', 'leiden_resolution', 'modality', 'norm_count', 'pca', 'pca_features', 'stdzn_max_value', 'stdzn_mean', 'stdzn_std'
obsm: 'X_pca', 'X_umap', 'pca_knn_distances', 'pca_knn_indices'
varm: 'de_res'
obsp: 'W_pca'

obsの箇所にmeta.dataが入っていて、何としても読み込みたい。
obsは以下のように$でアクセスするだけでmeta.dataスロットに合った形状のデータが得られた。

またXの箇所には発現マトリクスが入っていた。これをCreateSeuratObjectの引数に使えばSeuratオブジェクトに変換できそうである。

str(h5ad$X)

str(h5ad$X)
Formal class 'dgRMatrix' [package "Matrix"] with 6 slots
..@ p : int [1:9279] 0 596 1315 1919 2880 3929 5117 5958 6828 7539 ...
..@ j : int [1:8865149] 56 251 315 331 348 355 357 391 412 446 ...
..@ Dim : int [1:2] 9278 19651
..@ Dimnames:List of 2
.. ..$ : chr [1:9278] "AAACCTGAGACGCACA" "AAACCTGAGGAATGGA" "AAACCTGAGTAATCCC" "AAACCTGAGTCGTTTG" ...
.. .. ..- attr(, "name")= chr "barcodekey"
.. ..$ : chr [1:19651] "MIR1302-2HG" "AL627309.1" "AL669831.2" "AL669831.5" ...
.. .. ..- attr(, "name")= chr "featurekey"
..@ x : num [1:8865149] 4.51 4.51 4.51 4.51 4.51 ...
..@ factors : list()

ただ行に細胞バーコード、列に遺伝子となっており、転置する必要がある。

data.frame(h5ad$X)[1:5,1:5]

通常のt関数は使えず、Matrixパッケージのt関数を使うと転置が可能であった。

data.frame(Matrix::t(h5ad$X))[1:5,1:5]

最終コマンド

最終的に以下の流れで読み込むことができた。

library(Seurat)
library(SeuratDisk)
library(anndata)

# Convertでh5adからh5seuratファイルを作成
SeuratDisk::Convert(source = "BL_hashing.h5ad", 
		    dest="h5seurat",
		    overwrite=TRUE)
		    
# meta.data引数とmisc引数をFalseにしてSeuratオブジェクトに変換		   
seuratobj <- SeuratDisk::LoadH5Seurat("BL_hashing.h5seurat",
				       meta.data = F,
				       misc = F)

# h5adファイルを読み込み
h5ad <- anndata::read_h5ad(filename = "BL_hashing.h5ad")

# h5adファイルのobsの箇所をSeuratオブジェクトのmeta.dataスロットに登録
seuratobj <- AddMetaData(seuratobj, metadata = h5ad$obs)

ちなみにdepositされていた8つのh5adファイルを以下のlapplyですべて読み込んだ例も載せておく。

 h5adのファイルパス一覧を使ってすべてSeuratオブジェクトに変換してリストにまとめる。

paths <- dir(pattern = "h5ad")
sclist <- lapply(paths, function(x){
  SeuratDisk::Convert(source = x, 
		    dest="h5seurat",
		    overwrite=TRUE)
  h5ad <- anndata::read_h5ad(filename = x)
  path_h5seurat <- gsub(pattern = "h5ad", replacement = "h5seurat", x = x)
  x <- SeuratDisk::LoadH5Seurat(file = path_h5seurat, meta.data = F, misc = F)
  x <- AddMetaData(x, metadata = h5ad$obs)
  return(x)
})

Discussion

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