アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
誰でもわかる(と思う!)LTR解析 その2
LTR解析について連載するとか言っておきながら何ヶ月たったことか・・・
どうもmisawatです。大丈夫です。生きてます。

さて、期間が大分空いてしまいましたが、第二回として
「LTRの研究事例」
について紹介していきたいと思います。

やはり「どんなことがわかるのか!」ってのが知りたいとこですよね。

さて、少しおさらいですがLTRはレトロトランスポゾンのRNAの両端に存在する繰り返し配列のことでした。
さらに、この領域はプロモーターとしての機能を有するというのも前回お話したとおりです。
これらの前知識を持った上で、実際にpublishされた論文を紹介していこうと思います。

まず一報目。
CAGE profiling of ncRNAs in hepatocellular carcinoma reveals widespread activation of retroviral LTR promoters in virus-induced tumors
Kosuke Hashimoto et al., Genome Res, 2015, 25, 1812-1824

主にnon-coading RNA(以下ncRNA)の研究をされているグループが2015年にGenome Research誌に出した論文です。
背景としてncRNAとガンとの関係性が日々とりざだされていますが、ガンとは多様なもので、その部位、癌腫等によって全く違うメカニズムを有していたりします。
よって、「〜ガンとncRNAの関係」の「〜」の部分が大変重要なんですね。

この研究では肝細胞癌とncRNAの関係性について調べています。
肝細胞癌はB型肝炎ウイルスやC型肝炎ウイルスによる肝炎からの発症が多い病気です。
その他の原因としてはアルコールとかですかね。

聡明な皆様ならもうLTRとの関連がわかったかもしれませんね。
そう、この論文のリザルトの中には
「肝細胞癌においてLTRプロモーターが活性化されていた」
「LTRに由来するncRNA発現増加量は通常の10倍以上であった」
というのがあります。

LTRの研究は癌研究にもつながっていくんですねぇ。
ちなみにこの論文はCAGE(Cap Analysis of Gene Expression)法を用いた肝細胞癌の転写開始地点及び発現増強部位の特定を行うところがからスタートしています。
手法としてもとてもためになる論文なので一度読んでみると良いかもしれないですね。



次に二報目。
MMTV/LTR Promoter-Driven Transgenic Expression of EpCAM Leads to the Development of Large Pancreatic Islets
Jeffrey R. Vercollone et al., J Histochem Cytochem August 2015 vol. 63 no. 8 613-625

先ほどは肝臓、今回は膵臓のお話ですね。
ここではガンの話ではなく、膵臓ランゲルハンス細胞の接着、分化、成長について調べています。
膵臓ランゲルハンス細胞の接着、分化、成長はEpCAMという遺伝子によって調節されているらしいです。

この遺伝子はマウスにおいてMMTV(Mouse mammary tumor virus)のLTRの支配下にあるんだそうです。
ガンとかではなく、細胞の成長分化にまで関与してしまうLTR・・・こわいですね!!

論文の内容としてはin vivoでEpCAMのValidationを行いました!という論文のようです。
in vivo・・・大変ですよねぇ。



最後に三報目。
結構長くなってきているので巻きで行きます。
Quantification of Total and 2-LTR (Long terminal repeat) HIV DNA, HIV RNA and Herpesvirus DNA in PBMCs
Massanella, M. et al., BioProtoc., 2015, 5(11), e1492.

これは論文というか、手法紹介に近い論文?ですね。

HIVにかかると、日和見感染と呼ばれる「普段はなんともないようなウイルスや細菌に感染する」現象が起こります。
日和見感染の代表的なウイルスにCMV(Cytomegalovirus)やEBV(Epstein‐Barr virus)などがあげられますね。
医学系の教育を受けた方なら「EBV」=「伝染性単核症」といった関連付けがされるのではないでしょうか。

この論文では、HIVの潜在的な貯留量を計測し、CMVやEBVとの関係性を見る手法を紹介しています。
定量自体はdroplet digital PCRという手法を用いています。
ターゲットはもちろん末梢血の単核球分画です。

「LTRの有無」ではなく、定量情報で議論できるのはいいですね。


さてさて、今回はLTRに関する論文を三報程紹介してみました。
やはりウイルスが絡むと研究の色が一気に医学よりになりますね。
それだけ応用が期待される研究分野だということなのかもしれません。


次回は・・・、いよいよ分析とかの話にはいりたいですね〜。
果たしていい公開データや紹介できそうな優良な分析ツールはあるのか!
乞うご期待です。。。




いや、やっぱりそんな期待しないでくだs・・・






〜つぶやき〜
最近、愛車にアルミテープを貼りました。
こういうオカルトチック(?)なテクノロジーが大好きです(生命科学においても)。
| misawat | バイオインフォマティクス | 17:54 | comments(0) | - |
データ解析トレンド
現在、アメリエフはNGSデータを中心に解析を行っています。
しかし、もともと起業した社長は遺伝統計学のバックグラウンドがあり、最初はSNPアレイデータ解析の会社でした(そもそも、創業当時はNGSが今ほど普及していませんでした)。

とはいえ、ここ数年ずっとNGSデータ解析のお問い合わせが主でした。
ところが、今年に入ったころくらいから、GWAS解析、CNV解析などのSNPデータ解析のお問い合わせを立て続けにいただきました。
NGSは網羅性、データ量では勝りますが、安価さと精度ではSNPアレイのほうがまだ勝っていることもあり、大規模な集団の解析ではSNPアレイの需要は高いということでしょうか。
東北大ToMMoと東芝が開発した日本人SNPを反映したSNPアレイも販売されており、SNPアレイの活用はますます進むのではないかと思います。

アメリエフでは、SNPデータ解析受託、SNPアレイデータ解析で一般的に使われるソフトウェアをプリインストールしたサーバの販売、ソフトウェアを使用するためのLinuxコマンド基礎を含めたトレーニングの提供を行っております。
ご興味のある方は、弊社HPお問い合わせフォームよりお問い合わせください。
| kubo | バイオインフォマティクス | 12:12 | comments(0) | - |
誰でもわかる(と思う!)LTR解析 その1
〜前置き〜

だんまりを決め込んでいました。
どうもmisawatです。

たまにはInputだけではなくOutputせねば!と謎の義務感にかられたので記事を書きにきました。

テーマは「LTR (Long Terminal Repeat)」の解析です。

・・・。

はい。


「繰り返し〜」とか「リピート〜」とか言うと遺伝子解析屋さんは「むぇ〜(´・ω・`)」ってなるかもしれませんね。


ただこのLTR、解析出来ると色々おもしろいことがわかってきます。
そこで、「LTRとは何か」から始まり「LTRの解析方法」までを解析する記事を連載していこうと思います。


今回は第一回として「LTRとは何か」を書いていきます。





〜本編〜

さて、LTRとは上述したとおりLong Terminal Repeatの頭文字をとった単語です。

LTRはレトロウイルスやLTR型レトロトランスポゾンのRNAの両端に存在する数百〜数千の長い長〜い繰り返し配列です。特徴的ですね。



・・・特徴的ですね!!


そう"特徴的"なんです。

つまり宿主への感染を感会えた時に、このLTRを目印に感染が確認出来ます。

更にこのLTRにはプロモータとしての機能があったりもするのでウイルス関連の研究をされている方には魅力的な研究対象かもしれませんね!


次回はLTRに関する研究について書こうと思います。


以上、misawatがお送りしました!





〜つぶやき〜
「頭文字」って単語を見ると無性に峠を攻めたくなるのは私だけじゃないはず・・・。
| misawat | バイオインフォマティクス | 14:44 | comments(0) | - |
ヒトデータのセキュリティ
今夏は慎重に行動していたおかげでほとんど蚊に刺されていなかったのですが、お盆休みで実家に帰ったところ、寝ている間に10箇所も刺されてしまいました。
せっかく気を付けていたのに、台無しです。

無理やりこじつけますが、気を付けるといえば、ヒトのシーケンスデータを扱うときに、セキュリティはどれくらい気にするものなのか、ご質問をたびたびいただきます。
弊社がヒトシーケンスデータを扱う際のセキュリティルールは、NBDCヒトデータ取扱いセキュリティガイドラインを参考にしています。

NBDCヒトデータ取扱いセキュリティガイドライン
こちらは、NBDCで提供しているNBDCヒトデータベースの共有データを利用する際のセキュリティルールです。
セキュリティレベルには標準レベル(Type-I)とハイレベル(Type-II)があり、弊社では標準レベルの基準と同等のセキュリティを適用したサーバでヒトデータの解析を行っています。

ひとつの指針として、ご参考になりましたら幸いです。
| kubo | バイオインフォマティクス | 15:14 | comments(0) | - |
HISAT2
RNA-seqのアライメントツールHISAT2についてご紹介します。

HISAT2は、当ブログ内でもたびたびご紹介してるTopHat2開発グループが作った、TopHat2の後継ソフトです。
TopHat2より速いしメモリもそんなに使わないとアピールしています。

ゲノムに対して小さなインデックスを細かく、たくさん作成することで効率的なマッピングをするということだそうです。

私はLinux環境で使うので、Linux用バイナリをダウンロードしてきて使ってみました。バージョンは2.0.4でした。
zipファイルを解凍したら、新しく作成されたディレクトリ内の"hisat2"、"hisat2-build"、"hisat2-inspect"をパスの通ったところに追加します。Bowtie2と同じ手順です。私の場合は/usr/local/binです。
$ cd hisat2-2.0.4
$ sudo ln -s /path/to/hisat2 /usr/local/bin/
(以下同)

本体の次はゲノムとインデックスを準備します。
手元のリファレンスゲノムからhisat2-buildで自分でビルドすることもできるようですが、GRCh38やhg19、その他メジャーなゲノムについては、インデックスを一緒に配布してくれているので、そこからダウンロードすればいいですね。既知のtranscriptやSNPを使用してビルドしたインデックスもあり、そちらを使用したほうがうまくマッピングできるようです。

あとはマッピングするだけです。

ところで、名前が紳士服と関係なくなってしまったのが少し残念ですね。
最初のころ、RNA-seq関連のソフトウェアが紳士服と関係あることにしばらく気づいていなかった社員のつぶやきでした。
| kubo | バイオインフォマティクス | 14:44 | comments(0) | - |
10x genomicsのGemCode技術

去年、The Scientist誌が選んだ革新的技術トップ10(Top 10 Innovations 2015)の1位を飾った技術、それが、スタートアップ企業10x genomicsが開発した、「GemCode」と呼ばれる次世代シーケンサーの新しい前処理技術です。

AAEAAQAAAAAAAAPZAAAAJGM4NzJlYjM0LTU3NDktNDg2MS05ZTI4LWE2ZDY4OGVmYWE1Zg.png (7.9 kB)

GemCodeの技術を簡単に説明すると、 Illuminaシーケンサー等によって得られるショートリードから、合成的にロングリード(Linked readsと呼ばれる)を生成するというものです。従来のIlluminaの次世代シーケンサーが苦手としていた「大きな欠失・挿入部位の解析」、「転座・融合部位の解析」や「HLAなどの高度な多型性を持つDNA領域の解析」を精度よく行うことができます。


方法としては、
(1)それぞれ異なるバーコード配列が付加されたゲルビーズのライブラリの液滴
(2)試薬と混合させたナノリットル(nl)のスケールのDNA分子の油液
をマイクロ流路を使って混合させます。(1)と(2)が合わさることで、10-100kbほどのDNA分子が断片化され、各DNA配列に同一のバーコード配列がライゲーションされます。このようにして、バーコード配列をライゲーションさせたDNA配列をIlluminaのシーケンサーで読みます。


このとき、同一のバーコードを持つDNA配列(シーケンスリード)群は、
(1)ゲノム上で近接している
(2)(母方もしくは父方の)同一ハプロタイプに由来している。
ことが想定されます。そのため、シーケンスされたリードをゲノムにマッピングし、local de novo assemblyを行うことで、合成的にに長いDNA配列(10-100kb前後)を構築することができるというわけです。 GemCodeの技術について、こちらの動画がわかりやすいので見てみてください。

Changing the Definition of Sequencing from 10X Genomics on Vimeo.


「でもお高いんでしょう?」


そんなことありません! 10x Genomicsによると、従来のIlluminaのWGSやExomeと同様のリード数(Coverage)で解析が可能とのことです。つまり、ランニングコストとしては通常のIlluminaのシーケンスと同じぐらいと考えてよさそうです。


さらに、ゲノムDNAは「1ng」からシーケンスできます。 現在のIlluminaの手法と比較しても、かなり少量でも試せます。
比較対象:http://dnatech.genomecenter.ucdavis.edu/sample-requirements/


Illuminaのショートリードの弱点を補完するGemCode技術。 様々な解析での活用が期待できると思います。


アメリエフでも、10x genomicsの技術により得られたシークエンスデータの解析を試しているところです。興味はあるけどシーケンスされたデータをどうやって解析したらいいかわからないという方は、一度アメリエフに相談してみてはいかがでしょうか。


執筆者自己紹介・インターンシップ感想

はじめまして、アメリエフにインターンに来ているimamachi-nです。 アメリエフでは、融合遺伝子の解析レポートや10x Genomicsの調査などを行ってきました。


私は大学では実験を行いつつ、Linux、PythonやRなどを利用してNGS解析も行っています。ただ今まで、PythonやRのコーディングに関しては独学でやっていたので、恥ずかしながら適当なコードを書いていました。アメリエフでは、書いたスクリプトを社員の方に繰り返し確認していただき、修正すべき点をチェックしてもらいました。また、Pythonでコーディング規約に即した書き方をしているかflake8を使いチェックを行う方法も教えてもらいました。今まで適当に書いていたスクリプトを見なおすことで、無駄のない効率的なコードや、あとで見直した時に可読性の高いコードを書く方法を学ぶことができました。


社内で使用していたSlackは、今では私の研究室でも試験的に導入しています。NGS解析について結果やデータの共有、ディスカッションを行う際にSlackを活用しています。今まではメールでやり取りしていたのですが、Slackを使うことで今まで議論してきたことなどがログとして残るので便利です。


社員の方と一緒に仕事ができるという点で、非常に充実したインターンシップだと感じました。 NGSの解析に興味がある方にとっては、アメリエフのインターンシップは最適だと思います。もちろん、アメリエフで働きたいという方も会社の雰囲気をつかむ上でインターンシップに参加するメリットは大きいと思います。



  • 社員一人ひとりのモチベーションが高く、新しいことにどんどんチャレンジしていく風土がある。

  • 新しいNGSの技術が出てくる中で、迅速にそれらに対応し、バイオインフォマティクスの専門家として適切な解析方法を提案できる。

  • ウェットの研究者が抱えるデータ解析上での課題に対して、綿密に対応してくれる。


それがアメリエフという会社だとインターンシップを通して感じました。

| 管理者 | バイオインフォマティクス | 14:46 | comments(0) | - |
Vagrantを用いた仮想環境構築

バイオインフォマティクス初心者が勉強を始める際に、最初に立ちはだかる壁は、Linuxという未知のOSとの相対だと思います。 最近はVirtualBoxなどの仮想化ソフトもネットにたくさんありますので、導入すること自体は皆さんできるかと思います。しかし、もし...

  • パソコンを新しくしたから、昔と同じ仮想環境を新しいパソコンでも構築したい
  • 複数の端末で同じ仮想環境を使用したい
  • 仮想環境のバックアップをとっておきたい

ということを思った時に、仮想化ソフト単品ではなかなか難しいと思います。だからといって、仮想環境内で色々なアプリを入れた後に、全く同じ環境を他の端末で作るのは、普通にやっていたらかなり時間を取られてしまいますし、ツールのアップデートが知らない間にされていて、微妙に挙動が合わなかったりします。

そのようなときに便利なのが仮想化ソフトのラッパーであるVagrantです。

Vagrantでできること

VagrantはVirtualBoxなどの仮想化ソフトに追加機能を与えてくれるアプリ、と思って頂けたら良いと思います。主な特徴として

  1. 構築環境を記述した設定ファイルを元に仮想環境の構築から設定を自動化できる
  2. 一度構築した仮想環境をboxというパッケージにエクスポートすることができる

があるのですが、要は一度設定した構築環境を複数の端末で再現するのにうってつけというわけです。 もちろんCLIも有りますので、サーバーのターミナル上で仮想環境を構築したい場合にも使えます。開発ツールのテストとかにも使えるので、初心者以外でも有用です。

使い方は色々な方がネットで説明されていますので割愛・・・。 今まで仮想化ソフト単品で仮想環境を構築再構築に明け暮れていた方は、Vagrantで自動環境構築にチャレンジしてみましょう!

執筆者自己紹介・インターンシップ感想

はじめまして、アメリエフにインターンに来ているnomatです。 アメリエフでは、ツールのテスト・デバッグ・マニュアル作成や、研究調査に携わっております。

私は大学では実験メインで活動してきた学生なものですから、インターン開始直後はTerminalで操作することすら覚束ない初心者丸出しの状態でした。しかし、アメリエフの解析担当者の皆様にご指導いただき、今ではバリバリコマンドを打って・・・いる気がします(多分)。 この手の操作方法やテクニックはウェットメインの学生でも有用なところが多いですし、最近は解析用OSSも豊富にありますので、是非とも学部時代とかに教えるようにしていただきたいですね。

アメリエフでは2ヶ月勤務させていただきましたが、とてもオープンな雰囲気で、充実した時間を過ごすことができました。勤務内容に関してはもちろん、日頃疑問に思っていたことを解析担当者の方に質問すると、スパっと答えてくださるので、色々と勉強させていただきました。 バイオインフォマティクスを始めてみようかなと考えるウェット系のラボの学生の多くは、周りにインフォの知識を持つ人が少ないため、独学でやっていると思います(私もそうでした)。しかし、やはり知識と技術を持つ人に聞きながら勉強すると、上達が早いです。 何より、情報系の人が当たり前のようにやっている「常識的な事」というのは、案外ネットの情報だけではつかみにくいものが多いです(ディレクトリの配置はこうしたら効率的、など)。これらを勉強できたのは私にとって大きな収穫でした。

最近のバイオ研究ではインフォマティクスは必須になりつつあります。これから研究に携わる学部生も、ウェットしかやってこなかった院生も、一度バイオインフォマティクスに触れてみましょう! その際にアメリエフのインターン、おすすめですよ!

| 管理者 | バイオインフォマティクス | 14:15 | comments(0) | - |
Cytoscapeによるネットワーク図示

Cytoscapeは、複雑なネットワークおよびその属性の図示、統合、分析に用いられるオープンソースのソフトウェアです。 遺伝子ネットワーク、ソーシャルネットワーク、路線図など、点(node)と線(edge)で構成されるデータセットを可視化することができます。 データの可視化によって、全体像や何らかの傾向が把握でき、そのデータが意味するところを理解する手助けとなる可能性があります。 Cytospaceには、プラグインを追加することによって、化合物を扱う(ChemViz)、外部のパスウェイデータベースであるReactomeを利用する(ReactomeFIViz)というように、機能を大幅に拡張できるという特徴があります。

Cytoscapeは様々な種類のネットワークを記述したフォーマットを読み込むことができます。 一番単純なものは、線の起点と到着点を示す”source”、"target"、それから点と点を結ぶ線の種類を示す”interaction type”の3つの列から構成されるSIF(Simple Interaction File)フォーマットです。 下は、ソフトウェアのサンプルデータ”galFilterd.sif”を読み込んだものです。 galfiltered_1.png (46.0 kB)

また、レイアウトやスタイルを変えることで、図の印象を変えることができます。 galfiltered_2.png (25.5 kB) galfiltered_3.png (68.6 kB)

点(node)の持つ情報(attribute)をもとにして、スタイルを書き換えることもできます。 以下は、サンプルデータ中の”galFiltered.cys”を図示したものです。 ノードの発現量が"色"に、ノードの持つ離散値の属性が"サイズ"に対応しています。 galfiltered_4.png (110.7 kB)

作成したネットワーク図を、ウェブブラウザから扱えるインタラクティブな図として出力することができます。 f8402978-66f4-4dc8-9ed1-eeb4059b7cc3.png (137.1 kB)

執筆者自己紹介・インターンシップ感想

はじめまして。 3月からインターンシップで勤務しております、nakamurahと申します。 アメリエフでは2か月間、BEDフォーマットチェックスクリプトの作成、Cytoscapeによるネットワーク図示、Python+NetworkXを用いたネットワーク図示について、関わらせて頂きました。 お蔭様で、普段の研究室生活ではめったに使わないLinuxやPythonに、慣れてきたように思います。 この知識を研究室の先生や学生に教えるなど、研究室の文明化に貢献でき始めているところです。

アメリエフでの勤務の感想として、ミーティングの時に見られるように、皆さん「目的」と「手法」、「利点と欠点」をはっきりさせて、論理的に議論される印象を受けました。 一方で、バイオインフォマティクスの分野は、次々と新しい手法が開発され、既存の手法でもバージョンが更新されて対応が必要になるように、常に最前線に立ち続けることが非常な困難を伴う分野であると感じました。

私が普段扱っている非モデル植物の研究は、最先端の技術が使われる生物から1回りも2回りも遅れています。 この状況を改善するために、これからも技術の向上に邁進していきたいと思います。

| 管理者 | バイオインフォマティクス | 14:35 | comments(0) | - |
10x Genomics Long Rangerのご紹介
前回のブログでも少し触れましたが、10x Genomicsが開発したGemCodeシステムは、ショートリードから擬似的にロングリードを生成する革新的な技術です。

今回はそのGemCodeシステムに対応したゲノム解析パイプラインLong Rangerのご紹介をしたいと思います。

今年の2月にGemCodeシステムの新機種Chromiumが発表され、先月Chromiumに対応したLong Rangerのバージョン2がリリースされました!
Long Rangerは、Whole GenomeおよびExome Sequencingに対応しています。
GemCodeシステムではilluminaのシーケンサーを使うため、BCLファイルが出力されます。
BCLファイルからFASTQへの変換は以下のコマンドを実行します。
$ longranger demux --run=/path/to/BCL/output

内部では、illuminaのbal2fastqが動いているため、あらかじめbal2fastqをインストールしておく必要があります。
出力は、リードとバーコードのFASTQファイルが分かれて出力されます。

次にこのFASTQファイルを入力として、以下のコマンドを実行すると、マッピングから変異検出、ハプロタイプフェージング、構造変異の検出を行ってくれます。
$ longranger run --id=sampleID --sex=female --fastqs=/path/to/fastqs --reference=/refdata-hg19-2.0.0

Long Rangerの変異検出はデフォルトでは、freebaysというソフトウェアを用いていますが、--vcmodeオプションでGATK(v2.4以上)で実行することができます。

Chromiumを使って調整したサンプルのシーケンスデータは、10xGenomicsのサイトで公開されています。

実際にWGSのデータ(FASTQファイル)を実行してみたところ、40時間ほどで実行が終了しました。解析には、Intelの16-coreのプロセッサー、256GBのメモリを搭載したマシンを用いました。

出力結果には、フェージング後のBAM、VCFファイルと構造変異のBEDファイルに加えて、専用のゲノムブラウザLoupeに用いる.loupeという形式のファイルも出力されます。

ゲノム解析パイプラインのLong Ranger以外にも、シングルセルRNA-seq解析パイプラインのCell Rangerやde novoアセンブリ用のSupernovaが公開されています。
| onouek | バイオインフォマティクス | 12:59 | comments(0) | - |
Yet Another Bioinformatics Library
バイオ分野で使われているプログラミング言語には、バイオインフォマティクス向けのライブラリが用意されており、BioPerlBioPythonなどをご存知の方は多いかもしれません。

最近注目されているGoogleによって開発されたGO言語にも、バイオインフォマティクス用のライブラリであるbiogoというものがあります。

今回はそのbiogoを使用している例として、Lariatというアライメントツールのご紹介をしたいと思います。

Lariatは10x Genomicsが開発したGemCodeシステムに対応したアライメントツールです。GemCodeシステムはバーコード配列を利用してショートリードから合成的にロングリードを生成する革新的なシステムです。

入力はFASTQ-likeなフォーマットで、入力ファイルに以下の情報が必要になります。
read header
read1 sequence
read1 quals
read2 sequence
read2 quals
10X barcode string
10X barcode quals
sample index sequence
sample index quals

Lariatの最初のステップでは、BWAのAPIを使ってアライメントを行っていきます。

次のステップでバーコードの情報を使って反復領域などのマッピングが難しい領域のリードを繋げていき、最終的なマッピングのポジションなどを決定します。
これにより、segmental duplicationなどの反復領域などへのリードのマッピングの正確性が向上すると考えられています。

ちなみに、biogoはBAMやSAMファイルのハンドリングに使われています。

Lariatの入力ファイルはFASTQ-likeなフォーマットが必要と書きましたが、Lariat自体はLong Rangerというパイプラインに組みこまれているので、自分で用意することはないと思います。
| onouek | バイオインフォマティクス | 14:02 | comments(0) | - |
STARとCufflinks
RNA-seq解析におけるマッピングソフトウェアの選択肢は、TopHat一強から、STARもずいぶん多く使われるようになってきました。
性能や速度から、STAR一択! と言い切る人もいますが、まだ根強くTopHatユーザもいるという印象です。

マッピングは解析の中でも、非常に重要ですが、実行時間がボトルネックということも多いので、マッピングソフトウェアの速度や精度は非常に気になるところですね。

STARでマッピングした後は、遺伝子発現解析を行うと思いますが、Cufflinksを使用される際は、「--outSAMstrandField intronMotif」オプションでXSタグを付ける必要があります。
※STARのバージョンは2.5.1bです。

■コマンド例

$ STAR --genomeDir /path/to/STARIndex/ --runThreadN 3 --sjdbGTFfile genes.gtf --outSAMtype BAM SortedByCoordinate --outFileNamePrefix samplename --outSAMstrandField intronMotif --readFilesCommand zcat --readFilesIn sample_1.fastq.gz sample_2.fastq.gz

■その他のオプションの説明


--outSAMtype BAM SortedByCoordinate
デフォルトでは出力ファイル形式がSAMなので、BAMファイルで出力します。SortedByCoordinateと指定しているので、染色体のポジションでソートされて出力されます。
--readFilesCommand zcat
入力fastqファイルをSTARで処理する前に実行するコマンドを指定できます。今回使用しているfastqファイルはgz圧縮しているので、zcatでファイルを解凍した結果をSTARで処理しています。
| kubo | バイオインフォマティクス | 14:26 | comments(0) | - |
AWKでフィルタリング
 以前の記事(TopHat-Fusionの結果の見方)で、
「Linuxコマンドのawkを使ってfusions.outをフィルタリングする方法もあります」と書いたので、
今回はawkコマンドのご紹介をしたいと思います。

AWKコマンドは主に、テキストファイルから要素を抜き出したりするのに便利なコマンドです。

基本形

$ awk '条件文{実行文}' ファイル名

では、実際にfusions.outの7列目までを用いて、フィルタリングの例を示してみたいと思います。
chr20-chr17 49411707 59430946 ff 9 3 9
chr20-chr17 49411707 59445685 ff 106 116 167
chr20-chr20 47538545 46365686 fr 17 11 9
chr17-chr17 57992061 57917126 ff 4 2 1
chr20-chr20 46415145 52210297 rf 22 18 27
chr2-chr17 142237963 37265642 rr 2 3 2

(1) 融合遺伝子の染色体番号
(2) 融合遺伝子の左側の遺伝子上のポジション
(3) 融合遺伝子の右側の遺伝子上のポジション
(4) 左側と右側の遺伝子の向き(f:forward, r:reverse)
(5) breakpoint上のリード数
(6) breakpointを挟むペアエンドのリード数
(7) 片側のリードがbreakpoint上にあるペアエンドのリード数

たとえば、5列目のbreakpoint上のリード数でフィルタリングを行いたい場合、
$ awk '$5>100{print}' fusions.out > fusions_filtered.out

とすると、fusions_filtered.outにbreakpoint上のリード数が100より多い行が出力されます。
chr20-chr17 49411707 59445685 ff 106 116 167

条件文ではなく、実行文にifを使っても同様の処理が可能です。
$ awk '{if($5>100) print}' fusion.out > fusions_filtered.out

このように、awkはフィルタリングにも使える便利なLinuxコマンドですので、出力結果の絞り込みを行うときなどさまざまな場面で役に立ちます。

よく使うコマンドですが、意外にもこのブログで書かれていなかったので書いてみました。
| onouek | バイオインフォマティクス | 14:08 | comments(0) | - |
snpEffのアノテーション書式
たびたびブログでもご紹介している、アノテーションソフトsnpEff、一年前の9月にバージョン4.0になったとご紹介しましたが、2015年1月現在のバージョンは4.2です。

snpEffで付与されるアノテーションの書式は、以前は下記のような「EFF=...」でした。
##INFO=<ID=EFF,Number=.,Type=String,Description="Predicted effects for this variant.Format: 'Effect ( Effect_Impact | Functional_Class | Codon_Change | Amino_Acid_change| Amino_Acid_length | Gene_Name | Gene_BioType | Coding | Transcript | Exon | GenotypeNum [ | ERRORS | WARNINGS ] )' ">
現在のバージョンでは、デフォルトでは「ANN=...」に変更しています。
##INFO=<ID=ANN,Number=.,Type=String,Description="Functional annotations: 'Allele | Annotation | Annotation_Impact | Gene_Name | Gene_ID | Feature_Type | Feature_ID | Transcript_BioType | Rank | HGVS.c | HGVS.p | cDNA.pos / cDNA.length | CDS.pos / CDS.length | AA.pos / AA.length | Distance | ERRORS / WARNINGS / INFO' ">
バージョン4.1でも、実行時に「-classic」オプションを使用することで、「EFF=...」の書式でアノテーションすることも可能です。古いバージョンと互換性を残してくれるのはうれしいですね。

ANN書式が旧書式に比べ発展している点として、連続する変異を合わせた影響を考慮するなど、便利なアノテーションが増えているようなので、今後はぜひともANN書式を使いこなしていきたいです。

書式の詳細については公式の説明をご覧ください。
| kubo | バイオインフォマティクス | 14:26 | comments(0) | - |
ファイルのコピーに失敗するときは転送速度を落としてみる
 こんにちは、久保(kubor)です。
 少し前からオフィスの近くで大阪風のお好み焼きを提供するお店を探しているのですが、どうも広島風の赤いお店が多い気がします。

外付けHDDからのコピーに失敗することがある

 ところで、弊社では日々お客様からのデータを受け取り、解析をしています。その際、データを受け取ると、まずは解析サーバにデータをコピーします。サイズはまちまちですが、ゲノムデータなどの大きいデータのコピーはよく失敗してしまいます。
 特に、USB3.0が使えるなど、高速転送できるときに遭遇する気がします。
 そんな時は、せっかくの高速転送ですが、諦めて転送速度を落としてあげると正常にコピーできることが多いです。

rsyncで転送速度を落としてコピー

 Linux上でのデータのコピーは、cpコマンドを使うことが多いと思います。ですが、少し複雑なことをしようとするとrsyncを使うことになります。基本的な使い方はcpと同じですが、オプションが豊富に用意されています。
例えば、データコピー時に権限や所有者、タイムスタンプなどをそのままにしてコピーするときに使う-aオプションは、-rlptgoDとオプション指定した場合と同義になります。つまりこれらを分割することもできるということです。僕は、コピー状況を確認するために-vをつけたり、-hを付けて読みやすい単位でファイルサイズを表示させたりして使用しています。
 さて、表題の転送速度を落とすオプションは以下です。
--bwlimit:転送速度(KB/s)を指定
 コピーに失敗するときには、このオプションで調節します。例えば、--bwlimit 10000とすれば、10MB/sで転送されます。
 /mnt/hdd/sampleA_R1.fastqがあったとしてカレントディレクトリへのコピーは以下のようにします。
$ rsync --avh --bwlimit 10000 /mnt/hdd/sampleA_R1.fastq ./
 コピーに失敗する際はお試し下さい。
| kubor | バイオインフォマティクス | 14:48 | comments(0) | - |
Chromeで開いてる全てのタブのタイトルとURLをMarkdown形式で取得する
こんにちは、雪のようなくちどけを経験できなかった青春時代は、アニメや漫画で上書きすればいいと思っている久保(kubor)です。
メルティーキッス くちどけラム&レーズンすごく美味しい!

いつのまにかChromeのタブが増える

さて、バイオインフォマティシャンにとって日頃使用するソフトたちの情報収集は欠かせない業務の1つです。様々な研究所、研究者たちがGitHubや彼らのWEBサイトにて情報を公開していますので、一重に情報が集約されていることは大変少ないです。
そういった理由もあり、調べ物をしているといつの間にかChromのタブが大量に開かれることになりがちです。

調査をある程度進めた後は、開いたページを参考文献として保存するのですが、 各タブのURLをコピーするのはかなり骨が折れます。
そんな時はChrome拡張機能の「Copy All URLs」を使うと幸せになれます。

Copy All URLsで便利にURLをコピーする

Copy All URLsのダウンロードとインストールは以下から行えます。
Copy All URLs | Chrome ウェブストア



インストールしたらブラウザの右上に傘のマークが出てきますので、こちらからCopyやPasteを選ぶことができます。
デフォルトではURLを取得するだけですが、Optionsから設定を加える事ができます。

Markdownでページタイトルと一緒に取ってくる

弊社ではテキストのフォーマットに制約を設けてはいませんので、各自好きなフォーマットを利用しています。僕は、Markdown派なので、URLもMarkdownで取得したいです。



ついでにページタイトルも取得しましよう。
Options > Format > Customを選択すると入力フォームが出てきます。
こちらに以下のテキストを入力します。
※空白行が無いと複数のURLが改行無しで取得されてしまいます。
[$title]($url)
これで完了です。



それでは早速タブをたくさん開いてみたので、URLを取得してみます。
[samtools/samtools](https://github.com/samtools/samtools)
[arq5x/bedtools2](https://github.com/arq5x/bedtools2)
[taoliu/MACS](https://github.com/taoliu/MACS/)
[apple/swift](https://github.com/apple/swift)
[pydata/pandas](https://github.com/pydata/pandas)
[hadley/ggplot2](https://github.com/hadley/ggplot2)
[tensorflow/tensorflow](https://github.com/tensorflow/tensorflow)
[アメリエフのブログ](http://blog.amelieff.jp/)
[ゲノム解析ならアメリエフ株式会社](http://amelieff.jp/)
上手く取得できましたね。
他にも好きなフォーマットにすることができますので、良ければお試し下さい。
| kubor | バイオインフォマティクス | 16:54 | comments(0) | - |
バイオインフォマティシャンでも、pecoがしたい!
先日、ついに日比谷線を利用しました、久保(kubor)です。
先日pecoが話題になっていた傍ら、社内ではなかなか広まらずに、寂しい思いをしているので、バイオインフォマティシャンにもおすすめの使い方を紹介いたします。

0. pecoってなに

標準入力で受け取ったテキストを、インクリメンタルサーチすることができるコマンドラインツールです。
更に、サーチ結果は、標準出力されます。

すなわち、テキストフィルタリング出力となります。

ちなみにpecoはpercolというソフトをGolangで書いたものです。
同様の機能を持つソフトとして、他にもzawなども有名ですが、数年前にpecoが高速だと話題になってからというもの、pecoの人気が突出しているように感じます。
以下が、各リポジトリです。

1. pecoのインストール

バイナリを落とすのが早いです。
以下からOSにあったものをダウンロードして、~/binディレクトリなどのPATHが通っているディレクトリにシンボリックリンクを作成します。
https://github.com/peco/peco/releases

2. 基本的な使い方

例えば、cat sample.txt | peco として、パイプでpecoを使えば、ウィンドウが開きsample.txtの中身をインクリメンタルサーチすることができます。
サーチ結果は、カーソルキーで選択することもできます。

これだけだとインクリメンタルサーチができるgrepですが、pecoの入出力をパイプすることで様々なことが可能になります。

3. スニペットを呼び出す

peco-snippet-loader

シンプルで、わかりやすい使い方として、スニペットをpecoで呼び出す方法を紹介します。zshを基準として書いていますので、bashの場合は適宜変更して関数を用意して下さい。

まず、解析によく使う便利コマンドや、ワンライナーをスニペットとして保存しておきます。
下記のように、ホームディレクトリにドットファイル(ファイル名の先頭にドットを入れると隠しファイルとなる)を作成するのがおすすめです。

~/.peco.snippet.bioinfo
#
# file: .peco.snippet.bioinfo
#

# grep
[remove header] grep -v "^#"
# perl
[print TSV header]perl -e 'for(split(/¥t/, )){$i++; print "$i $_¥n"}'
# samtools
[bam to sam] samtools view
[sam to bam] samtools view -bS
[sort bam] samtools sort
[make fasta index] samtools faidx
# sort
[sort BED file] sort -V -k1 -k2n

次に、.zshrcに以下を追記します。
保存後、`source ~/.zshrc`を実行した後、Ctrl+xと入力するとスニペットを選択する画面が起動します。
function peco-snippets-loader() {
if ls ~/.peco.snippet* >/dev/null 2>&1; then
snippet=`cat ~/.peco.snippet* | grep -v "^#" | peco`
BUFFER="$(echo $snippet | sed -e 's/^¥[.*¥] *//') "
CURSOR=$#BUFFER
else
echo "~/.peco.snippet* is not found."
fi
zle reset-prompt
}
zle -N peco-snippets-loader
if type peco >/dev/null 2>&1; then
bindkey '^x' peco-snippets-loader
fi
ポイントは`cat ~/.peco.snippet* | grep -v "^#" | peco`で、スニペットからコメント行(#から始まる行)を除外して、pecoに渡しています。
そして最後に選択したスニペット行頭の説明文をsedで取り除いています。

とっつきにくいですが使ってみると楽しいので、是非お試し下さい。
僕は、この他にもよく使うディレクトリへのショートカットや、ヒストリー検索をpecoで選択できるようにしています。
| kubor | バイオインフォマティクス | 16:06 | comments(0) | - |
TopHat-Fusionの結果の見方
TopHat-Fusionは、マッピングソフトのTopHatから派生した機能で、融合遺伝子を検出するのに広く使われています。
TopHat-Fusionを実行すると、fusions.outという候補となるポジションやリード数などの情報が出力されます。fusion.outの説明は、NGS Surfer's Wikiに書かれていますので、詳しくはこちらをご覧ください。
出力されたファイルに対してtophat-fusion-postというコマンドを実行すると、フィルタリングされた結果(result.txt)が作られます。

このresult.txtの見方を、BCAS4-BCAS3の検出を例に説明したいと思います。

MCF7 | BCAS4 | chr20 | 49411707 | BCAS3 | chr17 | 59445685 | 106 | 116 | 167

(1) MCF7 : 融合遺伝子が検出されたサンプル名
(2) BCAS4 : 融合遺伝子の左側の遺伝子名
(3) chr20 : 左側の遺伝子がある染色体番号
(4) 49411707 : 左側の遺伝子のポジション
(5) BCAS3 : 融合遺伝子の右側の遺伝子名
(6) chr17 : 右側の遺伝子がある染色体番号
(7) 59445685 : 右側の遺伝子のポジション
(8) 106 : breakpoint上のリード数
(9) 116 : breakpointを挟むペアエンドのリード数
(10) 167 : 片側のリードがbreakpoint上にあるペアエンドのリード数

以上が、フィルタリング後のresult.txtの列の説明でした。
参考:http://ccb.jhu.edu/software/tophat/data/result.html

Tophat-fusion-postによるフィルタリングは便利ですが、Linuxコマンドのawkを使ってfusions.outをフィルタリングする方法もあります。
| onouek | バイオインフォマティクス | 14:20 | comments(0) | - |
Pythonを書くときに、タブじゃなくてスペースでインデントする
こんにちは。
先日、金さん(kimk)がお土産で買ってきてくださった栗羊羹があまりに美味しかったので、いっそ羊羹をまとめ買いするかどうか、悩んでいる久保(kubor)です。

羊羹のまとめ買いで悩む方が読者にいるかどうかはわかりませんが、プログラムを書くときにスペースインデントかタブインデントか悩むことはあるかと思います。

例えば、Python2では問題になりませんが、Python3だとスクリプト中のインデントに使うスペースとタブの混在は許されません。
したがってどちらかに統一する必要があります。
この場合、タブで書くのか、スペースで書くのか、悩むならスペースで書くのがおすすめです。

タブは環境依存で幅が変わる

多くのテキストエディタにおいて、タブはスペース8個分で表示されますが、
テキストエディタや、設定によって表示幅が変わります。

これによって大きく可読性が損なわれることもあります。
スペースでインデントするメリットは、この幅が変わらないという点です。

スペースでインデントするためにスペースキーの連打は不要

pep8ではインデントのスペース幅は4個分を推奨しています。
とは言え、毎回スペースを4回もタイプしていると気が狂いそうです。

テキストエディタ側で、タブキーをタイプした際にスペースが4個入力されるように設定していまいましょう。

例えば、Vimをお使いの方は、「~/.vimrc」に以下を記述します。
set expandtab
set shitwidth=4
set softtabstop=4
各項目の意味は以下の通りです。
  • expandtab:タブキーをタイプしたときに、タブの代わりにスペースを入力
  • shiftwidth:Vimが自動的にインデントするときのスペースの数
  • softtabstop:タブキーを入力した時のスペースの数

これで気が狂わずに済みます。
この設定を有効にしている時に、本当にタブを入力したいときは、
Ctrl+V [TAB]
で入力可能です。

自分が入力したのがタブなのかスペースなのか見えないからわからない時には、不可視文字を表示させる設定が有用です。
「Vimで不可視文字を表示させる方法」
| kubor | バイオインフォマティクス | 16:03 | comments(0) | - |
Pythonでタブ区切りテキストの読み込み
世の中には、2種類の人間がいる。
タブ区切りテキストを使う人と、使わない人だ。

こんにちは、タブ区切りテキストを使う人、久保(kubor)です。
Pythonでタブ区切りテキストを読み込むときは、csvライブラリが便利です。
標準ライブラリなので環境依存をそれほど考えなくて良いです。

以下にBEDファイルのヘッダー以降について、「end」と「行番号」を出力する例を示しました。

5行目のdelimiter='¥t'としているのがポイントです。
ここを変えればどんな区切り文字でも対応できます。

6行目のnext(reader)で、sample.tsvのヘッダーを飛ばしています。

行番号の取得はline_num()メソッドを使います。

sample.tsv
chr    start    end
chr1 18900 19356
chr2 1133 35211
chr3 21234 552312

サンプルスクリプト
import csv
tsv_file = sample.tsv
with open(tsv_file, 'r') as f:
reader = csv.reader(f, delimiter='¥t')
next(reader)
for row in reader:
print reader.line_num, row[2]
2 19356
3 35211
4 552312
ファイルハンドルは、with文を使うとスッキリするので好きです。
| kubor | バイオインフォマティクス | 14:04 | comments(0) | - |
Excelの1セルにコンマ区切り文字列を入れる
解析結果を出力する際、複数の数値をコンマで結合して羅列したいことがあります。

例えばBEDフォーマットでは、エクソンの長さや開始位置を示すblockSizesやblockStartsには「468,69,147,159」や「0,608,1434,2245」のような値を入れることになっています。
※BEDフォーマットって何?という方は、過去記事「BEDフォーマット完全解説」をご覧ください。

このようなファイルをExcelで開くと、コンマが桁区切りと認識され、1つの数値に変換されてしまうことがあります。

下の図の赤で囲ったセルには「152,159,198,136,456」という文字列が入っています。
「152」「159」「198」「136」「456」という5つの数値をコンマで結合した文字列を入れたかったのに、「152159198136456」という数値として認識されてしまいました。


対策として、Excelにファイルを読み込むときにデータ型を明示的に文字列として指定すれば良いのですが、毎回やるのは少し面倒くさいですね。

もっと簡単な対策として「152, 159, 198, 136, 456」のようにコンマの後に半角空白を入れるのがおすすめです。これで勝手に変換されなくなります。

もしくは「152,159,198,136,456,」のように、最後にもう一つコンマをつけてもいいです。UCSC GenomeBrowserのBED12のblockSizesやblockStartsには最後にコンマがついています。
| hat | バイオインフォマティクス | 14:53 | comments(0) | - |
Vimで不可視文字を表示させる方法
こんにちは、まだ都営三田線に乗ったことがない久保(kubor)です。

先日CentOS6.7がリリースされましたね。
個人的に、Vim7.4に対応したのが嬉しいです。
嬉しいので、VimのTipsを書かせてください。

「Vimで不可視文字を表示させる方法」です。

Vimに限らず、テキストエディタでコーディングしているとタブやスペースが見えず、不安になると思います。怖いですよね。

Vimでは、~/.vimrcに以下の2行を書いて、表示させることが出来ます。
set list
set listchars=tab:»-,trail:-,nbsp:%,eol:↲

set listで、不可視文字を表示させます。
set listcharsで、表示させる(置き換える)文字を設定します。
表示文字の設定はお好みでどうぞ。

listcharsについて、項目の対応は以下です。
他にも有りますが基本的にはこの4つで十分だと思います。

  • tab:タブ文字を表示
  • trail:行末のスペースを表示
  • nbsp:ノーブレークスペースを表示
  • eol:改行を表示



余分なスペースもわかるし、タブとスペースの区別も出来て便利です。

ちなみに、一時的にこれらを非表示にしたいときは、
コマンドで対応するといいと思います。
:set nolist
| kubor | バイオインフォマティクス | 14:43 | comments(0) | - |
フォーマットもいろいろ
バイオインフォマティクスでは多くのファイル形式が使われますが、GFFとGTFは名前も似ていてややこしいですね。ということで、今回はGFFとGTFの違いに触れながら、フォーマットの説明をしたいと思います。

GFF(General Feature Format)はゲノムの配列に関連した特徴を示した9列からなるタブ区切りのフォーマットです。
chr1    hg19_rmsk       transcript      10001   10468   .       +       .       ID=(CCCTAA)n;geneID=(CCCTAA)n
chr1 hg19_rmsk exon 10001 10468 1504.00 + . Parent=(CCCTAA)n
chr1 hg19_rmsk transcript 10469 11447 . - . ID=TAR1;geneID=TAR1

それぞれの列の説明です。
1. seqname:染色体またはスキャフォールドの名前
2. source:データベースやプロジェクト名
3. feature:配列のタイプ(例:repeat, exon, promoter, etc)
4. start:配列の開始位置
5. end:配列の終了位置
6. score:任意のスコアが入る。ない場合は '.'
7. strand:+(forward)、-(reverse)または '.'
8. frame:0から2の数値で、翻訳を開始する塩基の位置を示す(0なら配列の1番目がコドンの1番目の塩基になる)。コーディング領域でない場合は '.'
9. attribute:グループ

GTF(Gene Transfer Format)はGFFと8列目までは同じ項目ですが、GTFでは9列目にgene id とtranscript idを含み、セミコロンとスペース区切りで連続した複数の情報が付加されています。
chr1    hg19_rmsk       exon    16777161        16777470        2147.000000     +       .       gene_id "AluSp"; transcript_id "AluSp";
chr1 hg19_rmsk exon 25165801 25166089 2626.000000 - . gene_id "AluY"; transcript_id "AluY";
chr1 hg19_rmsk exon 33553607 33554646 626.000000 + . gene_id "L2b"; transcript_id "L2b";


ここまで簡単に2つのフォーマットの説明をしましたが、もっと詳しく知りたい方は、下記のリンクをご参照ください。
UCSC FAQ
Sanger Institute
The Sequence Ontology Project
Brent Lab Homepage
| onouek | バイオインフォマティクス | 14:33 | comments(0) | - |
Python デバッグ
Pythonに限った話ではありませんが、一生懸命プログラムを書いて、いざ実行すると、いくつものバグに遭遇します。エラーが出る、あるいは、エラーは出ないのに、プログラムが意図したとおりに動いてくれないということも。

こういう場合、きっとどこかで条件文を書き間違えたとか、整数同士で割り算をしたら結果の小数点が切り捨てられて整数になってしまったとか、そんな些細なミスを犯しているに違いありません。
しかし、どこで何がおかしくなっているのか、何行にもわたる長いプログラムを眺めていても、たいていの場合は解決しません。
そんなとき助けてくれるのが、デバッガです。

以前「Perl デバッグ」という記事にて、Perlのデバッガについてご紹介しています。
今回の記事では、同じように動くPythonの標準のデバッガをご紹介します。
大体同じような内容ですので、Perlの記事と併せて読んでいただけましたら幸いです。

Pythonを実行するときに -m pdb と付けることで利用できます。
(pdbって、python debugの略でしょうか?)

$ python -m pdb sample.py
> /home/work/directory/sample.py(9)()
-> """
(Pdb) コマンドを入力します

コマンドはPerlとほとんど同じです。

list(l):実行が止まっている行のソースが前後あわせて10行表示されます。

next(n):一行実行します。関数やモジュールは飛ばします。

step(s):nextと違い、関数やモジュールの中まで呼び出して一行ずつ実行します。

print(p):「print 変数」と実行することで、実行時点での変数の中身を表示できます。変数に意図した値が入っていない、といった問題を確認できます。

continue(c):ブレークポイントまでプログラム実行を進めます。ブレークポイントが設定されていないとき(あるいは、ブレークポイントが条件文の中に設定されていて、その条件が一切満たされないとき)は、プログラムの最後まで実行します。

break(b):ブレークポイント(continueコマンドの実行を止める場所。長いコードの、目的の場所の前や、ループの中などに設定すると便利です)を設定します。ブレークポイントは複数設定することができます。

clear(cl):設定したブレークポイントを削除します。

quit(q):途中で終了します。あるいは「「ctrl+D」でも終了できます。

他にもいろんなコマンドがあるので、気になった方は調べてみてください。
Pythonを書いてみたはいいけど、なぜかうまく動かない……とお悩みの方、これでバグをつぶしていきましょう!
| kubo | バイオインフォマティクス | 15:25 | comments(0) | - |
無害なメッセージと有害なメッセージ
hatです。私は昔から自分の体調を気にしすぎてしまうところがあります。

頭が痛いといっては脳腫瘍を疑い、お腹が痛いといっては盲腸を疑い、指のささくれが膿んだといっては指が壊死するのではないかと悶々とし、しかし結果的には何でもなくてここまで元気に生きてきました。

10年ほど前でしょうか、足の親指にある黒い塊を皮膚がんではと気に病み、死んだ後見られて困るものなどを身辺整理した後、勇気を出して皮膚科を診察した結果、「大きなホクロです」と診断されたことがあります。
せっかくなので「でも変な形をしているんです」などと少し食い下がってみましたが、「ただの変な形のホクロです」ときっぱりと言われました。

そう言われると確かにそれはホクロにしか見えず、安堵するやら恥ずかしいやらの複雑な心境で帰ってきました。そして、一目で良性か悪性か判断できるとは、お医者さんはさすがプロだなあと思いました。

ソフトウェアを実行していて警告やエラーが出た時に、気にするべきものか無視していいものか判断が迷うことがあります。

GATK(v1.6)を実行していて、次のような警告が出ることがあります。
WARN 09:58:00,492 RestStorageService - Error Response: PUT '/GATK_Run_Reports/xxx.report.xml.gz' -- ResponseCode: 403, ResponseStatus: Forbidden,

GATKは実行時にログをBroadInstituteまたはAmazonクラウドへ送るようにしているのですが、なんらかの理由でログ送信が失敗した場合にこの警告が出るようです。
他にエラーが無く、結果ファイルが出力されていれば、実行自体はうまくいっているので気にしなくて大丈夫です。

ちなみにログを送らないようにするには、https://www.broadinstitute.org/gatk/request-key からGATKのキーファイルを申請し、GATK実行コマンドに「-et NO_ET -k キーファイル」オプションを付けて実行する必要があるそうです。

参考
http://gatkforums.broadinstitute.org/discussion/1250/what-is-phone-home-and-how-does-it-affect-me
http://gatkforums.broadinstitute.org/discussion/1721/warnings-and-errors-during-score-calibration

| hat | バイオインフォマティクス | 16:12 | comments(0) | - |
BED、VCFをスッキリと染色体番号順にソートする方法
こんにちは、あんドーナツ好きで有名な久保(kubor)です。
ブログのネタを日頃から探しているのですが、なかには、過去に先輩が取り上げていそうで、意外にも紹介していない話題が結構あります。

例えば、sortコマンド。
このコマンドは、言わずと知れた基本的なUnixコマンドですが、
社内でもあまり知られていない、素敵なオプションがあります。

このオプションがあればBEDファイルやVCFファイルを綺麗にソートすることが出来ますよ。

では、実際にBEDファイルをソートしてみます。
sort -k 1,1 target.bed
chr1 172113674 172113784 MIR199A2
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr2 37869024 37899342 CDC42EP3
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 161581163 161583014 AGPAT4-IT1
chr6 39282473 39290330 KCNK16
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
-kオプションで一列目を指定していますが並び順が変ですね。
ポジションもソートされていませんね。
まずはポジションをソートしましょう。
sort -k 1,1 -k 2,2n target.bed
chr1 172113674 172113784 MIR199A2
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
-k 2,2n オプションで2列目を数値扱いでソートしました。
ポジションはバッチリですね。
あとは、染色体です。「-V」オプションを追加します。
sort -V -k 1,1 -k 2,2n target.bed
chr1 172113674 172113784 MIR199A2
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
綺麗にソートされましたね。
文句なしです。

本来「-V」オプションは、ソフトウェアのバージョン番号をソートするためのオプションですが、染色体番号にも使用可能です。

※ 残念ながらOS Xの場合は、付属のBSD版のsortコマンドに、「-V」オプションがありませんので、GNU版のsortが入っているCoreutils(http://www.gnu.org/software/coreutils/coreutils.html)をインストールする必要があります。

ちなみに、VCFファイルも1列目が染色体番号、2列目がポジションと、BEDファイルの例と同じですので、まったく同じオプションでソート可能です。
ですが、オプションを毎回書いてると面倒なので、エイリアスしておくといいですよ。

sortBedはbedtoolsですでにあるコマンド(http://bedtools.readthedocs.org/en/latest/content/tools/sort.html)ですし、sortFineで、登録しちゃいましょう。

bashをお使いの方は、~/.bashrc
zshをお使いの方は、~/.zshrc
などへ以下を記述します。
alias sortFine="sort -V -k 1,1 -k 2,2n"
必要であれば設定ファイルを再読み込みします。
source ~/.bashrc
or
source ~/.zshrc

それでは、使ってみましょう。
sortFine target.bed
chr1 172113674 172113784 MIR199A2
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
うーん、これは相当素敵です。相当ファイン、sortFineですね。
| kubor | バイオインフォマティクス | 14:24 | comments(0) | - |
prefetch すらっと落とす SRA
こんにちは、久保(kubor)です。
先日転んだ時に頭を打ったので、それがきっかけで、僕に異常な言動がないかhatさんに確認してもらったのですが、
「打つ前から変わらずに変なことを言ってるよ。」
とのことでした。ひどい。

そんなhatさんが「すら(SRA)っとクイック(Q)に変換」を公開していたので、僕はこれに対抗して、
「NCBIのSRAから複数ファイルをまとめてダウンロードする方法」を紹介いたします。
組み合わせると”すらすら”とfastqを手に入れることが叶います。
今回は柑橘のシーケンスデータ「DRR008634 - DRR008641」8ファイルをダウンロードします。

まずは、NCBIのSRAからファイルリストを手に入れましょう。
どれでもいいのでアクセッション番号で検索して、「Run Selector」を確認します。


目的の8ファイルに辿りつけましたので、これらのアクセッション番号を手に入れます。
「Accession List」からテキストファイルをダウンロードします。



このファイルには、改行区切りでアクセッション番号が記述されています。
DRR008634
DRR008635
DRR008636
DRR008637
DRR008638
DRR008639
DRR008640
DRR008641

これを使用してSRA Toolkitでまとめてダウンロードしましょう。
使うコマンドは、「prefetch」です。SRA Toolkitに含まれていますので、ご確認ください。
prefetch --option-file SRR_Acc_List.txt

「--option-file」の引数で指定したファイルの中身のアクセッション番号を、順にSRAデータベースからダウンロードしてくれます。

ダウンロードされたファイルは、以下に保存されます。
$HOME/ncbi/public/sra

fastq形式への変換は「すら(SRA)っとクイック(Q)に変換」を御覧ください。
| kubor | バイオインフォマティクス | 17:17 | comments(0) | - |
ChemmineRを使ってみよう【4】
ChemmineRの紹介連載4回目です。
前回は、ChemmineOBというパッケージを使って化合物のPubChem fingerprintを取得しました。
今回では、そのfingerprintを使って、類似比較・クラスタリング解析を行います。

fingerprintによる類似検索

fpset[1]と類似したものをfpset全体から検索します。
fpSim(fpset[1], fpset, method="Tanimoto", cutoff=0.4, top=4)
650001
1
化合物のID類似度が出力されます。今回解析した10サンプル中に、cutoff以内で類似している化合物がなかったようです。

fingerprintを用いたクラスタリング解析
simMA <- sapply(cid(fpset), function(x) fpSim(fpset[x], fpset, sorted=FALSE))
hc <- hclust(as.dist(1-simMA), method="single")
plot(hc,hang=-1)


ヒートマップの作成
クラスタリングまでしたら、ヒートマップだって描きたいですよね。
library(gplots)
heatmap.2((1-simMA), Rowv=as.dendrogram(hc), Colv=as.dendrogram(hc), col=bluered(256), density.info="none", trace="none")



これで化合物間の類似関係が視覚的にわかりますね!
(10サンプルなので少々つまらない結果になりましたが……)


非常に簡単な紹介となりましたが、以上でChemmimeRに関する連載は終了です。
私にとって、あまり触れたことのないケモインフォマティクス解析の良い勉強となりましたが、皆様の研究の一助にもなれば幸いです。
| kubo | バイオインフォマティクス | 16:11 | comments(0) | - |
ChemmineRを使ってみよう【3】
前回までで、SDFデータの読み込みと、データのvalidationをご紹介しました。
今回は、読み込んだ化合物のfingerprint/atom pair descriptorを取得する方法です。

SDFsetからfingerprint/atom pair descriptorへの変換には、ChemmineOBという新しいパッケージをインストールします。
Open Babelというツールを扱うためのRインターフェースです。
biocLite("ChemmineOB")
library(ChemmineOB)
※ Open Babelのデータベースの関係で、OSによって挙動が違うようです。詳細はChemminOBのREADMEをご覧ください。

化合物のFingerprintの取得
fingerprintを取得します。
fpset<-fingerprintOB(sdfset,"FP2")
fpset
An instance of a 1024 bit "FPset" of type "FP2" with 10 molecules
fingerprint以外にもatom pair descriptorでも解析できます。
apset <- sdf2ap(sdfset)
apset

An instance of "APset" with 10 molecules

次回では、fingerprintを使った類似検索・比較を行います。
| kubo | バイオインフォマティクス | 16:25 | comments(0) | - |
ChemmineRを使ってみよう【2】
ChemmineRの紹介記事2回目、前回はパッケージのインストールとSDFファイルの取得(読み込み)の方法をご紹介しました。
今回は、取得した読み込んだデータのアクセス方法と、SDFデータの確認方法をご紹介します。

前回は、CIDからgetIds関数で、またはSDFファイルの読み込みによってSDF情報を読み込んで、sdfsetというSDFsetオブジェクトを作成しました。
sdfset
An instance of "SDFset" with 10 molecules
class(sdfset)
[1] "SDFset"
attr(,"package")
[1] "ChemmineR"
SDFsetオブジェクトの各要素には、以下のようにしてアクセスします。
sdfset[1] # 1番目の要素にアクセス
sdfset[5:10] # 5〜10番目の要素にアクセス
SDFsetオブジェクトのヘッダを見るときはheader関数を使います。1つめの化合物のヘッダは次のようにしてみることができます。
header(sdfset[1])
$CMP1
Molecule_Name
"650001"
Source
" -OEChem-07071010512D"
Comment
""
Counts_Line
" 61 64 0 0 0 0 0 0 0999 V2000"

続きから、読み込んだSDFデータの確認をしていきます。
続きを読む >>
| kubo | バイオインフォマティクス | 16:00 | comments(0) | - |
すら(SRA)っとクイック(Q)に変換
NCBI SRAからダウンロードしたファイルがsraフォーマットの場合、以下のコマンドでまとめてfastqに変換すると便利です。

$ find . -name '*.sra' -exec fastq-dump {} ¥;


fastq-dumpについては
NCBI SRA Toolkitの使い方
も、findコマンドについては
findの-execオプション
もご覧ください。

| hat | バイオインフォマティクス | 15:08 | comments(0) | - |
ChemmineRを使ってみよう【1】
ChemmineRというパッケージの使い方をご紹介します。
オンラインでケモインフォマティクス解析を行えるChemMine ToolsというツールのRパッケージです。
ぱっと見て「けみなー」と読みましたが、もともとのツールがChemMineなので、「けむ・まいなー」と読むのでしょう。

連載で、公式のマニュアルにある「Five Minute Tutorial」から一部と、PubChem fingerprintを使った解析の一部を簡単に説明いたします。


インストールと読み込み
Bioconductorに登録されています。
source("http://bioconductor.org/biocLite.R")
biocLite("ChemmineR")
library(ChemmineR)

SDF(Structure data format)ファイルの取得
解析を始める前に、解析したい化合物の情報が必要です。化合物のCID(PubChem ID)か、SDFファイルを使って解析します。

解析したい化合物のCIDのリストから始める場合は、getIds関数でSDFの情報を取得します。getIdsには用意したCIDのベクトルを与えます。少し時間がかかると思います。
sdfset <- getIds(c(650001,650002,650003,650004,650005,650006,650007,650008,650009,650010))
取得したSDFデータは念のためファイルに書き出しておくと、同じ化合物を繰り返し解析するときに便利だと思います。
write.SDF(sdfset, file="sub.sdf", sig=TRUE, cid=TRUE, db=NULL)
書き出したファイルの読み込みは以下の通りです。CIDではなく、SDFファイルを用意した人も、同様にSDFファイルを読み込むことができます。
sdfset <- read.SDFset("your_file.sdf")
データがない場合はパッケージのサンプルデータを使います(今回例で用いているものと同じです)。
data(sdfsample)
sdfset<-sdfsample[1:10]
次は読み込んだSDFファイルの重複や有効性の確認方法をご紹介します。
| kubo | バイオインフォマティクス | 14:49 | comments(0) | - |
あとはまかせない
以前、atコマンドを使うと時間差でコマンドが実行できて便利だという記事(あとはまかせた!)を書きました。

atコマンドはとても便利で普段から活用しているのですが、仕掛けたタスクをとりけしたい場合もあると思います。

そんな時は、「at -l」で仕掛けたタスクを確認して、「at -d タスク番号」でタスクを除去できます。

$ at -l
4 2015-04-08 20:01 a hat

$ at -d 4

まかせるつもりだったけどやっぱりやめた!」という時にぜひお試しください。
| hat | バイオインフォマティクス | 15:09 | comments(0) | - |
お花見メタゲノムに参加します!
東京に開花予想が出た頃、お花見メタゲノムプロジェクトの採取キットが届きました!


個人で参加するので解析まではお手伝いできないかもしれませんが、これを機にメタゲノム解析について少しでも勉強できればいいなと思っています。

サンプルの採取場所は、自宅近くの目黒川沿いにする予定です。
あの綺麗な満開の花に、どんな細菌がうようよいるのでしょう。ゾクゾクします。

最後に、私の好きな都々逸をご紹介します。

お酒呑む人 花なら桜
今日もさけさけ 明日もさけ


いよっ!粋ですね。それでは皆様も良いお花見を。
| hat | バイオインフォマティクス | 14:21 | comments(0) | - |
ReactomePAを使ってみた
パスウェイデータベースReactomeをRから操作できるBioconductorパッケージ、ReactomePAを使ってみました。

インストール
> source("http://bioconductor.org/biocLite.R")
> biocLite("ReactomePA")

テスト実行
> library(ReactomePA)
> library(DOSE)
> data(geneList)
> de <- names(geneList)[abs(geneList)>1]
> require(ReactomePA)
> x <- enrichPathway(gene=de, pvalueCutoff=0.05, readable=T)
> barplot(x, showCategory=8)
遺伝子のパスウェイ解析をする時に便利そうです。
| hat | バイオインフォマティクス | 17:23 | comments(0) | - |
findの-execオプション
findコマンドの-execオプションを使うと、findの結果を他のコマンドで実行することができます。

例:現在のディレクトリにある *.gz ファイルを全て解凍する
$ find . -name '*.gz' -exec gunzip {} ¥;

・findの結果が{}に入ります。
・末尾の「¥;」の前には半角空白1つ入れてください。

これは、先日お客様にトレーニングを行っている時に、お客様から教えていただきました。
今までこのような場合はfindの結果をパイプしてxargsに食わせていたのですが、このほうが簡単に書けそうです。

よく使うfindコマンドにこんなオプションがあったとは。
勉強になりました!
| hat | バイオインフォマティクス | 13:43 | comments(0) | - |
アメリクのご紹介
アメリクとは……
アメリク

バイオ研究の解析に使用するソフトや解析手法について、無償で調査するサービスです。
調査結果はアメリエフのブログでご紹介いたします→ソフトの紹介例
「学会で●●というソフトについて聞いたけど、どんなソフトだろう? メリットやデメリットは何だろう?」のような疑問の解決などに、ぜひご利用ください。

申込みフォームはこちら
◇ ブログの公開時期などは、ご希望にそえない場合がございますのでご了承ください。
◇ ご依頼者の許可なく個人情報を第三者に開示することはございません。
| kubo | バイオインフォマティクス | 14:35 | comments(0) | - |
SAM/BAMのフラグの意味を教えてくれるサイト
SAM/BAMのFLAG列には、そのリードまたはペアとなるリードの
マッピング状態に関する情報がビットの和で示されています。

例えばあるリードがFLAG=12だった場合、12は2進数で1100となり、
SAM Specによると「マップされなかったリード」は0x4で
「相方がマップされなかったリード」は0x8ですので、
「自分も相方もマップされなかったリード」であることがわかります。

仕事中、FLAGの意味を確認したいことがしばしばあるのですが、
この「FLAGを2進数変換して解釈する」ステップが地味に面倒くさく
感じていました。

この件に関して、いつも頼りにさせていただいている
NGS Surfer's Wikiに、非常に便利なサイトの紹介があるのを見つけました。
続きを読む >>
| hat | バイオインフォマティクス | 15:06 | comments(0) | - |
FastQCの新機能
FastQCにv0.11.1から次のチェック項目が追加されています。

・Per tile sequence quality
タイル単位のクオリティが出せるようになりました。


・Adapter Content
これまではOverrepresented sequencesでアダプタも確認していましたが、アダプタ混入率を示す専用の項目ができました。


Duplication moduleやK-mer moduleも改良されたそうです。
また、まだ試していませんが、OK/Warn/NGの閾値がカスタマイズできるようになったそうです。

後はレポートHTMLが、-oで指定したディレクトリ直下に単体で出力されるようになったのが、個人的には嬉しいです。
前はディレクトリの中に他のファイルと一緒にHTMLが入っていたので、サンプル数が多いとディレクトリをポチポチ開いてレポートを見て回るのが地味に面倒くさかったんですよね。

オープンソースのソフトなのに高機能でどんどんブラッシュアップもされていて、ユーザとしては本当にありがたいです。
| hat | バイオインフォマティクス | 14:18 | comments(0) | - |
OR条件でgrep
Rでgrep

続・Rでgrep

grepのこんなオプション

grepを数えるオプション

普段解析でよく使うコマンドのためか、このブログはgrepネタが多い気がします。
ということでまたgrepネタです。

grepで特定の文字列を含む行を抜き出す場合、文字列を¥|でつなげることで「OR」にすることができます。

例:「ABC」または「DEF」を含む行を抜き出す
$ grep 'ABC¥|DEF' ファイル

例:SAMファイルから「ヘッダー行」または「XT:A:U」を含む行を抜き出す
$ grep '^@¥|XT:A:U' SAM | samtools view -Sb - > BAM

複数のgrepコマンドをパイプでつないで「AND」にするのはよくやっていたのですが、「OR」ができることは今日ふと気になって調べて初めて知りました!
| hat | バイオインフォマティクス | 14:19 | comments(0) | - |
snpEffバージョンアップ
アノテーションソフトsnpEffが1.4.0にバージョンアップしました。
GRCh38/hg38に対応しています。待ってました。

まだ安定ではないようで、頻繁にマイナーバージョンアップを繰り返しています(2014-09-25現在、最新バージョンは1.4.0 Eです)が、早く新しいゲノムバージョンのアノテーションを行ってみたいものです。

その他、西アフリカで大流行を引き起こしているエボラウイルスに対応しているのが時事を反映していますね。
早く収まるといいのですが……。
| kubo | バイオインフォマティクス | 15:58 | comments(0) | - |
初めてのClinVar
だんだん肌寒くなってきましたね。kitanoです。
私が入社して3週間になろうとしていますが、今回はその中で私が初めて知った言葉をひとつピックアップして解説していこうかと思います。

その言葉とは「ClinVar」です。
ヒトの研究に携わっていらっしゃる方々には馴染みのある言葉だと存じます。しかし、植物を扱っていた私にとっては初めての言葉でした。

「ClinVar」(http://www.ncbi.nlm.nih.gov/clinvar/)とはNCBIでサービス運営されているデータベースです。
データベースの内容としては、NCBI、National Library of Medicine(NLM)、National Institutes of Health(NIH)が主体となって作った、「医学的に重要なヒトの変異と表現型に関連した報告を提供する無料で入手可能なアーカイブ」です。
 データソースは、ヒトの変異位置情報としてdbSNP、dbVarと密接に連携し、表現型の記述としてはMedGenをベースにしています。
 ClinVarに投稿されたそれぞれの報告には投稿者・変異・表現型、またSCV000000000.0のフォーマットで割り当てられたアクセッション番号が記載され記述されています。しかし、臨床的な解釈が相反する場合には、個別にRCV000000000.0のフォーマットでアクセッション番号が割り当てられ、各変異の医学的な重要性を評価するために、同じ変異と表現型の組み合わせの報告、さらに他のNCBIのデータベースを集約し評価しています。
 また、ClinVar内のデータはXML、VCF、タブ切り替えをセットとしてHTMLやダウンロード含む複数のフォーマットで利用できます。

 このようにClinVarはヒトの変異と表現型に焦点を当てたデータベースであり、臨床医学の発展のためにも重要なデータベースの一つであると考えられます。


【参照】
CliVar:http://www.ncbi.nlm.nih.gov/clinvar/
NCBI:http://www.ncbi.nlm.nih.gov/
National Library of Medicine(NLM):http://www.nlm.nih.gov/
National Institutes of Health(NIH):http://www.nih.gov/
dbSNP:http://www.ncbi.nlm.nih.gov/snp/
dbVar:http://www.ncbi.nlm.nih.gov/dbvar/
MedGen:http://www.ncbi.nlm.nih.gov/medgen
| kitanog | バイオインフォマティクス | 15:28 | comments(0) | - |
STAR(1)
RNA-seqデータを高速にマッピングするソフトウェアSTARについて
ご紹介します。

A. Dobin et al, Bioinformatics 2012;
doi: 10.1093/bioinformatics/bts635
"STAR: ultrafast universal RNA-seq aligner"

◆STARのインストール
$ wget http://STAR.googlecode.com/files/STAR_2.3.0e.Linux_x86_64.tgz
$ tar zxvf STAR_2.3.0e.Linux_x86_64.tgz


◆ゲノムのダウンロード
$ wget ftp://ftp2.cshl.edu/gingeraslab/tracks/STARrelease/STARgenomes/hg19_Gencode19.tgz
$ tar zxvf hg19_Gencode19.tgz


◆STAR実行
$ STAR --genomeDir hg19_Gencode19 --readFilesIn FASTQ1 FASTQ2 --runThreadsN 3


ヒトRNA-seqデータ(※1)をTophat(※2)でマッピングした結果と
比較してみました。
続きを読む >>
| hat | バイオインフォマティクス | 10:24 | comments(0) | - |
BEDtoolsのオプション
BEDtoolsは何かと便利なツールです。当ブログでも、intersectBedの記事のページがアクセスランキング上位になっています。。
つい最近までBEDtoolsのバージョンはずっと前にインストールした2.13を使っていましたが、先日より新しいバージョンのものを使ったところ、intersectBed、windowBedなどのツールに、2.13にはなかったオプションを発見しました。

"-header" です。

旧バージョンのBEDtoolsでは、BED、VCFなどのフォーマットのヘッダ行が無視され、データ行のみが出力されていました。特にVCFはヘッダがないと不便なので、いつもヘッダのみを最初に出力し、次にデータ行を追記していましたが、面倒だなと思っていました。
具体的にどのバージョンから実装されているのかは現在わかりませんが、この -header オプションで、データ行だけでなくヘッダも出力できるようになったので、その手間が省けるようになりました。

ざっとよく使うツールだけ確認しましたが、intersectBed、mergeBed、flankBed、windowBed、sortBedあたりは-headerオプションが使えます。

このオプション追加はツールの性能自体の変化ではありませんが、こうして便利な機能が増えると嬉しいです。
| kubo | バイオインフォマティクス | 14:37 | comments(0) | - |
ヒートマップ図の色指定
Rでヒートマップ図を描く時、gplotsパッケージに用意されている色のセットを使うと、色指定が楽です。

> install.packages('gplots')
> library(gplots)


論文などを見ていると、最近は青→白→赤が多いような気がするのは私だけでしょうか。
青→白→赤は、「bluered(色数)」で指定できます。色数を多くするとグラデーションを細かくできます。
> heatmap(as.matrix(mtcars), col=bluered(256))


従来の緑→黒→赤は「greenred(色数)」で指定できます。
> heatmap(as.matrix(mtcars), col=greenred(256))


「もう少し深みのある緑がいい」といったこだわりのある方は、以下のようにして自分で色を作るのもいいかと思います。
> g2r <- c(rgb(0, seq(0.5,0,-0.01), 0), rgb(seq(0,1,0.01), 0, 0))
> heatmap(as.matrix(mtcars), col=g2r)

| hat | バイオインフォマティクス | 14:38 | comments(0) | - |
grepを数えるオプション
特定の文字列を含む行を抜き出す便利なコマンドgrep。
最近便利に使っているオプションをご紹介します。

それは grep -c
grepコマンドで取得できた行数を出力してくれるオプションです。cはcountのcと覚えています。
これを知らなかったときは、grepコマンドの結果をwcコマンドで数えていましたが、-cオプションの方が便利です。

例として下記のようなVCFファイル sample.vcfがあったとします。


ヘッダ行を除いたデータ行の出力コマンドは以下になります。
$ grep -v "^#" sample.vcf

※ -v は指定文字列を除いた行を出力するオプション、^は行頭を示す記号です。

データ行を数えるときは以下のコマンドです。
$ grep -c -v "^#" sample.vcf
3

ちなみにこのオプションを知らなかったときは、こう数えていました。もちろん結果は同じですよ。
$ grep -v "^#" sample.vcf | wc -l
3
| kubo | バイオインフォマティクス | 15:52 | comments(0) | - |
hg38調査(1)
hg20/GRCh38が出る と浮かれていた、去年の夏。

やっと出たhg38(hg20という名前では無かった)をダウンロードして
染色体の長さやNの数を調べた 今年の正月。

それで満足してその後hg38のことはすっかり忘れていましたが、
気が付くと8月も後半。小学生は夏休みの宿題が気になりだす時期です。
私も放置していたhg38のことが気になってきました。

ということで、RNA-seqとChIP-seqの公開データを同じコマンドで
hg19とhg38にマッピングし、マッピング率の違いを調べました。
ゲノムの精度が上がっているなら、マッピング率が上がることが
期待されます。

以下、結果です。左がhg19、右がhg38のマッピング率です。

RNA-seq Paired-end(tophatでマッピング)
サンプルA:85.28% → 85.79%
サンプルB:91.28% → 91.76%

ChIP-seq Single-end(bowtieでマッピング)
サンプルC:69.83% → 70.34%
サンプルD:98.19% → 99.04%

今回選んだ4サンプルは、いずれもhg38でマッピング率がわずかに向上しました。
いいですね!

次回、Reseqデータでも検証してみたいと思います。
| hat | バイオインフォマティクス | 14:47 | comments(0) | - |
名前が2つ?
学生時代に研究していた植物のタンパク質には、2つの名前がありました。
そのタンパク質をつくる遺伝子をノックアウトすると起きる現象由来の名前と、タンパク質の機能由来の名前です。そのタンパク質の研究は、ノックアウトにより起きる現象から始まったため、前者のほうが歴史的に使われている名前なのですが、後者の方がタンパク質の性質を表す適切な名前、ということで2つの派閥があったようです。
2つも名前があるなんて、ややこしいなと思っていました。

ゲノムの世界でも似たようなことがありました。
主に使われているリファレンスゲノムには、Genome Reference Consortiumが配布しているものとUCSC Genome Browserが配布しているものがありますが、これまでは、両者のbuild番号はたとえばヒトではGRCh36/hg18やGRCh37/hg19と、大きく隔たっていました。
ところが、昨年末に公表されたゲノムのbuild番号はGRCh38/hg38。UCSC Genome Browserの番号が大きく跳んで、Genome Reference Consortiumにそろえられて、互換性のあるバージョンがわかりやすくなっています。

詳細については公式の発表をご覧ください。
| kubo | バイオインフォマティクス | 16:16 | comments(0) | - |
Pythonで計算する時の注意点
Pythonで計算結果を小数点以下まで得たい時、単純に

val = 3 / 10

としてしまうと、結果が0になります。

val = float(3) / 10

のように、どちらかをfloat型にして計算すると0.3が返ってきます。

Perlだと $val = 3 / 10; で0.3が返ってくるので、Pythonでもうっかりやってしまいがちです。

気を付けたいと思います。
| hat | バイオインフォマティクス | 13:20 | comments(0) | - |
headとtail
会社(神田)の近くにおいしい鯛焼屋さんが二軒あります。
どちらも餡がたくさん詰まっていて、熱々の焼き立てをほうばると
小麦粉と重曹の香ばしい香りが鼻に抜け、たまらぬおいしさです。
会社が神田に移転して良かったと思うひとときです。

ところで、みなさん鯛焼は頭から食べる派でしょうか?
尾から食べる派でしょうか?

ということで、本日はLinuxのheadコマンドとtailコマンドの話です。
続きを読む >>
| hat | バイオインフォマティクス | 16:54 | comments(0) | - |
VCFのアノテーション
SnpSiftを使うと、VCFにdbSNPや1000Genomesのアノテーションをつけることができます。

(1)アノテーション用のデータを以下からダウンロードして解凍します。
ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz

(2)SnpSiftを実行します。

$ java -jar SnpSift.jar annotate 00-All.vcf X.vcf > X_anno.vcf
※「X.vcf」はアノテーションをつけたいVCFファイルを指定します。

デフォルトではたくさんの項目がつきますので、例えばdbSNPのバージョンと1000Genomesのアリル頻度だけつけたい場合は、以下のように指定します。
$ java -jar SnpSift.jar annotate -info dbSNPBuildID,CAF 00-All.vcf X.vcf > X_anno.vcf

デフォルトではVCFのID列とINFO列に情報が入りますが、例えばdbSNPのIDだけをID列だけにつけたい場合は、以下のように指定します。
$ java -jar SnpSift.jar annotate -id 00-All.vcf X.vcf > X_anno.vcf

この他、VCFであれば(おそらく)なんでもアノテーションに使うことができます。
遺伝子アノテーションソフトのSnpEffは弊社内でもよく使っているのですが、一緒にこんなに便利なソフトが提供されていたとは...。
| hat | バイオインフォマティクス | 14:33 | comments(0) | - |
HGVDについて (3) 平均depth編
HGVDについて (1)
HGVDについて (2) サンプル数編

※※ 前回の記事まで、データベースの名称を誤っておりました。ブログをご覧くださった皆様、並びに関係者様各位にご迷惑をおかけしてしまい申し訳ございません。

前の記事に引き続き、Human Genetic Variation Databaseで公開している日本人ゲノムデータの、データの確からしさの判断に使えそうな項目を見てみます。
データベースに含まれている項目の詳細については、前の記事をご参照ください。

今回は、Mean_depthについて調べてみます。
まず全SNPのMean_depthの分布のグラフを書いてみます。
HGVDのdepth(all)

depthが1〜200xのSNPが多いようです。
基本統計量を求めてみると、Mean_depthが浅いSNPもあります。
もう少し、Mean_depthの浅いSNPに注目してみるため、Mean_depthが50以下のSNPでMean_depthの分布を調べてみました。

HGVDのdepth(50以下)

全体における割合は少ないですが(2.81%)、Mean_depthが10未満のSNPが13,534個ありました。
もちろん、これはそのSNPをシーケンスしたサンプル全体における平均なので、ただdepthが浅いと判断するのではなく、サンプル数や、depthの標準偏差も見ながら判断するする必要があると思います。


以上でHGVD(Human Genetic Variation Database)についての紹介を終えたいと思います。
連載の最初にも書きましたが、人種特異的な遺伝情報のデータベースは、病気の遺伝子研究に役立つと考えられます。
これからも様々なデータベースが作られ、発達していくのが待ち遠しいですね。
| kubo | バイオインフォマティクス | 15:09 | comments(0) | - |
GRCh38
ヒトの新しいリファレンスゲノムGRCh38が公開されました。
GRCh37(≒hg19)と比較してどれくらい変わったのでしょうか?

まずは、各染色体の長さや、N(不明塩基)の数がGRCh37→GRCh38で
どのように変わったか調べてみました。



染色体長の変化
・各染色体につき最大30Mbpくらい増えたり減ったりしている
・ゲノム全体では7.4Mbp減少

N数の変化
・ゲノム全体で83Mb減少
・ゲノムに占めるNの割合は7.6%から4.9%に減少

GenomeRef:Announcing GRCh38によると、GRCh37で
セントロメア領域にあった大量のNが、GRCh38では減らせている
そうです。
また、1000ゲノムプロジェクトと一致しない領域やIndelが
修正されているそうです(1q21、10q11、chr9など)。

新しいリファレンスゲノムで新しい知見がたくさん見つかることを
期待しています。
| hat | バイオインフォマティクス | 16:47 | comments(0) | - |
HGVDについて (2) サンプル数編
前の記事
次の記事

Human Genetic Variation Databaseで公開している日本人ゲノムデータの、データの確からしさの判断に使えそうな項目を見てみます。

配布されている圧縮ファイルをダウンロードし、解凍しますと、READMEとタブ区切りテキストの2つのファイルができます。

READMEによると、日本人から検出された478,228個のSNPについて、1行につき1SNPの以下の情報が記載されています。
hgvdのカラム説明
■#Sample(サンプル数):データベースに含まれる1,208人のうち、N人でそのSNPのポジションをシーケンスできたという意味です。Altのアリルが検出されたサンプルの数ではありません。

このうち、#Sample、Mean_depthが、確からしさを判断するのに使えそうです。

まず#Sampleの分布を見てみます。
HGVDのサンプル数ヒストグラム
(※クリックで大きくなります)

まず左の図をご覧ください。
横軸が各SNPをシーケンスできたサンプル数、縦軸がその登場回数です。
200〜500サンプル、あるいはそれ以上でシーケンスできているSNPが多いようですが、案外サンプル数が少ないSNPも登録されているようです。
試しに、サンプル数が100より少ないSNPでヒストグラムを書いたのが右図です。40サンプル以下、その中でも5サンプルより少ないサンプルでしか読まれていないSNPが多数あるようです。

ちなみに、データベースのうち、100サンプル未満でしかシーケンスされていないSNPをすべて除くと、約95.0%のSNPが残ります。

次はMean_depthの分布を見てみます。



データベースの名前を「Human Genome Variation Database」と記載しておりましたが、これは別のデータベースの名前でした。
正しくは「Human Genetic Variation Database」でした。
ブログをご覧くださった皆様、並びに関係者様各位にご迷惑をおかけしてしまい申し訳ございません。
| kubo | バイオインフォマティクス | 15:58 | comments(2) | - |
bwa memの-Mオプション
bwaバージョン0.7が去年の2月にリリースされて、
ほぼ1年が経とうとしています。
もう0.7に切り替えた方も多いのではないでしょうか。

bwa 0.7では、memというコマンドが使えるようになりました。

memは従来のbwaswの後継で、ロングリード(〜1Mbp)に使える
コマンドですが、bwaswよりも速度・精度とも向上しているとのことです。

また、ショートリード(70-100bp)に対する精度もよいそうです。
ショートリードはこれまで「fastqごとにaln」→「samse/sampe」で
実行するのが主流で、中間ファイルができたりして面倒くさかったので、
mem一発でやれると便利ですね。

このように素晴らしいmemですが、1点注意があります。
memには-Mというオプションがありますが、bwaの結果を
PicardやGATKで処理する場合は、この-Mオプションを必ずつけて
実行してください。

$ bwa mem -M GENOME.fa 1.FASTQ 2.FASTQ > BAM

何故かと言うと、memはローカルアライメントを行い、
キメラリードを検出することができるようになりました。
キメラリードは、BAMのflagに0x800(supplementary alignment)が
立ちます。
ですが、PicardのMarkDuplicateはこの仕様に対応していないため、
このBAMを入力するとエラーになってしまうのです。

-Mをつけて実行したBAMでは、flagに0x800の代わりに従来の0x100
(secondary alignment)が立つため、Picardで動かすことができます。

試しに、同じデータを-Mありと-Mなしで実行した結果のBAMファイルを
比較したところ、本当にflagだけが1792(0x800=2048、0x100=256、2048-256=1792)違っていました。

キメラリードを検出できるのは良いので、早くPicardが0x800に
対応してくれるといいなあと思っています。
| hat | バイオインフォマティクス | 14:32 | comments(0) | - |
CircosでSelfChainを描く(3)
CircosでSelfChainを描く(1)
CircosでSelfChainを描く(2)
のつづきです。

手順
(1)データファイル用意
(2)設定ファイル作成
(3)Circos実行


データファイルも設定ファイルも用意できましたので
いよいよ「(3)Circos実行」を行います。

次のコマンドを実行すると、circos.pngというファイルができます。
$ circos -conf test.conf



Circosはこれ以外にもいろいろな形状のトラックを描くことが
できますので、またご紹介したいと思います。
| hat | バイオインフォマティクス | 14:36 | comments(0) | - |
CircosでSelfChainを描く(2)
CircosでSelfChainを描く(1)のつづきです。

手順
(1)データファイル用意
(2)設定ファイル作成
(3)Circos実行


今回は「(2)設定ファイル作成」を行います。
以下のような内容のtest.confという名前のファイルを作成します。


<<include colors_fonts_patterns.conf>>

<image>
<<include etc/image.conf>>
</image>

karyotype = karyotype.human.txt

<ideogram>
<spacing>
default = 30u
<pairwise /hs/ /hs/ >
spacing = 0.5r
</pairwise>
</spacing>
thickness = 25p
stroke_thickness = 1
stroke_color = black
fill = yes
fill_color = black
radius = 0.3r
show_label = yes
label_font = default
label_radius = dims(ideogram,radius_outer) + 10p
label_size = 28p
band_stroke_thickness = 0
band_stroke_color= black
band_transparency= 4
</ideogram>

<links>
radius = 0.7r
crest = 0.5
thickness = 1
color = red
bezier_radius = 0r
bezier_radius_purity = 0.3
<link>
file = data.txt
<rules>
<rule>
condition = var(intrachr)
color = blue
radius = 0.7r
bezier_radius = 1.7r
</rule>
</rules>
</link>
</links>

<<include etc/housekeeping.conf>>


以下のような設定を行っています。

1. 外円にヒトの染色体領域を描画する。
2. その内側にSelfChainのデータファイルを読み込んで描画する。
_a. 同一染色体上の関係は青線で外側に描画する。
_b. 異なる染色体間の関係は赤線で内側に描画する。


詳しい説明はCircosのチュートリアルをご覧ください。

次回は描画を行います。
| hat | バイオインフォマティクス | 14:44 | comments(0) | - |
CircosでSelfChainを描く(1)
Circosを使ってみようでご紹介したCircosを使って、
以下のようなヒトのChainSelfの図を描いてみたいと思います。



手順
(1)データファイル用意
(2)設定ファイル作成
(3)Circos実行

今回は「(1)データファイル用意」を行います。

ChainSelfは、ヒトゲノムを自分自身にアライメントして、
類似領域を探したデータです。

データはUCSCから入手可能です。
UCSC hgdownloadから、chainSelf.txt.gzをダウンロードして
解凍します。

ChainSelfのように、ゲノム上の2領域間をつなぐデータは
CircosではLinkトラックという形式で表します。

Linkトラックとして読み込むために、ダウンロードしたファイルから
染色体1、スタート位置1、エンド位置1、
染色体2、スタート位置2、エンド位置2
だけ抜き出したファイルを作成します。

また、Circosの表記に合わせ、染色体名中のchrをhsに置換します。
さらに、全データだと数が多いため、スコアが10Mより大きい箇所にだけ
絞り込みました。

以下のコマンドを実行すると、上記の変換を行うことができます。
$ awk '$2>10000000{print $3, $5, $6, $7, $10, $11}' chainSelf.txt | sed -e "s/chr/hs/g" > data.txt

これでデータファイルの作成は完了です。

次回は、Circosの設定ファイルを作成したいと思います。
| hat | バイオインフォマティクス | 15:09 | comments(0) | - |
HGVDについて (1)
先月の11月12日、日本人のゲノム情報データベースHuman Genetic Variation Databaseが公開されました。このことは一般紙でも報道されました。
以前の連載でも、ゲノム情報のデータベース、特に人種や民族別のデータベースが、高スループットの遺伝子解析機器を用いた解析に有用だと少し触れました。ぜひ、今後このデータベースが日本人の遺伝子研究に活用され、様々な病気の治療に役立って欲しいと思います。

このデータベース上のSNPの半分以上は、dbSNPなどこれまでの既知変異データベースに含まれていない日本人特異的なものだということで、日本人の遺伝性疾患の研究に、日本人特有のデータベースがいかに重要かということがわかります。
一方で、アレル頻度が非常に低いSNPも大量に含まれているとも言われています。


次回では、せっかくなので、実際に中身を覗いてみます。
ホームページによると、このデータベースは日本人1208名のエクソームシーケンシングデータから発見されたアリルやジェノタイプの頻度が載っています。
他のゲノム情報データベースと比べてみますと、HapMapではサンプル270人、1000人ゲノムプロジェクトではサンプル2577人ということですので、日本人集団だけで1208人というのは、大規模に感じられます。
そのほか、データベースには信頼性の指標となりそうな、SNPをカバーできたサンプル人数や平均カバレージなどの情報が含まれていると説明されています。
これらの数字を見てみたいと思います。


HGVDについて (2) サンプル数編
HGVDについて (3) 平均depth編


データベースの名前を「Human Genome Variation Database」と記載しておりましたが、これは別のデータベースの名前でした。
正しくは「Human Genetic Variation Database」でした。
ブログをご覧くださった皆様、並びに関係者様各位にご迷惑をおかけしてしまい申し訳ございません。
| kubo | バイオインフォマティクス | 14:22 | comments(0) | - |
DeNovoGear
tokunagaです。

有言実行です。
第2回使ってみたシリーズです。

今回ご紹介するのは

DeNovoGear

サイト
http://denovogear.weebly.com/index.html
論文
http://www.nature.com/nmeth/journal/v10/n10/full/nmeth.2611.html


trioデータからDeNovoのmutationやIndelを検出するツールです。
likelihood-based errorというモデルを使用し、
false positiveを減らしているようです。

Inputファイルはpedファイルとbcfファイルです。

pedファイルについては過去のブログ記事をご参照ください。

bcfファイルは簡単に言うとVCFファイルのバイナリ版で、
検出したvariantの情報が記載されています。
詳しい説明はこちらに。
http://samtools.sourceforge.net/mpileup.shtml

このDeNovoGearでは

1.trioデータから新規変異を検出する(X染色体については別解析)
2.pairデータから新規変異を検出する
3.Phasing denovo mutations

といったパターンの解析が出来るようです。

テストデータもありましたので、1のtrioデータで実行してみました。
実行コマンドは以下です。
$ ./denovogear dnm auto --bcf sample_CEU.bcf --ped ¥
 sample_CEU.ped

上記を実行しますとCHILD_IDのサンプルに特異的な変異を検出してくれます。
outputの指定なしだと以下のように結果が出力されました。


Denovogear0.5.3
Created SNP lookup table
First mrate: 1 last: 1
First code: 6 last: 6
First target string: AA/AA/AA last: TT/TT/TT
First tref: 0.0002388 last: 0.99301

Created indel lookup table First code: 6 last: 6
First target string: RR/RR/RR last: DD/DD/DD
First prior: 0.05 last: 0.114

Created paired lookup table
First target string: AA/AA last: TT/TT
First prior 1 last: 1
DENOVO-SNP CHILD_ID: NA12878_vald-sorted.bam.bam chr: 2 pos: 214668360 ref: G alt: A maxlike_null: 3.95324e-12 pp_null: 0.000399465 tgt_null(child/mom/dad): GG/GG/GG snpcode: 1 code: 6 maxlike_dnm: 9.9301e-09 pp_dnm: 0.999601 tgt_dnm(child/mom/dad): AG/GG/GG lookup: 4 flag: 0 READ_DEPTH child: 48 dad: 76 mom: 34 MAPPING_QUALITY child: 59 dad: 59 mom: 59

Total number of SNP sites interrogated: 36
Total number of SNP sites passing read-depth filters: 36
Total number of INDEL sites interrogated: 1
Total number of INDEL sites passing read-depth filters: 1
Total number of Paired sample sites interrogated: 0
Total number of Paired sample sites passing read-depth filters: 0
Done !


--out_vcfオプションで
vcfファイルを出力することも可能です。

ご興味のある方は試されてみてはいかがでしょうか?
| tokunaga | バイオインフォマティクス | 17:46 | comments(0) | - |
MEDIPS
ご無沙汰しております。
tokunagaです。

ここ最近、
このような工夫を施したツールやアルゴリズムを開発しました等、
バイオインフォマティクスに関する論文が数多く出てきてますね。
全てを把握することはとてもじゃないですが難しいことです。

そこで折角なので、
実際に最近論文で報告された解析ツールを
試しに使用いたしまして、
この場で簡単にご紹介したいと思っております。

今後も使ってみたシリーズで記事書いていきたいと
密かに思ってます。

宜しくお願いいたします!

まず本日ご紹介するのは
MEDIPS

Bioconductorのサイト
http://www.bioconductor.org/packages/release/bioc/html/MEDIPS.html
論文
http://www.ncbi.nlm.nih.gov/pubmed/24227674


Methylated DNA immunoprecipitation sequencing(MeDIP-seq)の
解析ツールで、RのBioconductorのパッケージです。

MeDIP-seqはメチル化したDNA断片を5-methylcytosine (5mC)という
抗体を用いて免疫沈降させて濃縮し、シーケンシングを行うという
技術ですが、その解析に特化されているツールです。

まず、入力形式はbamファイルです。
パッケージのテストデータはヒトのES細胞で、
bowtieでマッピングしたbamファイルを使用してます。

各サンプルごとのクオリティーチェックとして、
Saturation分析、シーケンスリードのカバレッジパターンが
グラフ化されます。




そしてCase Studyでは複数サンプル間でメチル化パターンが
異なる領域を同定することが出来ます。

chr start stop ID
1 chr22 19136001 19136200 ID_1
2 chr22 19753401 19753500 ID_2
3 chr22 20785001 20785200 ID_3
4 chr22 20785301 20785600 ID_4
5 chr22 24820701 24820800 ID_5
6 chr22 24890501 24891100 ID_6
7 chr22 28035101 28035200 ID_7
8 chr22 28193101 28193200 ID_8
9 chr22 28193301 28193400 ID_9
10 chr22 29707001 29707200 ID_10


今回ご紹介するのはここまでですが、
他にChIP-seqにも活用できるとも書いてありました。

また、論文では他にも色んな図を描くことができるようなので、
ご興味のある方はチェックしてみてください。
| tokunaga | バイオインフォマティクス | 16:46 | comments(0) | - |
NBDCヒトデータベース運用開始
独立行政法人科学技術振興機構(JST)バイオサイエンスデータベースセンター(NBDC)が、ヒトデータに関するデータベースを共有するためのプラットフォームとして、「NBDCヒトデータベース」の運用を開始したそうです。

http://humandbs.biosciencedbc.jp/

公的資金を用いたプロジェクトで、ヒトを対象とした研究の成果として産生されたデータ(遺伝情報、臨床情報、画像情報など)を受け入れるそうです。
NGSデータは、DDBJのDRA(DDBJ Sequence Read Archive) / JGA(Japanese Genotype-phenotype Archive)へ格納することになると記載がありました。

データは大きく4つに分類されます。
1. オープンデータ
 すでに発表された論文の参照データなど。
2. 制限公開データ
 関連研究に従事したことのある研究者が使用できるデータなど。
3. 公開待機データ
 データ提供者による成果の公開の後、1あるいは2として公開される予定のデータなど。
4. 匿名化前・公開留保データ
 匿名化前のデータ。

データ提供者と利用者が守るべきセキュリティレベルが、データごとに異なるの点がポイントです。



利用可能な研究データはまだないようですが、順次公開されるそうです。
表現が異なる点もありますので、正確な内容はこちらをご覧ください。
| きむ | バイオインフォマティクス | 17:54 | comments(0) | - |
WebLogoを使ってみよう
WebLogoを使うと、複数配列の塩基の構成比を描画することができます。

以下に示すのは、miRBaseから取得した、ヒト、チンパンジー、ゴリラ、
マウス、ラットのmir-25配列です。

>hsa-mir-25 MI0000082
GGCCAGUGUUGAGAGGCGGAGACUUGGGCAAUUGCUGGACGCUGCCCUGGGC
AUUGCACUUGUCUCGGUCUGACAGUGCCGGCC
>ptr-mir-25 MI0003072
GGCCAGUGUUGAGAGGCGGAGACUUGGGCAAUUGCUGGACGCUGCCCUGGGC
AUUGCACUUGUCUCGGCUGAGACGCGCCCGCC
>ggo-mir-25 MI0003061
GGCCAGUGUUGAGAGGCGGAGACUUGGGCAAUUGCUGGACGCUGCCCUGGGC
AUUGCACUUGUCUCGGUCUGACAGUGCCGGCC
>mmu-mir-25 MI0000689
GGCCAGUGUUGAGAGGCGGAGACUUGGGCAAUUGCUGGACGCUGCCCUGGGC
AUUGCACUUGUCUCGGUCUGACAGUGCCGGCC
>rno-mir-25 MI0000856
GGCCAGUGUUGAGAGGCGGAGACACGGGCAAUUGCUGGACGCUGCCCUGGGC
AUUGCACUUGUCUCGGUCUGACAGUGCCGGCC


これをWebLogoのサイトに入力すると、以下のような画像が出力されます。
よく類似していますが、69-77塩基目あたりなどが少し異なるようです。



ローカル環境で実行したい場合は、ソースコードをダウンロードして
インストールし、次のように実行します。
$ ./weblogo -F png < 入力.fa > 出力.png

WebLogoは、塩基配列だけでなくアミノ酸配列でも実行できます。
| hat | バイオインフォマティクス | 14:01 | comments(0) | - |
Circosを使ってみよう
構造変異の論文などでよく見かけるこのような図は
Circosというソフトウェアで描くことができます。

Circosのサイトより

◆Circosのインストール
http://circos.ca/software/download/circos/から適当なバージョンのアーカイブをダウンロードして解凍します。

◆Circosに必要なPerlモジュールのインストール
解凍してできたcircos-0.xxディレクトリ下のbinディレクトリに
チェック用のスクリプトが用意されています。

・必要なPerlモジュール一覧を確認する
$ ./list.modules

・それらのモジュールがインストールされているかどうか確認する
$ ./test.modules

→足りないモジュールをCPANなどでインストールします。
私の環境では、CPANで force install でないとインストールできないモジュールがいくつかありました
また、先にyumでgd-develを入れておかないとGDがインストールできませんでした

◆Circosの実行
circos-0.xxディレクトリ下のexampleディレクトリに移動します。
プログラムについてきたテストデータと設定ファイルで、
Circosを実行してみました。

$ ../bin/circos -conf etc/circos.conf -outputfile mycircos.png

結果画像はPNG形式とSVG形式で出力されます。
実行ディレクトリにmycircos.png、mycircos.svgという名前で
このような画像ファイルが生成されました。


クリックで拡大します

CircosのWebサイトのサンプルページを見ていると、様々な図が
描けるようです。これは使いこなせるようになっておきたいですね。
| hat | バイオインフォマティクス | 15:08 | comments(0) | - |
Bio-SamTools
Bio-SamToolsを使ってBAMから情報を抜き出す方法をご紹介します。

例えば、あるBAMファイルのchr21:19,660,000-19,660,600の領域に
マップされたリードの「ID」と「アライメント位置」と
「ペアワイズアライメント結果」を表示したい場合には、
次のようなPerlスクリプトを実行します。

#!/usr/bin/perl
use strict;
use warnings;
use Bio::DB::Sam;

my $bam = BAMファイル;
my $obj = Bio::DB::Sam->new(-bam => $bam);

my @ali = $obj->get_features_by_location(
__-seq_id => 'chr21', -start => 19660000, -end => 19660600);

foreach my $a(@ali){
__my $id = $a->query->name;
__my ($chr, $sta, $end) = ($a->seq_id, $a->start - 1, $a->end);
__my ($ref, $mat, $que) = $a->padded_alignment;

__print "ID: $id¥n";
__print "Position: $chr:$sta-$end¥n";
__print join("¥n", $ref, $mat, $que), "¥n¥n";
}


このような結果が返ってきます。



アライメントの詳細を見たい時に便利ですね。
この他にも、CIGARやクオリティスコアや「ペアリードがマップ
できているか」など、さまざまな情報が取得できます。
| hat | バイオインフォマティクス | 15:10 | comments(0) | - |
気軽にbamファイルの中を見る
入社して四か月になり、私も簡単な解析ならできるようになりました。
とはいえ、解析中に起きたエラーについて先輩に相談したところ
「ちょっとBAMファイルをless(コマンド)で確認してみようか」
と言われ、
「BAMってバイナリデータでは……?」
と狼狽したので、まだ初心者のようです……。
(バイナリデータとは、人が読める文字で書かれたテキストデータではない、コンピュータが読みやすい二進数で記されたデータのことです。テキストデータと同じように開こうとするとエラーが出るか、無理やり普通のデータと同じように開いてみても、文字化けのようになります。下図はBAMファイルをWindows上のテキストエディタで開いてみたものです)
バイナリの文字化け

BAMファイルの中身を閲覧したいときには、samtools viewでBAMをSAMに変換します。もちろん、SAMファイルとして出力してもいいでのすが、少し中を確認したいだけのときは、パイプ(|)でlessコマンドを続け、扱いやすいテキストとして気軽に読むようにできます。

$ samtools view example.bam | less
| kubo | バイオインフォマティクス | 15:32 | comments(0) | - |
CELファイルの中身
tokunagaです。
本日はCELファイルについてお話ししたいと思います。

CELとは、AffymetrixのArrayで出力される遺伝子発現やジェノタイピングデータの含まれたファイルのフォーマットです。
CELファイルには、DATファイル(Affymetrixのスキャナーから出力されるフォーマット)のピクセル値から計算されたintensityの結果が含まれています。その他に標準偏差値、ピクセル値、外れ値のフラグ等が入っています。細かいコンテンツはArrayの種類やバージョン、使用したアルゴリズムによってそれぞれです。
これらのデータはAffimetrixが提供しているツールや、Rのパッケージを使って情報を取り出すことが可能なようです。

また中身の情報を見るだけではなく、CELファイルから様々な解析結果の図を描くことが出来るようなので、機会がありましたらその方法についてもご紹介もしたいと思います。
| tokunaga | バイオインフォマティクス | 14:43 | comments(0) | trackbacks(0) |
各リードが1回だけ登場するBAMを作るには
マッピング率を計算する時など、BAMファイルから、マッピングされたリード数を取得したい場合があると思います。

samtools idxstatsでもマッピング結果の統計を出すことができるのですが、マルチマッピングのBAMの場合は重複してカウントされてしまうようで、正確なリード数が得られません。

samtools viewに-Fオプションをつけて実行すると各リードが1回だけ登場するBAMを作成することができます。

【実行例】
samtools view -bh -F 256 -o hoge_uniq.bam hoge.bam

【解説】
・BAMの2列目は「フラグ」列です。
 同じリードが2回目以降マッピングされた場合、フラグに0x100(=256)が立ちます。
・samtools viewの-Fオプションは、「そのフラグを持つデータを除く」の意味なので、-F 256とすると「1回目のマッピング箇所のみ残す」ことができます。

このBAMに対してsamtools idxstatsをかけると、正確なマッピング数やアンマッピング数を得られます。
| hat | バイオインフォマティクス | 13:51 | comments(0) | trackbacks(0) |
GATKの「known sites」
今日は、次世代シーケンサの解析に用いられるGATKについてお話します。
GATK実行時に「known variants/sites」を指定しなければいけないコマンドが多々あります。
GATKのマニュアルにかかれたコマンド例や、論文を参考にしていましたが、網羅的にまとめているページがありました。
http://www.broadinstitute.org/gatk/guide/article?id=1247


「UnifiedGenotyper」「HaplotypeCaller」では、アノテーションに使用しているだけなので、ご自由時どうぞ
それ以外の「RealignerTargetCreator」「IndelRealigner」「BaseRecalibrator」「VariantRecalibrator」「VariantEval」は、解析の結果が変わってきてしまうので、推奨の「known variants/sites」を使用するように
と、書かれています。

特に「VariantEval」では、1000人ゲノムプロジェクトの影響を除いたほうがいいので、dbsnp129を推奨していました。
| きむ | バイオインフォマティクス | 15:18 | comments(2) | - |
BioHackathon2013
6月23〜28日に開催されたBioHackathon2013に参加してきました。

初日のシンポジウム@スカイツリーでは、日本SGI様と共同で
ドリンクコーナーのスポンサーをさせていただきました。

二日目以降はライフサイエンス統合DBセンターで、
ひたすらハッキングします。
私はH-invitationalデータベースのグループに入って、
RfamデータをH-invitationalにリンクするための調査や
データ作成を行いました。

休憩時間やバンケットではたくさんの研究者の方とたくさんお話を
させていただくことができ、とても楽しかったです。

★★★
____早く世界中の有益なデータがつながって
____簡単に検索できるようになりますように!
________________________________★★★
| hat | バイオインフォマティクス | 14:58 | comments(0) | - |
テスト用Fastqファイルを作る
マッピングソフトウェアの性能比較をする際に、
正解がわかっているデータを使いたい場合があります。

そんな時に便利な、ゲノム配列からFastqファイルを作成できる
ArtificialFastqGeneratorというソフトウェアをご紹介します。

例えばヒトゲノムhg19から150bp Paired-EndのFastqを作成するには
以下のようなコマンドを実行します。
$ java -jar ArtificialFastqGenerator.jar -R hg19.fa -O OUTNAME -RCNF 2 -S ">chr1" -E ">" -RL 150

パラメータの意味です。
・-R:ゲノム配列のMultiFastaファイル
・-O:出力ファイル名
(この例では、OUTNANE.1.fastq, OUTNANE.2.fastqができます)
・-RCNF:N入り配列を許すか
(0:N入りを許す、1:全部Nの配列は除く、
 2:1つでもNが入っている配列は除く)
・-S、-E:-Rで指定したFastaのどの配列を使うか(この例ではchr1)
・-RL:リード長

エラーを入れたり、実際のFastqを参照してリアルなクオリティスコアを
入れたりすることもできます。

とても便利なソフトウェアです。
| hat | バイオインフォマティクス | 14:14 | comments(0) | - |
hg20
あまり梅雨らしくないこの6月、いかがお過ごしでしょうか。
今年も「あっ、数日前から明けてました」と後付けで梅雨明け宣言されて、なし崩し的に夏になる予感がしています。

私は四季の中で夏が一番好きなので、夏に向かいつつあるこの時期は毎年とてもわくわくします。
野山を駆けまわって遊んで、麦茶を飲んで汗だくで昼寝した、子供のころの夏休みを思い出すからでしょうか。
ひぐらしが鳴きだすと「夏が終わってしまう...」と寂しくなってしまいます。

今年の夏も、実家に帰って野山を駆け回ったり昼寝したりしたいと思います。

ですが、今年の夏の一番の楽しみといえば、なんといっても
hg20/GRCh38が出る
ということではないでしょうか。

hg19/GRCh37と比べてどれくらい変わっているのでしょうか?
ランダム配列がどのくらい減っているのでしょうか?!
楽しみです!
| hat | バイオインフォマティクス | 16:01 | comments(0) | - |
オープンソースでパスウェイ解析
家系解析でもがん解析でも、原因遺伝子候補を絞り込んだ後のアノテーションが重要になります。
アノテーション方法の一つに、変異があった遺伝子がどのようなパスウェイ上にあるかを調べるパスウェイ解析があります。
今回、オープンソースのソフトウェアとデータを使ったパスウェイ解析の手順をご紹介します。

CytoscapeはCytoscape Consortiumが提供している、オープンソースの分子間相互作用ビューアです。
Windows版、Mac版、Linux版があります。

CytoscapeにReactome用プラグインを入れると、オープンソースのパスウェイデータベースであるReactomeを使った解析ができるようになります。

1. Cytoscapeのインストール
Cytoscapeのサイトから、お使いのマシン用のインストーラをダウンロードしてインストールします。

2. Reactome用プラグインのインストール
Cytoscape プラグインページの「Download and Launch the Reactome FI plugin」からcaBigR3.jarをダウンロードして
Cytoscapeのプラグインフォルダ(Windowsの場合デフォルトで C:¥Program Files¥Cytoscape_vx.x.x¥plugins)に置き、
Cytoscapeを再起動します。

3.テストデータのダウンロード
Cytoscape プラグインページの「Use the Reactome FI plugin」→「Gene Set/Mutation Analysis」からGeneSampleNumber.txtをダウンロードして適当な場所に置きます。
これは、遺伝子とサンプル数が2列1行になったデータです。

4. テストデータ読み込み
Cytoscapeのメニューバーの「Plugins」→「Gene Set/Mutation Analysis」を選択します。
Choose data fileでダウンロードしたテストデータ、Specify file formatは「Gene/sample number pair」を選択し、「OK」を押します。

Reactomeのパスウェイデータを使用して、遺伝子の関係が図示されます。
変異のあったサンプル数が多い遺伝子ほど大きく表示されます。


原因遺伝子候補がどのようなパスウェイ上にあるか、どのような遺伝子と関係があるか、見ることができます。
| hat | バイオインフォマティクス | 15:45 | comments(0) | trackbacks(0) |
ソフトウェア結果比較【BAMソート編・1】
NGS解析では、同じ処理を行うのにいくつものソフトウェアがあって、どれを使ったらいいのか迷うことがあります。

例えばBAMファイルをゲノムポジションでソートする場合、以下の選択肢があります。
1. samtools sortを使う
2. picard toolsのSortSam.jarを使う

同じBAMファイルを異なるソフトウェアでソートしても、結果は同じになるのでしょうか?

テストデータとして、NCBI SRAのSRR077861の先頭1000リードをhg19にBWAでマッピングした結果のBAMを使いました。

・テストデータ
sample.bam(222,146 bytes)


【samtools sortでソート】
・実行コマンド
samtools sort sample.bam sort_samt

・結果BAM
sort_samt.bam(219,312 bytes)


【SortSam.jarでソート】
・実行コマンド
java -jar /usr/local/src/java/picard-tools-1.75/SortSam.jar I=sample.bam O=sort_pica.bam SO=coordinate

・結果BAM
sort_pica.bam(221,293 bytes)


結果BAMが異なりました!
いったいどこが異なるのでしょうか。

比較結果は次回で。
| hat | バイオインフォマティクス | 10:42 | comments(0) | trackbacks(0) |
Mendeleyの紹介
本日は、文献管理ソフトの一つであるMendeleyを簡単にご紹介したいと思います。

最近、論文を読む機会が増えてきたので論文管理ソフトを探していました。論文管理ソフトと言えばEndNoteが有名ですが、値段が高いことがネックです。そこで、フリーの論文管理ソフトを探してネットをさまよっていたところ、このMendeleyを発見しました。
基本的には、ダウンロードしてきたPDFを登録して一覧で管理できるソフトなのですが、私が便利と思っている、以下の機能があります。


・PDFの登録が簡単
自分のPC内にあるPDFフォーマットの論文をドラッグ&ドロップするだけで、タイトルや著者名、発行年や巻数などを自動でimportしリスト化してくれます。

・クラウドベースである
異なる場所にある端末からログインしても、別の端末で集めた論文を読むことができます。Evernoteと同じ感覚でしょうか。

・ウェブ上の論文を簡単に取り込める
Web importerというプラグインを導入すると、ウェブブラウザで表示した論文をクリック一つで取り込めます。主要な論文誌、またGoogle scholarやPubMed、CiNii等の検索結果からも取り込むことができます。

・Citationの書式を選べる
Mendeleyでは論文のcitationの書式がいくつか選べます(まだ数は少ないですが)。これをそのままwordに張り付けることもできるので、論文執筆の際には大きな助けになるでしょう。

・論文にマーカーやメモが張り付けられる
Mendeley上で表示したPDFには、マーカーでラインを引くことや付箋のようなメモを貼り付けることができます。

・オンラインコミュニティがある
MendeleyのページからMendeleyのオンラインコミュニティに参加することができます。研究者同士の交流や、情報交換が可能になるでしょう。また自分の興味のある分野を登録しておくと、その分野の論文を自動的に検索してくれる機能もあります。


以上、ざっと書いてきましたが、他にも便利な機能がたくさんあるようです。Mendeleyの使い方を纏めたページなどもありますので、ご興味がある方は是非ご覧ください。
| deda | バイオインフォマティクス | 17:22 | comments(0) | trackbacks(0) |
coverageBedの使い方(2)
coverageBedの使い方(1)のつづきです。

■カバレッジ計算
マッピング結果とゲノム、2つのBEDが用意できたら、以下のコマンドでカバレッジを計算します。

$ coverageBed -a map.bed -b genome.bed

以下のような結果が出力されます。

chr1___0__249250621__1__100__249250621__0.0000004
chr10__0__135534747__0____0__135534747__0.0000000
__:
chr2___0__243199373__1__100__243199373__0.0000004
chr20__0___63025520__0____0___63025520__0.0000000
__:

出力項目は、左から
1列目「genome.bedの1列目」
2列目「genome.bedの2列目」
3列目「genome.bedの3列目」
4列目「genome.bedにオーバーラップしたmap.bedの領域数」
5列目「カバレッジが0より大きいgenome.bedの塩基数」
6列目「genome.bedの領域長(=3列目と同じ値)」
7列目「map.bedによりオーバーラップされたgenome.bedの割合」
です。

7列目を見れば、染色体ごとのマッピング率がわかります。

ゲノム全体のマッピング率は、以下のコマンドで計算できます。
(今回のテストデータでは、6.46059e-08になりました)

$ coverageBed -a map.bed -b genome.bed | awk '{a+=$3;b+=$5} END {print b/a}' -

-bで指定するBEDを特定の領域(Exomeのキャプチャ領域など)にすれば、その領域に対するカバレッジを計算できますし、
以下のようなゲノム領域を10kb刻みで区切ったBEDにすれば、10kbごとのカバレッジを計算することができます。

chr1____0__100
chr1__100__200
chr1__200__300
__:

また、-aの代わりに-abamを使うと、BAMファイルを指定することもできます。
$ coverageBed -abam map.bam -b genome.bed
| hat | バイオインフォマティクス | 19:39 | comments(0) | trackbacks(0) |
coverageBedの使い方(1)
以前bedtoolsの一つであるintersectBedの使い方についてご紹介しましたが、今回はcoverageBedについてご紹介します。

マッピング結果がどのくらいゲノム全体をカバーできているか知りたい時、coverageBedを使うとカバレッジが簡単に計算できます。

■用意するもの
・マッピング結果のBEDファイル
__--マッピング領域の染色体名・スタート・エンド

BEDは1行に最低3列、オプションで12列まで書くことができるフォーマットですが、coverageBedでは3列のBEDしか受け付けないようです。

列数が3より多い場合は、以下のコマンドで列数を3列だけにしておきます。
$ awk '{print $1"¥t"$2"¥t"$3}' original_map.bed > map.bed

今回、マッピング結果のBEDファイルとして以下のようなファイルを作成しました。
chr1__0__100
chr2__0__100

・ゲノムのBEDファイル
__--各染色体の染色体名・スタート・エンド

ゲノムのBEDファイルは、ゲノムのfastaファイル(ここではgenome.faという名前のものを使いました)から以下のコマンドで作ることができます。
染色体の長さが分かれば手作業で作ってもよいです。

$ samtools faidx genome.fa
__→genome.fa.faiというファイルができます

$ awk '{print $1"¥t0¥t"$2}' genome.fa.fai > genome.bed
__→以下のような内容のgenome.bedができます

chr1__0__249250621
chr2__0__243199373
__:
chrY__0__59373566
chrM__0__16571

以上で準備は完了です。次回でカバレッジの計算を行います。
| hat | バイオインフォマティクス | 19:02 | comments(0) | trackbacks(0) |
鎖鋸
皆様こんにちは。detです。

本日は、とあるデータベースをご紹介したいと思います。

タイトルにもある鎖鋸とは、DRA/SRAなどのリードアーカイブの中から、論文が発表済みのリードデータのみを収集したデータベースです。DRAやSRAには論文が発表されていないリードデータも含まれるため、リードの質をチェックするために、こちらのデータベースで確認してみるのもいいかもしれません。
では、トップページにアクセスしてみましょう。 

何とも言えないキャラクターが登場します。



キャラクターの横に、検索窓がありますので、私がDRAで適当に検索した
SRA038332 と入力してみます。
もし論文があれば、タイトルとアブスト等の情報が出てきます。
また下の方には、リードに対してFastQCをかけた結果も表示されています。



以上のように、IDを入れるだけで、論文情報とリードクオリティを確認することができます。
まだ、アップされていない情報もあるようですので、今後の発展に期待したいですね。

それでは、また。
| deda | バイオインフォマティクス | 18:05 | comments(0) | trackbacks(0) |
統合データベース講習会に参加してきました
皆様こんにちは。detです。

先週の土曜と日曜に、NBDC主催の「統合データベース講習会(AJACS駿河)」に参加してきました。

今回は生命科学系のデータベースに関する講演だけでなく、超高速シーケンサーのデータ解析パイプラインに関する講演もあり、様々なお話が伺えて非常に勉強になりました。

今日は、その講習会で知ったKNApSAcKという面白いデータベースについてご紹介したいと思います。これは、奈良先端科学技術大学院大学の金谷先生が開発されているもので、薬、生物、食品、ヒトの様々な情報に関するデータベースの集合体として公開されています。このKNApSAcKのHPにアクセスすると、ポップな柔らかい感じのトップページが現れます。特に私が面白いと思ったのは、左にあるDietNavi(病気予防データベース)です。その中の各病気をクリックすると、その病気を予防する効果がある食物の詳細なリストが表示されます。このデータベースを活用すれば、様々な病気の予防に一役買ってくれるかもしれません。
他にも、食べ合わせリスト等のユニークなデータベースがあり、いろいろ見てみるのも楽しいでしょう。

まだまだ構築途中のようですので、今後の発展に期待したいと思います。
| deda | バイオインフォマティクス | 17:23 | comments(0) | trackbacks(0) |
intersectBedの使い方
bedtoolsBEDフォーマットのファイルを扱うのに便利なツール群です。

今回はその中の1つ、intersectBedについてご紹介します。
intersectBedを使うと、複数BED間で重複している領域を簡単に抽出することができます。

テストデータとして、2つのBEDファイルを用意しました。

==> A.bed <==
chr1 10 30 A_1
chr1 40 50 A_2
chr1 70 80 A_3

==> B.bed <==
chr1 10 15 B_1
chr1 20 25 B_2
chr1 30 40 B_3
chr1 60 90 B_4

A.bedとB.bedの位置関係は次のようになります。


-aで指定したBED中の要素のうち、-bで指定したBED中の要素にオーバーラップするものだけを抽出するには、-waオプションをつけます。

$ intersectBed -a A.bed -b B.bed -wa
chr1 10 30 A_1
chr1 10 30 A_1
chr1 70 80 A_3

→A_1とA_3だけが出力されました。
 A_1が2回出力されるのは、2回オーバーラップしているためです。

-uオプションをつけると、重複した結果がマージされます。

$ intersectBed -a A.bed -b B.bed -wa -u
chr1 10 30 A_1
chr1 70 80 A_3

→A_1の出力が1回だけになりました。

-waの代わりに-wbを使うと、bの情報も出力されます。

$ intersectBed -a A.bed -b B.bed -wb
chr1 10 15 A_1 chr1 10 15 B_1
chr1 20 25 A_1 chr1 20 25 B_2
chr1 70 80 A_3 chr1 60 90 B_4

→左4列は「a中のオーバーラップ領域」、
 右4列は「オーバーラップしたbの情報」です。

-fオプションを使うと、オーバーラップの割合を指定できます。
 例えば、全長の3割以上がbにオーバーラップされるものだけを抽出するには、-f 0.3と指定します。


$ intersectBed -a A.bed -b B.bed -wa -f 0.3
chr1 70 80 A_3

→A_3だけが出力されました。
 A_1(20bp)は、B_1によって5bp、B_2によって5bpで計10bpオーバーラップしているのですが、出力されません。おそらく積算ではなく、1要素によるオーバーラップ領域長で見ているのだと思われます。


-vオプションをつけると、オーバーラップ「しなかった」ものが出力されます。
$ intersectBed -a A.bed -b B.bed -wa -v
chr1 40 50 A_2

-cオプションをつけると、オーバーラップしたbの数も出力されます。
$ intersectBed -a A.bed -b B.bed -wa -c
chr1 10 30 A_1 2
chr1 40 50 A_2 0
chr1 70 80 A_3 1

-waoオプションをつけると、オーバーラップしたbの情報と、オーバーラップした塩基数が出力されます。
$ intersectBed -a A.bed -b B.bed -wao
chr1 10 30 A_1 chr1 10 15 B_1 5
chr1 10 30 A_1 chr1 20 25 B_2 5
chr1 40 50 A_2 . -1 -1 . 0
chr1 70 80 A_3 chr1 60 90 B_4 10

以上、A.bedを主体にした結果を出しましたが、Bを主体にした結果を出したい場合は、-a B.bed -b A.bedのように、-aと-bの指定をひっくり返せばOKです。

また、今回のテストBEDにはストランド情報がないので実行していませんが、-sオプションで「同じ向きに限定」、-Sで「逆の向きに限定」することができます。

すごいのはこのintersectBed、BAM・VCF・GFFにも使えることです(※BAMの場合は-aの代わりに-abamでファイルを指定してください)。
マッピング結果への簡単なアノテーション付けくらいならこれだけで行えてしまうくらいなので、まだ使ったことがない方はぜひお試しください。

| hat | バイオインフォマティクス | 11:23 | comments(0) | trackbacks(0) |
これ、なんて読みますか
ゲノム解析の定番ソフトウェアに「Picard」というのがあります。
私はずっと「ピ・カード」と呼んでいましたが
学会などでは「ピカール」と呼ばれるのをよく聞きます。
正式には「ピカール」なんでしょうか?

マイクロRNAのデータベース「miRBase」も、
私はずっと「ミルベース」と呼んでいたのですが
「ミベース」と呼ばれることが多いことを最近知りました。

この他にも私が勝手な呼び方をしているものがありそうです。
いくつかのバイオインフォマティクス用語の私の読み方を書いてみました。
「それは間違っている」「こう呼ぶのが公式」
「自分はこう呼んでいる」というのがありましたら教えてください。

FASTA「ファスタ」
FASTQ「ファスト・キュー」
VCF「ブイ・シー・エフ」
BWA「ビー・ダブリュ・エー」
Bowtie「ボタイ」
SOAP「ソプ」
GATK「ジー・エー・ティー・ケー」
IGV「アイー・ジー・ブイ」
Exon「エソン」

この中で、一番意見が分かれ、結論が出ていないのが
Exon(エキソン/エクソン)だと思います。

ちなみに私が「エキソン」派なのは、
「キ」のほうがmRNAがシャープにスプライシングされる
感じが伝わると思うからです。

みなさんはどちら派でしょうか?
| hat | バイオインフォマティクス | 15:04 | comments(0) | trackbacks(0) |
Rでアノテーション付け
tokunagaです。

これまでにVCFtoolsを使ったVCFファイルの加工や比較の方法をご紹介しました。
VCFtools
VCFtools
今回はVCFファイルに関するRパッケージをご紹介いたします。

VariantAnnotation

Bioconductorのパッケージです。
VCFフォーマットのファイルを入力することができ、変異の概要やアミノ酸の変化、またはSIFTやPolypenスコアなどのアノテーション情報を見ることができます。
マニュアルには比較的分かりやすい使用手順が記載されています。

機会がありましたらまた詳しい情報を載せたいと思います。

まだR初心者ですが、Rには色んな便利なパッケージがありますので、しっかり勉強してご紹介いていきたいと思っております。
| tokunaga | バイオインフォマティクス | 07:50 | comments(0) | trackbacks(0) |
PBSIM
tokunagaです。
本日はBioinfomaticsで気になる記事を見つけましたのでご紹介いたします。

Bioinformatics. 2012 Nov 4
PBSIM: PacBio reads simulator–toward accurate genome assembly
Ono Y, Asai K, Hamada M


PacBioのシーケンサーから出力されるリードの特徴としてCLR(長いがエラー率の高いリード)とCCS(短いがエラー率の低いリード)があるそうです。その特徴を考慮し、model-basedとsampling-based methodというアルゴリズムを使ってゲノムのアセンブルを行うシュミレーターについての記事です。

以下のサイトからLinux用にコンパイル済みのツールとソースコードがダウンロードできます。
http://code.google.com/p/pbsim/downloads/list

もしご興味のある方は試してみてはいかがでしょうか?
また、気になるツールを見つけましたらご紹介いたします。
| tokunaga | バイオインフォマティクス | 16:11 | comments(0) | trackbacks(0) |
VCFtools
tokunagaです。

本日は以前ご紹介したVCFtoolsでちょっと気になっていたコマンドを調べましたのでご紹介したいと思います。
vcf-compareというVCFファイル同士の簡単な比較を行ってくれるコマンドです。

前回ご紹介したように前処理としてVCFファイルをbgzipで圧縮して、tabixでインデックスを付けた後、以下のコマンドを実行します。
vcf-compare A.vcf.gz B.vcf.gz

実行すると以下のような情報が標準出力で出てきます。
# This file was generated by vcf-compare.
# The command line was: vcf-compare(r731) A.vcf.gz B.vcf.gz
#
#VN 'Venn-Diagram Numbers'. Use `grep ^VN | cut -f 2-` to extract this part.
#VN The columns are:
#VN 1 .. number of sites unique to this particular combination of files
#VN 2- .. combination of files and space-separated number, a fraction of sites in the file
#2サンプル間で位置情報が共通していたSNV/Indel数と割合
VN 106845 A.vcf.gz (39.2%) B.vcf.gz (44.3%) 
#Bに特有なSNV/Indel数と割合
VN 134310 B.vcf.gz (55.7%) 
#Aに特有なSNV/Indel数と割合
VN 165380 A.vcf.gz (60.8%) 
#SN Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
#位置情報が共通し、且つREFの塩基が一致した数
SN Number of REF matches: 106339 
#ALTの塩基が一致した数
SN Number of ALT matches: 105653 
#REFの塩基が一致しなかった数
SN Number of REF mismatches: 506 
#ALTの塩基が一致しなかった数
SN Number of ALT mismatches: 686 
#遺伝子型を比較したサンプル数
SN Number of samples in GT comparison: 0 
# Number of sites lost due to grouping (e.g. duplicate sites): lost, %lost, read, reported, file
#位置情報が重複し除外されたSNV/Indelや割合
SN Number of lost sites: 157 0.1% 272382 272225 A.vcf.gz 
#位置情報が重複し除外されたSNV/Indelや割合
SN Number of lost sites: 167 0.1% 241322 241155 B.vcf.gz 

2サンプル以上でも実行することが可能です。
今回はデフォルトで実行しましたがオプションも結構な数あるようです。
とりあえずVCFファイルを大まかに比較したい場合には便利かもしれません。

また、機会がありましたら他のコマンドもご紹介したいと思います。
| tokunaga | バイオインフォマティクス | 12:26 | comments(0) | trackbacks(0) |
VCFtools
tokunagaです。

今日はVCFtoolsについてご紹介したいと思います。
URL:http://vcftools.sourceforge.net/
VCFToolsは、NGSのデータ解析で出力されたVCFファイルを加工するのに役に立つツールです。

VCFtoolsを使用する際にはbgzip、tabixの使用が必要となります。

例えばマージする場合
まずbgzipで固めます。
bgzip $FILE1.vcf
bgzip $FILE2.vcf


次にtabixでインデックスを付けます。
tabix $FILE1.vcf.gz
tabix $FILE2.vcf.gz


これでVCFToolsを実行することが出来ます。
vcf-merge $FILE1.hetero.gz $FILE2.hetero.gz|bgzip -c > $OUT.vcf.gz


前処理はちょっと手間はかかりますが、複数のVCFファイルのマージや比較、共通部分の抜き出しなど様々なコマンドがあります。

また機会がありましたら便利なコマンドをご紹介したいと思います。
| tokunaga | バイオインフォマティクス | 18:42 | comments(0) | trackbacks(0) |
便利なオプション
akbです。
今日はLinuxコマンドの便利なオプションを
ご紹介したいと思います。


.侫.ぅ襪筌妊レクトリの検索
$ find [パス][オプション][ファイルまたはディレクトリ名]

【optionの説明】
-name <パターン>: ファイル名がパターンと同じファイルを検索する

【実行例:/home/以下を指定してhoge.txtを検索する】

$ find /home -name hoge.txt

▲妊レクトリの作成
$ mkdir [オプション][ディレクトリ名]

【optionの説明】
-p: ディレクトリを含むディレクトリを作成する

【実行例:Tokyoディレクトリを含むJapanディレクトリを作成する】

$ mkdir -p Japan/Tokyo

ファイルの末尾(デフォルト:10行)を表示する
$ tail [オプション][ファイル名]

【optionの説明】
-f: 更新中のファイルの末尾をリアルタイムに表示する

ぅレンダーを表示する
$ cal
【実行例】


テ時を表示する
$ date
【実行例】



恥ずかしながら最近"cal"コマンドを知りました。ついつい使いたくなるコマンドです。
是非、お試しください。
| akb | バイオインフォマティクス | 16:08 | comments(0) | trackbacks(0) |
Galaxy による QC
こんにちは。detです。

今日はGalaxyを用いたQCについてご紹介いたします。

これまで、このブログで「QCの道」というタイトルでFASTX-Toolkitの使い方をご紹介してきました。

今回は、このQC機能をGalaxy上で実行しつつ、Galaxyの基本的な使い方を紹介したいと思います。

それでは、まずGalaxyのパブリックページを開きましょう。



上記のようなページが開いたら、まずは右上の[user]をクリックして、登録を選択します。



上記のようなページが開きますのでユーザー登録しましょう。



ユーザー登録が終わったら、左パネルの一番上にある、Get Data の中の Upload File を選択します。すると、真ん中の画面に上記のようなページが出現します。今回はfastq 形式のデータをQCしますので、File Formatには、fastqsangerを選択してください。次に[ファイルを選択]をクリックして、処理したいfastqファイルを選択します。選択できたら、下の方にある、青い[Execute]ボタンをクリックします。



すると、上記のように、右パネルにデータがアップロードされます。



次に、左パネルに沢山あるメニューの中から、NGS: QC and manipulation
という項目をクリックします。するとさらにメニューが表示されますので、その中から、GENERIC FASTQ MANIPULATION の FASTQ Trimmer を選択します。



すると、上記のようなページが現れます。
これは、以前の記事でご紹介した FASTX-Toolkit の fastx_trimmer とほぼ同じ機能を持ちます。指定した塩基の数だけ、末端からトリミングされます。では、上記のように、ファイルと、Offsetの数値を指定して、Execute してみましょう。



すると、上記のように右パネルに処理後のデータを示す項目が出現します。ファイル名すぐ横の目玉マークをクリックすると、下の図のように、トリミングされた結果が表示されます。



また、左下のフロッピーアイコンをクリックすると処理後のデータをダウンロードできます。

Galaxyは基本的には以上の操作を繰り返すことでグラフィカルに処理を行うことができます。

今日はこのくらいで。detでした。
| deda | バイオインフォマティクス | 20:38 | comments(0) | trackbacks(0) |
Galaxy の紹介 その3
こんにちは。detです。

本日は先日ご紹介したGalaxyに関する記事の続きです。

Galaxyは、基本的には、公開のパブリックサーバーにジョブを投げて利用することになります。しかし、重いデータや、外部に出したくないデータを解析したいこともあるでしょう。特に次世代シーケンサーから得られる出力は膨大になりがちです。
解決策として、Galaxyを自分の手元にあるローカルなサーバーにインストールして使うことができます。

Galaxyの動作メカニズムを簡単に説明します。まずユーザーがデータをGalaxyサーバーに投げると、apacheやnginx等のwwwサーバー(proxy)がそのデータをGalaxy本体へ渡します。GalaxyはそれらのデータをMySQLやPostgreSQLなどのSQLデータベースで保管し、適宜データを取り出しながらジョブを実行します。
つまり、Galaxy本体のインストールとは別に、SQLサーバーとproxyのインストールが必要になってきます。
(このあたりの概要はこちらの資料をご参照ください)

Galaxyの開発者はSQLサーバーにはPostgreSQL、proxyにはnginxを用いることをお勧めしているため、これらの組み合わせでGalaxyのローカル環境を構築するのが、安定稼働への近道でしょう。

Galaxy Wiki 内の以下のページにローカルへのインストール方法について詳しく載っています。これらを参考にGalaxyサーバーを構築することができます。

・Galaxyの入手と起動まで
・詳細設定について

構築途中に上手くいかず困ったときは、Galaxyのメーリングリストで質問するか、過去の質問を検索するとヒントが得られるかもしれません。
| deda | バイオインフォマティクス | 11:28 | comments(0) | trackbacks(0) |
Galaxy の紹介 その2
こんにちは。detです。

本日は以前ご紹介したGalaxyに関する記事の続きです。

先月下旬にシカゴでGalaxyに関する国際会議(GCC2012)が開かれました。全世界から数百人の参加者が集まり、活発なディスカッションが繰り広げられたようです。次世代シーケンサーの登場による、解析基盤としてのGalaxyの盛り上がりは今後ますますホットになっていくと思います。

この国際会議の前日にはGalaxyのトレーニングワークショップ(WS)が開催されました。本日は、この中から最も初心者向けのWSと思われるWS2: Introduction to Galaxy の内容をご紹介したかったのですが、残念ながらまだ当日のスライドが公開されていません。

そこで、本日は、GCC2012関連のスライドではありませんが、Galaxyを用いたNGSデータ解析に関する資料をご紹介したいと思います。

まずは、こちらの資料。Galaxyの概要について大まかに知ることができます。次に、Galaxyを用いたNGSデータ解析の資料。この資料では、以下のデータ解析項目についてGalaxyの画面キャプチャを豊富に取り入れながら紹介しています。(これらの資料が開けない場合は時間をおいて試してみてください。)

・Prepare, Quality Check and Manipulate FASTQ reads
・Mapping
・SAMTools
・SNP & INDEL analysis
・Peak Calling / ChIP-seq
・RNA-seq analysis

Galaxyを利用することで、Linuxに関する知識がなくてもこれらのNGSデータ解析が可能になります。非常に便利なツールと言えるでしょう。

次回もGalaxyについて紹介していきたいと思います。
今日はこのくらいで。detでした。
| deda | バイオインフォマティクス | 10:57 | comments(0) | trackbacks(0) |
Galaxy の紹介 その1
こんにちは。detです。

本日はゲノムデータ解析インターフェイスツールであるGalaxyをご紹介したいと思います。



Galaxyはペンシルベニア州立大学のNekrutenko labとエモリー大学のTaylor labの共同で開発されているシステムであり、ウェブ上で利用できます。このGalaxyは様々なフリーのツールを内包しており、例えば raw data のクオリティコントロールから、bwaやbowtieでのマッピング、samtoolsやGATKを用いたSNV/InDel検出、MACSを用いたピーク検出などをGUIで利用することができます。
Linuxのコマンドラインに親しみが薄い方にも敷居の低い、便利なシステムと言えるでしょう。

また、このGalaxyはスタンドアロン版もあるため、自分の手元にあるサーバー上で動かすこともできます。そのため、外部に出せないデータやアップロードするには大きすぎるデータの解析も可能です。

Galaxyはチュートリアルが非常に充実しており(沢山あって困ってしまうぐらい)それらを有効利用することで誰でも解析が出来るように工夫されています。

今日は、Galaxyビギナー向けのチュートリアルを一つご紹介します。それは、Galaxyページ内のGalaxy 101というチュートリアルです。このチュートリアルでは初歩的(基本的)なGalaxyの使い方を学ぶことができます。このチュートリアルは動画でも見ることができます。

今後Galaxyについては、少しずつこのブログで紹介していきたいと思っています。
| deda | バイオインフォマティクス | 15:18 | comments(0) | trackbacks(0) |
QCの道 その7
こんにちは。detです。
今日はQCの道 その6の続きです。

FASTX-Toolkitの使い方について、引き続き紹介いたします。

・fastx_collapser
FASTA/Q ファイルの中で、同じ配列のリードが重複して存在していた場合に1つを除いてすべて削除します。入力がFASTQ形式の場合は強制的にFASTA形式で出力され、ID行は連番で振り直されます。重複していたリードのID行には[連番-重複数]の形式で重複数が記載されます。

【optionの説明】
-h: ヘルプを表示します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。
-v: 処理前後でのリード数などを出力してくれます。

【実行例: 重複リードを削除】
$ fastx_collapser -i input.fastq -o out.fasta -v -Q 33


・fastx_reverse_complement
FASTA/Q ファイルの各リードの塩基配列をその相補鎖で置き換えます。FASTQファイルを入力に選んだ場合はクオリティ行も対応するように並べ替えられます。

【optionの説明】
-h: ヘルプを表示します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。
-v: 処理前後でのリード数などを出力してくれます。

【実行例: 相補鎖の作成】
$ fastx_reverse_complement -i input.fastq -o out.fastq -v -Q 33


---これまでの記事へのリンク---
QCの道 その1
QCの道 その2
QCの道 その3
QCの道 その4
QCの道 その5
QCの道 その6
| deda | バイオインフォマティクス | 17:34 | comments(0) | - |
cmpfastq と cmpfastq_pe
以前の記事cmpfastqというfastqファイルのペアエンドリードを揃えるツールをご紹介いたしました。本日は、cmpfastqの改良版である
cmpfastq_peとの相違についてご紹介いたします。

cmpfastq_peはcmpfastqから以下の点を改良したバージョンのようです。

1. Illumina CASAVA v1.8の出力に対応
2. メモリ使用量を削減

ペアエンドリードを揃えるアルゴリズム自体は、cmpfastqおよびcmpfastq_peで同一であり、片方のfastqファイルの各リードのIDをハッシュに保存し、そのハッシュに保存されたIDがもう片方のfastqファイルに含まれるかどうかでペアエンドリードを揃えています。
cmpfastq_peでは、そのIDを抜き出す際の正規表現を工夫することで、CASAVA v1.7 以前と v1.8 に対応しています。
また、cmpfastq_peでは、ハッシュの処理等に工夫があり、メモリ使用量・処理時間などのコストがかなり削減されています。

cmpfastqが重くて動かないときは、cmpfastq_peを使用してみるのもいいかもしれません。
| deda | バイオインフォマティクス | 15:37 | comments(0) | trackbacks(0) |
シーケンサーとしての私
こんにちは、hatです。
趣味で1年ほど前から三味線教室に通っています。

先日お稽古をしていて、楽器の演奏はシーケンシングに通じるものがあるなあと思いました。
シーケンサーが塩基を読んで結果データを出すように、私は楽譜を読んでメロディーを出している!

そう考えると、シーケンサーとしての私の性能は酷いもので、

・シーケンスエラーが多い(間違える、すっとばす)
・シーケンス速度が遅い(速く弾けない)
・ロングリードは読めない(楽譜をめくる時止まる)
・Depthが薄い(稽古不足)

と、とても使い物にならないレベルです。

今度複数名で弾く演奏会があるのですが、一曲通して弾くのは難しそうです。
そこで、せめてお気に入りのフレーズだけはちゃんとに弾こう!と目的をターゲットリシーケンシングに切り替えることにしました。

普段からシーケンサーの性能について好き勝手な文句ばかり言っている私ですが、シーケンサーの中の人もきっと大変なんだろうなあと思います。

| hat | バイオインフォマティクス | 13:42 | comments(0) | trackbacks(0) |
QCの道 その6
こんにちは。detです。
今日は前回のQCの道 その5の続きです。

FASTX-Toolkitの使い方について、引き続き紹介いたします。

・fastx_quality_stats
FASTA/Q ファイルのリードに含まれる塩基のポジション毎の統計量を算出し、表形式で出力してくれます。FASTA形式の入力を与えたときは、クオリティに関する項目は出力されません。

【出力される各項目の説明】
column: 5'末端から数えた塩基の位置
count: そのポジションの塩基数
min: そのポジションで最も低いqv
max: そのポジションで最も高いqv
sum: そのポジションのqvを全て足した値
mean: そのポジションのqvの平均値
Q1: そのポジションのqvの第1四分位数
med: そのポジションのqvの中央値
Q3: そのポジションのqvの第3四分位数
IQR: そのポジションのqvの四分位数範囲
lW: そのポジションのqvの統計量を箱ひげ図で表した時の最小値
RW: そのポジションのqvの統計量を箱ひげ図で表した時の最大値
A_Count: そのポジションに含まれるA塩基の数
T_Count: そのポジションに含まれるT塩基の数
G_Count: そのポジションに含まれるG塩基の数
C_Count: そのポジションに含まれるC塩基の数
N_Count: そのポジションに含まれるN塩基の数
max-count: 全てのポジションで最も塩基数が多いポジションの塩基数

【optionの説明】
-h: ヘルプを表示します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。

【実行例: FASTA/Qファイルの統計量計算】
$ fastx_quality_stats -i test.fastq -o out.text -Q33


・fastq_quality_boxplot_graph.sh
FASTQファイルのリードに含まれる塩基のポジション毎のクオリティ図を箱ひげ図の形式で作成します。inputに、上記 fastx_quality_stats の出力を必要とします。

【optionの説明】
-i: 入力ファイルには fastx_quality_stats の出力ファイルを指定します。
-o: 出力ファイルを指定します。デフォルトでpng形式で出力されます。指定しない場合は、バイナリを無理やり標準出力しようとするので、必ず指定してください。
-p: 出力ファイル形式を PostScript に変更します。
-t: 出力される図に記載されるタイトルを指定できます。

【実行例: FASTQファイルのクオリティ図作成】
$ fastq_quality_boxplot_graph.sh -i out.txt -t hogehoge -o boxplot.png

【実行結果例】



・fastx_nucleotide_distribution_graph.sh
FASTQファイルのリードに含まれる塩基のポジション毎に、塩基の種類の分布図を作成します。inputに、上記 fastx_quality_stats の出力を必要とします。

【optionの説明】
-i: 入力ファイルには fastx_quality_stats の出力ファイルを指定します。
-o: 出力ファイルを指定します。デフォルトでpng形式で出力されます。指定しない場合は、バイナリを無理やり標準出力しようとするので、必ず指定してください。
-p: 出力ファイル形式を PostScript に変更します。
-t: 出力される図に記載されるタイトルを指定できます。

【実行例: FASTQファイルのクオリティ図作成】
$ fastx_nucleotide_distribution_graph.sh -i out.txt -t hogehoge -o boxplot.png

【実行結果例】


それでは、また次回に。


---これまでの記事へのリンク---
QCの道 その1
QCの道 その2
QCの道 その3
QCの道 その4
QCの道 その5
| deda | バイオインフォマティクス | 17:29 | comments(0) | trackbacks(0) |
QCの道 その5
こんにちは。detです。
今日は前回のQCの道 その4の続きです。

FASTX-Toolkitが持つ機能について、引き続き紹介いたします。

・fastx_artifacts_filter
FASTA/Q の各リードにおいて、塩基が特定の種類に偏っている場合にそのリードを除去してくれます。他の塩基を四個以上含んでいる場合は、許容されます。

【optionの説明】
-h: ヘルプを表示します。
-z: 出力をgzip形式で圧縮します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。
-v: 処理前後でのリード数などを出力してくれます。

【実行例: 塩基の偏りをチェック】
$ fastx_artifacts_filter -v -i test.fastq -o out.fastq -Q 33


・fasta_nucleotide_changer
FASTA/FASTQ ファイルの、U塩基をT塩基へ、もしくはT塩基をU塩基へ置換してくれるツールです。出力は全てFASTA形式へ変換されます。

【optionの説明】
-h: ヘルプを表示します。
-r: T塩基をU塩基へ置換します。
-d: U塩基をT塩基へ置換します。
-z: 出力をgzip形式で圧縮します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。

【実行例: T塩基からU塩基へ置換】
$ fasta_nucleotide_changer -d -i test.fastq -o out.fastq -Q 33

それでは、また次回に。


---これまでの記事へのリンク---
QCの道 その1
QCの道 その2
QCの道 その3
QCの道 その4
| deda | バイオインフォマティクス | 09:37 | comments(1) | trackbacks(0) |
Perl デバッグ
さて、あなたがPerlでプログラムを作成し、use strict や use warnings で怒られないところまで進んだとします。
次に、やることはプログラムを実際に動かすところですね。
ところが、実は、ここでも上手くいかないことが多々あります。
出るはずのない数字、結果、エラー・・・。
そのプログラムにはバグが潜んでいるのです。

バグとは、コンピュータプログラムの製造(コーディング)上の誤りや欠陥を表します。(wikipediaより)

プログラムが長くなればなるほど、バグを含む確率は上がります。
どこでおかしな値が出ているのか、どこで動作がおかしくなっているのか、それを突き止めるのは大変な時もあります。

そこで、お勧めするのが、デバッガと呼ばれる、バグ修正を助けてくれるツールです。
perlを実行するときに、以下のように -d をつけて実行してみましょう。そうすると、そのプログラムはデバッグモードで実行することができます。

$ perl -d hoge.pl

すると、実際の命令行の冒頭が表示されます。そこで、以下のコマンドを用いながら、プログラムの動作を確認していきます。

n: 一行実行
n と一回打ち込んで enter を押すと一行ずつ実行されます。ソースの内容がそのまま表示されるので、どの行を実行しているのかもよくわかります。ちなみに、n の場合はサブルーチンは飛ばして進みます。

s: ステップイン
基本は上の n と同じで一行ずつ実行していきますが、s を使った場合は、サブルーチンの中も辿っていきます。

enter: 直前の s か n を繰り返す
一回 s か n を実行した後は、enterキーを押すだけでその処理を繰り返すことができます。

r: サブルーチンから抜け出す
s でサブルーチンに入ったのはいいけど、同じ繰り返しが何回もあって、早く元のルーチンに戻りたい・・・という時に、r コマンドでサブルーチンを脱出できます。

c: 指定行まで実行
例えば、c 300 とすれば、300行目まで一気に実行し、そこで止まってくれます。また、プログラム中にブレークポイントを設定しておけば、c と打つだけでそのブレークポイントの次まで一気に実行できます。

ちなみにブレークポイントはソースの中に次の一文を追加すると設定できます。
$DB::single = 1;

p: 変数を表示
p $hoge とすれば、$hoge の中身を表示してくれます。これはとても便利で、450行目の時の $hoge の中身が知りたい時は、c 450 と打った後、p $hoge と打つとその中身がわかります。

q: 終了
q と打てば、デバッガは終了します。

他にも色々なコマンドがありますので、興味がある方は是非ネットでいろいろ調べてみてください。

あなたのバグが早くとれますように・・・。
| deda | バイオインフォマティクス | 17:21 | comments(0) | trackbacks(0) |
fastq format
今日から次世代シーケンサー解析で使われているformatについて書いていきたいと思います。
まず、今回はfastq formatについてです。

・構成
塩基配列と各塩基に対するquality valueの情報が書かれているテキストファイルです。各リードは4行で構成されています。


[例]
 @SRR022885.1 BI:080102_SL-XAR_0001_
  FC201E9AAXX:6:1:752:593
◆CGTACCAATTATTCAACGTCGCCAGTTGCTTCATGT
 +
ぁIIIIIIIIII>IIIIIII@IIII.I+I>35I0I&+/

次回に続きます。
| tokunaga | バイオインフォマティクス | 14:49 | comments(0) | trackbacks(0) |
Perl プログラムの性能解析(NYTProf)
こんにちは。detです。

少々複雑なPerlプログラムを組んだ時に、
どの部分にどれだけ時間がかかっているか、
知りたいことがあると思います。

そんな時は、プログラムの各処理ごとに時間や実行回数などを計測し、
出力してくれるプロファイラを使うと便利です。

いくつか種類がありますが、今回はソースコードを見ながら、
html形式のGUIでプロファイル結果が確認できるDevel::NYTProfという
Perlプロファイラを簡単にご紹介いたします。

まず、CPANから、Devel-NYTProf-4.06.tar.gz をダウンロードします。
http://search.cpan.org/~timb/Devel-NYTProf-4.06/lib/Devel/NYTProf.pm

次は、READMEの指示通りに、

$ perl Makefile.PL
$ make
$ make test
$ make install


として、インストールします。

インストールできれば、あとは簡単です。
hoge.plのプロファイリングをしたい時は、以下のように実行します。

$ perl -d:NYTProf hoge.pl
$ nytprofhtml

とすれば、nytprof/index.html ができるので、
そのhtmlをウェブブラウザで開くと、以下のようになります。



それぞれどのサブルーチンが何回呼ばれて、
どのくらい時間がかかっているかが表示されます。
さらにそのサブルーチンをクリックすると、以下のように、
ソースコードを見ながら確認ができます。



かなり見やすいので、うまく使えば、
かなり便利なのではないでしょうか?

ちなみに、NYTProfのNYTはアメリカの有名な新聞社、
NewYorkTimes社のことです。社内のエンジニアの方が、
このプロファイラを作成したとか。
新聞社からこういうものが産まれてくるのが面白いですね。

それでは、detでした。
| deda | バイオインフォマティクス | 09:38 | comments(0) | trackbacks(0) |
Microarray解析
tokunagaです。
今回は、遺伝子発現変動を網羅的に調べたいときに使われているmicroarrayのデータ解析について流れを簡単にご紹介します。

生データ
 発光強度をスキャナーの画像から数値化する
   ↓
バックグラウンド補正
 出力された時点で補正されている場合もある
   ↓
正規化
 サンプル別のバラつきを揃える
   ↓
フィルタリングによる絞込み
 群間の有意差、p-value、FC、シグナル強度etc.の情報を基に絞り込む
   ↓   
クラスタリング、ヒートマップ作成
   ↓
アノーテーション付け
 遺伝子、GO etc.の情報を付与する

後半は目的によって変化してくると思います。
今はRのパッケージが色々とそろっているようですね!
| tokunaga | バイオインフォマティクス | 18:21 | comments(0) | trackbacks(0) |
perlに挑戦中
perl初心者も初心者のtokunagaです。
まだまだ勉強中です。

先ほどまでvcfファイルをフィルタリングするプログラムを作成していました。
vcfとはVariant Call Formatの略で、次世代シーケンサーのデータから検出された多型を記述する一般的な形式です。
そして無事作成し終わり、ちゃんと実行できたことに感動しました。
ただ、絞り込みすぎて最後まで残る多型の数が少ないので、絞り込み条件を引数で設定できるようにしようと思います・・・・

メタ文字も最初は、なんだろうこの文字列は・・・・と思っておりましたが、使ってみると便利ですね。

近々今度はRにも挑戦予定です。
| tokunaga | バイオインフォマティクス | 12:14 | comments(0) | trackbacks(0) |
perlを他の言語から見つめてみる
こんにちは。detです。

今回は、他言語出身者(特にC言語出身者)がperlに触れるときに、
最初に躓きやすい点について簡単に記事にしてみました。

1.条件分岐において
Cでifの次にもう一度条件分岐するときは、else if と書きます。
perlも同じような制御構造をもつので、つい、else ifと書いてしまいますが、これは怒られます。else if ではなく elsif と書きます。
if(...){..}elsif(...){...}
こんな感じですね。

2.条件分岐において2
条件分岐の後の、実行文が一つだけだと、
つい以下のように書いてしまいます。
if($i<2) $i++; 
これ、perlでは動きません。ちゃんと{}でくくってあげる必要があります。
if($i<2){ $i++; }

3.出力について
printf ではなく print です。
ついprintfと書いてしまいます。
perlにもフォーマット指定出力用にprintfがありますが、
通常は print を使います。

4.変数について
まず、変数の頭に、$を必ず付けること。
これを案外忘れることがあります。
また、myなどの宣言は変数の型を表すのではなく、
その変数の有効範囲を表します。

今日はこのくらいで。detでした。



| deda | バイオインフォマティクス | 13:09 | comments(0) | trackbacks(0) |
Perl
こんにちは。detです。

先週から、仕事で使用するスクリプトをPerlで書き続けています。
少しずつPerlの記法というものが、理解できて来ました。

なにより、Perlは自由な面が多く、これがメリットでもあり、
デメリットでもあると思います。

たとえば、以下の点などは、私のようなCしか経験のない人間には、
驚きでした。

1:型宣言が無くても変数が使える
型宣言をしなくても、文字でも数字でも何でも格納できてしまいます。
最初はどの変数が何に対応しているのかわかりにくくなるのでは、と思いましたが、これが慣れるとなかなか楽。(スコープの関係で、myはつけます)

2:if修飾子などの省略記法が色々ある
Perlは省略記法が優れていて、どんどんスクリプトを省略できます。
例えばif修飾子は、命令の後ろにつけることでその命令の条件分岐ができます。
my $i = 1;
print "hello!¥n" if $i == 1;


3:多彩な正規表現
なにより、Perlを特徴付けているのが、正規表現でしょう。
多彩な正規表現により、文字列処理がかなり容易に行えます。

ほかにも色々ありますが、今日はこのくらいで…。
| deda | バイオインフォマティクス | 09:52 | comments(0) | trackbacks(0) |
間違えていたのは最初でした・・・・
tokunagaです。

前回の記事で私は、Linux基礎の課題で多くのエラーが出て苦しんでおりましたが、
CentOSをインストールしなおしたところ拍子抜けするほどうまくいきました。
どうやら最初のセットアップの手順を間違えていたようで・・・・
色々なツールが入ってなかったようです。
それまでの苦労が嘘のようにあっさりと最後まで課題をやり遂げられました。
おかげでエラーには強くなれたような気はします。

皆さんもお気をつけください・・・
| tokunaga | バイオインフォマティクス | 18:00 | comments(0) | trackbacks(0) |
英語の論文を読む際に便利なツール
akbです。
ブログ更新の順番で、1週間の早さを実感しております。

さて、今日は英語の論文を読む際に便利なツールを紹介します。

1.準備するもの
下記の 銑をインストールします。
.屮薀Ε供Fire fox
▲▲疋ン:Greasmonkey
(https://addons.mozilla.org/ja/firefox/addon/greasemonkey/)
スクリプト:『Fast look up JP and EN』
(http://userscripts.org/scripts/show/12901)

2.インストール後の確認
ブラウザ画面の右上か右下にサル(monkey)のアイコンが表示されていればOKです。
サルのアイコンに色が付いてない場合は、アドオンが有効になっていないので、アイコンをクリックして"有効"にチェックをいれます。
サルのアイコンの隣にある下矢印を押すと『Fast look up JP and EN』にもチェックが入っていることが確認できます。
ここで、『ユーザースクリプトコマンド』を選択し、『Fast look up JP and EN-reset setting』の項目をクリックします。

3.起動確認
画面上で、『Alt+y』を押します。
すると、ポップアップ画面が表示されるので使いたい辞書にチェックをいれます。チェックは数字で表記されますので、使いたい順番にクリックしていきます。(画像を参照)


4.いざ実践
web上の英語ページにいきましょう。
検索したい単語の上にマウスのカーソルを持っていき、Altキーを押しながらダブルクリックします。
辞書が表示されれば成功です。
(画像はRNAという単語の上でダブルクリックをした結果です)


これで、論文講読のスピードも上がると思いますので、ぜひお試しください。
| akb | バイオインフォマティクス | 18:24 | comments(0) | trackbacks(0) |
遺伝子ビジネス
本日は、気になる記事を御紹介します。

遺伝子ビジネスの最前線

国家プロジェクトである「オーダーメイド医療実現化プロジェクト」が始り、早9年。
(弊社もそのうちの1つとなるでしょうが)次世代シーケンサーの普及と共に
様々なバイオベンチャーが旗揚げしているように思います。

バイオインフォマティシャンのニーズも今後ますます増えそうな予感です。
| | バイオインフォマティクス | 17:30 | comments(0) | trackbacks(0) |
創薬分野とバイオインフォマティクス
本日は、創薬・診断分野とバイオインフォマティクスの関係について記述します。
現在、創薬・診断分野に関する基礎研究は大きく進展してきています。ゲノム、エピゲノム、RNA、タンパク質などのさまざまな切り口や統合的なアプローチによって疾患のメカニズムや生命現象が明らかとなりつつあります。
しかし、これらの基礎研究が直ちに市場に出ない現状があります。その理由の一つには、創薬のプロセスが関係しています。


図に示すように、1つの医薬品が製品化されるまでには、9-17年程度の期間が必要です(開発費も500億ほど必要です)。また臨床試験の段階で、統計学的な調査に時間を要する事も創薬のプロセスを長引かせる原因となっています。
このような状況の中バイオインフォマティクスの分野は、プロセスの効率化に向けた技術として注目されています。例えば経済産業省では、生体分子の機能・構造・ネットワーク解析やそれら研究を強力に推進するためのバイオツールやバイオインフォマティクスの開発などに向けた取り組みを行ってきており、これらを活用することで、日本発の革新的な診断・治療技術を発信していくことを目標としています。
弊社でもそのお手伝いができればと日々邁進しています。
| | バイオインフォマティクス | 19:12 | comments(0) | trackbacks(0) |
アジアとバイオインフォマティクス
本日より、友人の結婚式に参加するため、シンガポール経由でマレーシアに行ってきます。初マレーシアです!とても興奮しております。
今日は、アジアとバイオインフォマティクスに関して記述したいと思います。
よろしくお願いいたします。

バイオインフォマティクスがゲノム解析において注目されていますが、アジアではシンガポールやインドがいち早く、バイオインフォマティクスの研究と実用化に乗り出しています。こちらのサイトは、バイオインフォマティクス事業を行っている世界の企業一覧です。インドは圧倒的に多く、シンガポールも代表社が挙がっています。
またシンガポールは上記の企業以外にも多くのバイオインフォマティクス事業を行っている企業や団体が存在します(参考HP)。

そのため、私のマレーシアにいる友人たちの多くが、バイオインフォマティクス関連の企業に勤めています。
日本ではまだ新しい事業ではありますが、世界ではメジャーな事業になってきているのではないかと感じています。
世界で勝負をするためには、多くの情報を得て、技術を高める必要があります。
私も努力していきます。
| | バイオインフォマティクス | 17:50 | comments(0) | trackbacks(0) |
GEN2PHEN
キムです。
最近、GEN2PHEN ConsortiumのHPをのぞく機会が増えてきました。
http://www.gen2phen.org/

GEN2PHENプロジェクトは、ヒトゲノム計画から怒涛の勢いで増大している‘parts list ’ ( the gene sequences ) と ‘toolkit ’ ( technologies )を、「Genotype -To -Phenotype Databases」作ることで、効果的に管理して利用しましょうというプロジェクトです。
NGS解析ツールの一覧
http://www.gen2phen.org/wiki/prediction-tools
は目的ごとに整理されています。
NGSデータ解析の後半で、必要になってくるツールが満載です。
| きむ | バイオインフォマティクス | 20:33 | comments(0) | trackbacks(0) |
第2世代シーケンサー
7/14(水)・15(木)と弊社のキムが、アプライドバイオシステムズの第2世代シーケンサーである「SOLiD」の講習会に参加させていただきました。

SOLiDの詳細

詳細は上記リンクや他の資料をご参照いただくとして、他の第2世代シーケンサーと比べてリード長が短い(50bp)反面、安価に大量にシーケンシングすることができるので、「De Novo」や「ターゲットシーケンシング」に適していることや、ChIP-Seqやメチル化などエピゲノム解析にも利用できるなど、幅広い場面で適用できるシステムであることが理解できました。

一次・二次解析(配列のアライメント等)は付属やフリーのソフトでできますが、その次の遺伝子ごとに変位をカウントするなど、独自の解析には弊社のようなバイオインフォマティクスの専門家が必要であることに確信を持ちました。

やはり、民間企業によるバイオインフォマティクスは益々 必要&重要になってきています。
| 社長 | バイオインフォマティクス | 14:05 | comments(0) | trackbacks(0) |
医療情報のデータマイニング
今日は、統計解析パッケージソフト「SPSS」の年に一度のセミナー「SPSS Data Mining Day 2010」が渋谷であったので参加してきました。

仕事の都合もあって、午前のセッションしか聞けなかったのですが、その午前のセッションの一つに「データマイニングの医療への応用」(武蔵野赤十字病院 消化器科 部長 黒崎雅之先生)という講演がありました。


内容は、まさに私たちが取り組んでいるエビデンスにもとづく個別化医療(Personalized Medicine、オーダーメイド医療)の実施がテーマで、
1.病気の情報
2.患者さんの情報
3.治療の方法
について、それぞれデータマイニングした結果を発表されていました。

具体的には、1.の病気の情報では、ある病気の進行度合いは年齢や性別で異なることが知られているが、生化学データやその他の臨床情報も加味したデータマイニングを行うことで、どういう場合に症状が進行しやすいのかエビデンスが得られている、という発表内容でした。

SNPタイピングによるGWAS解析で、C型肝炎ウイルス除去におけるインターフェロンの効果の度合いに、ある遺伝子が関わっていることが昨年発表されたので、その情報も含めたデータマイニングをされていました。

まだまだ国内の医療は、IT化やデータ解析によるエビデンス形成に遅れをとっている分野ですので、私たちもより一層がんばって社会に貢献していきたいと、気持ちを新たにしました。


【参考】
SPSS
Genome-wide association of IL28B with response to pegylated interferon-alpha and ribavirin therapy for chronic hepatitis C
| 社長 | バイオインフォマティクス | 15:46 | comments(0) | trackbacks(0) |
   1234
567891011
12131415161718
19202122232425
262728293031 
<< March 2017 >>

このページの先頭へ