Scanpyの真の力を僕達はまだ知らない

これまで、老化細胞データを読み込むために、ScanpyというPythonパッケージを利用していました。read_h5ad関数でH5ADファイルからデータを読み込んだ後は、他のツールを使って主成分分析（PCA）やt-SNE/UMAPを実行しました。

実は、データ読み込みだけでなく、PCAによる次元削減やt-SNE/UMAPでの可視化までをScanpyだけで完結させることができます。それだけでなく、不要なサンプルを削除したり正規化をかけたりする「前処理」もScanpyで行うことができます。

本記事では、読み込んだデータに対して前処理から可視化までをScanpyで実行する方法について紹介します。

データの読み込み

まずは、今回使用するデータを準備します。ファイルはwgetコマンドでダウンロードし、Scanpyのread_h5ad関数でデータを読み込みます。

import os
import scanpy as sc

os.system('wget https://storage.googleapis.com/calico-website-mca-storage/lung.h5ad')
adata = sc.read_h5ad('lung.h5ad')

これで、20,755細胞×15,106遺伝子のカウントデータが読み込まれ、変数adataに格納されました。

AnnData object with n_obs × n_vars = 20755 × 15106
    obs: 'age', 'animal', 'batch', 'n_counts', 'n_genes', 'percent_Rn45s', 'percent_mito', 'replicate'
    var: 'feature_types-0', 'feature_types-1', 'feature_types-10', 'feature_types-11', 'feature_types-12', 'feature_types-2', 'feature_types-3', 'feature_types-4', 'feature_types-5', 'feature_types-6', 'feature_types-7', 'feature_types-8', 'feature_types-9', 'gene_ids-0', 'gene_ids-1', 'gene_ids-10', 'gene_ids-11', 'gene_ids-12', 'gene_ids-2', 'gene_ids-3', 'gene_ids-4', 'gene_ids-5', 'gene_ids-6', 'gene_ids-7', 'gene_ids-8', 'gene_ids-9', 'genome-0', 'genome-1', 'genome-10', 'genome-11', 'genome-12', 'genome-2', 'genome-3', 'genome-4', 'genome-5', 'genome-6', 'genome-7', 'genome-8', 'genome-9', 'n_cells-0', 'n_cells-1', 'n_cells-10', 'n_cells-11', 'n_cells-12', 'n_cells-2', 'n_cells-3', 'n_cells-4', 'n_cells-5', 'n_cells-6', 'n_cells-7', 'n_cells-8', 'n_cells-9'
    obsm: 'X_pca', 'X_umap', 'X_tsne'

細胞ごとの合計カウントを100万に統一する

細胞によっては、mRNA量が多すぎるものや少なすぎるものがあります。この量はmRNAを取得する段階で増減してしまうものなので、実験的なバイアスを減らすという意味でも合計カウントを統一することが必要となります。

いくつに統一するかは自由ですが、この分野では合計カウントを100万に統一することが一般的です。これをCount Per Million（CPM）と読んだりします。

Scanpyでは、pp.normalize_totalを使うことでCPMを簡単に実現できます。

sc.pp.normalize_total(adata, target_sum=1e6)

パラメータ	初期値	説明
adata		Anndata型の注釈付きデータセット
target_sum	None	いくつに統一するか、Noneの場合は中央値が指定される
inplace	True	Trueならadataを更新、Falseならコピーを出力

データを対数変換する

この分野では、mRNAのカウントデータを対数変換することがよくあります。こうすることで、データが正規分布に従うようになることが理由のようです。

Scanpyでは、データの対数変換はpp.log1pで実現できます。

sc.pp.log1p(adata)

パラメータ	初期値	説明
X		Anndata型などのデータセット
base	None	対数変換の基底、Noneの場合は自然対数となる
copy	False	Trueならコピーを出力、FalseならXを更新（Anndataを渡したときのみ）

閾値未満の遺伝子を除去する

カウントデータは、最終的には有意差分析などの二次解析に使われることになります。なので、どの細胞においても全く活性化されていない遺伝子は、解析に不要なので除去した方が良いです。

Scanpyでは、mRNAカウントが0のときに、その遺伝子が不活性であると定義します。pp.filter_genesを使うと、指定した値未満の細胞でしか活性化していない遺伝子を除去できます。

sc.pp.filter_genes(adata, min_cells=10)

パラメータ	初期値	説明
data		Anndata型の注釈付きデータセット
min_counts	None	カウントの許容最小値
min_cells	None	遺伝子活性あり細胞の許容最小数
max_counts	None	カウントの許容最大値
max_cells	None	遺伝子活性あり細胞の許容最大数
inplace	True	Trueならadataを更新、Falseならコピーを出力

※ min_counts、min_cells、max_counts、max_cellsはどれか一つのみ指定可能

他にも、pp.highly_variable_genesを使えば、分散が一定以上の遺伝子のみを抽出することもできます。

sc. pp.highly_variable_genes(adata)
adata = adata[:, adata.var.highly_variable]

他にも、前処理のための関数は色々用意されているので、マニュアルを参考に試してみると良いと思います。

次元削減と可視化

最後に、これまでは別のライブラリで実装してきた次元削減（PCA）と可視化（t-SNEとUMAP）を、Scanpyで再現します。

sc.pp.pca(adata, n_comps=50)
sc.pl.pca(adata, components=['1,2', '1,3', '2,3'], color='age')

まずは、PCAからです。こちらは、tl.pcaで実行できます。

ScanpyではPCA結果を可視化する関数pl.pcaも用意されており、主成分を指定して二次元のプロットを表示させることができます。

sc.tl.tsne(adata)
sc.pl.tsne(adata, color='age')

sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.pl.umap(adata, color='age')

こちらは、t-SNEとUMAPをScanpyで実行しました。tl.tsneまたはtl.umapを使用します。

これらのプロットは、adata.obsのカテゴリデータで色付けすることができます。例えば、若年細胞と老化細胞で遺伝子活性の偏りが無いかをここで調べることができそうです。

ちなみに、別のライブラリで出力したt-SNEとUMAPの結果は以下の通りでした。

これを見る限りでは、明らかにScanpyの出力結果と異なるプロットが生成されているようです。二者の違いは現状では分かっていませんが、今後も調査を続けてt-SNEとUMAPの使い方をマスターしていきたいと思います。

まとめ

本記事では、Scanpyを使ってデータの前処理と次元削減・可視化を実行する方法について紹介しました。

Scanpyには、他にもt検定で遺伝子の有意差検定を行ったり、細胞の分化経路推定を行うための機能も有しているようです。また、Scanpyの開発に参加することもできるようで、コントリビューターになりたい人のためのページも公開されています。もし、Scanpyに実装したい解析方法があれば、ここから追加のリクエストを送ることができそうです。

このように、すでに様々な機能が存在するうえに、今後もさらなる発展が期待できるので、シングルセル解析を行いたい人はScanpyを一度は試してみると良いと思います。

以上です。