RとPythonどちらを使う？

10x GenomicsのChromiumが販売され始めてからシングルセル解析は身近なツールになり多くの研究者に使われるようになっています*1。

それに伴い解析のためのツールも続々と開発されており、誰でもデータを利用しようと思えばできる環境がすでにできつつあります*2。そんな贅沢な状況の中で少し厄介なのは、開発に使われる言語がRとPythonで二分されていることです。

基本的な解析*3は、RではSeurat、PythonではScanpyですべて行うことが出来ます。CITEseq等のマルチオミクスデータへの対応はSeuratの方が優れている一方で、Scanpyで利用できるよう実装されている解析手法もあり *4、どちらもうまく活用したいのが本音です。

また、Trajectory inference・バッチ補正・ダブレット推定・細胞種アノテーションなどの工程は独立したR/Pythonライブラリとして公開されていることが多いです。そして、これらのライブラリを試そうと思う時に実装されている言語が自分が使っているものと違うと、少し面倒くさい訳です。

そんな中で最近はreticulateパッケージを使ってRとPythonをシームレスに使う方法が楽だと感じているのでメモします。

reticulateでRとPythonを両方使う

RパッケージreticulateはRとPythonをシームレスに実行するための機能が実装されています。具体的には以下のようなことが出来ます。(※ Pythonにもrpy2という類似した機能を持つライブラリが存在していますが、後述の理由によりscRNA-seq解析においてはreticulateの方が問題がなさそうでした)

RとPythonの間で変数を共有する
R/Rstudio上でPythonコードを走らせる

library(reticulate)
use_python("/Users/maruyamatooru/opt/anaconda3/envs/single_cell/bin/python") # 使いたいpythonのパスを指定する

a = list(c(1,2), c(3,4))

np = import("numpy")
mat = np$array(a)
mat

class(mat)

     [,1] [,2]
[1,]    1    2
[2,]    3    4

[1] "matrix"

これと同様の方法でR上でScanpyなどのPythonライブラリを操ることができるようになります。Seuratで前処理・クラスタリングした結果をScanpy等のPythonライブラリ上で使ってみましょう。

library(Seurat)

# Seuratでデータを読み込む
m = Read10X("public_dataset/10x/Parent_NGSC3_DI_PBMC_filtered_feature_bc_matrix/filtered_feature_bc_matrix/")
sc_seurat = CreateSeuratObject(m, assay="RNA")

# ミトコンドリアリード率・検出遺伝子数で細胞を選出
sc_seurat[["percent.mt"]] = PercentageFeatureSet(sc_seurat, pattern = "^MT-")
sc_seurat = subset(sc_seurat, subset = nFeature_RNA > 200 & percent.mt < 10)

# 正規化
sc_seurat = NormalizeData(sc_seurat, normalization.method = "LogNormalize", scale.factor = 10000)

# PCA & UMAP & Clustering
sc_seurat = FindVariableFeatures(sc_seurat, selection.method = "vst", nfeatures = 2000)
sc_seurat = ScaleData(sc_seurat)
sc_seurat = RunPCA(sc_seurat, features = VariableFeatures(sc_seurat))
sc_seurat = FindNeighbors(sc_seurat, dims = 1:30)
sc_seurat = FindClusters(sc_seurat, resolution = 0.5)
sc_seurat = RunUMAP(sc_seurat, dims = 1:30)
sc_seurat

An object of class Seurat 
36601 features across 8295 samples within 1 assay 
Active assay: RNA (36601 features, 2000 variable features)
 2 dimensional reductions calculated: pca, umap

以上の工程で出来上がったSeuratオブジェクトをPythonのAnnDataクラスへと変換してみます。

# AnnDataへ変換
sc = import("scanpy")
anndata = import("anndata")
sc_scanpy = anndata$AnnData(t(sc_seurat[["RNA"]]@data), 
                            obs=sc_seurat@meta.data,
                            var=data.frame(Symbol=rownames(sc_seurat)))
rownames(sc_scanpy$var) = rownames(sc_seurat)
sc_scanpy$obsm = list(pca=sc_seurat@reductions$pca@cell.embeddings[,1:2],
                      umap=sc_seurat@reductions$umap@cell.embeddings[,1:2])
sc_scanpy$var["highly_variable"] = sapply(sc_scanpy$var_names$tolist(), function(x) x %in% sc_seurat[["RNA"]]@var.features)
sc_scanpy

AnnData object with n_obs × n_vars = 8295 × 36601
    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.5', 'seurat_clusters'
    var: 'Symbol', 'highly_variable'
    obsm: 'pca', 'umap'

このような流れでscanpy.external API *5やその他のPythonライブラリと連携すれば一つのスクリプトで一括して解析を行うことが可能です。