シングルセル解析ソフトScanpyを試してみる

この記事は創薬 Advent Calendar 2018 17日目の記事です。

シングルセル解析ソフトScanpyを試してみる

PythonのシングルセルRNA-seq解析ツールであるところのScanpyを阪大医学部Python会の@yyoshiakiさんに教えてもらったので、試してみました。

RだとSeuratというパッケージがいいらしいですが、Pythonの方をプッシュされたのでそちらを行いました。

筆者の知識レベル

  • in vitroの実験はやったことない。
  • FACSは同期がやってるので概要はなんとなく知ってる
  • RNA-seqはほぼ知識ゼロ

チュートリアル

最初のチュートリアルのZheng 2017をやりました。

論文の概要

ドロップ型のscRNA-seqシステムの開発

サンプル中の細胞をひとつひとつ高い確率でキャプチャでき、高速にmRNAを数えることができる。

  • droplet内に一つずつ細胞を入れ、その細胞のmRNAをエマルジョン(a Gel bead in EMulsion, GEM)内でまるっと数えることができる(一細胞RNA-seq)
  • cDNAへの逆転写はGEM内でバルクで進行する
  • GEMごとのプライマーにバーコード(14bp)がついてる
    • どのRNAがどの細胞由来かわかる
  • GEMごとのバーコードに加え、さらにプライマーごとにランダムなUMI(unique molecular identifier, 10 bp)がついてる
    • UMIの数を数えることで、その遺伝子の数がわかる。
    • UMIにはサンプルに対応するIDも入ってるので、複数サンプルを平行に処理することが可能
  • タグ(バーコード、UMI)がついてるので、PCR増幅はGEMを壊してまとめて行える

元文献では色々なサンプルを使って測定したデータを使って解析してましたが、今回はperipheral blood mononuclear cells (PBMCs) のクラスタリングをチュートリアルで行っていました。

チュートリアルのコード(改変あり)

元データ

1細胞データ

遺伝子発現ラベル(バルク)

全ファイル同じフォルダに入れてください

本家(scanpy == 1.0.0)とscanpyのバージョンが違い(scanpy == 1.3.6)、できないものがあったので、できたやつだけ書きました。

各種パッケージのインポート

import numpy as np
import pandas as pd
import matplotlib.pyplot as pl
import scanpy.api as sc

sc.settings.verbosity = 1  # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.set_figure_params(dpi=80)  # low dpi (dots per inch) yields small inline figures
sc.logging.print_versions()

scanpy==1.3.6 anndata==0.6.16 numpy==1.14.2 scipy==1.2.0 pandas==0.23.4 scikit-learn==0.20.2 statsmodels==0.9.0

データの読み込み

%%time
adata = sc.read('matrix.mtx', cache=True).T  # transpose the data
adata.var_names = pd.read_csv('genes.tsv', header=None, sep='\t')[1]
adata.obs_names = pd.read_csv('barcodes.tsv', header=None)[0]

Variable names are not unique. To make them unique, call `.var_names_make_unique`.


CPU times: user 1.16 s, sys: 386 ms, total: 1.54 s
Wall time: 1.65 s

adata

AnnData object with n_obs × n_vars = 68579 × 32738

重複した名前をユニークにする

adata.var_names_make_unique() #最初からユニークでは？
adata

AnnData object with n_obs × n_vars = 68579 × 32738

バルク(1細胞じゃない)での遺伝子発現ラベルの読み込み

adata.obs['bulk_labels'] = pd.read_csv('./zheng17_bulk_lables.txt', header=None)[0].values

変数を減らして、正規化後、Cell Rangerでフィルタリング

記事筆者はCell Rangerについて何も知らない(力尽きた)

%%time
sc.pp.filter_genes(adata, min_counts=1)  # only consider genes with more than 1 count
sc.pp.normalize_per_cell(adata)          # normalize with total UMI count per cell
filter_result = sc.pp.filter_genes_dispersion(
    adata.X, flavor='cell_ranger', n_top_genes=1000, log=False)

CPU times: user 2.07 s, sys: 979 ms, total: 3.04 s
Wall time: 3.08 s

プロット

対数プロットになるらしいが……

sc.pl.filter_genes_dispersion(filter_result, log=True)

f:id:kyusque:20181220163745p:plain

PCAによる次元圧縮

%%time
adata = adata[:, filter_result.gene_subset]  # filter genes
sc.pp.normalize_per_cell(adata)  # need to redo normalization after filtering
sc.pp.log1p(adata)  # log transform: X = log(X + 1)
sc.pp.scale(adata)
# the PCA is *not* contained in the recipe sc.pp.recipe_zheng17(adata)
sc.tl.pca(adata, n_comps=50)

CPU times: user 5.54 s, sys: 668 ms, total: 6.21 s
Wall time: 3.18 s

PCAのloadingのプロット

便利。

記事筆者が知らない遺伝子ばっかりだ…

sc.pl.pca_loadings(adata)

f:id:kyusque:20181220163837p:plain

K meansによるクラスタリング

1細胞のmRNAからの情報のみによる分類でも、バルクのグループに分かれてますね。

%%time
sc.pp.neighbors(adata)

CPU times: user 20.5 s, sys: 4.07 s, total: 24.6 s
Wall time: 21.7 s

%%time
sc.tl.umap(adata)

CPU times: user 56 s, sys: 615 ms, total: 56.6 s
Wall time: 49.4 s

sc.pl.umap(adata, color='bulk_labels')

... storing 'bulk_labels' as categorical

f:id:kyusque:20181220163901p:plain

メモ

  • 他にもマーカー遺伝子を取り出す手法などありましたが、現バージョンでそのままできなかったので今回は諦めました。
  • 移植後の患者からの細胞群の解析(こういう実験はin vivoに入るのでしょうか？in vitro？)など元データの文献自体も面白そうだったので、輪読に使えるのかなーと思います。
  • 元文献を読んでると知らない知識が多く、質問事項がえらいことになったので、おいおい勉強していこうと思います。