PythonによるscRNA-seq解析 その1¶

使用データ¶

マウス15.5日胚の眼球から取り出した網膜の細胞。10X Chromium 3′ v2で解析。2つのReplicatesでそれぞれ約2,600細胞。

Mouse developing retina from Lo Giudice et al. https://doi.org/10.1242/dev.178103

Lo Giudice, Q., Leleu, M., La Manno, G. & Fabre, P. J. Single-cell transcriptional logic of cell-fate specification and axon guidance in early-born retinal neurons. Development 146, dev178103 (2019).

Pythonに関係ない部分の解析¶

データのダウンロード¶

好きな方法で公共データベースからダウンロードする。レプリケイトがふたつあって、それぞれ以下のSRR IDに対応。

fastq-dump SRR8181428 --split-files --gzip
fastq-dump SRR8181429 --split-files --gzip

SRR8181428 (Rub browser)(メタデータ )

SRR8181428_1.fastq.gz ：Read 1  
SRR8181428_2.fastq.gz ：Read 2  
SRR8181428_3.fastq.gz ：Index

SRR8181429 (Rub browser)(メタデータ )

SRR8181428_1.fastq.gz ：Read 1  
SRR8181428_2.fastq.gz ：Read 2  
SRR8181428_3.fastq.gz ：Index

Cell Ranger は、以下の命名ルールに沿った fastq file 名が必要

'[サンプル名]_S[サンプル番号]_L00[レーン番号]_[リードタイプ]_001.fastq.gz'

今回は、以下のようにファイル名を変更した。

Read 1: SRR8181428_S1_L001_R1_001.fastq.gz
Read 2: SRR8181428_S1_L001_R2_001.fastq.gz
Index : SRR8181428_S1_L001_I1_001.fastq.gz

Read 1: SRR8181429_S1_L001_R1_001.fastq.gz
Read 2: SRR8181429_S1_L001_R1_001.fastq.gz
Index : SRR8181429_S1_L001_I1_001.fastq.gz

Cell Rangerによる解析¶

Cell Ranger v.9.0.0でマウスゲノム（mm10）へのマッピング、遺伝子ごとのカウントを実行。

レプリケイトは（論文でそう名付けられてたので）E2, F2という名前のバッチに設定。

なお、最新版のcellrangerでは "--create-bam" オプションが必須になっている点に注意。マッピング結果のBAM形式ファイルを利用する予定がなければ、計算量とディスク節約のためにfalseに設定すれば良いが、ツールによって（後半で扱うvelocytoなど）BAMファイルが必須の場合がある。

cellranger count --id=RetinalBatchE2 \
   --fastqs=/path/to/SRR8181428\
   --sample=SRR8181428 \
   --create-bam=true \
   --transcriptome=/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A

cellranger count --id=RetinalBatchF2 \
   --fastqs=/path/to/SRR8181429\
   --sample=SRR8181429 \
   --create-bam=true \
   --transcriptome=/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A

講習では、このCell Ranger解析の結果ディレクトリ（の一部、講習で使う部分だけを抜き出したもの）を配布している。

こんな感じの構成になっていて、Cell Rangerのアウトプットが見えるはず。

./data の中身

 data -- RetinalBatchE2 - outs - filtered_feature_bc_matrix - barcodes.tsv.gz  # cell index  
                                                            - features.tsv.gz  # 遺伝子 Id  
                                                             - matrix.mtx.gz    # count データ  

      -- RetinalBatchF2 - outs - filtered_feature_bc_matrix - barcodes.tsv.gz  # cell index
                                                            - features.tsv.gz  # 遺伝子 Id
                                                            - matrix.mtx.gz    # count データ

ライブラリのimport¶

定番のライブラリのimportに加えて、scanpyをscの別名でimportする。

Figのちょっとした設定。各ライブラリのバージョンは以下の通り。(インストールした際の python のバージョンなどで ライブラリのバージョンが異なっている可能性はあるかもしれません。)

In [1]:
import warnings
warnings.filterwarnings('ignore')

import os
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
import seaborn as sns

sc.logging.print_header()
sc.settings.set_figure_params(dpi=100, facecolor='white')

scanpy==1.10.4 anndata==0.11.3 umap==0.5.7 numpy==2.1.3 scipy==1.15.1 pandas==2.2.2 scikit-learn==1.6.1 statsmodels==0.14.4 igraph==0.11.8 pynndescent==0.5.13

データの読み込み¶

Scanpyの特徴は、データフレームを拡張した anndata と呼ばれる オブジェクトを使う点にある。anndataを使うと、ひとつのオブジェクトに遺伝子発現量のデータ、サンプルや細胞のアノテーション、遺伝子の情報などをまとめて格納できる。 anndataを使うことで、実験の情報が詰まったひとつのオブジェクトに対して処理を次々に実行し、さらに処理結果をそこに蓄積していくことができる。

Scanpyでは、10x Genomicsのデータは、結果のディレクトリを指定することでそのままロードが可能な関数が用意されている。

ディレクトリには、3つのファイルが書き出されている。

barcodes.tsv: 個別の細胞を識別するバーコードが記述されたテキストファイル。各行ごとにひとつのバーコードが記述されている。
genes.tsv: 計測された遺伝子が記述されたタブ区切りテキストファイル。1列目に遺伝子のENSEMBL Gene ID、2列目の遺伝子のシンボルが記述されている。
matrix.mtx: 各細胞（バーコード）について各遺伝子の発現がいくつ観測されたかのカウント情報をまとめたファイル。このファイルはMatrix Market Exchange Formatsという形式で記述されており、疎行列（含まれる値の多くがゼロであるような行列）を比較的コンパクトに記述するための形式となっている。

Scanpyの10xデータ用読み込み関数を使うと、これらを同時に読み込んで、適切に構造化されたオブジェクトができあがる。

In [2]:
adata_E2 = sc.read_10x_mtx(path='./data/RetinalBatchE2/outs/filtered_feature_bc_matrix/')
adata_F2 = sc.read_10x_mtx(path='./data/RetinalBatchF2/outs/filtered_feature_bc_matrix/')
In [3]:
adata_E2
Out[3]:
AnnData object with n_obs × n_vars = 3392 × 32285
    var: 'gene_ids', 'feature_types'

3392 細胞 × 32285 遺伝子

In [4]:
adata_F2
Out[4]:
AnnData object with n_obs × n_vars = 3611 × 32285
    var: 'gene_ids', 'feature_types'

3611 細胞 × 32285 遺伝子

2つのバッチをひとつのanndataオブジェクトに統合する。バッチのラベルもここで設定。

In [5]:
adata = adata_E2.concatenate(adata_F2, batch_categories=['E2', 'F2'])

観測値（細胞）に関するデータは、obs でアクセスする。 observationsの略。このデータフレームに細胞に関するデータを追加していくことができる。

In [6]:
adata.obs
Out[6]:
batch
AAACCTGAGATGTCGG-1-E2 E2
AAACCTGCAATCCAAC-1-E2 E2
AAACCTGCACATGGGA-1-E2 E2
AAACCTGCAGCAGTTT-1-E2 E2
AAACCTGGTTCCTCCA-1-E2 E2
... ...
TTTGTCACATCGATGT-1-F2 F2
TTTGTCAGTAGAGGAA-1-F2 F2
TTTGTCAGTCATCCCT-1-F2 F2
TTTGTCATCAGTTCGA-1-F2 F2
TTTGTCATCCGAATGT-1-F2 F2

7003 rows × 1 columns

変数（遺伝子）に関するデータは、var でアクセスする。variable の略。現在は、遺伝子のシンボル（名前）と、Ensembl Gene IDなどが登録されている。

In [7]:
adata.var
Out[7]:
gene_ids feature_types
Xkr4 ENSMUSG00000051951 Gene Expression
Gm1992 ENSMUSG00000089699 Gene Expression
Gm19938 ENSMUSG00000102331 Gene Expression
Gm37381 ENSMUSG00000102343 Gene Expression
Rp1 ENSMUSG00000025900 Gene Expression
... ... ...
AC124606.1 ENSMUSG00000095523 Gene Expression
AC133095.2 ENSMUSG00000095475 Gene Expression
AC133095.1 ENSMUSG00000094855 Gene Expression
AC234645.1 ENSMUSG00000095019 Gene Expression
AC149090.1 ENSMUSG00000095041 Gene Expression

32285 rows × 2 columns

カウントテーブルには、以下のように X でアクセスできる。疎行列として格納されている。

In [8]:
adata.X
Out[8]:
<Compressed Sparse Row sparse matrix of dtype 'float32'
	with 14123833 stored elements and shape (7003, 32285)>
In [9]:
adata
Out[9]:
AnnData object with n_obs × n_vars = 7003 × 32285
    obs: 'batch'
    var: 'gene_ids', 'feature_types'

E2: 3611 細胞 × 32285 遺伝子
F2: 3392 細胞 × 32285 遺伝子 を足したデータ 7003 細胞 × 32285 遺伝子 が格納されている。

前処理¶

観測された遺伝子が極端に少ない細胞、割り当てられた細胞が極端に少ない遺伝子をフィルタリング。

In [10]:
# 観測された遺伝子が極端に少ない細胞
sc.pp.filter_cells(adata, min_genes=200)
# 割り当てられた細胞が極端に少ない遺伝子
sc.pp.filter_genes(adata, min_cells=20)
In [11]:
adata
Out[11]:
AnnData object with n_obs × n_vars = 7001 × 15001
    obs: 'batch', 'n_genes'
    var: 'gene_ids', 'feature_types', 'n_cells'

7003 細胞 × 32285 遺伝子 -> 7001 細胞 × 15001 遺伝子

ミトコンドリア遺伝子の発現割合でフィルタリングする予定なので、遺伝子側のメタデータにミトコンドリア遺伝子か否かの情報を付与。

In [12]:
adata.var['mt'] = adata.var_names.str.startswith('mt-')
In [13]:
adata.var['mt']
Out[13]:
Xkr4              False
Gm1992            False
Gm19938           False
Rp1               False
Mrpl15            False
                  ...  
CAAA01118383.1    False
Vamp7             False
Tmlhe             False
4933409K07Rik     False
AC149090.1        False
Name: mt, Length: 15001, dtype: bool

sc.pp.calculate_qc_metrics は、クオリティコントロールに関連したいくつかの指標を一挙に計算してくれる便利な関数。

qc_varsに遺伝子側のメタデータラベルを設定すると、それについてもカウントとパーセンテージを勝手に計算してくれる。

In [14]:
# inplace オプション: 計算した結果を .obs と .var に書き込むか？
# percent_top オプション: Which proportions of top genes to cover. If empty or None don’t calculate. Values are considered 1-indexed, percent_top=[50] finds cumulative proportion to the 50th most expressed gene.
# log1p オプション: Set to False to skip computing log1p transformed annotations.
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=False, log1p=False, inplace=True)

それぞれの分布をプロット。

In [15]:
sc.pl.violin(adata, ["n_genes_by_counts", "total_counts", "pct_counts_mt"], jitter=0.4, multi_panel=True)
No description has been provided for this image

n_genes_by_counts: 各細胞ごとにカウントがあった遺伝子数 #total_counts: 各細胞ごとの合計カウント数 #pct_mt_counts: ミトコンドリアのカウントのパーセント¶

In [16]:
sc.pl.scatter(adata, 'total_counts', 'n_genes_by_counts', color='pct_counts_mt', size=40)
No description has been provided for this image
In [17]:
fig = plt.figure()
sns.displot(adata.obs['pct_counts_mt'], kde=False)
plt.show()

<Figure size 400x400 with 0 Axes>
No description has been provided for this image
In [18]:
fig = plt.figure()
sns.displot(adata.obs['pct_counts_mt'][adata.obs['pct_counts_mt'] < 10], kde=False)
plt.show()

<Figure size 400x400 with 0 Axes>
No description has been provided for this image

分布を見て方針を決めて、フィルタリングを実行。

In [19]:
print('Total number of cells: {:d}'.format(adata.n_obs))

sc.pp.filter_cells(adata, min_counts = 2000)
print('Number of cells after min count filter: {:d}'.format(adata.n_obs))

sc.pp.filter_cells(adata, max_counts = 13000)
print('Number of cells after max count filter: {:d}'.format(adata.n_obs))

sc.pp.filter_cells(adata, min_genes = 1000)
print('Number of cells after gene filter: {:d}'.format(adata.n_obs))

adata = adata[adata.obs['pct_counts_mt'] < 6]
print('Number of cells after MT filter: {:d}'.format(adata.n_obs))

Total number of cells: 7001
Number of cells after min count filter: 5174
Number of cells after max count filter: 4981
Number of cells after gene filter: 4976
Number of cells after MT filter: 4773
In [20]:
adata
Out[20]:
View of AnnData object with n_obs × n_vars = 4773 × 15001
    obs: 'batch', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'n_counts'
    var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'

7003 細胞 × 32285 遺伝子 -> 7001 細胞 × 15001 遺伝子 -> 4773 細胞 × 15001 遺伝子

正規化や対数変換などの前に、この段階のカウントマトリックスを別のレイヤーに退避させておく。あとで確率モデル推定のときに必要になる。レイヤーに保管したデータは基本的に以降の操作の影響を受けないが、細胞や遺伝子のslicingは適用されるので注意。（常にマトリックスのshapeは一致する）

In [21]:
adata.layers['counts'] = adata.X.copy()

正規化¶

ライブラリサイズによる正規化、対数変換などの前処理は、scanpy.pp以下にいくつか便利な関数がある。

ここでは、細胞ごとのカウントの和が10,000になるように正規化してから、対数変換。

In [22]:
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
sc.pp.log1p(adata)

この段階のデータ（いろいろとややこしい変換をしていない「ピュア」なデータ）も、あとでプロットや結果解釈のときに使いたいので退避させておく。rawは特別なレイヤーで、数値だけでなくobsやvarのメタデータも含めてフリーズしてくれる。また遺伝子側のsliceの影響を受けない（細胞側のsliceは適用される）

In [23]:
adata.raw = adata

特徴量選択（発現量の変動が大きい遺伝子）¶

発現量変動の大きい遺伝子のみを抽出して、データのサイズを小さくする。

発現量変動の大きい遺伝子とは、細胞間の発現量の分散が指定した値以上に大きい遺伝子を取得すればいいが、遺伝子の平均発現量毎に分散も異なるので、正規化分散をしてから特徴量選択をする。 seurat, cell_ranger, seurat_v3 のアルゴリズムが選択できる。 本講習ではデフォルトの seurat を使用する。

In [24]:
sc.pp.highly_variable_genes?

Signature:
sc.pp.highly_variable_genes(
    adata: 'AnnData',
    *,
    layer: 'str | None' = None,
    n_top_genes: 'int | None' = None,
    min_disp: 'float' = 0.5,
    max_disp: 'float' = inf,
    min_mean: 'float' = 0.0125,
    max_mean: 'float' = 3,
    span: 'float' = 0.3,
    n_bins: 'int' = 20,
    flavor: "Literal['seurat', 'cell_ranger', 'seurat_v3', 'seurat_v3_paper']" = 'seurat',
    subset: 'bool' = False,
    inplace: 'bool' = True,
    batch_key: 'str | None' = None,
    check_values: 'bool' = True,
) -> 'pd.DataFrame | None'
Docstring:
Annotate highly variable genes :cite:p:`Satija2015,Zheng2017,Stuart2019`.

Expects logarithmized data, except when `flavor='seurat_v3'`/`'seurat_v3_paper'`, in which count
data is expected.

Depending on `flavor`, this reproduces the R-implementations of Seurat
:cite:p:`Satija2015`, Cell Ranger :cite:p:`Zheng2017`, and Seurat v3 :cite:p:`Stuart2019`.

`'seurat_v3'`/`'seurat_v3_paper'` requires `scikit-misc` package. If you plan to use this flavor, consider
installing `scanpy` with this optional dependency: `scanpy[skmisc]`.

For the dispersion-based methods (`flavor='seurat'` :cite:t:`Satija2015` and
`flavor='cell_ranger'` :cite:t:`Zheng2017`), the normalized dispersion is obtained
by scaling with the mean and standard deviation of the dispersions for genes
falling into a given bin for mean expression of genes. This means that for each
bin of mean expression, highly variable genes are selected.

For `flavor='seurat_v3'`/`'seurat_v3_paper'` :cite:p:`Stuart2019`, a normalized variance for each gene
is computed. First, the data are standardized (i.e., z-score normalization
per feature) with a regularized standard deviation. Next, the normalized variance
is computed as the variance of each gene after the transformation. Genes are ranked
by the normalized variance.
Only if `batch_key` is not `None`, the two flavors differ: For `flavor='seurat_v3'`, genes are first sorted by the median (across batches) rank, with ties broken by the number of batches a gene is a HVG.
For `flavor='seurat_v3_paper'`, genes are first sorted by the number of batches a gene is a HVG, with ties broken by the median (across batches) rank.

The following may help when comparing to Seurat's naming:
If `batch_key=None` and `flavor='seurat'`, this mimics Seurat's `FindVariableFeatures(…, method='mean.var.plot')`.
If `batch_key=None` and `flavor='seurat_v3'`/`flavor='seurat_v3_paper'`, this mimics Seurat's `FindVariableFeatures(..., method='vst')`.
If `batch_key` is not `None` and `flavor='seurat_v3_paper'`, this mimics Seurat's `SelectIntegrationFeatures`.

See also `scanpy.experimental.pp._highly_variable_genes` for additional flavors
(e.g. Pearson residuals).

Parameters
----------
adata : 'AnnData'
    The annotated data matrix of shape `n_obs` × `n_vars`. Rows correspond
    to cells and columns to genes.
layer : 'str | None', optional (default: None)
    If provided, use `adata.layers[layer]` for expression values instead of `adata.X`.
n_top_genes : 'int | None', optional (default: None)
    Number of highly-variable genes to keep. Mandatory if `flavor='seurat_v3'`.
min_mean : 'float', optional (default: 0.0125)
    If `n_top_genes` unequals `None`, this and all other cutoffs for the means and the
    normalized dispersions are ignored. Ignored if `flavor='seurat_v3'`.
max_mean : 'float', optional (default: 3)
    If `n_top_genes` unequals `None`, this and all other cutoffs for the means and the
    normalized dispersions are ignored. Ignored if `flavor='seurat_v3'`.
min_disp : 'float', optional (default: 0.5)
    If `n_top_genes` unequals `None`, this and all other cutoffs for the means and the
    normalized dispersions are ignored. Ignored if `flavor='seurat_v3'`.
max_disp : 'float', optional (default: inf)
    If `n_top_genes` unequals `None`, this and all other cutoffs for the means and the
    normalized dispersions are ignored. Ignored if `flavor='seurat_v3'`.
span : 'float', optional (default: 0.3)
    The fraction of the data (cells) used when estimating the variance in the loess
    model fit if `flavor='seurat_v3'`.
n_bins : 'int', optional (default: 20)
    Number of bins for binning the mean gene expression. Normalization is
    done with respect to each bin. If just a single gene falls into a bin,
    the normalized dispersion is artificially set to 1. You'll be informed
    about this if you set `settings.verbosity = 4`.
flavor : "Literal['seurat', 'cell_ranger', 'seurat_v3', 'seurat_v3_paper']", optional (default: 'seurat')
    Choose the flavor for identifying highly variable genes. For the dispersion
    based methods in their default workflows, Seurat passes the cutoffs whereas
    Cell Ranger passes `n_top_genes`.
subset : 'bool', optional (default: False)
    Inplace subset to highly-variable genes if `True` otherwise merely indicate
    highly variable genes.
inplace : 'bool', optional (default: True)
    Whether to place calculated metrics in `.var` or return them.
batch_key : 'str | None', optional (default: None)
    If specified, highly-variable genes are selected within each batch separately and merged.
    This simple process avoids the selection of batch-specific genes and acts as a
    lightweight batch correction method. For all flavors, except `seurat_v3`, genes are first sorted
    by how many batches they are a HVG. For dispersion-based flavors ties are broken
    by normalized dispersion. For `flavor = 'seurat_v3_paper'`, ties are broken by the median
    (across batches) rank based on within-batch normalized variance.
check_values : 'bool', optional (default: True)
    Check if counts in selected layer are integers. A Warning is returned if set to True.
    Only used if `flavor='seurat_v3'`/`'seurat_v3_paper'`.

Returns
-------
Returns a :class:`pandas.DataFrame` with calculated metrics if `inplace=False`, else returns an `AnnData` object where it sets the following field:

`adata.var['highly_variable']` : :class:`pandas.Series` (dtype `bool`)
    boolean indicator of highly-variable genes
`adata.var['means']` : :class:`pandas.Series` (dtype `float`)
    means per gene
`adata.var['dispersions']` : :class:`pandas.Series` (dtype `float`)
    For dispersion-based flavors, dispersions per gene
`adata.var['dispersions_norm']` : :class:`pandas.Series` (dtype `float`)
    For dispersion-based flavors, normalized dispersions per gene
`adata.var['variances']` : :class:`pandas.Series` (dtype `float`)
    For `flavor='seurat_v3'`/`'seurat_v3_paper'`, variance per gene
`adata.var['variances_norm']`/`'seurat_v3_paper'` : :class:`pandas.Series` (dtype `float`)
    For `flavor='seurat_v3'`/`'seurat_v3_paper'`, normalized variance per gene, averaged in
    the case of multiple batches
`adata.var['highly_variable_rank']` : :class:`pandas.Series` (dtype `float`)
    For `flavor='seurat_v3'`/`'seurat_v3_paper'`, rank of the gene according to normalized
    variance, in case of multiple batches description above
`adata.var['highly_variable_nbatches']` : :class:`pandas.Series` (dtype `int`)
    If `batch_key` is given, this denotes in how many batches genes are detected as HVG
`adata.var['highly_variable_intersection']` : :class:`pandas.Series` (dtype `bool`)
    If `batch_key` is given, this denotes the genes that are highly variable in all batches

Notes
-----
This function replaces :func:`~scanpy.pp.filter_genes_dispersion`.
File:      ~/miniforge3/envs/pags2024/lib/python3.11/site-packages/scanpy/preprocessing/_highly_variable_genes.py
Type:      function
In [25]:
# top 2000genes のみを抽出する
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor='seurat')
print('\n','Number of highly variable genes: {:d}'.format(np.sum(adata.var['highly_variable'])))

 Number of highly variable genes: 2000

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
In [26]:
sc.pl.highly_variable_genes(adata)
No description has been provided for this image

計算結果は自動的に adata.var 、つまり、遺伝子に関するメタデータを格納したオブジェクトに追加される。 highly_variable の項目が True の遺伝子が高発現変動遺伝子。

In [27]:
adata.var
Out[27]:
gene_ids feature_types n_cells mt n_cells_by_counts mean_counts pct_dropout_by_counts total_counts highly_variable means dispersions dispersions_norm
Xkr4 ENSMUSG00000051951 Gene Expression 2817 False 2817 1.286673 59.762891 9008.0 True 1.159619 1.738898 1.326996
Gm1992 ENSMUSG00000089699 Gene Expression 533 False 533 0.088416 92.386802 619.0 False 0.115807 0.803469 -0.190706
Gm19938 ENSMUSG00000102331 Gene Expression 588 False 588 0.096843 91.601200 678.0 False 0.165617 0.838521 -0.025215
Rp1 ENSMUSG00000025900 Gene Expression 21 False 21 0.004142 99.700043 29.0 True 0.014691 1.690833 3.998795
Mrpl15 ENSMUSG00000033845 Gene Expression 1963 False 1963 0.359234 71.961148 2515.0 False 0.617765 0.783958 -0.605477
... ... ... ... ... ... ... ... ... ... ... ... ...
CAAA01118383.1 ENSMUSG00000063897 Gene Expression 205 False 205 0.029853 97.071847 209.0 False 0.066917 0.785096 -0.277454
Vamp7 ENSMUSG00000051412 Gene Expression 1735 False 1735 0.315241 75.217826 2207.0 False 0.534669 0.744384 -0.595847
Tmlhe ENSMUSG00000079834 Gene Expression 738 False 738 0.117269 89.458649 821.0 False 0.200372 0.883025 0.184900
4933409K07Rik ENSMUSG00000095552 Gene Expression 20 False 20 0.002857 99.714327 20.0 False 0.003138 0.025211 -3.865090
AC149090.1 ENSMUSG00000095041 Gene Expression 2249 False 2249 0.533210 67.876018 3733.0 False 0.692310 1.040684 0.181685

15001 rows × 12 columns

次元削減¶

主成分分析 (PCA) は、データ点全体の分布の大局的 (グローバル) な構造をよりよく表現する。 それに対して、t-SNE, UMAP は、どちらかというと局所的 (ローカル) な構造を重視する。高次元空間で近くに配置されている点を低次元空間でも近くに配置していく。その結果、どの細胞とどの細胞が似ているか、似ている細胞集団から構成されるクラスタがいくつあるのかといった生物学的な解釈に結びつけやすいので、scRNA-seq解析でよく使われている。しかしながら、クラスタ間の距離関係、このクラスターは、あっちのクラスターとどのくらい離れている、、とかはわからない。クラスタ間の距離の比較はできない。あくまで隣接している細胞同士の距離関係が最適化するように学習されている。

t-SNE のアルゴリズム
高次元空間と低次元空間の距離関係を一致させて、どうやって比較するか？

  1. 低次元側のデータ点の座標をランダムに (あるいは主成分分析やなんらかの高速な手法で暫定的に) 配置する
  2. 高次元空間上の「距離」と低次元空間上の「距離」を計測するための関数を設定する
  3. 高次元空間上の「距離」と低次元空間上の「距離」を比較して、誤差を表現する (コスト関数) を設定する
  4. コスト関数をなんらかの最適化計算手法で最小化する

tSNE と比較した UMAP の利点

  1. 計算の実行時間が tSNE より早い (が、実装によっては t-SNE も十分速い)
  2. 大規模なデータセットでも t-SNE ほどメモリを消費せずにすむ
  3. データ全体のグローバルな構造を tSNE よりもよく保存する (と言われている)
  4. アルゴリズムが代数的位相幾何学の強固な理論的基盤に支えられている (?)

　(独習Pythonバイオ情報解析　第11章 シングルセル解析② 次元削減 より)

t-SNE, UMAP ともに、計算効率のために、データの前処理として PCA を行って、特徴量の数を削減したセットを使用する。

主成分分析（PCA）¶

PCAはscanpyの前処理関数で簡単に実行できる。とりあえず、50次元まで落としてみる。 use_highly_variableのフラグをオンにすると、遺伝子全体ではなく前項で決定した高発現変動遺伝子のみを使って次元削減をする。

In [28]:
sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')

PCAで落とした座標は、観察値のメタデータを格納する obsm に自動的に入る。obsm は Multi-dimensional annotation of observations の略。複数の次元でひとつの何かを表現するような観測値のメタデータがここに格納される。

In [29]:
print(adata.obsm['X_pca'].shape)

(4773, 50)

7003 細胞 × 32285 遺伝子 -> 7001 細胞 × 15001 遺伝子 -> 4773 細胞 × 15001 遺伝子 -> 4773 細胞 × 50 次元

プロットしてみる。Scanpyでは基本的に「前処理」（Preprocessing）に関わる関数がscanpy.ppに、「プロット」（Plot）に関わる関数がscanpy.plに入っている。

いまのところ特に色つけるようなデータも無いので、適当に細胞ごとのカウントで色分けをした。遺伝子名をcolorに指定するとその遺伝子の発現量によるプロットも可能。

In [30]:
sc.pl.pca(adata, color='total_counts')
No description has been provided for this image

遺伝子 Xkr4 の発現量による色付けの例

In [31]:
sc.pl.pca(adata, color='Xkr4')
No description has been provided for this image
In [32]:
sc.pl.pca(adata, color='total_counts', 
          components=['1,2', '2,3', '1,3'])
No description has been provided for this image
In [33]:
# projection='3d'
# を指定すると3次元表示（見づらい）
sc.pl.pca(adata, color='total_counts', projection='3d')
No description has been provided for this image
In [34]:
sc.pl.pca(adata, color='batch', projection='3d')
No description has been provided for this image

寄与率をグラフで見る

In [35]:
sc.pl.pca_variance_ratio(adata, n_pcs=50)
No description has been provided for this image

t分布型確率的近傍埋め込み（t-SNE）¶

まず、t-SNE, UMAP共通のステップとして、データから「近傍グラフ」（neighborhood graph）の構築が必要。

In [36]:
sc.pp.neighbors(adata)

データ点間の接続関係（細胞の近傍関係）は、全細胞vs.全細胞のペアの情報を記録する obsp （Pairwise annotation of observations の略）に格納される。

In [37]:
adata.obsp['connectivities']
Out[37]:
<Compressed Sparse Row sparse matrix of dtype 'float32'
	with 96074 stored elements and shape (4773, 4773)>

接続関係を元にしてt-SNEを実行。

In [38]:
sc.tl.tsne(adata)

t-SNE結果のプロット。

In [39]:
sc.pl.tsne(adata, color='total_counts')
No description has been provided for this image

UMAP¶

近傍グラフはすでに計算しているので、scanpy.tlのumap関数を使えばオーケー。

In [40]:
sc.tl.umap(adata)

プロットはt-SNEとほとんど同じ。scanpy.pl以下にUMAP用の描画関数がある。

In [41]:
sc.pl.umap(adata, color='total_counts')
No description has been provided for this image

あきらかに似た形状の、ふたつの構造がある。。。

各細胞が由来するバッチで色分けしてみると・・・

In [42]:
sc.pl.umap(adata, color='batch')
No description has been provided for this image

ということで、このふたつの構造はもろにバッチ効果の結果だった。なんらかの方法で補正が必要だが、それは後述。

クラスタリング¶

グラフベースのクラスタリング Leiden クラスタリング¶

高次空間のデータ点からK近傍グラフ構造を作成して、データ点のクラスタリングを、グラフのコミュニティ検出問題として扱う。グラフ構造から「コミュニティ (点のカタマリ、クラスタ)」を見つける。モジュラリティ:コミュニティへの割り当ての良さの指標。このモジュラリティが最大になるようなコミュニティ割り当てを目指す。

Leidenクラスタリングを実行してみる。モジュラリティの計算に影響を与える "resolution" パラメータが存在する。
"resolution" は値が大きいほどクラスタが多くなる。default=1

In [43]:
sc.tl.leiden(adata, resolution=0.5, key_added='leiden_r0.5')

クラスタリングの結果は観測値のメタデータ（obs）の中に格納される。

In [44]:
adata.obs
Out[44]:
batch n_genes n_genes_by_counts total_counts total_counts_mt pct_counts_mt n_counts leiden_r0.5
AAACCTGAGATGTCGG-1-E2 E2 1776 1776 3776.0 112.0 2.966102 3776.0 1
AAACCTGCAATCCAAC-1-E2 E2 1697 1697 3113.0 97.0 3.115965 3113.0 1
AAACCTGTCCAATGGT-1-E2 E2 2102 2102 3925.0 91.0 2.318471 3925.0 8
AAACGGGAGGCAATTA-1-E2 E2 3440 3440 10695.0 224.0 2.094437 10695.0 4
AAACGGGCACGGCGTT-1-E2 E2 2019 2019 4322.0 228.0 5.275335 4322.0 4
... ... ... ... ... ... ... ... ...
TTTGTCAAGGGCTCTC-1-F2 F2 3009 3009 7465.0 143.0 1.915606 7465.0 9
TTTGTCACATCGATGT-1-F2 F2 2507 2507 5026.0 135.0 2.686033 5026.0 3
TTTGTCAGTAGAGGAA-1-F2 F2 2573 2573 6615.0 176.0 2.660620 6615.0 0
TTTGTCAGTCATCCCT-1-F2 F2 2389 2389 4923.0 121.0 2.457851 4923.0 5
TTTGTCATCAGTTCGA-1-F2 F2 1610 1610 2835.0 54.0 1.904762 2835.0 3

4773 rows × 8 columns

プロットする点の座標じたいはUMAPの結果なので、プロットはUMAP版の関数を使い、色分けだけをクラスタリング結果で指定すればいい。

In [45]:
sc.pl.umap(adata,
           color=['batch', 'leiden_r0.5'],
           ncols=2,
           frameon=False)
No description has been provided for this image

深層生成モデルの利用¶

scVIが仮定しているデータの生成プロセスは以下の確率モデルに基づく。
詳細は、2023年の講習会で東さんが説明されているのでそちらの動画をご参照ください。(https://www.genome-sci.jp/lecture202212)

全細胞数がN, 細胞のインデックスが $n \in \{1..N\}$、全遺伝子数がG、遺伝子のインデックスが $g \in \{1..G\}$、細胞nのバッチIDが$s_n$であるとしたとき、細胞nの遺伝子gのカウントは、

$$z_n \sim Normal(0, I)$$

$$l_n \sim logNormal(l_\mu, l_\sigma^2)$$

$$\rho_n = f_w (z_n, s_n)$$

$$w_{ng} = Gamma(\rho_n^g, \theta)$$

$$y_{ng} = Poisson(l_n w_{ng})$$

$$h_{ng} = Bernoulli( f_h^g(z_n, s_n) )$$

$$ \begin{equation} x_{ng} = \begin{cases} y_{ng} & \text{if}\ h_{ng}=0 \\ 0 & \text{otherwise} \end{cases} \end{equation} $$

$l_\mu, l_\sigma$は、スケーリングファクターでバッチの数をBとしたとき$l_\mu,l_\sigma \in \mathbb{R}^B$で、実際のバッチごとの平均と分散に設定しておく。

ガンマ分布のパラメータ、ベルヌーイ分布のパラメータは、正規分布からサンプリングした$z_n$を使ってニューラルネットワークで表現する（VAEにおけるreparametrization trick）。

最終的なカウントはガンマとポアソンの複合分布なので、負の二項分布でモデル化してることになる。さらにベルヌーイ分布で、scRNA-seqにありがちなdrop-outイベントを表現。これらの組み合わせで、ゼロ過剰負の二項分布(Zero-inflated negative binomial distribution; ZINB)を表現している。

事後分布はVAEをSGDで最適化して変分推論。

学習後、$z_n$にアクセスすれば、バッチの効果が除去された細胞の潜在ベクトルが得られるし、$\rho_{ng}$にアクセスすればdrop-outなどのイベントが生じなかった場合の遺伝子発現の期待値が得られる。

"drop-out" とは

scRNA-seqでは、各細胞のシーケンスの量が少ないため、特定の遺伝子の発現が検出されないことがある。遺伝子の発現データが欠落していることをドロップアウトという。ドロップアウトが発生すると、データセット内の遺伝子の発現データが不完全になり、遺伝子間の相関や差異が正確に評価されない。この問題を克服するために、統計的手法を用いて欠落したデータを推定 (imputation) または埋めることで、遺伝子発現プロファイルの正確性を向上させる必要がある。

In [46]:
import scvi

データを格納したanndataオブジェクトと、バッチのIDを指定したobsのカラムのラベルを指定して、モデルをセットアップする。カウントデータの確率モデルなので、（正規化や対数変換したマトリックスではなく）カウントデータを格納したレイヤーを指定する必要があることに注意。

In [47]:
scvi.model.SCVI.setup_anndata(
    adata,
    layer='counts',
    batch_key='batch',
)

モデルを作る。ここで中間層のノード数とか潜在ベクトルの次元サイズなど設定できるが、とくに変更しなくていいと思う。

In [48]:
model = scvi.model.SCVI(adata)

scVIモデルのトレーニング（バッチ補正）¶

学習を実行する。GPU使わない場合はめっちゃ時間かかるので今回の講習では割愛。学習済みのモデルパラメータを用意したのでそれをロードして学習できたことにする。

とはいえ、深層学習としてはたいして規模の大きいデータセットでもないので、iMac（3.6GHz 10コア Intel Core i9, メモリ128GB）でCPUのみでトレーニングしても30分程度。

In [49]:
#model.train()
In [50]:
#model.save('./models/scVI_model', overwrite=True)
In [51]:
model = scvi.model.SCVI.load('./models/scVI_model', adata=adata)

INFO     File ./models/scVI_model/model.pt already downloaded

学習し細胞ごとのバッチの効果が除去された細胞の潜在表現$z_n$は以下の関数で取得できる。PCAやUMAPの座標と同じように、obsmに格納しておく。

In [52]:
adata.obsm['X_scVI'] = model.get_latent_representation()

発現量期待値 $\rho_{ng}$ 、つまりdenoiseされた発現量テーブルは以下の関数で取得できる。適当なライブラリサイズで規格化してくれる。この数値はあとで使いたいので別のレイヤーに格納しておく。

In [53]:
adata.layers['scvi_normalized'] = model.get_normalized_expression(library_size=1e4)
In [54]:
adata
Out[54]:
AnnData object with n_obs × n_vars = 4773 × 15001
    obs: 'batch', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'n_counts', 'leiden_r0.5', '_scvi_batch', '_scvi_labels'
    var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'log1p', 'hvg', 'pca', 'batch_colors', 'neighbors', 'tsne', 'umap', 'leiden_r0.5', 'leiden_r0.5_colors', '_scvi_uuid', '_scvi_manager_uuid'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_scVI'
    varm: 'PCs'
    layers: 'counts', 'scvi_normalized'
    obsp: 'distances', 'connectivities'

scVI潜在表現+UMAPによる次元削減¶

潜在表現でちゃんとバッチ効果が補正されたかどうか確認してみる。

近傍グラフ計算の scanpy.pp.neighbors は、デフォルトではPCAで計算した主成分（adata.obsm['X_pca']）を元にグラフを構築する。ここではPCA結果の座標ではなく、scVIで推定された細胞の潜在ベクトルを指定してグラフ構築してみる。バッチ効果が補正されたグラフが構築されるはず。

In [55]:
sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.umap(adata)
sc.pl.umap(adata, color='batch')
No description has been provided for this image

scVI潜在表現+Leidenによるクラスタリング¶

同じように、scVI潜在表現で構成された近傍グラフを元にLeidenクラスタリングを実行してみる。上のセルで scanpy.pp.neighbors を実行したことにより、adata.obspにはscVI潜在表現で構築された近傍グラフが入っているので、ここでは普通に scanpy.tl.leiden を実行すればいい。

In [56]:
sc.tl.leiden(adata, key_added="leiden_scVI", resolution=1.0)
sc.pl.umap(adata,
           color=['batch', 'leiden_scVI'],
           ncols=2,
           frameon=False)
No description has been provided for this image
In [57]:
# 遺伝子 Top2a, Mcm6の比較でdrop-outが補正されたか確認。
sc.pl.scatter(adata, x='Top2a', y='Mcm6', color='leiden_scVI')
No description has been provided for this image
In [58]:
sc.pl.scatter(adata, x='Top2a', y='Mcm6', use_raw=False, layers='scvi_normalized', color='leiden_scVI')
No description has been provided for this image

Doublet検出¶

Doublet とは: 2つ以上の細胞が単一細胞としてシーケンスされている状態。

SOLOによるDoublet検出。

Bernstein, Nicholas J., et al. "Solo: doublet identification in single-cell RNA-Seq via semi-supervised deep learning." Cell systems 11.1 (2020): 95-101.

時間に余裕があればやる。

In [59]:
results = []
for batch in ['F2', 'E2']:
    tmp_solo_model = scvi.external.SOLO.from_scvi_model(
        model, restrict_to_batch=batch)
    tmp_solo_model.train()
    result = tmp_solo_model.predict(soft=False)
    # 意図は不明だがSOLO予測のindexになぜか "-0" が付加されるので消しておく
    result.index = result.index.str.replace("-0$", "", regex=True)
    results.append(result)
results = pd.concat(results)

INFO     Creating doublets, preparing SOLO model.                                                                  

GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs

Training:   0%|          | 0/400 [00:00<?, ?it/s]
Monitored metric validation_loss did not improve in the last 30 records. Best score: 0.120. Signaling Trainer to stop.
INFO     Creating doublets, preparing SOLO model.                                                                  

GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs

Training:   0%|          | 0/400 [00:00<?, ?it/s]
`Trainer.fit` stopped: `max_epochs=400` reached.
In [60]:
results[adata.obs.index]
Out[60]:
AAACCTGAGATGTCGG-1-E2    singlet
AAACCTGCAATCCAAC-1-E2    singlet
AAACCTGTCCAATGGT-1-E2    singlet
AAACGGGAGGCAATTA-1-E2    doublet
AAACGGGCACGGCGTT-1-E2    singlet
                          ...   
TTTGTCAAGGGCTCTC-1-F2    singlet
TTTGTCACATCGATGT-1-F2    singlet
TTTGTCAGTAGAGGAA-1-F2    doublet
TTTGTCAGTCATCCCT-1-F2    singlet
TTTGTCATCAGTTCGA-1-F2    singlet
Length: 4773, dtype: object
In [61]:
adata.obs['SOLO_prediction'] = results[adata.obs.index]
In [62]:
adata.obs['SOLO_prediction'].value_counts()
Out[62]:
SOLO_prediction
singlet    4426
doublet     347
Name: count, dtype: int64
In [63]:
sc.pl.umap(adata, color='SOLO_prediction', frameon=False)
No description has been provided for this image

DEG解析¶

仮説検定ではなく、ベイズ統計の枠組みでモデル比較を行う。基本的にはある集団と別の集団それぞれの遺伝子発現量期待値（VAEでdenoiseした発現量$\rho_{ng}$）のlog fold changeの値が、ある閾値以下となるか、それ以上の変化があるか、のふたつのモデルで比較する。

groupbyにラベルを指定すると、そのクラス（この場合はleidenクラスタリングのひとつのクラスタ）とその他全部で自動的に比較してくれる。比較のグループを一部に制限したり、ペアを指定したりいろいろと複雑に組み合わせは指定できる。詳細はAPIのドキュメントに。

In [64]:
DEs = model.differential_expression(groupby='leiden_scVI')
DEs.head()

DE...:   0%|          | 0/13 [00:00<?, ?it/s]
Out[64]:
proba_de proba_not_de bayes_factor scale1 scale2 pseudocounts delta lfc_mean lfc_median lfc_std ... raw_mean1 raw_mean2 non_zeros_proportion1 non_zeros_proportion2 raw_normalized_mean1 raw_normalized_mean2 is_de_fdr_0.05 comparison group1 group2
Prc1 0.9806 0.0194 3.922891 0.000584 0.000093 0.0 0.25 3.560691 3.684415 1.862878 ... 3.777792 0.459458 0.849850 0.185050 5.998869 0.816749 True 0 vs Rest 0 Rest
Pif1 0.9804 0.0196 3.912431 0.000042 0.000006 0.0 0.25 4.036641 4.243926 2.424492 ... 0.201201 0.033114 0.148649 0.025810 0.314090 0.054497 True 0 vs Rest 0 Rest
Kif23 0.9796 0.0204 3.871609 0.000229 0.000041 0.0 0.25 3.767542 3.817360 2.176383 ... 1.487993 0.194058 0.668168 0.110299 2.423138 0.367420 True 0 vs Rest 0 Rest
Sapcd2 0.9792 0.0208 3.851782 0.000012 0.000002 0.0 0.25 4.814266 4.931459 3.254749 ... 0.039039 0.007305 0.036036 0.007305 0.053764 0.014499 True 0 vs Rest 0 Rest
Ccnb1 0.9792 0.0208 3.851782 0.000462 0.000095 0.0 0.25 3.523082 3.689767 2.410103 ... 3.084096 0.465302 0.728228 0.154371 4.645929 0.852238 True 0 vs Rest 0 Rest

5 rows × 22 columns

proba_de: DEG の確率
proba_not_de: DEG ではない確率
balyes_factor: ベイズ統計の枠組みでモデル比較した時の指標のひとつ。
comparison: どのグループとどのグループで比較をした結果か。

それぞれのクラスタ「特異的」遺伝子、DEG確率の高い順に上から3つくらい見てみる。

In [65]:
markers = {}
for i, c in enumerate(adata.obs['leiden_scVI'].unique()):
    comparison_label = "{} vs Rest".format(c)
    cluster_df = DEs.loc[DEs['comparison'] == comparison_label]
    cluster_df = cluster_df[cluster_df['lfc_mean'] > 0]  # lfc_mean > 0: そのクラスター特異的に発現している遺伝子
    cluster_df = cluster_df[cluster_df['bayes_factor'] > 3]
    cluster_df = cluster_df[cluster_df['non_zeros_proportion1'] > 0.1]
    markers[c] = cluster_df.index.tolist()[:3]
In [66]:
markers
Out[66]:
{'3': ['Ccn1', 'Top2a', 'Gpx3'],
 '2': ['Nell1', 'Myo16', 'St18'],
 '0': ['Prc1', 'Pif1', 'Kif23'],
 '5': ['Kif19a', 'Egflam', 'Lhx4'],
 '8': ['Ccn1', 'Ddit4l', 'Top2a'],
 '6': ['Ccn1', 'Slfn9', 'Ccne2'],
 '11': ['Chst9', 'Pmfbp1', 'Gm32647'],
 '1': ['Thsd7b', 'D130079A08Rik', 'Cntn2'],
 '9': ['Sp8', 'Pou4f2', 'Unc5d'],
 '10': ['Lrfn5', 'Synm', 'Chl1'],
 '7': ['Gm32629', 'Krt73', 'Rgs16'],
 '4': ['Penk', 'Mybl1', 'Neurog2'],
 '12': ['Grm1', 'Fxyd7', 'Car10']}

図でLeidenクラスタのデンドログラムを表示したいので事前に計算しておく。

In [67]:
sc.tl.dendrogram(adata, groupby='leiden_scVI', use_rep='X_scVI')

クラスタ間での遺伝子発現量の比較は、scanpy.pl.dotplotが見やすくて便利。上で抽出したマーカー遺伝子を指定して、表示する発現量の数値としてはadata.rawに格納したデータで描画する。

In [68]:
sc.pl.dotplot(
    adata,
    markers,
    groupby='leiden_scVI',
    dendrogram=True,
    color_map="Blues",
    swap_axes=True,
    use_raw=True,
    standard_scale='var')
No description has been provided for this image

クラスタ平均値ではなくクラスタごと細胞ごとの全部の（scVIノーマライズされた）発現量でヒートマップを描くこともできる。

In [69]:
sc.pl.heatmap(
    adata,
    markers,
    groupby='leiden_scVI',
    layer='scvi_normalized',
    standard_scale="var",
    dendrogram=True,
    figsize=(8, 12)
)
No description has been provided for this image

データの保存¶

In [70]:
adata.write(filename='./data/retinal_1.h5ad')
In [ ]:
In [ ]: