#ajacs

NBDCと東大農アグリバイオインフォマティクス教育研究ユニットによる「次世代シークエンサ(NGS)ハンズオン講習会」にアシスタントとして参加してきました。

http://biosciencedbc.jp/human/human-resources/workshop/h27

BioLinuxのインストールにはじまり、Linuxでのコマンドラインの使い方、スクリプト言語やNGS解析の基礎、代表的なNGSアプリケーションについてそれぞれの解析方法など、約2週間かけて行われるハードコアな講習会です。ハードコアですが、普段からデータ解析にチャレンジしてみたい、でもどこから取りかかればよいか分からないという人にとっては、えいやと時間を取ってみっちりとトレーニングを受ける、絶好の機会ではないかと思います。「全日程参加できる方から優先的に参加を受け付ける」としていたこともあり、他のトレーニングコースに比べてモチベーションのより高い参加者の方が多い印象があります。

しかしこのようなハードな講習会に来られている方の多くは、普段データ解析をするインフォマティシャンやそのテクニックに触れる機会が少ないのではないかと推察されます。そうであればなおのこと、折角時間を取って覚えたことも使わないと忘れてしまいますし、本当に大変なのはトレーニングを受けたあと、いかに学んだことを日々の暮らしに組み込んでいけるかという点だと思います。

その点で、アシスタントといえど軽い気持ちでサポートをするわけにはいきません。たとえば「講師の言う通りにコマンドを打ったけれども動きません」とヘルプを求められたとき「こうすれば動きますよ」と単に直してあげるだけなら簡単です。しかし「なぜ動かなかったのかを推測する」「推測される原因から解決策を探す」その方法を含めて示さないことには、別の場面で応用ができませんし、そうなってしまうと折角のトレーニングの機会が無駄になってしまいます。自分のラボに戻れば講師もアシスタントもいないわけですから、2週間のトレーニングコースを普段からコンピュータを使ってデータを解析している人の習慣や考え方に気付くチャンスとして活用してもらうことで、長く役に立つスキルを持って帰ってもらえるといいのかもしれない、アシスタントをしながらそんなことを考えていました。一言「エラーメッセージをそのままGoogleにコピペして上から全部ドキュメントを読め」というだけといえば、まあそうなんですが。

アシスタントが大変なら講師はもっと大変で、大勢の人を相手にハンズオンをやるというのは、僕自身の経験からも、それはもう本当に大変です。しかし僕が参加したアメリエフ株式会社 服部さんによるPythonのコースと理研森岡さんによるChIP-Seqデータ解析のコースはどちらも充実したよい内容でした。参加者の多いハンズオンは一度に多くの人に対してコースを提供できるメリットと同時に、参加者の興味や技術レベル、リテラシレベルの分散によってはいまいち響かないものになってしまうリスクもあるのですが、お二人のコースに関してはとてもうまくいっていたと思います。興味のある方は、 こちらのページ に資料がアップされていますので、是非チェックしてみてください。

講師のテンションがおかしい

さて久々にエントリを書いたかと思えばお利口さん気取りの感想文だけを置きにきたわけではもちろんありません。ではなんなのかというと、まずは @suimye こと理研森岡さんによる愉快なトークでお届けしたChIP-Seqデータ解析コース、その講習用ページをご覧ください。https://github.com/suimye/NGS_handson2015

ChIP-Seqの原理からはじまりプロットまで、講習の流れが示されているところをざっと見て、ページ一番下の「中級者以上コース（一匹狼たちへの課題）」というリンクに注目してください。括弧内の文言の意味がよくわかりませんがクリックするとこちらのページに飛びます。ページを開いてみてください。https://github.com/suimye/NGS_handson2015/wiki/NGS_senior

控えめに言ってもファンキー過ぎるページですが、僕はこういうの好きなので気にしません。いずれにせよ、1日のトレーニングは朝から夕までおおむねこのテンションで続きました。参加者の方々は始終「こんなときどんな顔をすればいいか分からない」というような表情でいらっしゃったのが傍から見ている分には大変愉快でした。笑えばいいと思うよ。

ともあれ、ウェブページです。何かヤバいクスリでもキメてるんじゃないかというテンションのあたりはとりあえず無視して見ていきますと「問題」という欄があります。

ファイル kaeru.nazo.fastq.gzは、ヒトの細胞を由来とするChIP-seqデータである。実習および講義で教わる操作および独自のアイディアで、このChIP-seqデータの転写因子名を調べよ

つまりChIP-Seqのトレーニングコースは自分でさっさとやってしまった！簡単過ぎるぜ、もっとタフなタスクをよこせ！というマッドな参加者を想定してこのような課題を用意してくださったようです。気が利くんだか、カエルの写真も脈絡がないし単に貼りたかっただけなんじゃねえーかとか、まあ、いいんですけど、こういうの、みなさん好きですよね？

やってみよう

というわけで、この課題にチャレンジしてみます。「素敵なプレゼント」とか、嫌な予感しかしないですけど。どうせカエルのぬいぐるみです。

正攻法でいくと、fqをヒトゲノムにマップしてピークコール、結合領域のモチーフから転写因子を探るなどのやり方がありますが、それだとあまり面白くないので、データベース屋として、データベースだけを駆使して転写因子を当ててみましょう。

まず、fastqをダウンロードして中身を見てみます。gzipで圧縮されていますが、気にせずlessで中身を見ます。

$ less kaeru_nazo.fastq.gz @KAERU.13652678 HWUSI-EAS366_145:6:26:14947:20927/1 TGGCTGGCAACAATAGATACTGGGGACTACTAGACA + #################################### @KAERU.46814293 HWUSI-EAS366_145:6:95:12328:13457/1 AATGGAATTGAATGGAATGGAATTGAATGGAATGGG + ,7>>21?############################# @KAERU.9734561 HWUSI-EAS366_145:6:19:12517:8892/1 TTGAGATGGAGTCTTGCTCTGTCGCCCAGGCTGGAG

おそらく答えとなる転写因子はあまりマニアックすぎないもので、クセのないデータを用意しているのでは、そして自身の研究データはコラボレータとの同意などの手間を考えると使いづらいので、公共DBで公開されたデータを使っているのではないか、公開データだとしたらfastqのリードIDの文字列中にSRAのIDが入っているのではないか…という予想からファイルの中身を見てみましたが、なんか KAERU とか書いてある。芸が細かいというかなんというか。

SRA IDがあればそのIDで検索すればすぐ分かるのですが、今回はそれを許してはくれないようです。残念。しかし一歩前進です。@で始まるID列、”HWUSI..."から始まるこの部分はシーケンサの機械固有のID、レーン番号、タイル番号、クラスタの座標と続く、イルミナのシーケンサから出力されたものだと思われます。そしてこの機種名の前に文字列が入っており、そこをKAERUとわざわざ書き換えているということは、これは元々SRA Run ID (SRR ID)があったものと考えられます。

IDはありませんでしたが、諦めるのはまだ早い。シーケンサの機械固有のID、これをキーにして公共DBから同じ機械によって出された配列データをリストアップすれば何か分かるかもしれません。

このようなシチュエーションにはこれまで遭遇したことがなかったので、各公開データとシーケンサの機械固有IDを紐付けたDBは残念ながら手元にありません。イルミナのシーケンサから得られた全ての公開fastqデータからシーケンサの固有IDを攫うにしても、なんせHiSeq2000だけで56万実験セットとかあるので、ちょっと時間がかかりそう。まあ、ダメ元でgoogleに突っ込んでみましょう。

なんか出た!!

2件ですがヒットしました。どうやら同じことを考えたのか、 @ma_ko というインターネットユーザがこのファイルを見てうんざりしたらしいという情報が2件目に見つかったのも収穫ですが、1件目にヒットしたこれはまさしく公開されたfastqのリードIDの文字列です。しかもURLから見るに、今回と同じNGS講習会の、昨年のテストデータのようです。これは怪しい。「実はこのデータは昨年使ったものでした〜！ワハハ！」みたいな展開、めっちゃありそう。

しかし、検索にヒットしたものはレーン番号が3番で、課題のfqは6番でした。わざわざレーン番号に細工をするというのはちょっと考えにくいので、「同じシーケンサから出た別のデータ」を疑って調べていきます。プロジェクト単位で見るために、この SRR445816 っていうやつを http://sra.dbcls.jp/search で見てみましょう。

手前味噌なんですけど、これはうちで作ってる公開データのレポートを出力する検索システムです。SRAの新しいスキーマにまだ対応していなくて、しばらく新しいデータが更新されておらず、現在システムのアップデート中なのですが。ともあれ、試しに検索してみましょう。

なんかそれっぽいの出た！！

転写因子 Oct4, Sox2, Klf4, c-Myc をゲノム上にマップするプロジェクトにおける一連のシーケンスデータのセットのようです。これは正解に近づいた予感がする。講師も「論文のピークを信じるな、ピークコールのスレッショルドは目的に応じて恣意的に決められるのだから自分の研究に利用したければ元データを再計算しろ」と言っていたし、この論文の再解析を行った可能性はありそう。

同じプロジェクトで得られた Run は全部で7つ。最初のSRR445816がgoogleでヒットした、昨年の講習会でテストデータに使われたデータですね。この下に並んでる中からレーン番号6番でシーケンスされたデータが見つかればそれが一番怪しい。 “FTP” をクリックしてDDBJのFTPに接続して、データをダウンロードして、それぞれのデータを順番に見てみます。

$ less SRR445817.fastq.bz2 @SRR445817.1114 HWUSI-EAS366_145:4:1:5805:1023 length=36 NNNNNNNNNTNNNNNNNNNNNNNNNCNGNNNNNNAN +SRR445817.1114 HWUSI-EAS366_145:4:1:5805:1023 length=36 #################################### @SRR445817.1117 HWUSI-EAS366_145:4:1:5974:1021 length=36 NNNNNNNNNANNNNNNNNNNNNNNNANCNNNNNNAN +SRR445817.1117 HWUSI-EAS366_145:4:1:5974:1021 length=36 #################################### @SRR445817.1118 HWUSI-EAS366_145:4:1:6008:1017 length=36 NNNNNNNNNCNNNNNNNNNNNNNNNTNTNNNNNNCN

4番レーン。ちがう

$ less SRR445818.fastq.bz2 @SRR445818.781 HWUSI-EAS366_145:5:1:1107:1023 length=36 NGTCCNNNNNNNGNNNNGGNAAGAACGNNAAANNNG +SRR445818.781 HWUSI-EAS366_145:5:1:1107:1023 length=36 #################################### @SRR445818.782 HWUSI-EAS366_145:5:1:1156:1021 length=36 NGTATNNNNNNNNNNNNGNNNACTTTGNNNANNNNC +SRR445818.782 HWUSI-EAS366_145:5:1:1156:1021 length=36 #################################### @SRR445818.784 HWUSI-EAS366_145:5:1:1335:1022 length=36 NAAGANNNNNNNNNNNNATNNGGGGGCNNTTTNNNA

5番レーン。ちがう

$ less SRR445819.fastq.bz2 @SRR445819.912 HWUSI-EAS366_145:6:1:1164:1034 length=36 GAGGCNNNNNNANNNNNNNACTGGCACAANATTNAA +SRR445819.912 HWUSI-EAS366_145:6:1:1164:1034 length=36 #################################### @SRR445819.913 HWUSI-EAS366_145:6:1:1323:1025 length=36 NGNNNNNNNNNNNNNNNNNNNNNAANNCNNNNNNNN +SRR445819.913 HWUSI-EAS366_145:6:1:1323:1025 length=36 #################################### @SRR445819.914 HWUSI-EAS366_145:6:1:1420:1030 length=36 NAGGNNNNNNNNNNNNNNNANCGAGGGCNNGANNNN

6番レーンだ!!!

ということでこの SRR445819 が同じシーケンサの機械の同じレーンから得られたデータだと分かりました。もちろん、まだ同一 Run である保証はないので、確かめてみましょう。

まず展開後の kaeru_nazo.fastq が 31,236,448行 (7,809,112 reads), 946MB で SRR445819.fastq が 240,113,360行 (60,028,340 reads), 12GB 。サイズは全然違いますが、ヒト細胞のChIP-Seqで700万リードというのは少ないので、これは実習の難易度 (というか計算時間) を下げるためにreduceしたと考える方が自然です。そうすると、 kaeru_nazo.fastq 中のリードが SRR445819.fastq からも出てくればまあ当たりだと考えていいんじゃないでしょうか。

$ head kaeru_nazo.fastq @KAERU.13652678 HWUSI-EAS366_145:6:26:14947:20927/1 TGGCTGGCAACAATAGATACTGGGGACTACTAGACA + ####################################

$ grep -A 4 'HWUSI-EAS366_145:6:26:14947:20927' SRR445819.fastq @SRR445819.13652678 HWUSI-EAS366_145:6:26:14947:20927 length=36 TGGCTGGCAACAATAGATACTGGGGACTACTAGACA +SRR445819.13652678 HWUSI-EAS366_145:6:26:14947:20927 length=36 ####################################

当たった。

もうちょっとちゃんと確かめた方がいいんだろうけど、まあ、もうこれでいいんじゃないかな…。

kaeru_nazo.fastq は SRR445819 由来ということにします。じゃあこの SRR445819 は先に挙げた論文のシーケンスのうち、どの転写因子の抗体を使ったものなのかを調べれば答えが分かります。通常、Sampleの情報に使った抗体や細胞株を記述してあって欲しいのですが…

書いてなーい。イェーイ。

ご安心ください。よくある話です。別のどこかのDBにこのサンプルのデータが登録/記述されていないか見てみましょう。

SRR445819に対応するサンプルのIDはSRS300774です。今回のデータに関連する一連のIDはSから始まっており、これはNCBIに登録されたデータであることを示しています（ちなみにEから始まるのがEBIが受け付けたデータ、Dから始まるのがDDBJ）。NCBI, EBI, DDBJではさまざまな登録データのサンプル情報を別に管理するBioSampleというDBを運用しており、ここに情報がある可能性が高い。おもむろにNCBI BioSampleにSRSのIDを突っ込んで検索してみます。

雑な検索だけど気にしない。

イエーイ！

c-Mycの抗体を使ったChIP-Seqに使われたSampleの情報が表示されました。つまりこれが先ほどの配列データの元のサンプルの正体です。

無事それっぽい答えにデータベースだけを使って辿りつくことができました。辿りつけましたが、なんかあんまり達成感、ないですね…。fastqを投げ込むと「これだよ」と答えを出してくれるシステムを作ろうかな。

講習終了後、参加者がいなくなってからこっそり講師に「あれc-Mycでしょ」と伝えると「正解―！！！」正解だそうです。「プレゼントあげるよー！！2つあるから2つあげる！！！」

2つももらってしまった。インチキなのに…。まあいいか。何が入っているんでしょう。

なんだこれ……………………………。

couple of slightly rushed but overall successful prism cards! they’ll be slipped into purchases and webstore packages. they won’t be here in time for supanova but i’ll have em around at smash! ^o^

doodle of the OC on break from homework