この記事は
Blue Collar Bioinformatics; Automated retrieval of expression data with Python and Rの翻訳です。
多少古いですが、自身の勉強もかねて。
誤字、誤訳の指摘大歓迎です。
以下翻訳
この2月、バイオロジーデータの統合を主題としたバイオハッカソン2010に参加するため日本に行きます。で、その準備としてユースケースを考えていました。
あるデータセットと格闘しているときに、どのような追加のデータがあれば問題の解決に役立つでしょうか？典型的な例として、遺伝子の発現量データの情報があります。興味のある遺伝子のリストが手元にあったとして、その発現量を知りたいという場合です。
セマンティックwebの夢をかなえた世界でならば、以下のようなRDFトリプルのセットにクエリを投げることになると思います。
- 発現量セット -> has 生物種 -> Mus musculus
- Mus musculus -> has 細胞腫 -> proB
- proB -> has biological state -> wild type
- 発現量セット -> has 識別子 -> a RefGene ID
- a RefGene ID -> has 発現量 -> 7.65
すると、興味のある細胞腫(proB)の通常状態での発現量が手にはいり、次いで手元の遺伝子の発現量が取得できます。理想としては発現量セットにメタデータがついていて欲しいところです。そうすれば7.65という数値の意味するものが実験のコンテキストからわかりますからね。
セマンティックウェブのチャンネルからはこのクエリを投げることはできませんが(訳注: 2010年の時点で生物学データのRDF化が十分でなかったため。今はわからない。)、同じような問題を自動で解くための道具は豊富に存在します。
NCBIの Gene Expression Omnibus(GEO)は発現量データをオンラインでホスティングしています。データセットへのクエリはBiopythonのEutilsAPIで可能です。
興味のあるデータセットが特定できたならば、BioconductorのGEOQueryでデータを取得し、UCSC RefGeneの識別子と対応させることができます。これらはすべてRpy2でpythonからRを実行することで連結させることができます。
訳注: 今なら多分
pipeRの方が良い
まずはトップレベルの目標を定めましょう。ここではマウスのproB細胞の野生型の発現量データを求めることです。
def main():
organism = "Mus musculus"
cell_types = ["proB", "ProB", "pro-B"]
email = "chapmanb@50mail.com"
save_dir = os.getcwd()
exp_data = get_geo_data(organism, cell_types, email, save_dir)
def _is_wild_type(result):
"""サンプルが野生型であるか否かをタイトルから判断"""
return result.samples[0][0].startswith("WT")
次いで、生物種と細胞腫に基づいて、入手可能なデータを取得するためのクエリをGEOに投げます。実際の仕事では、あなたのプライベートな_is_wild_type関数を、対話的な機能を果たすものに置き換え、ユーザーである生物科学者に実験手法に応じて選ばせるものとなるかもしれません。
最終的な結果をpickle化して手元に保存していることにも注目してください。こうして1度だけクエリを行うことで外部サーバーの負荷を軽減しています。
def get_geo_data(organism, cell_types, email, save_dir, is_desired_result):
save_file = os.path.join(save_dir, "%s-results.pkl" % cell_types[0])
if not os.path.exists(save_file):
results = cell_type_gsms(organism, cell_types, email)
for result in results:
if is_desired_result(result):
with open(save_file, "w") as out_handle:
cPickle.dump(result, out_handle)
break
with open(save_file) as save_handle:
result = cPickle.load(save_handle)
GEOにクエリを送信して結果を取得するという難題はBiopythonのEntrezインターフェイスによって行われています。クエリをくみ上げた後、結果は説明付きのシンプルなオブジェクトにパースされ、タイトルとペアになった発現量とサンプルごとのIDを保持します。
def cell_type_gsms(organism, cell_types, email):
"""Entrezを使用して特定の生物種と細胞腫に関してGEOを取得する。
"""
Entrez.email = email
search_term = "%s[ORGN] %s" % (organism, " OR ".join(cell_types))
print "Searching GEO and retrieving results: %s" % search_term
hits = []
handle = Entrez.esearch(db="gds", term=search_term)
results = Entrez.read(handle)
for geo_id in results['IdList']:
handle = Entrez.esummary(db="gds", id=geo_id)
summary = Entrez.read(handle)
samples = []
for sample in summary[0]['Samples']:
for cell_type in cell_types:
if sample['Title'].find(cell_type) >= 0:
samples.append((sample['Title'], sample['Accession']))
break
if len(samples) > 0:
hits.append(GEOResult(summary[0]['summary'], samples))
return hits
これで、実験の種類を選べば、第一関門突破です。
次は、発現量データの値を取得します。これはRefGene IDを発現量にマップした辞書を作成することで行われます。つまり、結果は以下のようになります。
WT ProB, biological rep1
[('NM_177327', [7.7430266269999999, 6.4795213670000003, 8.8766985500000004]),
('NM_001008785', [7.4671954649999996, 5.4761453329999998]),
('NM_177325', [7.3312364040000002, 11.831475960000001]),
('NM_177323', [6.9779868059999997, 6.3926399939999996]),
('NM_177322', [5.0833683379999997])]
実際に面倒なところは、発現量の取得と遺伝子IDとのマッピングですが、これはBioconductorのGEOqueryライブラリが行ってくれます。遺伝子と発現量の対応表、およびより高次のメタデータはローカルのファイル中に保存されます。前者はRスタイルのタブ区切りテーブルで、後者はJSONで保持されています。両方ともローカルに読みやすい形で保存しましょう。
def write_gsm_map(self, gsm_id, meta_file, table_file):
"""Bioconductorを用いてGEOの発現量データを取得し、表に保存する。
"""
robjects.r.assign("gsm.id", gsm_id)
robjects.r.assign("table.file", table_file)
robjects.r.assign("meta.file", meta_file)
robjects.r('''
library(GEOquery)
library(rjson)
gsm <- getGEO(gsm.id)
write.table(Table(gsm), file = table.file, sep = "\t", row.names = FALSE,
col.names = TRUE)
cat(toJSON(Meta(gsm)), file = meta.file)
''')
この関数によってプローブ名と発現量との対応表ができます。より便利なものとするため、expression マイクロアレイのプローブを生物学的な遺伝子のIDに対応させる必要があります。これはGEOライブラリをもう一度呼び出すことによって行われます。データは再びRのデータフレームの形で取得し、興味のある要素のみ、タブ区切りのファイルで保存します。
def _write_gpl_map(self, gpl_id, gpl_file):
"""Rを用いてGEOプラットフォームのデータを取得して表形式で保存する"""
robjects.r.assign("gpl.id", gpl_id)
robjects.r.assign("gpl.file", gpl_file)
robjects.r('''
library(GEOquery)
gpl <- getGEO(gpl.id)
gpl.map <- subset(Table(gpl), select=c("ID", "RefSeq.Transcript.ID"))
write.table(gpl.map, file = gpl.file, sep = "\t", row.names = FALSE,
col.names = TRUE)
''')
最後に、今までの処理をまとめます。それぞれの対応表ファイルをダウンロード、パースし、発現量 -> プローブID -> Transcript ID と対応関係をつなげていきます。最終的な辞書は遺伝子のIDを発現量に対応付けたものになるので、これを返します。
def get_gsm_tx_values(self, gsm_id, save_dir):
"""「転写産物」 => 「GEO, GSMファイルから取得した発現量」の対応表を作成する
"""
gsm_meta_file = os.path.join(save_dir, "%s-meta.txt" % gsm_id)
gsm_table_file = os.path.join(save_dir, "%s-table.txt" % gsm_id)
if (not os.path.exists(gsm_meta_file) or \
not os.path.exists(gsm_table_file)):
self._write_gsm_map(gsm_id, gsm_meta_file, gsm_table_file)
with open(gsm_meta_file) as in_handle:
gsm_meta = json.load(in_handle)
id_to_tx = self.get_gpl_map(gsm_meta['platform_id'], save_dir)
tx_to_vals = collections.defaultdict(list)
with open(gsm_table_file) as in_handle:
reader = csv.reader(in_handle, dialect='excel-tab')
reader.next() # header
for probe_id, probe_val in reader:
for tx_id in id_to_tx.get(probe_id, []):
tx_to_vals[tx_id].append(float(probe_val))
return tx_to_vals
完全なスクリプトは今までのスニペットを統合したものになります。クエリからはじめて、最終的には必要な発現量データを取得するところまで実行しました。
とはいえ、RDFを用いた理想的なものに比べるとだいぶ回り道になってしまいました。このように、疑問ドリブンなデータへのアクセスをできる限りシンプルなものにしたい、というのはバイオインフォマティクスのツールを作成する者が繰り返し直面する課題です。
訳はここまで
スクリプト自体はhttps://github.com/chapmanb/bcbb/blob/022c50b3132bc3d009ed1d533a3954bd564a9ed3/rest_apis/find_geo_data.py
にあります。
筆者のbcbbというリポジトリには他にもいろいろ役立つバイオインフォマティクス用のプログラムがあります。