🎉

【シングルセル解析】[Seurat] scRNAseq公開データ読み込み例 ~ h5adファイル ~

2023/05/19に公開

bioinformatics

前置き

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

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

h5adファイル

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

SeuratDiskのインストール

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

例) zenodo.3969338

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

読み込むにはSeuratDiskパッケージのConvert()機能と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を使うには

おそらくRのanndataパッケージを使うには、Pythonのanndataパッケージも必要となる。

PCにPython環境が作られていなければ、R anndataの機能を使おうとした際にMinicondaを入れるか多分聞いてくるのでそれに従うとよい。
何もエラーが無ければ問題なし。
既にAnaconda3などを入れておりパスが通っていればその環境を使おうとする。

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

私の場合はreticulateパッケージの紹介のようにr-reticulate環境を作成して、そこにPython 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)
})

例) zenodo.7813150

https://zenodo.org/record/7813151

deposit dataの「scRNA_BCR_TCR.h5ad」ファイルを読み込む際に講じたエラー対策を記録しておく。

まずはスタンダードな読み込みを試す。

library(Seurat)
library(SeuratDisk)

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

LoadH5Seurat()での読み込み時にまた以下のエラーが出た。

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

上記で紹介したようにmeta.data = FALSE、misc = FALSEオプションを入れてみたが、エラーが解決しない。

sc <- LoadH5Seurat(file = "scRNA_BCR_TCR.h5seurat",
                   meta.data = FALSE,
                   misc = FALSE)

read_h5ad()でRに読み込み

仕方ないので、まずはanndataパッケージのread_h5ad()でRに読み込んで、そこからSeuratオブジェクトへ変換することにする。

reticulateパッケージのuse_condaenv()機能を使って、仮想環境を指定しておく。

anndataの使用準備

library(anndata)
library(reticulate)

use_condaenv("r-reticulate")

read_h5ad()機能でdepositデータのh5adファイルを読み込む。

h5adファイルの読み込み

h5ad <- read_h5ad("scRNA_BCR_TCR.h5ad")

中身はごちゃごちゃしてわかりにくいが、obsがmeta.dataスロットの箇所、obsmにUMAPの結果があることが推測できる。

h5adの中身

Seuratオブジェクト作成

発現マトリクスはXの箇所にsparse matrixが保存されており、行に細胞バーコード、列に遺伝子名の行列になっている。

CreateSeuratObject()機能に、転置した発現マトリクスXとmeta.dataのobsを指定する。depositされた細胞数/遺伝子数が減少してしまわないようにmin.cells=引数、min.features=引数は指定しない。

Seuratオブジェクトを作成

sc <- CreateSeuratObject(counts = Matrix::t(h5ad$X), 
                         meta.data = h5ad$obs)

UMAP情報の追加

depositデータにはUMAPの値が保存されていた。これを読み込めれば報告通りのUMAPの値が使用できる。

UMAP情報はobsmに入っていた。(UMAPのマトリクスが入ったリストになっていた。)

h5adからUMAP値を取り出し

umap <- h5ad$obsm
umap <- umap$X_umap

SeuratオブジェクトのreductionsスロットにはSeuratObjectパッケージのDimReducクラスのオブジェクトが必要となる。次のようにCreateDimReducObject()機能を使ってreductionsスロットに保存用のオブジェクトを作成する。

UMAP値をSeuratオブジェクトに入れる

# 列名を変更
colnames(umap) <- c("UMAP_1", "UMAP_2")
# UMAPの行名に細胞バーコードを入れる
rownames(umap) <- colnames(sc)

# DimReducクラスのオブジェクト作成
umap <- SeuratObject::CreateDimReducObject(embeddings = umap, 
                                          assay = "RNA", 
					  key = "UMAP_")

# Seuratオブジェクトのreductionsスロットに追加
sc@reductions[["umap"]] <- umap

meta.dataの細胞名っぽいところをactive.identに設定して、DimPlot()で確認。

Idents(sc) <- "refined_clusters"
DimPlot(sc)