2016 BioHackathon 1日目
Written by bonohu in misc on 月 13 6月 2016.
本日6/13から6/17までBioHackathon。今回も引き続きAOE2.0に向けてhackします。メタボにならないように、今回はランチはホテルでとらない作戦で。夜にガッツリ方針を議論。自分のやるべき方向がさらに明瞭に。
Category: misc
2016 BioHackathon Symposium
Written by bonohu in misc on 日 12 6月 2016.
今年もBioHackathonは、シンポジウムからの参加。皆さん、進捗していますなあ。
2016 SPARC Japan ワーキンググループキックオフミーティング
Written by bonohu in misc on 木 09 6月 2016.
今年度(2016年度)から国際学術情報流通基盤整備事業(SPARC Japan)のセミナー企画ワーキンググループのメンバーとして関わることになり、そのキックオフミーティングで一ツ橋の国立情報学研究所へ。担当として2つのセミナーの企画運営に加わることに。 実は、3年前に第3回 SPARC Japan セミナー2013「オープンアクセス時代の研究成果のインパクトを再定義する:再利用とAltmetricsの現在」で話をしているこのセミナーだったのだが、流れが大分DBCLSで目指してきた文献やデータの流通を「交通整理」し、「知のめぐり」を良くする方向に生命科学分野だけでなく、全体的にそうなってきたから今回セミナー企画ワーキンググループメンバーとして呼ばれたのかなと。これまで蓄えてきた経験を活かして貢献していきたい。
Homebrewと暮らす
Written by bonohu in misc on 月 06 6月 2016.
Homebrew便利だが、しばらく経つとソフトウェアがアップデートされていってローカルに持っているインデックス情報が古くなり、いろいろと不都合が起きるようになる。そこで、以下のコマンドを日課(or週課?)のように打って欲しい。
[shell]
brew update -v
brew upgrade -v
[/shell]
これをすることでインデックスが更新され、古くなったソフトウェアは自動的にアップデートされる。
そして、たまには
[shell]
brew cleanup -v
[/shell]
するとよい。古くなったバージョンのソフトウェアを消去してくれる。今、したら
This operation has freed approximately 17.2G of disk space.
と出て、17.2Gbyteディスクスペースが空いたらしい。
localBLASTのblastdbcmdを使って配列エントリの特定の領域を取得する
Written by bonohu in misc on 木 02 6月 2016.
新規なモデル生物では、ゲノム配列が公開されていても多くの場合染色体ごとに配列が一本につながっていない。どうなっているかというと、scaffoldという単位でゲノム断片配列が記述されている。
そういった場合でも、コマンドラインで使用する際のNCBI BLASTに含まれる makeblastdb というプログラムを使ってインデックスを作成、blastdbcmd を使ってエントリ名(-entry)と領域(-range)を指定することで、その領域の塩基配列を取得することが可能である。
1 2 3 4 | |
このコマンドでは、一行目でhoge.faファイル(FASTA形式)に対してインデックスを作成、二行目でscaffold001の2000塩基から2500塩基までの配列を抽出する例を示している。通常の makeblastdb のオプションに加えて …
'Mishima.syk #8'
Written by bonohu in misc on 日 29 5月 2016.
今回が8回目のMishima.syk（三島創(製?)薬勉強会）。三島で開催されたSPARQLthon44の次の日に開催をはたらきかけて、異分野交流を目論んでみた。みなさん、いろいろと忙しい中、複数人の人がその流れで参加してくれて、私のねらいは当たったかなと。ただ、誤算はハンズオンの講師の方が緊急入院…。一日も早いご回復を祈るばかりです。
自分は、本業の仕事をもっと知ってもらいたいと思って、直前のEBI出張と絡めて発表してみた。DRY解析教本の人、というイメージ以外にも公共データベースのなかのひとと思っていただけるようなイメージも強くできたかな!?自分の発表スライドは、今回新たに用意されたMishima.sykのgithubのサイトにアップしてみた。
今回の勉強会は、沼津駅北口のプレザヴェルデでの開催で、ちょうどその日がぬまづみなとの街バルの日ということもあって、懇親会は盛り上がりました。一次会として、みんなで適当なグループに(自然発生的に)分かれて、まずバルに参加。バルチケット3枚で、2時間半という設定が絶妙でしたね。そして、その後、バスク料理のお店、うさぎの木に集合して、2次会。素晴らしい采配だったかと。とても楽しかったですね。またみんなで …
Trinityで転写産物を定量する
Written by bonohu in misc on 金 27 5月 2016.
Trinityの出力結果は、転写単位ごとの塩基配列のFASTA形式ファイルだが、付属のプログラム(align_and_estimate_abundance.pl)を使うと、転写量を見積もって定量、そして出力してくれる。TrinityのウェブサイトのTrinity Transcript Quantificationに詳細が書かれている。
homebrewでこれまで色々入れたMacOSX10.11.5の環境では、RSEM(RNA-Seq by Expectation-Maximization)を別途インストール(本家ウェブサイトからダウンロードの上、makeしてmake install)さえすれば以下のような感じのオプションで動いた。上記ウェブサイトに書かれている例ほとんどそのままだが、変更点は--aln_method bowtie2を指定してbowtie2を使うようにしたのと、--thread_countで使っていいthread数を環境に合わせて変えたぐらい。
[shell] /usr/local/Cellar/trinity/2.2.0/util/align_and_estimate_abundance.pl --thread_count 12 --transcripts trinity_out_dir/Trinity.fasta --seqType fq --left hoge_1.fq …
Trinity on MacOSX
Written by bonohu in misc on 月 09 5月 2016.
de novo transcriptome assemblyをするソフトウェアの定番のTrinityがhomebrewで入って動くようになりました(v2.2.0)。これまでhomebrewには入っていたものの、私の環境ではきちんのインストールできないでいました。必要なときはその都度DDBJのスパコンにお世話になっておりました。 それが手元のMacで動くようになったのは素晴らしいことです。
brew install -v trinity
でインストールして、
Trinity --seqType fq --left hoge_1.fq.gz --right hoge_2.fq.gz --max_memory 16G --CPU 4
というコマンドで実行するだけです(ペアエンドの場合)。入力としてはペアエンドのファイル2つ、hoge_1.fq.gzとhoge_2.fq.gzで、見ての通りfastqファイルはgzip圧縮を解かなくても実行可能です。
また、シングルエンドの場合は、
Trinity --seqType fq --single fuga.fq.bz2 --max_memory …データベースとは1
Written by bonohu in misc on 水 04 5月 2016.
昨今、インターネット上には利用可能なデータリソースが溢れている。 データベース(DB)とは、生命の持つ情報を再利用可能な形で電子化したデータとして集めたものである。インターネットの普及とともにウェブインターフェースでDB本体へアクセスできるようになってきたことから、ウェブサイトのことがDBだと勘違いしている人も居るぐらいである。 SwissProt有償にします事件以後、いつかのサイトはさまざまな事情からユーザー登録を必要とするものが出てきて、さらにはそれを有償で提供するものも出てきた。しかしながら、生命科学のDBは現在でもその多くが誰でも無償でアクセスできるようになっている。それは、DBはつかってもらってなんぼ、だからである。論文のオープンアクセスと同じである。
「データベース生物学」カテゴリーとして、DBを使った研究に関して多くの人に伝えたい事を、今後不定期に私見を主体に書き記していこうと思う。
圧縮したファイルで処理する
Written by bonohu in misc on 日 17 4月 2016.
NGSデータ解析の出発点は、シーケンサーやSRA (Sequence Read Archive)から取得したFASTQ形式のファイルなのは同じだろうが、これを圧縮しない(ファイル拡張子が.fqや.fastqのことが多いが、ようするにテキスト形式のファイル)まま、入力として使っているだろうか?配列解析プログラムの多くは、gzip圧縮やbzip2圧縮したファイル(ファイル拡張子は、それぞれ.gzと.bz2)を直接読み込めるようになっていることもあるので、まずファイル解凍するのではなく、圧縮形式から実行可能かどうか、確かめたほうがよい。解凍するプロセスがなくなり、作業効率が格段に向上するからだ。
- v2.2.0のTrinity、入力のFASTQファイルがgzippedでも対応しているのを確認(20160506)
また、データ解析する際に中間ファイルが多数出てくる。そして、そのファイル群は多くの場合、テキスト形式のファイルで、ファイル圧縮もかかっておらず、サイズがでかい。それらをそのまま置いておくとHDやSSDの肥やしになるし、その都度ファイルを解凍してから処理するのも面倒だし。かといって消してしまうと後で困ることも。というわけで、どうせ圧縮するのなら前にも触れた並列版bzip2のpbzip2を使って圧縮し、ディスクスペースをどんどん空けていこう。
それらの中間ファイルをテキスト処理するには、例えばPerlだと以下の様な感じでちょっとコードを変えるだけで済む。参考まで。
[perl] my …
行番号を付ける
Written by bonohu in misc on 月 11 4月 2016.
バイオの業界ではとくに、IDが重複していることがある。Rでの処理の場合、それでは入力が受け付けてもらえず実行できないことがよくある。そんな時には行番号を付け、それをuniqueなIDにしてしまうという解決策がある。その処理をRの内部ですることもあるが、今回はUNIXコマンドで行番号を付けてくれるnlコマンドを使う方法でやってみたので、書いておく。
[shell] nl -n ln -v 0 file.txt [/shell]
細かいオプションはここを参考にした。-n lnで左詰めに、-v 0で0からスタート(file.txtにヘッダ行が1行あったので)。
コマンド処理万歳。
joinコマンドの出力をタブ区切りに
Written by bonohu in misc on 水 06 4月 2016.
joinというUNIXコマンドは、2つのファイルの同じ値のある行をjoinして、一行にまとめてくれるコマンド。
join -j 1 file1 file2
でタブ区切りテキストのfile1とfile2の一番左側のカラムの値が同じ行が結合されて出力される。
このコマンドのやってくれる事自体は大変便利なのであるが、デフォルトではその後の出力はスペース区切りとなり、ちょっと不便だった。これをタブ区切りに変えるやり方を知ったので、こちらにも転記しておく。
join -j 1 -t "$(printf '\011')" file1 file2
結合する前のファイルはsortされていないといけない。そうなっていないと
join: file1:7: is not sorted
というエラーが出る。
sort -u file1 > file1s
sort -u file2 > file2s
のようにsortしたファイルで処理するようにする。-uオプションを付けることでsortしたうえに重複した行も除くように。
あとjoinするのに用いるword(キー)もユニークでなくてはならない。すなわち、二回以上出てきてはダメ。実は今日やっていたのはこちらの制約に引っかかり …
scikit-learnのインストール
Written by bonohu in misc on 水 30 3月 2016.
pythonの機械学習ライブラリ。ふとしたキッカケで入れてみようと、思い立ったが吉日。
pythonのバージョン3が入っていなければ、まずそれを。
[shell] brew install -v python3 [/shell]
これでインストールされるpip3を使って必要なライブラリを。
[shell] pip3 install scikit-learn [/shell]
の他に、scipy, numpyあたりが必要なようです。
このサイトにあるサンプルコードを動かして自分のやりたいデータ解析への応用を探っているところ…。
RNA-Seq実験ハンドブック
Written by bonohu in misc on 土 19 3月 2016.
RNA-Seq実験ハンドブックという本が出版されます。出版社のページによれば2016年03月23日発行予定とのこと(2016年3月発行予定に変わっていました)。「RNA-Seqはこうして誕生したーそしてデータ再利用へ」と題したコラムを書かせてもらいました。短いですが、遺伝子発現解析の歴史的なことを簡単に書かせていただき、本業で推進しているデータ再利用の重要性を説いております。そしてRefExの宣伝も。
[amazon template=thumbnail&asin=4758101949]
Amazonのページには2016年3月26日発売となっているようですが。手元に届くのが楽しみです。
第60回 人工知能学会 分子生物情報研究会(SIG-MBI)
Written by bonohu in misc on 金 18 3月 2016.
続けて同じ場所で開催されたSIG-MBIにかなり久しぶりに参加。SIG-MBIの開催が第60回にもなっていたとは。続けて開催されているのは大変素晴らしい。敬意を表します。なお、自最後に参加した2005年の第1回オープンバイオ研究会との共催の時には、「第3回人工知能学会生命知識研究会(SIG-BMK)」だった模様。時が経つのは早い。
研究会ですが、発表者がその第1回オープンバイオ研究会の発表者とメンツがあまり変わりなく。まあ、他にも発表する機会が今回の研究会クラスタであったからと思いたいが、参加するのなら発表するようにしたほうがいいんじゃないかなと老婆心ながら。
自分は、「公共遺伝子発現データの再利用に向けた取り組み」と題して、NBDCの共同研究としてDBCLSでやってきた遺伝子発現データを使いやすくするプロジェクトに関して20分ちょい話をさせていただいた。DBCLSやNBDCについてもよく聞く略語であっても知る機会が少ないと思い簡単に紹介させていただいたが、プチやはりだった年表形式のほうがよかったかなと反省。若手の方から質問がなかったのは興味がなかったのか、私の話し方が悪くてわからなかったのか…。わかりやすく伝わるよう、中身もプレゼンも、さらなる精進が必要。
BH15.15成果報告会と第20回オープンバイオ研究会
Written by bonohu in misc on 木 17 3月 2016.
朝からみんなでバスでJAISTに移動して成果報告会。遺伝子発現目次の進捗状況に関して報告。いつもどおり、みんなで議論して情報を共有したことが大きな成果。githubにアップしたmetadataからの抜き出しスクリプト群の作成が目に見える成果ですが、それ以外には、solrからのデータ抽出のバグというか、データが大きくてもうまくいくやり方を経験的に知り得たことが大きかったかと。
その後の第20回オープンバイオ研究会では、公共データベースの利用を促進する取り組みを、非モデル生物種に関して注力しています、という話を四コマプレゼンでしてみた。統合DBプロジェクトに関わっていても私自身が普段やっていることに関してDBCLS以外の方は知ることがない内容の話だろうとふとと思って、納品されたばかりのテンプレートを使ってプレゼンを急造して即席で。その後の第3回オープンサイエンスアワードでは関わっているプロジェクトが何個か表彰していただいた。たいへん励みになります。ありがとうございました。
BH15.15 3日目
Written by bonohu in misc on 水 16 3月 2016.
ArrayExpressのメタデータにBioProjectがなかったので、GEO由来のエントリだけでもBioProjectIDが付けられないか模索中…。すべてのデータをスクレイピングする以外に解決法はないのか?いろいろ調べた結果、結局bioprojectのXML(bioproject.xml)から生成するスクリプトを書いて、SRA由来でない発現データにもBioProjectIDを付与するための元データが得られた。GEO由来のエントリだけだが。ArrayExpressにしかないものはまた別にrescue方法を考える。
やり残していたExperiment(xRX)とRun(xRR)の対応表を1つの実験毎に重複なく作成するスクリプト作成。これで発現定量した結果とそれに相当するExperiment,さらにはStudy(xRS)を結び付けられるように。データ統合化まであと一歩!
それらをまとめて、Perlスクリプトだけだがgithubに新規のrepository(AOE)を作成してpush。これだけ1日に何回もgit push -u origin masterしたのは久しぶりだったり。
Transcriptome analysis (RNA-seq)だけでなく、すべてのstudyデータに対してもjoinしたデータを作成して、今後のデータ更新の元となるナニカが出来た感。
BH15.15 2日目
Written by bonohu in misc on 火 15 3月 2016.
2日目は抽出したメタデータの抜けをチェックして出来る限り取りこぼさないようparserを改良。1つのメタデータに複数のデータ単位が含まれていることをランチ前に発見して、ランチタイム遅延。対処する目処をつけてからカレー、その後本格実装。タグの多様性がみとめられたので、IDの正規表現で最悪補完できるようにした。その結果、ファイルに抜けが減り、TIBCO Spotfireでのjoinが急激に遅くなったため、join用にスクリプト作成。AOEに実装するのに必要なデータが大分揃ってきた。あとは、Series-Runの対応表と更新日時のデータ。そして、AOEにすでにあるArrayExpress(GEO)由来のデータとの統合。
BH15.15 1日目
Written by bonohu in misc on 月 14 3月 2016.
今日2016年3月14日から、国内版Biohackathon(BH15.15)。2015年15月ということで、2016年3月の開催となったが、このパターンでの開催もついに今回で最後(というかルール破綻)になるとのこと。「SRA/BioProject/BioSample etc の統合的な検索」の中の「遺伝子発現データ検索への応用」というサブテーマで取り組みます。
先月のSPARQLthonでSRAのメタデータからの抜き出しに問題があったのが、原因を究明。どうも、csv形式での抜き出しに問題があるようで、json形式に変更してきっちりとすべて抜き出してこれるように。これで必要なメタデータのうち、Submission date以外は取れたはず。1日目の進捗はとりあえずそれということで。
あれから5年
Written by bonohu in misc on 金 11 3月 2016.
これからどうなるのか。そう思ったあの日から5年経ち、活動本拠地は静岡県三島市の国立遺伝学研究所に移動したものの、相変わらず生命科学データベースの普及活動を中心とした研究生活を続けられている。「幅優先」だった活動は、「深さ優先」にシフトし、その頃からも実はやっていた個別の共同研究の割合が増えてきているが。
その前の年に始めた共同研究は、皮肉にも5年前の今日に取ったデータが決め手となり、紆余曲折あったものの約3年前に論文に。その共同研究は継続しており、今年のこの日にも打ち合わせが。研究を続けられていることは素晴らしい。関係者の皆さんに感謝の念がやまない。これまでに持ちかけられた多くの共同研究は途中で頓挫してしまっており、論文化していないものが多数であるからだ。
この先5年、どうなるかなんて全くわからない。現に5年前の今日は、今日の日のことを全く予測できなかった。ただ、やるべきと思うことを続けてやっていくことだけだ、と心に再確認した5年目の3.11だった。ありがとう。そして、これからもよろしく。
標準出力と標準エラー出力をまとめる
Written by bonohu in misc on 金 04 3月 2016.
UNIXのコマンドラインには標準出力と標準エラー出力があります。プログラムの出力結果は通常標準出力で、
[shell] sh run.sh > log.txt [/shell]
とすることでlog.txtにその出力結果が記録されます。このようにしてもまだ画面に何か表示されることがあります。それが標準エラー出力の結果です。
[shell] sh run.sh > log.txt 2> logerr.txt [/shell]
とすることで、標準出力はlog.txtに、標準エラー出力はlogerr.txtに記録されます。2つのファイルに分かれさせたくない場合は以下のように。そうするとlog.txtの1つのファイルに出力が記録されます。 [shell] sh run.sh > log.txt 2>&1 [/shell]
ちなみにこれは使っているシェルがbashやzshなどで、cshやtcshだとコマンドが異なりますのでご注意を。
タブをgrep
Written by bonohu in misc on 木 03 3月 2016.
grep便利ですね。タブ区切りのテキストで必要な情報を持つ行だけ取ってくるとか、多用されます。ですが、 [shell] grep chr1 hoge.gff [/shell] では1番染色体(chr1)のデータだけ抽出したいのにchr11やchr12のデータなども引っかかってきます。そこでワードがタブで区切られていることを利用して、chr1(tab)で検索すると1番染色体のデータだけが取ってこれそうですよね。
それを実現するには以下のようにすればそれが実現できます。 [shell] grep chr1$'t' hoge.gff [/shell]
samtools view も並列化
Written by bonohu in misc on 水 02 3月 2016.
ときすでにsamtoolsのバージョンは1.3。single threadで遅く、時にはそれが律速になっていたsambam変換も並列化されている。 samtools sortの時と同様に、-@フラグで使用するCPU数を指定。例えば、
[shell] samtools view -@ 8 -bS file.sam > file.bam [/shell]
で。
indexのほうは並列化オプションは現状ないらしい。*.bamとかやっても複数のファイルのindexは作れず、既存のbamファイルを壊してしまうので注意。
[shell] foreach f (*.bam) samtools index $f end [/shell]
とするしかないのかな。
'ファイルの移動方法: tar編'
Written by bonohu in misc on 火 01 3月 2016.
大量のファイルの移動が必要な季節になってきましたが、我々DB屋は年がら年中です。複数のディレクトリ階層構造を持ったファイルの同期にはrsyncが一番ですが、単にそういったデータを転送するのであれば大学院生の頃(約20年前)から使っているtarを使った方法があります。
[shell] cd /Volumes/old_HDD tar cf - |(cd /Volumes/new_HDD; tar xvf -) [/shell]
移したいディレクトリのあるところにcdして、その次の行のコマンドを実行します。/Volumes/new_HDDは移したい先のディレクトリ名です。サイズや複数のファイルを持つディレクトリを移すのに今でもよく使っています。
SPARQLthon41
Written by bonohu in misc on 水 17 2月 2016.
2016年2月16,17日に理研和光にてSPARQLthon41に参加。引き続き、DBCLSで開発中の遺伝子発現目次AOEに公共遺伝子発現データベースArrayExpress(含むGene Expression Omnibus)に無いデータ(=Transcriptome sequenceによる遺伝子発現定量)を補完するのに取り組む。 SRA(Sequence Read Archive)のデータをSolrで検索できるようにしてもらったのを利用して、加えるべきデータセットを生成するスクリプト作成とか。メタデータを精査していくと、NCBIのXMLにはsubmission dateのデータがないことが発覚。その部分以外はAOEに載せるべきデータはとりあえず作成。
-
studyから
-
PRJ (BioProject ID)
-
GSE (GEOのSeries ID)
-
xRP (SRAのProject: SRP,ERP,DRPから始まるID)
-
Title
-
-
experimentから
-
xRX (SRAのExperiment: SRX,ERX,DRXから始まるID)
-
PRJ
-
SAMN (BioSample ID)
-
Sequencer
-
-
sampleから
-
xRS (SRAのSample: SRS …
-
統合データベース講習会AJACS薩摩
Written by bonohu in misc on 日 31 1月 2016.
2016年1月26,27日に鹿児島大学桜ヶ丘キャンパスでの統合データベース講習会AJACS薩摩に講師として参りました。今年度、3回目。
前回参加のAJACS津軽では、遺伝子発現DBの使い方とR/Bioconductorを使ったデータ解析入門を担当し、遺伝子発現データ解析主体の内容でした。今回は、次世代シークエンサー(NGS)データ解析の入門編を合計4.5時間かけてガッツリとやって参りました。DRY解析教本も出版されており、説明しやすかったというのが正直なところ。
1日目のラストに前半「ゲノム配列とそれに付与されたデータの使い方」でゲノム配列データの探し方とゲノムブラウザの使い方をわりとしっかり目に。それに続く2日目の午前に、NGSデータ化をどう探すか、どう活用するか、の講習を。
とくに2日目の方は新作コンテンツで、3時間も割り当てられたため、余ってしまうんじゃないかと危惧していましたが、実際やってみると計ったかのように時間通りでした。みんなでアクセスするとDBCLS SRAがつながらなくなったり、USBメモリで配布したデータのやりとりで問題があったり、細部ではいろいろと不都合ありましたが、なんとかやり通せたかと。今後はもっと若い人に講師をやっていただけるよう、していきたいですね。
command lineで github に pull request を送る
Written by bonohu in misc on 木 21 1月 2016.
統合データベース講習会の資料はmarkdownで書かれており、それらをgithub上にアップすることで講習会資料のPDFになったり、統合TVのウェブサイトから見えるようになったり、systematicに処理されるようになっています。それをアップするにはAJACS-trainingさんにpull requestを送り、以前のコンテンツと置き換えて貰う必要があるのですが、これを今まではgithubのウェブインターフェース上で行っておりました。
今回、「コマンドラインでGitHubにプルリクエストを送る」というqiitaの記事を参考に(というかほとんどそのままですが)、コマンドラインのみでやってみたので備忘録的に書き記しておきます。 自分の書いたmarkdownのファイルは~/hoge.mdという想定で、AJACS58のbono1というところにpull requestを送る場合の例です。
[shell] git clone AJACS-training/AJACS58 cd AJACS58 git checkout -b bono1 cp ~/hoge.md bono1/README.md git add bono1/README.md git commit -m "draft version" git fork git push …
2016年の計
Written by bonohu in misc on 金 01 1月 2016.
2015年は、5月に領域融合レビュー「次世代シークエンサーにより得られたデータの解析」を公開、10月に「次世代シークエンサーDRY解析教本」を出版し、現場でDB利用を’method’の一つとした研究が増えていくよう実例作りに貢献することができた。また、非モデル改め次世代モデル生物のデータ解析手段の開発を共同研究ベースで複数進めることができたのは当初の予想を上回り、とても良かった。
一方で、公共DBの全レコードを対象としたデータ解析研究がおろそかになってしまったのが残念であった。時間は有限なので、その辺うまくバランスを取っていくことを2016年の課題としたい。できれば、一緒にやってくれる仲間を増やしていきたい。そのためにも、我々のやっていることや考えていることの情報発信を積極的に。統合牧場の閉鎖など、諸事情あって低調になってしまっている状況は反省すべき。
2015年年頭のブログエントリで課題としていた他の遺伝研の部門と融合的な活動、とりわけDDBJを補完する活動が低調だったのは反省すべき。遺伝研に来て3年目となる今年は、特に重点的に取り組んでいきたい。考えていることの体系化も含めた情報発信、今年も頑張ります!
sequenceserver
Written by bonohu in misc on 木 29 10月 2015.
先月のBioHackathonで参加者たちがそれぞれに持ってきた研究紹介のパンフレットを見ていて発見したのが、このsequenceserver。何の事はない、Ruby製のlocalBLASTのGUIなinterfaceなのだが、これがinteractiveに複数BLASTをかけてそのalignmentの結果をチェックしたい時になかなか便利。便利だと思っている点をまとめると
-
それはqueryに対してどこにヒットしたかという、かつてのpaintBLAST的な出力がトップに出る
-
queryをコピペしたら塩基配列かアミノ酸配列か自動判別してくれる
-
検索対象のDBが、塩基配列系とアミノ酸配列系でリストされ、チェックボックスで選べる
-
足したいパラメータを追加で指定できる
といったところ。繰り返して使う時の設定しないといけない点の少なさがとても良いかと。
インストールは、以下のとおり、gem installで。 [shell] sudo gem install sequenceserver [/shell]
だが、これがかなり時間がかかるので、
[shell] sudo gem install -V sequenceserver [/shell]
にしてverboseモードにしたほうが精神衛生上よいかと。インストールできたら
[shell] sequenceserver [/shell]
でサーバー起動、ウェブブラウザも起動してhttp://localhost:4567/へアクセスしに行く。そうするとlocalBLASTの画面が出て、利用可能なDBがリストされる。デフォルトだとcurrent directoryのため、なにも表示されないかもしれないので …
トーゴーの日シンポジウム2015
Written by bonohu in misc on 月 05 10月 2015.
今年も、10月5日にトーゴーの日シンポジウム。今年は2日にわたって、東京大学弥生講堂・一条ホールにて。全体的な印象としてはデータベースが生命科学研究のいろんな分野において使われ、重要な位置を占めるようになってきたな、と。 成果発表会的にRefExとAOEに関するポスター発表。意外に反響があって、もっと頑張らねば。そのためにはやるべきことの断捨離とアウトソーシングが必要かな。
BioHackathon2015
Written by bonohu in misc on 月 14 9月 2015.
2015年9月14日〜18日までの丸一週間、長崎にてBioHackathon2015に参加中。ハッシュタグは #biohack15。今回は、AOEにGEO/ArrayExpressに入っていないRNA-seqデータを含めて、名前の通りAll of Gene expressionデータが検索できるようにしようと目論んでおります。これまで同様、それらのスクリプト群は、githubからオープンにします。
SRA形式ファイルの料理法
Written by bonohu in misc on 火 25 8月 2015.
以下の様な声をたまに聞く。
DDBJのDRAにはFASTQ形式のファイルがなくてダウンロードできない
と。FASTQ形式はテキスト形式で、圧縮しない状態ではディスク容量を喰う。ブログエントリを書いている2015年8月末時点ですでに約2.5ペタバイトものデータ容量がSRAだけであり、その中間ファイルとか含めると…。そのためにFASTQ形式そのものでは保存されておらず、SRA形式という圧縮もかかった形式で利用可能となっている、というのがそれに対する答え。
ではSRA形式はどう処理すればよいのか?FASTQ形式のファイルに変換できるのか?それには、NCBI謹製のSRA Toolkitを使えば良い。様々なプラットフォーム向けのそれがダウンロードサイトから利用可能となっている。Macならhomebrewで
[shell] brew install -v homebrew/science/sratoolkit [/shell]
このSRA Toolkitの中のfastq-dumpというプログラムを使えばFASTQ形式のファイルが得られる。
[shell] fastq-dump hoge.sra [/shell]
ただ、ペアエンドのSRAエントリの場合には以下のようにしてファイルが分割されるように指示してやる必要がある。
[shell] fastq-dump --split-files hoge.sra [/shell]
統合データベース講習会AJACS津軽
Written by bonohu in misc on 月 24 8月 2015.
2015年9月3,4日に弘前大学での統合データベース講習会AJACS津軽に講師として参ります。今年度二回目の「代表招集」。
前回、大阪大学吹田キャンパスでのAJACS千里では、ゲノム配列データベースの使い方と遺伝子発現データベースの使い方の2コマ担当でしたが、今回は遺伝子発現DBの使い方とR/Bioconductorを使ったデータ解析入門ということで2日目のすべてを担当します。今回は、遺伝子発現解析の入門編を合計3時間かけてガッツリとやる予定です。
担当分の前半は、前回の「遺伝子発現データベースの使い方」のアップデート版、後半はR/Bioconductorを使ってデータ解析する入門編です。これまで出してきた論文でR/Bioconductorを使って図表を作ってきましたが、その実例を皆さんと手を動かしてやってみたいと思っていますので、ふるってご参加ください。当日の資料は、絶賛作成中ですが、以下のページですでに公開してしまっています :-p
blast2html
Written by bonohu in misc on 日 23 8月 2015.
paintBLASTとか知っているだろうか?かつてNCBIのウェブサイトでBLAST検索した時に上部に、queryのどの部分にDB中の配列がマッチしたかを可視化してくれるクリッカブルイメージがあったと思うのだが、それを生成するためのツールがそれだったと思う。 それを使ってBLASTの結果を可視化したいと思い、探したのだが、どうも維持されていない。いろいろ検索して行きあたったのがpythonのツールでblast2html。localBLASTの結果(xml形式)をHTMLに変換して出力し、paintBLASTのような結果も表示くれるものらしい。というわけで早速
[shell] blastn -query query.fa -db db.fa -outfmt 5 > hoge.xml [/shell]
としてXML形式(-outfmt 5はXML形式出力の暗号)の(local)BLASTの結果を用意して
[shell] python blast2html.py -i hoge.xml -o hoge.html [/shell]
で実行してみるとエラーが。
UnicodeEncodeError: 'ascii' codec can't encode …
圧縮を並列化
Written by bonohu in misc on 土 22 8月 2015.
データ量は大きくなるばかりで、それをファイル圧縮して保存しておくなり、転送するなりが当たり前になっている。かつてはCPUのクロック数が上がることでこの圧縮が早くなることが望めたが、最近はクロック数はそれほど伸びなくなってきた。その代わりに1つのコンピュータが持つコア数が増え、並列計算でそれを実現することが良い解法となってきている。ファイル圧縮のプログラムとして長らくgzipやbzip2を使ってきたが、これらは並列化されておらず、1つのファイルあたり数Gから数十GbyteあるFASTQファイルの圧縮／解凍に長い時間待たされることが多くなってきていたが、そんな我々に朗報。
pigzがgzipの、pbzip2がbzip2のthread版として利用可能。Macならhomebrewで一発で入ります。
[shell] brew install -v pigz brew install -v pbzip2 [/shell]
使い方はこんな感じで、使うthread数を指定する必要は特にないようだ。
[shell] pigz large.fastq pigz -d large.fastq.gz [/shell]
デフォルトでCPUを使えるだけ使うようなので、&つけてバックグラウンドで走らせないほうがいいかと。
SPARQLthon35
Written by bonohu in misc on 火 11 8月 2015.
2015/08/10-11は、DBCLS三島というか、DDBJのある国立遺伝学研究所にてSPARQLthon35回目。前回の続きで、SRAに登録されているものの遺伝子発現データベースに登録されていないエントリがどれぐらいあるかの調査。予想に反して多かったので、別ソースでIDの対応が付けられないか、SRAをやっている同僚の仲里さんに手伝ってもらって対応をとってみた。
前回、5234ものRNAseqのSRAエントリがBioProjectID単位でArrayExpress(GEO)に載っていないというショッキングな結果だったが、そのうちいくつかは今回の名寄せで補われたものの、やはりかなりの数ある。BioProjectIDの発行元別にその結果を以下に示す。
-
PRJN: 4728 -762 = 3966
-
PRJE: 337 - 48 = 289
-
PRJD: 169 - 0 = 169
その一方、ChIPseqのほうは、多くが今回の名寄せで対応が着いた模様。こちらもBioProjectIDの発行元別に。
-
PRJN: 492 - 354 = 138
-
PRJE: 45 - 19 = 26
-
PRJD: 55 - 0 = 55
byobu使ってみた
Written by bonohu in misc on 木 06 8月 2015.
screenの使い方の投稿の修正をポストしたついでにtwitterに流したら、byobuというものがあってこっちのほうがナウでヤングだという御意見を頂戴する。Linuxでの紹介記事も教えてもらったが、
[shell] brew install -v byobu [/shell]
したらやっぱりhomebrewにもあってサクッと入った。細かい使い方は全く一緒なんだが、Function Keyを使ってバッファの遷移ができるとか、エスケープシーケンスを使わなくても「使える」仕様はなかなか。しかし、エスケープシーケンスをこれまでのように使ってscreen互換的にも使えるので満足。一番下の行にバッファーリストとか出てて、このプログラムが起動していることが一目瞭然だし。 byobuを起動後、初めてCntl+aを入力した時にエスケープシーケンスとしてそれをつかうのか、はたまたEmacs使いなのか、的にinteractiveに訊かれ、自分の好みのエスケープシーケンスを設定できる。F9キーを押してもその設定ができるから便利。 MacBookProとかならFunction Keyがあるからよいのだが…HHKキーボードだとちと…。まあ、一度エスケープシーケンスを設定してしまえば関係ないのだけれども。
SPARQLthon34
Written by bonohu in misc on 金 17 7月 2015.
2015/07/16-17の両日は、SPARQLthonの34回目でDBCLS柏。前回からやっている、SRAに登録されているものGEOやArrayExpressといった遺伝子発現データベースに載っていないデータをリストアップしてみようの続き。比較するために、BioProjectのIDに変換して同じだたかどうかを判定。データゴニョゴニョは前回でやっていたので、それの解釈が主な今回の仕事。まずは、その重なりをVenn図で可視化。
ArrayExpress(GEO)にない5234エントリ(5226+8)のうち(RNAseqの方)、BioProjectIDの由来で分類。すなわち、NCBI由来のPRJN、EBI由来のPRJE、DDBJ由来のPRJD別に。RNAseqの方は、
-
PRJN: 4728
-
PRJE: 337
-
PRJD: 169
同様にChIPseqの592エントリでは
-
PRJN: 492
-
PRJE: 45
-
PRJD: 55
結構抜けているし、それはDDBJ由来というよりは…という結果に。ちょっと抜けていると判定されたデータが本当にそうなのか、大規模に調べてみる必要がありそう。RNAseqデータ、ArrayExpress(GEO)にはないけど、「SRAにはありまぁす」という結果になるかもしれない …
screenコマンドでプロセスは続くよどこまでも
Written by bonohu in misc on 水 24 6月 2015.
コマンドラインでのデータ処理はすぐに終わらないものもあります。解析すべきデータ量が多い昨今、とくにその傾向が強いようです。そういう場合に役立つのがこのscreenコマンドで、大学院生時代からかれこれ15年以上愛用しています。使い方は簡単で、まずはscreenを起動するサーバーで
[shell]
screen
[/shell]
と打ってscreenを起動します。そうすると起動画面が出たあとはコマンドプロンプトが出て普通にUNIXコマンドラインを使えば良いです。
ホームディレクトリ以下の.screenrcに
escape ^Z^Z
と書いて、defaultだとControl(Cntl)+aになっているescapeキーをCntl+zにしてつかっているので、以下それを想定して書きます。変更する理由はCntl+aはコマンドライン操作において一番先頭にカーソルを移動するショートカットで、直前に実行したコマンドを再度実行する際の編集などで多用するからです。
作業終了するときには、Cntl+zを押してからdをタイプしてdetachします。そうすると目の前の端末からscreenのセッションが「離れ」ます。これを実行する前にデータのダウンロードや時間のかかるプログラムを実行開始しておいてもそれが続けて実行されつづけます。
再び再開するのは、そのマシンにsshでログインして、 [shell] screen -r [/shell] でスクリーンのセッションを呼び戻します(resume)。detachするまえのセッションが再開され、たとえばデータのダウンロードを実行していたのが今はどうなっているか、わかります。また、 [shell] screen -x …
SPARQLthon33
Written by bonohu in misc on 火 23 6月 2015.
2015/06/22-23の両日は、SPARQLthonの33回目でDBCLS柏。以前より、この会で取り組み続けている「NGSデータ解析ツールのDocker化とその環境づくり」を続けてやっているところ。最近は主にRNAseqデータの定量をするためのそれに取り組んでいます。日本のSRA(=DRA)に登録されたRNAseqデータはGEOにもArrayExpressにも自動的には登録されないということが指摘され、そういったデータをリストアップしてみようと動き出したところです。
領域融合レビュー「次世代シーケンサーにより得られたデータの解析」公開
Written by bonohu in misc on 月 18 5月 2015.
ぼうのブログの「2010年代のバイオインフォマティクス」カテゴリを創設したときに
-
概念的なまとめを書いた「薄い本」パート
-
「データ解析プロトコル」パート
に分けて文書化を行っていくと書いたが、ついに、その「薄い本」パートの第一弾となる「次世代シーケンサーにより得られたデータの解析」が公開された。もはや次世代ではなくなりつつある次世代シーケンサーから得られたデータを、広く多くの方に知ってもらうという位置づけで、日々データ解析する人たちから話を聞いて無料コンサル的なことをしている私の頭のなかにあるこのデータ解析の全容を書かせてもらった。私自身も小規模ながらデータ解析を実際にやっており、そういう意味では現場の声かと。多くのそういったデータ解析をやっている方の目に触れ、「なんだ、こんな内容でも書いていいのか」と思われ、自身も情報発信されんことを願っています。
領域融合レビュー, 4, e008 (2015) DOI: 10.7875/leading.author.4.e008 Hidemasa Bono: Sequence data analysis in life science utilizing next generation …
ツキとあたしとfigshare
Written by bonohu in misc on 日 10 5月 2015.
外部での講演、講義、研究打ち合わせが続いた2015年春の大型連休の前後だった。並行していろいろ書き物。振り返ってみると悲しいことを忘れようと仕事に打ち込んでいたかのよう。
-
北海道大学遺伝子病制御研究所 (2015/04/21)
-
東京大学大学院農学生命科学研究科 (2015/04/27)
-
群馬大学大学院医学系研究科 (2015/05/01)
そして5月8日はDBCLS三島にてお客様が来て、やっていることを発表。その際使ってきたプレゼンテーション(略してプレゼン)は上記の東京大学大学院農学生命科学研究科「情報生命工学」の講義の一部を担当させててもらった時の講義スライドがベースとなっている。この2015年春バージョンのスライドをベースに各場面に合わせて改変して発表していることを白状する。このスライドはfigshareというサイトにアップされており、以下の様にciteできる。
Bono, Hidemasa (2015): データベースを活用した知のめぐりのよい生命科学. figshare.
http://dx.doi.org/10.6084/m9.figshare.1391958
Retrieved 02:04, May 10, 2015 (GMT …GFFでIGVにアノテーションを表示する
Written by bonohu in misc on 金 24 4月 2015.
IGVに表示するためのGFFを出力する捨てコード(Perl)書き。入力はBLASTの結果のタブ区切り出力(-outfmt 6指定して出てくるもの)。以下のコードをhoge.prlとして、 [shell] perl hoge.prl sample1 < BLASTout.txt > BLASTout.gff [/shell] のようにして使う。
9カラム目のattributesは任意だと思っていたのだが、ここでID=を指定しないとIGVにtrackとして表示されない模様。NAME=はそのtrackの該当箇所にマウスインした時に表示される詳細ウインドウに表示されるので、書いた方が良い。 また、6カラム目のscoreだが、0-1000の範囲で書けるということで、パーセント一致度(percent identity)を10倍して小数点以下切り捨てで。しかし、IGVのどこにその値が反映されて表示されているのだろう?
[perl]
my $prefix = shift(@ARGV);
while(
NCBI HomologeneのFASTA出力ヘッダを生物種名に変更
Written by bonohu in misc on 水 22 4月 2015.
NCBI Homologeneからいろんな生物種でのある遺伝子のホモログをリストアップした際、デフォルトではそれのmultiFASTA形式ファイルでgiとRefSeqのIDがラベルとなって出力される。それでは人間様にとってはぱっと見わかりにくいので、FASTAヘッダ中に[Homo sapiens]のように含まれている生物種名を抽出してラベルとして書き換える。
以下の捨てコードをhoge.prlとして、homologeneの出力を標準入力にして [shell] perl hoge.prl < homologene.txt > homologene2.txt [/shell] のように使う。
ただ、一つの生物種に複数該当エントリがあるとIDが重複し、多重配列アラインメントプログラムで文句言われるので注意。その場合は中身を見て重複しているIDは違う名前にすべし(例えば、Homo_sapiensをHomo_sapiens2にするなど)。
[perl]
while(
匍匐前進
Written by bonohu in misc on 月 13 4月 2015.
年度初めからショッキングな出来事が多く、仕事がなかなかはかどっている感じがしませんし、実際そうだと思います。匍匐前進、こなしていくしかないですね。
突然のreboot
Written by bonohu in misc on 土 04 4月 2015.
- Docker
前日に引き続き。同じようにO/Nでtophat走らせていたが、やっぱり突然virtualboxがreboot。なんなのだろう?
root@21ee049caf37:/data# tophat -o test_out -p 12 -r 100 hg19 1_1.fastq 1_2.fastq
[2015-04-03 01:46:13] Beginning TopHat run (v2.0.14)
-----------------------------------------------
(中略)
[2015-04-03 10:24:39] Searching for junctions via segment mapping
Broadcast …Tophat2 on docker in Mac
Written by bonohu in misc on 金 03 4月 2015.
- Docker
昨晩帰りがけにO/Nでのrunを仕込んだので、時間はかかるが一晩で終わるかと思いきや、途中でVirtual box自体が落ちてて終わらず。
仕方がないので、1万readに絞ったデータセットを作って、それで実行。するとすぐに実行終了してbamファイルが出来ていた。これでdocker上でもtophatはちゃんと動いて、大容量ボリュームに結果を記録できることが無事出来た。
前に作ったtophatのdockerfileとは異なり、必要なツールをインストールするだけで実際の実行は切り離し、docker上でinteractiveにできるように。これをいかに人手をかけずに実行できるようにするかが次の課題に。
SPARQLthon31
Written by bonohu in misc on 木 02 4月 2015.
- Docker
SPARQLthonが、職場のある三島での開催の今回。だが、自分は引き続きDockerでbioなツールを動かす系の構築に勤しむ。今回のキーワードは、「ゆくゆくはSPARQLを知らないといけない」。
前回のhackathon後、de novo transcriptome assembly系のTrinityをDocker上で動かすことを試みた。プログラムは走るものの、最後まで終了せず。実行に時間がかかるので、それと並行して以前にも取り組んでいて放置していたReference genome有り系のtophatのDockerfileとそれを動かす環境を検討。
このサイトを参考に固定バージョンのtophat(とbowtie)をインストールして動かせるようにするDockerfileを書いた。現状、Bowtieは2.2.5、tophatは2.0.14。そして、以下のコマンドでdockerを起動。
[shell]
coreos> cd share
coreos> docker run -it -v pwd:/data bonohu/ubuntu-tophat2
[/shell]
起動したらdocker上では
[shell]
cd /data …
ブログ回帰
Written by bonohu in misc on 水 01 4月 2015.
twitterによる情報発信に依存するようになってもう5年以上になるだろうか。引っ越しする前の職場では多くのアルバイト学生さんが居て、彼らに自分たちの仕事や状況を伝えるのに便利なツールとしてtwitterを使ってきた。また、不特定多数の人に自分のことを伝えるツールとしても活用させてもらってきた。
確かにtwitterによる短文は書くのが楽でいいのだが、相手に伝わる文章を書くことをしなくなってしまって、結果としてかなりの時間をさいて一生懸命になって書いても相手に伝わる文章にならなくなってしまっている状況に気がついた。そう感じることがここ一週間で複数回あり、もちろんtwitterだけのせいではないと思うのだが、相手に伝わる文章を書く努力を怠っていることには違いないと自戒。年度頭という良い機会でもあるので、ブログをメインに情報発信する生活に回帰したいと思う。主に「雑感」カテゴリになると思うが、今年度は出来る限りブログ発信していきたい。
Picard MarkDuplicates
Written by bonohu in misc on 月 23 3月 2015.
NGSな変異解析にはsamtoolsの他にGATKやPicardといったツールを使うらしい。Picardだが、「ピカール」と読むらしい実は「ピカード」という英語発音らしい。インストールはいつものhomebrew。
[shell] brew install -v picard-tools [/shell]
いろいろできるらしいが、まずは重複を見つける(MarkDuplicates)。
[shell] picard MarkDuplicates ASSUME_SORTED=true REMOVE_DUPLICATES=true INPUT="hoge.bam" OUTPUT="hoge_rem.bam" METRICS_FILE=duplicate [/shell]
初めて動かしてみて、バージョンが1.128に上がっていて「0から始めるエクソームデータ解析 v1.0」のとおりには動かなかったので、ここに今日の時点で動いたコマンドを書き記しておく。