シングルセル解析用のデータベース作成 ~SQL編~

前回の記事では、シングルセル解析の結果を格納するためのデータベース(RDB)のテーブル設計を行いました。

本記事では、実際にRDBにデータを格納して、そのデータをSQLクエリで取得するところまでを行います。


テーブル設計のおさらい


まずは、前回作成したテーブル設計図を一度確認しておきます。こちらは、使い勝手を考慮して一部変更されています。


f:id:emoriiin979:20210213105333p:plain


変更箇所としては、研究テーブルのリレーションを変更し、一つの研究で複数の検体が採用されているという形式にしました。こうすることで、例えばTabula Murisデータが複数の研究で利用されているという状況もカバーできます。

細胞データの階層構造としては、個体>組織>検体>細胞という関係となっています。新たに検体テーブルを追加したので、細胞と組織のリレーションは削除しました。一方で、動物種>組織>細胞タイプ>細胞サブタイプの階層関係も別に存在するので、組織テーブルは二つのリレーションを兼務するものとなっています。


テーブルの作成


まずは、定義したテーブルをDBMS上で作成します。今回使用したDBMSは、今後Google Colabに持っていきやすいものが良かったので、SQLiteを選択しました。

※ただし、結論を先に言ってしまうと、この選択はあまり良くありませんでした。シングルセル解析データのような大量データはGoogle Colabにアップロードしづらいので、MySQLなどでDBサーバーを立てて、そこからリモート接続してデータを取得するようにした方が良いです。


SQLite3をダウンロードして、コマンドプロンプトやターミナルなどで下記のコマンドを入力します。


> sqlite3 singlecell.db


すると、SQLiteのコマンド待ち状態となるので、そこへSQLクエリを入力していきます。


> create table projects(
>   project_id text,
>   title text,
>   primary key(project_id)
> );


MySQLなどと同じで、CREATE TABLEコマンドで新たなテーブルを作成できます。

SQLiteの注意点として、外部キーの設定ができない(できるけど面倒)という問題があります。絶対に付けたいという人でなければ問題ないかと思いますが、もし必要ならMySQLなど他のDBMSを使うと良いでしょう。


この作業を、残りのテーブルに対しても実行します。


> create table project_individuals(project_id text not null, individual_id text not null, primary key(project_id, individual_id));
> create table individuals(individual_id text not null, organism_id text not null, primary key(individual_id));
> create table samples(sample_id text not null, individual_id text not null, tissue_id text not null, primary key(sample_id));
> create table tissues(tissue_id text not null, organism_id text not null, tissue_name text not null, primary key(tissue_id));
> create table organisms(organism_id text not null, organism_name text not null, primary key(organism_id));
> create table cells(cell_id text not null, sample_id text, type_id text, subtype_id text, age int, primary key(cell_id));
> create table cell_types(type_id text not null, tissue_id text not null, type_name text not null, primary key(type_id));
> create table cell_subtypes(type_id text not null, subtype_id text not null, subtype_name text not null, primary key(type_id, subtype_id));
> create table gene_expressions(cell_id text not null, gene_id text not null, count int, primary key(cell_id, gene_id));
> create table genes(gene_id text not null, browser_id text not null, gene_name text, primary key(gene_id));
> create table gene_browsers(browser_id text not null, browser_name text, primary key(browser_id));


作成したテーブルは、.tablesコマンドで確認できます。


f:id:emoriiin979:20210212183605p:plain


CSVファイルの作成と取り込み


SQLiteでは、CSVファイル内のデータを読み込みテーブルへ格納することができます。ここからは、テーブルに取り込むためのCSVファイルをPythonで作成します。


まずは、データセットで使用されているIDや値をテーブル定義で作成したIDに変換するための辞書を作成します。


# old & animal -> individual_id
ia_dict = {
    'young': { 3: 'IDV000001', 4: 'IDV000002', 5: 'IDV000003' },
    'old': { 1: 'IDV000004', 2: 'IDV000005', 3: 'IDV000006' },
}
# replicate -> sample_id 
sr_dict = { 1: 'SMP001', 2: 'SMP002' }
# cell_type -> type_id (cell ontology)
tc_dict = {
    'leukocyte': 'CL:0000738',
    'lung endothelial cell': 'CL:1001567',
    'alveolar macrophage': 'CL:0000583',
    'myeloid cell': 'CL:0000763',
    'classical monocyte': 'CL:0000860',
    'stromal cell': 'CL:0000499',
    'B cell': 'CL:0000236',
    'T cell': 'CL:0000084',
    'non-classical monocyte': 'CL:0000875',
    'natural killer cell': 'CL:0000623',
    'mast cell': 'CL:0000097',
    'type II pneumocyte': 'CL:0002063',
    'ciliated columnar cell of tracheobronchial tree': 'CL:0002145',
}
# subtype -> type_id & subtype_id
in_dict = {
    'Npnt stromal cell': [tc_dict['stromal cell'], 'CLS001'],
    'Hhip stromal cell': [tc_dict['stromal cell'], 'CLS002'],
    'Gucy1a3 stromal cell': [tc_dict['stromal cell'], 'CLS003'],
    'Dcn stromal cell': [tc_dict['stromal cell'], 'CLS004'],
    'CD4 T cell': [tc_dict['T cell'], 'CLS001'],
    'CD8 T cell': [tc_dict['T cell'], 'CLS002'],
}


作成した辞書を使って、シングルセル解析データの中身をCSVファイルに変換します。今回使用するデータは、いつものようにAdata型データとして読み込みます。


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


次に、CSVファイルに変換するための関数を定義します。


def outCSVFile(filename, **kwargs):
    pd.DataFrame(kwargs).to_csv(filename + '.csv', header=False, index=False)


それでは、準備が完了したのでadataから必要なデータを取得・変換し、CSVファイルを作っていきます。


# projects
project_id = np.array(['PRJ000001'])
title = np.array(['Murine single-cell RNA-seq reveals cell-identity- and tissue-specific trajectories of aging'])
outCSVFile('projects', project_id=project_id, title=title)
# organisms
organism_id = np.array(['ORG000001', 'ORG000002'])
organism_name = np.array(['Homo Sapiens', 'Mus Musculus'])
outCSVFile('organisms', organism_id=organism_id, organism_name=organism_name)
# individuals
individual_id = np.array(['IDV000001', 'IDV000002', 'IDV000003', 'IDV000004', 'IDV000005', 'IDV000006'])
organism_id = np.array(['ORG000002'] * 6)
age = np.repeat(np.array(['young', 'old']), 3)
outCSVFile('individuals', individual_id=individual_id, organism_id=organism_id, age=age)
# project_individuals
project_id2 = np.array(['PRJ000001'] * 6)
individual_id2 = individual_id
outCSVFile('project_individuals', project_id=project_id2, individual_id=individual_id2)
# tissues
tissue_id = np.array(['TSS000001'])
organism_id = np.array(['ORG000002'])
tissue_name = np.array(['Lung'])
outCSVFile('tissues', tissue_id=tissue_id, organism_id=organism_id, tissue_name=tissue_name)
# samples
individual_id2 = np.repeat(individual_id, 2)[:-1]
sample_id = np.tile(np.array(['SMP001', 'SMP002']), 6)[:-1]
tissue_id2 = np.array(['TSS000001'] * 11)
outCSVFile('samples', individual_id=individual_id2, sample_id=sample_id, tissue_id=tissue_id2)
# gene_browsers
browser_id = np.array(['BRW000001'])
browser_name = np.array(['Ensembl'])
outCSVFile('gene_browsers', browser_id=browser_id, browser_name=browser_name)
# genes
gene_id = np.array(adata.var['gene_ids-0'].tolist())
browser_id2 = np.array(['BRW000001'] * len(gene_id))
gene_symbol = np.array(adata.var.index.tolist())
outCSVFile('genes', gene_id=gene_id, browser_id=browser_id2, gene_symbol=gene_symbol)
# cell_types
type_name = np.array(adata.obs['cell_type'].unique())
tissue_id2 = np.array(['TSS000001'] * len(type_name))
type_id = np.array([tc_dict[x] for x in type_name])
outCSVFile('cell_types', type_id=type_id, tissue_id=tissue_id2, type_name=type_name)
# cell_subtypes
d = obs[['cell_type', 'subtype']]
subtype_name = np.array(d[~d.duplicated()][d['cell_type'].cat.set_categories(d['subtype'].cat.categories) != d['subtype']]['subtype'].tolist())
type_id2 = np.array([in_dict[x][0] for x in subtype_name])
subtype_id = np.array([in_dict[x][1] for x in subtype_name])
outCSVFile('cell_subtypes', type_id=type_id2, subtype_id=subtype_id, subtype_name=subtype_name)
# cells
cell_id = np.array(['CLL{:0=12}'.format(i+1) for i in range(len(obs))])
individual_id2 = np.array([ia_dict[age][animal] for age, animal in zip(obs['age'], obs['animal'])])
sample_id2 = np.array([sr_dict[rep] for rep in obs['replicate']])
type_id2 = np.array([tc_dict[cell_type] for cell_type in obs['cell_type']])
subtype_id2 = np.array([in_dict[x][1] if x in in_dict else '' for x in obs['subtype']])
outCSVFile('cells', cell_id=cell_id, individual_id=individual_id2, sample_id=sample_id2, type_id=type_id2, subtype_id=subtype_id2)
# gene_expressions
df_X = pd.DataFrame(adata.X.todense())
df_X.columns = gene_id
df_X['cell_id'] = cell_id
df_X.melt(id_vars='cell_id', var_name='gene_id', value_name='count').to_csv('gene_expressions.csv', header=False, index=False)


これで、全てのCSVファイルの作成が完了しました。

これらのファイルを、.importコマンドを使って取り込みます。


.mode csv
.import cell_subtypes.csv cell_subtypes
.import cell_types.csv cell_types
.import cells.csv cells
.import gene_browsers.csv gene_browsers
.import gene_expressions.csv gene_expressions
.import genes.csv genes
.import individuals.csv individuals
.import organisms.csv organisms
.import project_individuals.csv project_individuals
.import projects.csv projects
.import samples.csv samples
.import tissues.csv tissues


データの取得


これで、全てのデータがテーブルに格納されました。最後にらSELECT文で必要なデータを取得できるか確認します。


select
  ct.type_name
  , g.gene_symbol
  , i.age
  , count(*) as n_data
from
  gene_expressions ge
  left outer join genes g
    on g.gene_id = ge.gene_id
  left outer join cells c
    on c.cell_id = ge.cell_id
  left outer join cell_types ct
    on ct.type_id = c.type_id
  left outer join samples s
    on s.individual_id = c.individual_id
    and s.sample_id = c.sample_id
  left outer join individuals i
    on i.individual_id = s.individual_id
where
  ge.gene_id = 'ENSMUSG00000025902'
  and ct.type_name = 'T cell'
group by
  ct.type_name
  , g.gene_symbol
  , i.age;


f:id:emoriiin979:20210213134734p:plain


このように、ファイルからでなくデータベースからデータの取得ができることが確認できました。あとは、使用しているプログラミング言語に応じて必要な処理を導入すれば完了です。


まとめ


本記事では、前回作成したテーブル設計図をもとに、SQLiteでデータベースを作成しました。

RDBにデータを格納しておけば大量データをメモリ外で捌くことができるようになりますし、複数人での運用も可能となるので、試せる人は是非試して欲しいと思います。

今回はSQL系のデータベースシステムを扱いましたが、他にもNoSQL系やHadoopのような分散システムなどが存在するので、それらについても調べてみると面白いかもしれません。


以上です。

シングルセル解析用のデータベース作成 ~テーブル定義編~

データサイエンティストに求められるスキルとして、データベースを活用したデータハンドリング技術が挙げられます。これはビッグデータ解析を行う上ではExcelでは行数不足が起こり得るということもありますし、何より全てのデータがメモリに収まりきらないことがあるので、データベースサーバーからデータを逐次取得する技術が必要だからだと考えられます。

これまでの老化解析で使ってきたシングルセルデータも、数が増えてくればH5ADファイル単位でのデータ管理では回らなくなってしまいます。これをデータベース化することができれば、データ自体はストレージに格納しつつ、クエリを入力することで必要なデータを検索してメモリに送り込むことができます。


本記事では、シングルセル解析のデータを格納するためのデータベース設計を行います。データベースシステムは色々存在しますが、今回は一番有名なリレーショナルデータベース(RDB)方式で設計を行いました。この記事のようにテーブル定義すれば、おそらく様々なシングルセルデータを格納できるのではないかと思います。

続きを読む

ツールをWeb化するためのサーバー探し!!

学術論文では、バイオインフォマティクスなどの解析ツールはPythonやRなどで作成され、そのソースコードがそのまま公開されることが多いです。このとき、自分がPythonなどを扱えるなら問題ないのですが、これまでプログラミング言語に触れたことが無い人にとっては、これらのツールを使うハードルが高すぎます。

一般的に、普段実験に携わっている人が、新たにプログラミングを習得して得られるメリットはそれほど多くありません。一方で、それらの解析ツール自体は自身の研究に必要であることが多いです。したがって、プログラミングを新たに習得することなく、Pythonなどで実装されたツールが利用できることが理想といえます。


そのための解決策の一つとして、「ツールのWeb化」があります。Webサービスは画面操作で実行できるため、プログラミングを知らない人でも扱うことができます。例えば、PythonでできたツールをDjangoなどのWebフレームワーク内で実行するようにすれば、ツールを実行するためのラッパーを簡単に作れそうです。

しかし、そのためにはWeb化したツールを置くためのWebサーバーが必要となります。そこで、Webサーバーの構築に使えそうなレンタルサーバーサービスの候補をいくつか調べてきました。本記事では、3つのレンタルサーバーサービスを紹介し、どれが今回の用途に適していそうかを検討します。

続きを読む

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

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

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


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

続きを読む

Matplotlibはややこしい

Pythonでグラフを描くなら、Matplotlibの使用は避けて通れません。「Python グラフ」などで調べると、ほぼ確実にplt.plotなどMatplotlibの関数を見ることになるでしょう。

しかし、いざMatplotlibを使ってみようとしても、思ったように動作しないことが多いです。ネットからサンプルコードを探して、同じようにスクリプトを書いてみても、変なグラフができてしまったり、グラフが表示されないなどのトラブルに見舞われることがあります。


いったいなぜ、このような事態が起こるのでしょうか。それはもしかすると、Matplotlibに潜む「ややこしさ」が原因かもしれません

本記事では、Matplotlibが持つ3つの「ややこしさ」について紹介します。Jupyter notebookでMatplotlibを動かすにあたって、これらを注意すれば上記のトラブルに遭遇する確率は減るかもしれません…。

続きを読む

Seleniumではてなブログのアクセス数を取得する

ブログがどれくらい読まれているかは、はてなブログの場合はアクセス解析ページを見ることで、確認できます。このアクセス解析ページでは、日ごとのアクセス数などがグラフ化されているので、何か分析ツールを用意しなくてもある程度の視覚化は済んでいます。

しかし、より複雑な分析をしたいという人もいるでしょう。その場合は、アクセス解析ページに載っているデータを収集して、アクセス情報のデータセットを用意する必要があります。


ここで、TwitterアナリティクスのようにCSVデータが用意されているなら簡単なのですが、残念ながらCSVデータが提供されていないアナリティクスの方が多数派のようです。はてなブログアクセス解析ページでもそのようなデータセットは用意されていませんでした。


それでは、毎日アクセス解析ページを確認して、必要なデータをコピペしていかなければならないのでしょうか。それを毎日やるのはとても大変な作業です。手作業でのデータ収集では、うっかり間違ったデータを取ってきてしまうこともあり得ます。

そこで、Seleniumを使ってデータ収集の自動化を図りますSeleniumWebサービス画面の自動テストツールですが、はてなブログアクセス解析データを取得するという、スクレイピングのような使い方もできます。本記事では、Seleniumの使い方を紹介して、実際にアクセス情報を取得してみます。

続きを読む

細胞種推定のための機械学習入門

若年組織と老化組織の違いとして、「細胞種比率の違い」があります。例えば、老化によって免疫細胞の割合が増減していたりすることがあります。

このような、細胞種の割合を調べるためには、当然ですが全てのデータで細胞種を同定しておく必要があります。データごとの細胞種が分かれば、あとはこれらの細胞数を集計することで、細胞種比率を算出することができます。


しかし、現状のRNA-Seq解析(細胞のmRNA全量を調べる解析のことです)では、解析している細胞の種類を特定することができません。そもそも、未発見の細胞が未だ多いこともあるため、100%の精度で細胞種を決定するのは不可能だと考えた方がいいでしょう。


そこで、これまでに様々な細胞種「推定」方法が開発されてきました。本記事では、これらの細胞種推定方法のうち、「深層学習」を活用したものを紹介します

続きを読む