Pandasで細胞数をカウントし、Matplotlibでグラフにする

老化組織と若年組織の違いとして、細胞数の比率の違いが挙げられます。例えば、老年マウスでは加齢に伴い脂肪組織マクロファージのサブタイプ組成が変化することが知られており、炎症を誘発することが分かっています。このように、老人と若者の生物学的違いを認識するうえで、身体組織の細胞組成を調べることはとても有効です。

しかし、細胞数を目視でひとつひとつ計数するのはとても手間がかかります。細胞数が数万あることも珍しくはありませんし、途中でカウントミスが発生しても気づくすべがありません。なので、これほど多くの細胞をカウントするにはプログラミングが不可欠です。

本記事では、Scanpyを用いて得られた老化・若年組織のシングルセルRNA-SeqデータをPandasのDataFrame型として読み込み、既存の'cell_type'情報を使って細胞数のカウントを行います。これらの作業は、ほぼ全てPandasとMatplotlibで用意されているメソッドで実現させました。ここでは、その方法について共有します。

データの準備

まずは、いつものようにシングルセルRNA-Seqデータを読み込みます。実行環境も、いつも通りGoogle Colabで行いました。

import os
os.system('wget https://storage.googleapis.com/calico-website-mca-storage/lung.h5ad')

import scanpy as sc
adata = sc.read_h5ad('lung.h5ad')
sc.pp.normalize_total(adata, target_sum=1e6)  # CPM
sc.pp.log1p(adata)  # 対数変換

このデータには、20755個の細胞に対して、細胞種が同定されています。この情報は、anndata.obsに格納されています。

本記事の目標

ところで、今使っているRNA-Seqには出典となる論文が存在します。こちらの論文でも同じデータを使っていて、老年マウスと若年マウスの細胞学的違いの調査結果が記載されています。

その中の一つとして、老年マウスと若年マウスの細胞比率の違いを示すグラフがあります。

J.C. Kimmel et al., 2019, Fig 1D

本記事では、肺（Lung）データの細胞数カウントを実施し、元論文のグラフ（真ん中）の再現を目指します。

Pandasで細胞数を数える

さっそく、Pandasで細胞数を数えてみます。

実は、細胞数の計数方法自体は論文著者が公開しているため、ほとんどはここのソースコード通りに実装しています。

state_var = 'subtype'
obs = adata.obs.copy()

cell_states = obs.groupby(
    ['age', state_var, 'animal', 'replicate']).count().loc[('old'), 'batch'].index.levels[0]
n_young = np.array(obs.groupby(
    ['age', state_var, 'animal', 'replicate']).count().loc['young', 'batch'])
n_old = np.array(obs.groupby(
    ['age', state_var, 'animal', 'replicate']).count().loc['old', 'batch'])
n_young_sum_replicate = np.concatenate([np.array(obs.groupby(
    ['age', 'animal', 'replicate']).count().loc['young', 'batch'])] * len(cell_states))
n_old_sum_replicate = np.concatenate([np.array(obs.groupby(
    ['age', 'animal', 'replicate']).count().loc['old', 'batch'])] * len(cell_states))

import pandas as pd
df = pd.DataFrame({
    'Cell State': np.concatenate([np.repeat(cell_states, len(np.unique(obs['animal']))*len(np.unique(obs['replicate'])))] * 2),
    'Number': np.concatenate([n_young, n_old]),
    'Proportion': np.concatenate([n_young/n_young_sum_replicate, n_old/n_old_sum_replicate]),
    'Age': pd.Categorical(['Young']*len(n_young)+['Old']*len(n_old))
})
df = df.loc[~np.isnan(df['Number']), :]  # 存在しないanimal idとreplicate idを削除

複雑そうに見えますが、やっていることはgroupbyとcountメソッドを使って「Old/Young」「細胞種」「サンプル」ごとに細胞数を数えているだけです。細胞数（Number）の他に、細胞比率（Proportion）を計算するために、全ての細胞種の合計数（例えばn_old_sum_replicate）もカウントしています。

なお、このデータでは若年・老年マウスが3匹ずつ使用されており、それぞれのマウスで2つのサンプルを採取しているため、合わせて12個のサンプルがあります。（ただし、老年マウスのサンプル一つは品質が低いと判断されたため、データセットからは除外されています。）

今回は元論文のグラフを再現するために、細胞種としては'subtype'を選択しました。肺のサブタイプは全部で17種類あるため、最終的に作成されるDataFrameは17細胞種×11サンプル=187行となります。

Matplotlibで有意差付き棒グラフを作る

カウントが完了したので、この結果を棒グラフとして描画します。上記のPandasは論文著者のGitHubを漁るだけでしたが、作図方法は1から調べたので、正直こちらが一番苦労しました。

作図のためのMatplotlibコードについては、Matplotlib公式が提供しているサンプルコードに基づいて作成しました。

まずは、作図するためのデータを整形します。

import numpy as np

labels = df['Cell State'].unique()
x = np.arange(len(labels))
height = 0.4

y_old = df[df['Age']=='Old'].groupby('Cell State').mean()['Proportion']
y_young = df[df['Age']=='Young'].groupby('Cell State').mean()['Proportion']

次に、t検定を使って老化細胞・若年細胞で細胞比率に有意差があるかどうかを調べます。こちらは元論文のグラフには描画されていませんが、論文内でt検定が実行されており、その結果をもとに細胞比率の考察が展開されているので、同じようにt検定を行います。

ちなみに、細胞比率（Proportion）をそのままt検定に適用するのは非推奨です。詳しくはこの論文やこちらの解説記事を参照してほしいのですが、簡単に説明すると、細胞比率のように「全ての値を足すと100%」になるという制約をもつ組成データには、t検定のような実空間データ向け解析方法は適用するべきではないというものです。

元論文では、相加対数比変換（ALR）を使って組成データを実空間データに変換することで、t検定が適用できるようにしています。ただ、いろいろ調べてみましたが、ALRの実装方法がよく分からず、これを実現するパッケージも見つからなかったので、本記事では代わりに有心対数比変換（CLR）を使って再現しました。

from skbio.stats.composition import clr
from scipy import stats

df['Prop CLR'] = clr(df['Proportion'])
p_vals = df.groupby('Cell State').apply(lambda df: stats.ttest_ind(
    df[df['Age']=='Young']['Prop CLR'],
    df[df['Age']=='Old']['Prop CLR']
)[1])

これで全てのデータが用意できたので、Matplotlibで棒グラフを作成します。今回は細胞種をラベルとして使いたいので、元論文と同じく横向きの棒グラフを作ります。

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10,8))
rects1 = ax.barh(x - height/2, y_old, height, label='Old')
rects2 = ax.barh(x + height/2, y_young, height, label='Young')
ax.invert_yaxis()

ax.set_yticks(x)
ax.set_yticklabels(labels)
ax.legend()

for i, (rect1, rect2) in enumerate(zip(rects1, rects2)):
    width = max(rect1.get_width(), rect2.get_width())
    
    text = ''
    if p_vals[i] < 0.05:
        text += '*'
    if p_vals[i] < 0.01:
        text += '*'

    ax.annotate(
       '{}'.format(text),
        xy=(width + 0.01, rect2.get_y() + rect2.get_height() / 2)
    )

    if text != '':
        ax.plot(
            [width + 0.004, width + 0.006, width + 0.006, width + 0.004],
            [rect1.get_y(), rect1.get_y(), rect2.get_y() + rect2.get_height(), rect2.get_y() + rect2.get_height()],
            c = 'black',
            linewidth = 0.5
        )

plt.show()

というわけで、上のコードを実行すると図左のようなグラフが作成されます。図右に元論文のグラフを置きましたが、元論文と同じ横向き棒グラフが作成できているのが分かると思います。

よって、上記のコードを実行すれば、同じデータから全く同じグラフを作成できることが証明されました。

まとめ

本記事では、PandasとMatplotlibを使って論文内の細胞比率グラフの再現に挑戦しました。

冒頭で述べた通り、細胞比率を比較するのは老化研究において有効な手段の一つです。この調査が、本記事のコードを実行することで計数から図示まで達成できることを示しました。

anndataの中身で使用したのは、RNA-Seqデータそのものといえるanndata.Xと、細胞種が定義されているanndata.obs['subtype']だけです。仮に細胞種データが無かったとしても、ScNymなどを使えば、細胞種を予測することができます。

もし、解析に使用できそうなデータが手元にあるなら、一度試してみてみるのもよいと思います。ぜひ、本記事の内容を活用してみてください。

以上です。