アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
ASHG2016でgnomADが公開

10/18~22はアメリカ人類遺伝学会(ASHG)2016がバンクーバーで開催されていました。


その中でBroad Institueから発表されたGenome Aggregation Database (gnomAD)について紹介したいと思います。


次世代シーケンサー(NGS)を使っている人にはお馴染みの図ですが、NGSの技術の発展と普及によりDNAのシーケンスコストがどんどん下がっています。




 そのため、NGSを使った大規模なヒトゲノムの解析が世界中で行われています。例えばイギリスでは今まさに1万人分のゲノムを解析するUK10Kプロジェクトが行われていたり、日本では東北メディカルメガバンクから日本人2,049人の全ゲノム解析の結果がiJGVDで公開されていたりします。  最も大規模なものとしては、Exome Aggregation Consortium (ExAC)による60,706人のエクソームのデータベースが公開されており、今年の8月にNature誌で発表されました。


 今回発表されたgnomADは、ExACのversion 1からエクソームのデータが約2倍の126,216人分になり、 15,136 人の全ゲノム解析の結果も追加されています!




gnomADのブラウザでは、WESとWGSを合わせたSNPとINDELの頻度が閲覧可能ですが、今後CNVの解析もして閲覧できるようになるようです。まだβ版でフィルタリングの条件を検討している段階ということです。また、変異詳細のページではPhenotypeの情報を報告できるようになっており、情報が集まってくれば変異との関連付けもされていきそうです。
ちなみに、gnomADは"ノマド"と発音するみたいです。


gnomADでは、全部で40テラバイトの変異データをパースすることで作成しているそうで、解析にはScalaで実装されたhailというソフトウェアを使っています。このソフトについてはまた別の機会に紹介できればと思います。


 gnomADは、疾患関連変異の探索に強力なツールとして使えると思いますが、インハウスデータの頻度情報と比較して閲覧ができると便利ですね。


| onouek | 次世代シーケンサー解析 | 13:43 | comments(0) | - |
Somatic SNV検出編
2014年6月21日に開催した、アメリエフ株式会社・第33回バイオインフォマティクス勉強会の「フリーソフトではじめるがん体細胞変異解析入門」のスライドをSlideShareにて公開いたしました。
主に、ブログでもご紹介したことがあるソフトウェアSomaticSniperをつかった解析のご紹介です。昨年の日本癌学会学術総会でもポスターで発表しました。



なぜわざわざ体細胞変異解析に特化したソフトウェアを使用するのか、どんな点に注意して実験・解析を行うのかについて説明しています。
また少しですがCNV検出にも触れています。
みなさまの解析のご参考になりましたら幸いです。
| kubo | 次世代シーケンサー解析 | 16:31 | comments(0) | - |
chimerascanで融合遺伝子を検出する
融合遺伝子検出ソフトウェアは数多くありますが、ベストなソフトがないのが現状だと思います。
TopHat-FusionやdeFuse(deFuseの記事)などが有名ですが、今回はchimerascanというソフトの使い方を紹介したいと思います。

1. アノテーションデータの準備


・UCSCからヒトリファレンスゲノム(hg19)をダウンロードして、解凍します。染色体のFASTAファイルのみ(chr1.fa〜chrY.fa )を1ファイルにマージします。

・UCSCからヒトトランスクリプトームのデータをダウンロードして、解凍します。chimerascanのダウンロードリストからダウンロードすることもできます。


2. ソフトウェアのインストール


chimerascanからソースコードをダウンロードして解凍します。

・chimerasanのディレクトリにあるsetup.pyを使ってコンパイルします。

・chimerascan_index.pyでchimerascanのインデックスを作成します。
インデックスの作成にはbowtie-buildを使うのでパスを通しておく必要があります。

これで、ソフトウェアの準備ができました。

3. ソフトウェアの実行


chimerascanのディレクトリには、ヒト前立腺癌のRNA-Seqデータ(Fastqファイル)があり、融合遺伝子 TMPRSS2-ERG が含まれています。
これを用いた場合、以下のように実行できます。
$ python /usr/local/bin/chimerascan_run.py -v --quals solexa /home/genome/hg19/chimerascan/ /usr/local/src/chimerascan-0.4.5/tests/vcap_pe_53bp/TMPRSS2-ERG_1.fq /usr/local/src/chimerascan-0.4.5/tests/vcap_pe_53bp/TMPRSS2-ERG_2.fq output_dir


outputに生成されたchimerascan.bedpeからhtmlを作成することもできます。
$ python /usr/local/src/chimerascan-0.4.5/chimerascan/tools/chimerascan_html_table.py --read-throughs -o chimeras.html chimeras.bedpe


chimerascanは実行までの準備が比較的簡単だと思うので、ぜひ試してみてください。
| onouek | 次世代シーケンサー解析 | 16:15 | comments(2) | - |
fastaの折り返し位置を変える
東京はすっかり夏の日差しです。
夏好きとしては今年も夏が来るのは嬉しいのですが、無防備に紫外線を浴び続けたつけが着実に肌に出てきています。いい美白ケア法がないか気になる今日この頃です。

fastaフォーマットの配列行は一般的に80文字未満で折り返すのがお作法のようですが、折り返す長さを変えたい場合もあると思います。

FastX-toolkitのfasta_formatterを使うと、折り返す長さを変えることができます。

例えばこのようなfastaがあったとして、


以下のコマンドで、80文字で折り返したfastaを作ることができます(-wに折り返す文字数を指定します)。
$ fasta_formatter -i IN.fa -o OUT.fa -w 80


-w 0とすると折り返しなし(ID行と配列行が1行ずつの繰り返し)となり、スクリプト等で処理しやすいファイルが作れます。

お肌の曲がり角もコマンド一発で折り返しなしにできたらいいですね!
| hat | 次世代シーケンサー解析 | 15:32 | comments(0) | - |
fastqのIDの書式の話
ある公開されているexomeデータのfastqファイルをダウンロードして解析しようとしたところ、うまくいきませんでした。
最初は何が何だかわからず困っていたのですが、fastqファイルを確認するとID行の書式がよく見かけるものと違いました。

例として、最初の1リードはこうなっていました。
@HWI-ST7001130:252:C258HACXX:5:1101:1723:2158¥R1
ACCTGTTGCTTCCCCTGGGACTGAAGGCAGAAGTGACTCCCGAAATTCTCTTTCTGGACTCAGAAGGAAACCAAAGCAATCGATGACACCGCATCNNNNAC
+
?7?=D;BDD4C?AEEC;313AEECE?+<CC;1?DD????BD:?60BB<4)B8)=)8==AC)7=))=CC;=?;??>6.;((-(,5(>>>A='3&05####+(
このfastqの書式はBWAでは受け付けられず、エラーとなります(他のソフトでは未確認です)。

このID行(1行目)には以下の問題があります。
1.末尾のforward/reverse readを区別する部分に「¥」を用いている
 → 「/」または半角スペースならBWAが動きました。
2.リードのforward/reverseを「1」「2」ではなく「R1」「R2」と表記している
 → 「/」のときは「1」「2」、半角スペースのときは「1」「2」または「R1」「R2」のどちらでも動きました。

この場合は、「¥R1」「¥R2」を以下のコマンドで「/1」「/2」に書き換えて対処しました。
$ zcat [sample]_R1.fastq.gz | perl -nle 'if(/¥¥R1$/){s/¥¥R1/¥/1/};print' | gzip -c > [sample_new]_R1.fastq.gz
$ zcat [sample]_R2.fastq.gz | perl -nle 'if(/¥¥R2$/){s/¥¥R2/¥/2/};print' | gzip -c > [sample_new]_R2.fastq.gz
([sample]、[sample_new]は任意のサンプル名をご使用ください)

もし同じようなデータがありましたらご参考までに。
| kubo | 次世代シーケンサー解析 | 14:37 | comments(0) | - |
bamにread groupを追記する
GATKは、BAMのフォーマットに厳しく(参照ページ)、たとえばヘッダにサンプル名を含むread groupのリストがあり、かつすべてのリードがそのread groupに属しているBAMしか受け付けません。

Read group(以下RG)は、たとえばBWAではマッピングのときに -R で指定することができます。
$ bwa mem -M -R "@RG¥tID:sample¥tSM:sample¥tPL:Illumina" genome.fa sample_R1.fastq.gz sample_R2.fastq.gz

マッピングの時にRGを付けていなくても、Picard toolkitのAddOrReplaceReadGroupsというツールで追記することができます。

inputのsample.bamにRGを追記してoutput.bamを作成する場合の例を見てみましょう。
$ java -jar AddOrReplaceReadGroups.jar I=sample.bam O=output.bam RGID=sample RGLB=sample RGPL=Illumina RGPU=run_barcode RGSM=sample

AddOrReplaceReadGroups.jarの必須オプションは以下の通りです。

1. INPUT=File
2. OUTPUT=File
3. RGLB=String
4. RGPL=String
5. RGPU=String
6. RGSM=String

Picardに比べるとだいぶ手間と時間がかかりますが、samtools mergeで追記する方法もあります(参照ページ)。
ご参考までに、手順は以下の通りです。
続きを読む >>
| kubo | 次世代シーケンサー解析 | 15:33 | comments(0) | - |
SeqCap Epi連載[5]|BSMAP methratio.pyでメチル化検出
先日鈴鹿サーキットに初めて行ってきた久保(kubor)です。レーシングカーのスピードもさることながら、時の流れも早いもので、僕がアメリエフに入社して2ヶ月が経ちました。これからもどうぞ、よろしくお願いいたします。

さて、SeqCap Epi連載第5回目の今回は、BSMAP*1でマッピングした結果から、メチル化塩基を検出するPythonスクリプト(methratio.py)を紹介します。BSMAPは次世代シーケンサー(NGS)でバイサルファイトシーケンスを行った結果得られるシーケンスリードをマッピングするためのツールソフトウェアです。このBSMAPについては、前回「SeqCap Epi連載[4]|BSMAPでバイサルファイトシーケンスのマッピング」こちらの記事で紹介致しました。ご興味のある方は併せてお読みいただくと、より解析手順を把握しやすいと思います。
続きを読む >>
| kubor | 次世代シーケンサー解析 | 17:04 | comments(0) | - |
VCFのフォーマットに関するブログ記事について
先日、VCF formatの記事に一部誤りがあるとご指摘を受け、記述を訂正いたしました。
VCFのGTに関する解釈についての部分です。

VCFのフォーマットのGTやPLの解釈については、こちらのページ、GTとPLの関係についてはGATKのガイドを参照しています。

当ブログをご覧の皆様にはご迷惑をおかけしまして申し訳ありませんでした。
ブログ記事に関しまして、弊社は極力正しい記載を心がけて参りますが、万一誤記があった場合は、ご指摘いただけましたら幸いです。
| kubo | 次世代シーケンサー解析 | 14:00 | comments(0) | - |
SeqCap Epi連載[4]|BSMAPでバイサルファイトシーケンスのマッピング
「Trick or treat!!シーケンシングデータをくれないといたずらしちゃうぞ!!」
ということで記事執筆時の今日はハロウィンですが、ジャック・オ・ランタンはいままで作ったことがない久保です(kubor)。新しい解析方法を試してみたいけど、最適な公開データが見つからない。そういう時に誰かがこの一言で、データを見つけ出してくださると嬉しいですよね。もちろん、お菓子をいただけるというのであれば、それもありがたく頂戴いたします。

さて、第4回目を迎えたSeqCap Epi連載ですが、今回はこのパイプラインの要とも言えるマッピングソフト(BSMAP)の紹介です。

|バイサルファイト変換を考慮したマッピングソフト

「BSMAP(Bisulfite Sequence Mapping Program)*1」は、バイサルファイトシーケンスで得られたショートリードをマッピングするためのソフトウェアです。バイサルファイト変換を行うと非メチル化シトシン(C)がウラシル(U)に変換されます。このUはシーケンサーではチミン(T)として検出されますので、バイサルファイト変換によりC -> Uと置換され、その結果シーケンシングデータではC -> Tの置換が起きます。すなわち、バイサルファイト変換が100%の確率で成功している場合、残ったCはすべてメチル化されたCであると言えます。

シーケンシングリードをマッピングするときには、このC -> Tの置換、および逆鎖におけるG -> Aの置換を考慮することが肝要です。BSMAPのようなバイサルファイトシーケンスに対応したマッピングソフトは、これらを考慮することが可能です。
続きを読む >>
| kubor | 次世代シーケンサー解析 | 14:50 | comments(0) | - |
SeqCap Epi連載[3]|Trimmomticでシーケンシング用アダプターを除去
学生時代ですと今頃の時期は、実験圃場のイネの刈り取り時期を気にしている頃ですが、今年は来月の学会の準備ばかりを気にしているバイオインフォマティクス事業部の久保(kubor)です。

さて、3回目のSeqCap Epi連載ですが、今回からはまだアメリエフブログで紹介したことのないソフトを紹介します。
まず今日紹介するのは「Trimmomatic*1」です。

このソフトはSeqCap Epi解析パイプラインにおいてかなり初めのステップで使用するツールです。では、どのようなソフトなのか紹介いたします。
続きを読む >>
| kubor | 次世代シーケンサー解析 | 14:11 | comments(2) | - |
SeqCap Epi連載[2]|解析パイプラインの概要
前回に引き続き、SeqCap Epi連載の第2回です。
今回は、SeqCap Epiによる実験データをどのように解析していくのか、ご紹介いたします。用いるソフトはすべてオープンソースのフリーウェアです。
このパイプラインでは、SeqCap Epi CpGiant Enrichment Kit(ロシュ・ダイアグノスティック)を使用して、ライブラリー調整を行い、その後、NGSによるシーケンスを実行して得られたFASTQファイルから解析を行います。

解析ワークフロー


解析のワークフローは次の通りです。

順序項目使用ソフト
1QC(クオリティコントロール)FastQC, Trimmomatic
2マッピングBSMAP
3フィルタリングbamtoolsなど
4基本データの算出picardなど
5メチル化解析BSMAP(methratio.py)
6メチル化程度の比較methylKit
7アノテーションの付加methylKit
8SNPを考慮したメチル化解析BisSNP
  1. FastQCを用いて、クオリティをチェックし、Trimmomatic*1による低クオリティリードのトリミングを行います。
  2. リファレンスゲノムに対して、リードをマッピングします。
  3. BSMAP*2でマッピングするとSAMファイルを得られますので、Picardを使用してBAMに変換します。その後、bamtoolsやBamUtilなどのソフトを使用して、適切な整形を行います。
  4. 次に、整形後のBAMファイルから、リード数、カバレージ、ハイクオリティリードの割合、ターゲット領域にマッピングされたリードの数などの基本的なデータを算出します。この際には、PicardやGATK toolkitを用います。
    また、この際に実験が上手く行えたかどうかを確認するため、バイサルファイト変換効率を算出します。
  5. BSMAPのmethratio.py*3というスクリプトを用いて、BAMファイル内のメチル化を検出します。
  6. methylKit*4を用いてサンプル間のCpGサイトにおけるメチル化程度を比較します。
  7. metylKitを使用して近傍の遺伝子の転写開始点との距離や、遺伝子名などのアノテーション情報を付加します。
  8. BisSNP*5によりSNPを検出します。
以上が全体的な流れになっています。

このパイプラインの特徴



最大の特徴は、メチル化解析にBSMAPとBisSNPの2つのソフトを使用している点にあります。本パイプラインにおけるこれらの違いは以下です。

BSMAP:すべてのシトシンのメチル化率を算出
BisSNP:すべてのSNPおよびCpGサイトのメチル化率を算出

すなわち、これらの使い分けをすることにより、BSMAPでは、より網羅的にメチル化率を検出し、BisSNPではSNPの可能性のあるシトシンを検出することが可能になっています。このように、多くのCpGサイトをターゲットとし、SNPも考慮することを実現したメチル化解析パイプラインです。

更に、アメリエフでは、この他にも機能を追加して、メチル化程度をIGVで可視化ができるよう、解析パイプラインをブラッシュアップ中です。どうぞご期待ください。



*1:usadellab - Trimmomatic: A flexible read trimming tool for Illumina NGS dataにて公開中。トリミングの際にパラメータを設定することで、指定したアダプタ配列のトリミングが可能。
*2:bsmap - Bisulfite Sequence Mapping Programにて公開中。バイサルファイトシーケンスデータ用のマッピングソフト。
*3:BSMAPに含まれているPythonスクリプト。BAMファイルからメチル化シトシンを検出する。
*4:methylkit - R package for DNA methylation analysisにて公開中。メチル化解析データをサンプル間で比較するためのRパッケージ。
*5:Bis-SNPにて公開中。GATK toolkitに基づくSNPの検出とメチル化シトシンの検出が可能。
| kubor | 次世代シーケンサー解析 | 18:12 | comments(0) | - |
SeqCap Epi連載[1]|NGSでメチル化解析
DNAのメチル化解析手法はこれまでマイクロアレイや、リアルタイムPCRを使用したものがありましたが、どれも解析可能なゲノム範囲に限りがありました。
しかし、NGS(次世代シーケンサー)を使用するSeqCap Epi CpGiant Enrichment Kit(ロシュ・ダイアグノスティック)なら、全ゲノム上の550万箇所以上におよぶCpGサイトのメチル化を検出することが可能です*1

例えば、代表的なメチル化解析ビーズアレイ(HumanMethylation450 DNA Analysis Kit)と比較した場合、こちらは45万箇所のCpGサイトに対応している*2ため、12倍以上のCpGサイトを解析することが可能です。

このように、網羅的なメチル化解析が可能ではありますが、その反面実験データが多くなり、解析が容易に行いにくいという見方もできます。そこで、弊社ではSeqCap Epiを用いた実験データの容易な解析を可能にするオリジナル解析パイプラインを開発中です。
この解析パイプラインは、ロシュ社が公開している解析パイプライン*3をベースに開発していることが大きな特徴です。ここへ更にアメリエフ独自のアレンジを加えていますので、ご期待ください。

次回は、解析パイプラインの概要を紹介いたします。

*1:SeqCap Epi CpGiant Enrichment Kit製品紹介ページ参照
*2:HumanMethylation450 DNA Analysis Kit製品紹介ページ参照
*3:「NimbleGen SeqCap Epi ターゲットエンリッチメントデータの評価方法」をベースに開発
| kubor | 次世代シーケンサー解析 | 13:42 | comments(0) | - |
samtools ゲノムのインデックスファイルの中身
次世代シーケンサ解析では、リファレンスゲノムやbamファイルなど、サイズの大きなファイルを扱う必要があります。
大きなファイルには扱う前にインデックス(目次)を作成することがしばしばあります。多くのソフトはインデックスファイルがないと動きません。
同じデータでも、ソフトごとに別のインデックスファイルが必要な場合もあります。

インデックスファイルには人が読めないバイナリデータも、プレーンテキストのものもありますが、Samtoolsが作成するゲノムファイルのインデックスファイルは、人にも読めるものです。
内容を知っていると便利なこともあるので、ご紹介いたします。
続きを読む >>
| kubo | 次世代シーケンサー解析 | 14:09 | comments(0) | - |
新しいタキシードはいかが?
RNA-Seqを何度も実行している弊社の某社員が、先日
「bowtie、tophat、cufflinksの名前が紳士関連なのに気づいた」
と言っていたので「気づくの遅っ」と思ったのですが、
私には二十歳を過ぎてからサザエさんの登場人物名が海産物関係である
ことに気付いた知り合いもいるので、意外と固有名詞の意味は考えない
ものなのかもしれないと思いました。

bowtie(マッピング)、
tophat(スプライシングを考慮したマッピング)、
cufflinks(発現比較)、
cummeRbund(cufflinksの結果を描画するRパッケージ)
などをまとめてtuxedo suiteと呼ぶそうです。

最近monocle(single cell用のRパッケージ)というのも
公開されました。
個人的には次は「ステッキ」だろうと予想していたので
外れました。

cufflinks 2.2.0以降から、tuxedoが新しくなったそうです。


cuffnorm、cuffquantという新しいコマンドを使うことにより
サンプルが多い場合の発現比較が行いやすくなるそうです。
single cell RNA-seq案件が増えることでサンプル数が多い解析も
多くなってくると思います。
早く検証を行ってお客様に情報提供できるようにしたいと思います。
| hat | 次世代シーケンサー解析 | 14:15 | comments(0) | - |
SomaticSniper (後)
SomaticSniperについての続きです。

SomaticSniperを実際に動かしてみます。
SomaticSniperはCentOSでも問題なく動作しますが、Ubuntuでの使用が推奨されています。

基本的な実行コマンドは下の通りになります。
必要なファイルは腫瘍サンプルと、ペアとなる正常サンプルのBAMファイルです。両方とも、リアライメントやduplicated readの削除など、補正されたものが推奨されています。
reference.faにはBAMファイルを作成する際に使用したゲノムを指定します。

$ bam-somaticsniper -f reference.fa tumor.bam normal.bam snp_output_file

オプションで、VCFやBED出力を指定することができます。
その他にも、検出するSNVのdepthやマッピング・ジェノタイピングのクオリティ、統計的なsomaticらしさ(somatic score)なども指定できますので、信頼性の低いSNVを除外する場合はこの時に指定できます。
論文で、設定を検討して信頼性が高い条件を挙げていますので、条件を厳しくして信頼性がきわめて高いSNVだけを検出したいときは、パラメータをデフォルトから変更して実行します。
実行後の出力ファイルに対しフィルタリングができる付属のスクリプトもあります。固まって検出されたSNVのフィルタリングができるなど、実行時に指定できるオプションよりフィルタリング条件が多いのでこちらも便利です。

注意点としましては、GATKの変異検出ツールUnifiedGenotyperの-Lオプションのように解析する領域を指定することができないので、
target sequence解析で使用する場合、SomaticSniperを実行する前のbamをintersectBedなどでターゲット領域に絞ってから解析する必要があります。
VCF形式で出力した場合は、解析後にターゲット領域に絞ることもできます。
もちろん前者の方がSomaticSniperの実行時間は短くなりますが、なんといってもSomaticSniperは早いので、あまりファイルサイズを増やしたくない事情がある場合は後者の方法を利用してもいいですね。

前の記事でも触れましたが、BED出力も注意が必要です。
UCSC定義のBEDのフォーマット(詳しくはBEDフォーマット完全解説の記事を読んでください)と異なり、SomaticSniperのBEDは以下のとおり4〜6列目がUCSCの定義と違っています。
1. #CHROM
2. START
3. STOP
4. REF/ALT
5. SOMATIC_SCORE
6. TUMOR_DEPTH
SomaticSniperのBEDを使うときは、気に留めたほうがいいでしょう。

まとめとして、実行が簡単な点や実行が速い点、VCFという扱いやすいフォーマットで出力できる点などが、便利なツールだと感じます。

簡単ですが、SomaticSniperの説明でした。
| kubo | 次世代シーケンサー解析 | 15:33 | comments(0) | - |
SomaticSniper (前)
以前、heshiさんがSomatic Mutationを検出するツールについてのブログ記事を書かれましたが、そのうちSomaticSniperについて少しご紹介します。
論文:(SomaticSniper: identification of somatic point mutations in whole genome sequencing data. Bioinformatics. 2012 Feb 1;28(3):311-7.)
ダウンロード
TumorサンプルとNormalサンプルのBAMを比較し、統計的に有意にsomaticであると判断したSNVを検出するソフトウェアです。
SNVのみに対応しています。

出力フォーマットはデフォルトではタブ区切り形式(少しsamtoolsのpileup形式に似ています)ですが、VCF出力も指定できますので、snpEffなどのVCFへのアノテーションを行うツールと併用したり、IGVで可視化することもできます。
VCFではLOHのフラグも立ち、自前でフィルタリングする時に扱いやすいです。
BEDファイルを使うツールを使いたいときは、BED出力も指定できます。ただし、汎用的なフォーマットではないので使用方法を選びます。

特徴的なのは、実行時間の短さです。数GBのデータを2つも読み込んでいるにも関わらず、あっという間に解析が終わります。

次の記事で、もう少し詳細にご紹介いたします。
| kubo | 次世代シーケンサー解析 | 15:09 | comments(0) | - |
insert sizeを求める
先日、構造多型を解析するツールについてご紹介しました。
(Paired-end/Split-read/Complex)
これらのツールの一部は、解析の際、シーケンスデータのinsert sizeの入力が必要です。
自ら実験して得たデータならinsert sizeはわかると思いますが、公共のデータベースの公開データはリード長は記載されていてもinsert sizeは不明な場合が多いです。

そこでinsert sizeを計算してくれるツールを調べていたところ、PicardsのCollectInsertSizeMetricsがありました。
BAM/SAMファイルからinsert sizeのヒストグラムを描画してくれるツールです。
様々なパラメータがありますが、必須のパラメータは3つだけです。

1. 入力するBAM/SAM
2. insert sizeを出力するファイル
3. ヒストグラムを出力するファイル

コマンドは以下のようになります(ファイル名・バージョンは任意です)。
$ java -jar [path]/picard-tools-1.75/CollectInsertSizeMetrics.jar INPUT=sample.bam OUTPUT=output.txt HISTOGRAM_FILE=hist.pdf

OUTPUTパラメータで指定したファイル(ここではoutput.txt)には、insert sizeの平均や中央値、最大値や標準偏差など、insert sizeに関する統計量が書かれています。

HISTOGRAM_FILEで指定したファイル(ここではhist.pdf)には、下図のようなヒストグラムが生成します。
insert size distribution

また、このヒストグラムを生成したデータでは、統計量はmedian insert sizeが200bp、mean insert sizeが213.9bpと計算されました。
ヒストグラムと統計量を見比べ、このデータのinsert sizeは200 bpだと判断しました。

詳細については公式のマニュアル(こちら)をご参照ください。
| kubo | 次世代シーケンサー解析 | 15:01 | comments(0) | - |
ヒト Exome データ解析受託キャンペーン
ヒト Exomeデータ解析の受託キャンペーンを開始いたします!

ご注文受付期間:2014.1.14〜3.31ご発注分まで

期間限定の特別価格でのご提供となりますので、Exomeデータ解析の外注をご検討されている方は、ぜひ一度弊社ホームページをご覧ください!

ご不明な点がございましたらお気軽にお問い合わせください。
たくさんのお問い合わせお待ちしております!
| akb | 次世代シーケンサー解析 | 14:46 | comments(0) | - |
Phased vs. Unphased
VCFファイルの一番右側には、以下のようにジェノタイプ情報が記されています。

----------------------------------------
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
----------------------------------------

このジェノタイプ情報が、パイプ(|)で区切られている場合とスラッシュ(/)で区切られている場合があります。例えば、1|0とか0/1のように。
今までこの違いを意識したことはなかったのですが、どうやらちゃんとしたルールがあるようです。

"|"で区切られているデータは、フェーズされている(Phasedな)データ
"/"で区切られているデータは、フェーズされていない(Unphasedな)データ
といいます。

詳細は1000人ゲノムプロジェクトのWebページで解説されています。
http://www.1000genomes.org/node/101
| heshi | 次世代シーケンサー解析 | 17:33 | comments(0) | - |
構造多型を解析するツール : paired-endとsplit-read mapping
これまで、ゲノムの構造多型を検出ツールのうち、split-readか、paired-endを用いているツールを紹介してきました。
どちらかの手法だけでなく、両方を合わせて使用するツールも存在します。
複数の検出方法を併用するメリットは、paired-end mappingによる感度の高さとsplit-read mappingによるbreakpointの精度の高さの両方が得られることの他に、complex SVを検出できることがあります。


SVMerge
複数のソフトをまとめて使って解析できるパイプラインです。
使用されているのは、BreakDancer/Pindel/SECluster/RetroSeq/RDXplorer/Exonerateで、それぞれpaired-end mapping、split-read mapping、対となるreadのクラスタ分布、可動性遺伝因子データベースへのマッピング、read depth、アライメントによってSVを検出し、その結果をマージしています。
検出後、さらにVelvetによるlocal de novoアセンブリでSVを確認し、偽陽性を減らしています。

DELLY
Deletion、tandem duplication、inversion、translocationを検出するツールです。それぞれのSV種が独立して検出されるため、結果のネストや重複の可能性があります。
論文中では、他のツールと比べて、depthが少ないデータにも強いと紹介されています。
ホームページを見ると、メンテナンスやBWAなどの依存しているツールのアップデートに合わせたアップデートがこまめに行われているようです。


3回にわたり、様々なSV検出ツールについてご紹介いたしましたが、現在も新しい検出ツールは登場し続けています。
もしかしたら、、ご紹介したものよりも精度も感度も抜群に優れた定番となる新たなツールが現れるかもしれません。
| kubo | 次世代シーケンサー解析 | 14:32 | comments(0) | - |
Somatic Mutation検出ツール
Somatic Mutationを検出する主なツールを紹介します。
腫瘍サンプルと正常サンプルのBAM/Pileupファイル、リファレンスゲノムのFASTAファイルが入力ファイルとして必要です。

・VarScan2(2012)
入力ファイルは、BAMファイルをSAMtoolsでPileupファイルに変換する必要があります。引用論文数が多く、広く浸透しています。
・SomaticSniper(2012)
計算時間が抜群に短いです。SNV検出用なので、Indelは検出しません。
・Strelka(2012)
アリル頻度を考慮することによって、腫瘍純度が低いサンプルであっても、ある程度の検出感度を保つアルゴリズムです。
・MuTect(2013)
アリル頻度をベイジアンモデルのパラメータにすることで、アリル頻度が低い場合やリードのカバレッジが低い場合でも比較的感度が高いと主張しています。
・EBCall(2013)
Somatic Mutationと誤認識しやすいシーケンシングエラーを見分けるため、患者以外から得られた複数の正常なリファレンスサンプルを用います。
・Seurat(2013)
腫瘍純度が低いサンプルであっても、正常サンプルのカバレッジを高めれば、ある程度の検出感度を保てるアルゴリズムです。
・Cake(2013)
VarScan2とその他変異を検出するツールを組み合わせるアルゴリズムです。網羅的に検出した変異の中から、8種のフィルタリングによってSomatic Mutationを検出します。

偽陰性を減らすには複数のツールを組み合わせて使用するべき、という意見が多いようです。
| heshi | 次世代シーケンサー解析 | 14:00 | comments(0) | - |
構造多型を解析するツール : pair-end mapping
前回は、次世代シーケンスデータを用いて構造多型を検出するツールのうち、split-read mappingを用いているものを紹介しました。
今回の記事ではpair-end mappingを用いているツールをご紹介します。


MoDIL
10-50 bpのindelの検出ツールです。
ライブラリのInsert sizeの分布とマップされたreadのinsert sizeを比較し、その変化からindelを検出します。そのため、検出されたindelサイズは直接観測したものではなく間接的に推測した値になりますが、論文では実際の値と非常に近い値となったようです。
十分なdepthのあるshort readのシーケンスデータに向いているツールです。

BreakDancer
10-100 bpのindelを検出するMini(開発停止)と、それより大きなSVを検出するMaxの2パッケージがあります。Indelやinversion、translocationは検出できますが、tandem duplicationの検出には対応していません。
ちなみに、このBreakDancerの出力ファイルは、そのまま先週ご紹介したPindelに入力することができます。その他のツールでも検出結果を一定のフォーマットに変換すれば入力できるようです。

GASV
アルゴリズムを改良し、depthまで考慮するようになったGASVProもあります。
array-CGHデータの解析にも使えるので、アレイとシーケンサの結果を同一ツールで解析し、比較することができます。
今も時々メンテナンスされているようです。


Split-readとpaired-endを用いてSVを検出するツールについて紹介させていただきました。
ところで、この2つの方法のメリットとデメリットについてはそれぞれご説明しましたが(Split-readPaired-end)、欠点と長所があるなら補い合えたらいいと思いませんか?
次回はsplit-readとpaired-endを併せて使用してSVを検出するツールについてご紹介いたします。
| kubo | 次世代シーケンサー解析 | 15:17 | comments(0) | - |
構造多型を解析するツール : split-read mapping
次世代シーケンスデータを用いて構造多型を検出するツールは、これまでに数多く登場しています。
そのうちのいくつかを簡単にご紹介いたします。

今回はSplit-read mappingを用いているツールです。

Pindel
論文が発表されたのは2009年ですが、最近でも頻繁にメンテナンスされているようです。
およそ10-80 bpの小さなSVの検出に推奨されています。100bp以上のSVを検出するのが得意なツールと組み合わせて使うといいようです。

SLOPE
Target sequenceデータの解析に特化したツールです。また、Target sequencingの強みであるread depthが十分にあるデータが必要なようです。
Pair-end readの解析はsam他マップされたreadを入力しますが、single-end readの解析も、fastq形式なら行えるそうです。

CREST
Split-read mappingではマップされなかったreadを解析に用いますが、CRESTは、soft-clipping read(完全にはマップされなかったread)のうち長75 bp以上のものを解析に使用するところが特徴です。
BWAやBowtie2など、local alignmentが行えるマッピングツールでないとsoft clipされないため、マッピングソフトがlocal alignmentに対応しているか、注意が必要です。


次回はpair-end mappingを用いてSVを検出するツールについてご紹介いたします。
| kubo | 次世代シーケンサー解析 | 16:24 | comments(0) | - |
構造多型の検出原理2
前回に引き続き、今回は次世代シーケンサーを用いて構造多型(SV)を検出するもう一つの方法、Paired-end mappingを用いた方法についてご紹介します。

マップされなかったreadを用いていたSplit-read mappingと異なり、Paired-end mappingでは、マップされたものの、位置や向きがリファレンスと一致しなかったread(anomalously mapped read)を利用します。

たとえば、サンプル配列にinsertionがあると、対となるreadは本来ゲノム配列にマッピングされる距離よりも近距離にマップされます。Deletionがあった配列では逆に離れてマッピングされます。
paired-end mapping INDEL

また、inversionならreadの方向が逆になりますし、translocationならまったく離れた場所にマッピングされることもあります。
paired-end mapping OTHERS

この方法の長所としては、染色体間のtranslocationなど、split-read mappingを用いる方法に比べ多様なSVの検出が得意であることがあります。
欠点としてはbreakpointの精確な検出ができないこと、single-end readの解析は検出原理上不可能であることがあります。
また、もともと対となるread同士の距離にばらつきがあるため、小さなindelは検出できないことがあります。


次世代シーケンサを用いた構造多型(SV)の検出原理の説明はここまでです。
次回からはそれぞれの方法について、どのような構造多型検出ツールがあるのかご紹介します。
| kubo | 次世代シーケンサー解析 | 12:42 | comments(0) | - |
構造多型の検出原理1
構造多型(SV)を検出するためには、主にSplit-read mappingを用いた方法とPaired-end mappingを用いた方法があります。

まずSplit-read mappingを用いる方法からご説明します。
SVが存在するサンプル配列から得られたreadは、リファレンス配列にマッピングしても、配列が一致しないためマッピングされません。Split-read mappingを用いる方法では、このマッピングされなかったread(unmapped read)を使用してSVを検出します。
SVのbreakpoint(境界点)上にあるreadは、一部はリファレンス配列と一致していますが、SVのせいで残りは一致せず、マッピングされません。そこで、いったんreadを分割(split)して、リファレンス配列と一致するところのみをマップします。
その後、残った部分をリファレンス配列にマップします。残った部分がどのようにマップされたかによって、SVの種類を判断します。

Indelの例
Split-read(Indel)

その他のSVの例
Split-read(others)


Split-read mappingを用いる方法では、readの分割点がbreakpointにあたるため、精度よくbreakpointを検出できます。
また、図では分割したreadをマップするためにリファレンスにマップされた対のread(anchor read)を参考にしていますが、これを必要としない原理ではsingle-end readの解析も行えます。
ただ、多様なSVを検出するのは不得意なツールが多いようです。

次回で、もう一つの検出方法であるPaired-end mappingを用いる方法をご説明いたします。
| kubo | 次世代シーケンサー解析 | 15:14 | comments(0) | - |
構造多型の検出
ゲノミクスで多様性をもたらす要因や疾患の原因として注目されることが多い変異はSNVやindelなどですが、その他にも注目されている変異のひとつに、構造多型(Structural Variation、SV)があります。
主なSVはinsertion(挿入)、deletion(欠失)、tandem duplication(直列重複)、inversion(逆位)、translocation(転座)です。
SV説明図
次世代シーケンサのデータ解析では、SNVやsmall indelとは変異の規模(数十bpからkb、あるいはMb単位)や機構も異なるこれらの変異を検出するためにはSNVなどとは違った方法が必要です。
大きく分けてSplit-read mappingによる検出とPair-end mappingによる検出の2種類がよく用いられています。
次はこの2つの方法について少しご説明しようと思います。
| kubo | 次世代シーケンサー解析 | 16:16 | comments(0) | - |
変異の絞り込み 【7】 罹患同胞を用いた絞り込み
変異の絞り込み 【1】論文紹介
変異の絞り込み 【2】変異検出
変異の絞り込み 【3】候補の絞り込み方
変異の絞り込み 【4】公開データベースを用いた候補の絞り込み
変異の絞り込み 【4.5】お詫びと訂正
変異の絞り込み 【5】変異のクオリティとインパクト
変異の絞り込み 【6】遺伝型による絞り込み



前回までは父、母、子(IV-1)のトリオのexome解析により、疾患関連変異候補を24個にまで絞り込みました。
今回は、子の同胞(兄弟)のexomeシーケンスデータを用いて、より絞り込みます。
パキスタン家系図

まず、同胞(IV-2、IV-3)のexomeシーケンスデータから変異を検出します。
次に、前回までで絞り込んだ24個の変異が存在する24遺伝子上にある3人の子どもの変異を抜き出します。
さらに、再びQmergeVCFを使い、3人の変異情報を統合します。トリオの解析ではLCA罹患者である子(IV-1)をcase、非罹患者である両親をcontrolと設定しましたが、今回は罹患者IV-1、IV-3をcase、非罹患者であるIV-2をcontrolと設定しました。
その結果、疾患関連候補の24遺伝子上に、case2人に共通する変異133個が検出されました。
このあとは、トリオとほぼ同じように、公開データや変異のクオリティによって絞り込みました。
最後の、遺伝型による絞り込みにおいては、LCA罹患者であるIV-1、IV-3はホモ、非罹患者のIV-2はヘテロまたは野生型の変異を抜き出しました。

なんと、疾患関連変異候補はたったの4個にまで絞り込まれました。
問題は、この4個にLCAを引き起こしている変異が含まれているかどうかです。

残った4つの疾患関連変異候補はこちらです。

4個の中に、NMNAT1の9番目のバリンをメチオニンに置換する変異(上から2つめ)が含まれており、元々このexomeシーケンシングデータを公開したFalkらの報告と一致しました。
この結果は前回の記事の総括を非常に強く支持するものになりました。


以上、全7回にわたり、Exomeシーケンスデータから検出した変異から、重大な疾患に関する変異の候補をどのように絞り込むかについて連載いたしました。
(一部自社製ソフトを使用していますが、方針としては汎用性があるものです)

本解析ではトリオだけでなく罹患同胞の情報も用いたことで、非常に精度の高い絞り込みが行えました。
他にも、今回用いているdbSNPや1000人ゲノムプロジェクトなどの公開データがより充実することや、人種や民族別のデータベースがあると、精度が上がると考えられます。

希少な遺伝性の疾患の原因探索や治療には、変異を検出した後のフィルタリングが重要ですが、今後データベースなどの発展や充実によって、遺伝性疾患の研究が進展していくといいですね。
| kubo | 次世代シーケンサー解析 | 14:27 | comments(0) | - |
変異の絞り込み 【6】 遺伝型による絞り込み
変異の絞り込み 【1】論文紹介
変異の絞り込み 【2】変異検出
変異の絞り込み 【3】候補の絞り込み方
変異の絞り込み 【4】公開データベースを用いた候補の絞り込み
変異の絞り込み 【4.5】お詫びと訂正
変異の絞り込み 【5】変異のクオリティとインパクト



前回の記事までで、変異の公開データベース、変異のクオリティ、さらにタンパク質に与える影響(インパクト)を用いて約100,000個の変異を920個に絞り込みました。
今回は最終段階である、遺伝型による絞り込みの話です。

パキスタン家系図-トリオ

この家系では、LCA罹患者の両親(III-4、III-5)はLCAを発症していないことから、LCA関連変異は劣性遺伝していると考えられます。そこで、920個のうち両親ではヘテロ、罹患者ではホモである変異を抜き出しました。
本解析で用いているQmergeVCFは遺伝型もアノテーションしていますので、引き続きawkコマンドで抜き出しました。
その結果、疾患関連変異候補をたったの24個にまで減らすことができました。

pakistan-filter

最終的に疾患関連変異候補に残った24個はこちらです。24遺伝子一覧
上から2つめの変異が第1回でご紹介した論文中でLCAの原因遺伝子と報告されているNMNAT1遺伝子上の変異です(NMNAT1 mutations cause Leber congenital amaurosis)。

5回に分けてご紹介してきましたが、総括すると
■ dbSNPや1000人ゲノムプロジェクトなど、既知の変異の公開データベース
■ 変異のクオリティ
■ 変異がタンパク質に与えると予測されるインパクト
■ 家系情報(トリオ)
などの情報は、Exome解析により得られた大量の変異から疾患関連候補を絞り込むのに有効であることが示されました。


……と、実は連載開始時点ではここで終わる予定だったのですが、さらにあと1回続きます。
ここまでの解析では父、母、子のトリオ(それぞれIII-5、III-4、IV-1)のExomeシーケンスデータを用いていますが、論文で用いられている子(IV-1)の兄弟(LCA罹患者と非罹患者が1人ずつ)のExomeシーケンスデータも公開されています。

パキスタン家系図


トリオに加え、兄弟のExomeシーケンスデータにより得られたより詳細な家系情報を用いて、さらに絞り込みを行いました。
結果を次回でお話しします。

変異の絞り込み 【7】罹患同胞を用いた絞り込み
| kubo | 次世代シーケンサー解析 | 14:46 | comments(0) | - |
変異の絞り込み 【5】 変異のクオリティとインパクト
変異の絞り込み 【1】論文紹介
変異の絞り込み 【2】変異検出
変異の絞り込み 【3】候補の絞り込み方
変異の絞り込み 【4】公開データベースを用いた候補の絞り込み
変異の絞り込み 【4.5】お詫びと訂正


さて、前回に引き続き変異を絞り込んでいきます。
第3回でご説明したとおり、次はGATKで設定した変異のクオリティと、snpEffでアノテーションされた変異が与えるインパクトによる絞り込みです。


本解析では、ExomeシーケンスデータからGATKを用いて変異を検出していますが、その変異検出の際の信頼性(正しい変異かどうか)は、変異のカバレージやクオリティによります。そこで、経験的に信頼性が高い条件を設定し、クオリティフィルタ(PASSなど)を付与しています。絞り込みでは、その条件に合致する変異を抜き出しています。

snpEffはよく使われるアノテーションソフトです。変異によりコドンが変化したことにより、遺伝子翻訳後のタンパク質がどう影響を受けるかについてアノテーションを行います。変異してもアミノ酸配列が変わらない同義置換やアミノ酸の一時配列が変化する非同義置換や停止コドンであるかどうか、また、その変化によるインパクトの大きさについてもアノテーションを付与します。
本解析では"HIGH"または"MODERATE"とされたものを抜き出しました。


これらの情報も、QmergeVCFでアノテーションされていますので、具体的なやり方は前回同様awkコマンドです。
クオリティで絞り込んだ結果、11,561→7,793個、さらにsnpEffで絞り込んだ結果、7,793→920個になりました。



変異のクオリティとインパクトに基づく絞り込みはここまでです。
病気の罹患者から検出された変異の総数は100,000個以上ありましたが、現在の段階で疾患関連変異だと考えられる変異は現段階で920個。かなり絞り込めてきました。

疾患関連候補遺伝子の絞り込みは次回で最後です。
次回は、家系情報にもとづく絞り込みを行い、最終的な候補変異を決定します。

変異の絞り込み 【6】遺伝型による絞り込み
変異の絞り込み 【7】罹患同胞を用いた絞り込み
| kubo | 次世代シーケンサー解析 | 13:55 | comments(0) | - |
変異の絞り込み 【4.5】 訂正とお詫び
変異の絞り込み 【1】論文紹介
変異の絞り込み 【2】変異検出
変異の絞り込み 【3】候補の絞り込み方
変異の絞り込み 【4】公開データベースを用いた候補の絞り込み

第4回の記事の訂正です。

公開データベースを用いた絞り込みの例として、QmergeVCFでdbSNPと1000人ゲノムプロジェクトのデータをアノテーションし、それをもとに絞り込みを行いました。
ところが、本解析でアノテーションに用いていたdbSNPのデータが、間違っていたことが判明いたしました。
正しいデータを用いてアノテーション付しなおしたところ、変異の数は
106,488

(dbSNPに未登録の変異の抜き出し)

14,575

(1000人ゲノムプロジェクトでのアレル頻度が5%以下)

11,561
となりました。
(第4回の記事では106,488→58,980→20,714と報告していました)
訂正してお詫び申し上げます。

変異の絞り込み 【5】変異のクオリティとインパクト
変異の絞り込み 【6】遺伝型による絞り込み
変異の絞り込み 【7】罹患同胞を用いた絞り込み
| kubo | 次世代シーケンサー解析 | 14:48 | comments(0) | - |
変異の絞り込み 【4】 公開データベースを用いた候補の絞り込み
変異の絞り込み 【1】論文紹介
変異の絞り込み 【2】変異検出
変異の絞り込み 【3】候補の絞り込み方


前回は、疾患関連変異候補の絞り込みの概要についてご説明しました。
今回の記事では、その絞り込みを実際にはどのように行うのか、ご紹介します。

絞り込む前に、本解析では、まず、弊社製のツールQmergeVCFを用いて、3サンプルのデータを統合しました。
QmergeVCFについては、後日別の記事で改めてご説明したいと思っていますが、要約すると、複数のVCFファイルを一つのファイル(タブ区切り)に統合し、アノテーションも行うツールです。
QmergeVCFの出力ファイルは、Caseから検出された変異について、以下の情報がアノテーション付されたものです。
■遺伝子名、ポジション
■公開データベースによるアノテーション
■変異のインパクト
■(非同義置換なら)コドンおよびアミノ酸の変化
■Caseにおける変異のステータス(遺伝子型、depth、GATKを用いて付与した変異のクオリティ…前回の記事参照、など)
■Controlの、Caseと同じポジションにおけるステータス
※Caseには存在しない、Control固有の変異の情報は含まれません。
※それぞれの項目について、変異がデータベースに登録されていなかったり、コドンを変えない場合、Controlでは変異が起きていなかった場合は"."と出力されます。

この複数サンプルについての情報が一つに統合されたファイルを用いて絞り込みを行いました。
(ここでは絞り込みをLinuxで行っているためawkコマンドを使いますが、Windows/Mac上ならExcelのデータフィルター機能も便利です)


まず、公開データベースに登録されている変異の除外から行いました。
LCA患者の総変異(SNV + Indel)106,488個のうち、dbSNP132に登録されている変異を除去します。(※最新のdbSNPは138ですが、本解析ではdbSNP132までに登録されている変異の除去のみを行いました
Qmergeの出力結果(output)では、11列目にdbSNPのバージョンの情報が書かれています。
qmerge出力イメージ
そこで、11列目の情報が132以上のもの、または登録されていないものだけを抜き出して出力します(filter1)。

$ awk '{FS="¥t"} NR>1 && $11>=132 || $11 == "." {print}' output > filter1
awkコマンドでは1列目を$1、2列目を$2…と指定します。FSは入力ファイルの要素の区切りの指定で、タブ(¥t)のみを要素の区切りとみなすように指示しています。NRは行の指定で、ここでは1行目以外の行に対し処理を行うように指定しています。

dbSNP132で絞り込みをかけた結果、候補変異が106,488個から58,980個、およそ6割に減りました。

次に、1000人ゲノムプロジェクトでアレル頻度が高い変異を除去します。ここでは、アレル頻度が5%以下の変異、もしくは登録されていないもの(“.”)だけを抽出します。QmergeVCFの出力結果では、1000人ゲノムプロジェクトのアレル頻度は13行目にあります。

$ awk '{FS="¥t"} $13 <= "0.05" || $13 == "." {print}' filter1 > filter2

この結果、58,980個の変異候補を20,714個に絞り込むことができました。


公開データベースを用いた絞り込みはここまでです。公開データベースによる絞り込みで、疾患関連候補のSNV/Indelを

106,488 → 20,714

と2割にまで減らすことができました。

次回は、変異のクオリティやインパクトによる絞り込みで、さらに変異候補の数を減らしていきます。

変異の絞り込み 【4.5】お詫びと訂正
変異の絞り込み 【5】変異のクオリティとインパクト
変異の絞り込み 【6】遺伝型による絞り込み
| kubo | 次世代シーケンサー解析 | 15:26 | comments(0) | - |
アダプタ除去ソフトの比較
シーケンシングデータをFastQCなどでチェックしていると
アダプタ配列が混入しているのを見つけることがあります。

アダプタ除去ソフトウェアはいろいろありますが
今回は以下の3ソフトウェアの使い方をご紹介します。

cutadapt
FastX-Toolkit(fastxclipper)
tagcleaner


丁度以下のアダプタが混入しているらしき公開データがあったので、
このアダプタを除去することにしました。

◆Illumina PCR Primer Index 1
CTACAGTCCGACGATCTCGTATGCCGTCTT

◆実行コマンド
各ソフトを以下のように実行しました。

$ cutadapt -b TACAGTCCGACGATCTCGTATGCCGTCTTC -m 10 -n 1 ¥
original.fastq 1> out_cutadapt.fastq

$ fastx_clipper -C -l 10 ¥
-a TACAGTCCGACGATCTCGTATGCCGTCTTC ¥
-i original -M 3 -o out_fastxclipper.fastq

$ tagcleaner.pl -fastq original.fastq -out out_tagcleaner ¥
-tag5 TACAGTCCGACGATCTCGTATGCCGTCTTC -minlen 10


パラメータが違うので厳密に同じ条件にはできませんでしたが
だいたい同じような条件で実行しています。

◆結果
「処理後にアダプタ断片がどれくらい残っているか」と
「処理前後でリード数がどのくらい減少したか」を確認しました。

・FastX-Toolkitは厳しめ
(アダプタ断片がよく取り除かれているが、リード数がかなり減って
おり、アダプタ以外も誤除去している可能性あり)
・tagcleanerは緩め
(リード数はあまり減らないが、アダプタ断片が残っている)
・cutadaptは両者の中間
(リード数はあまり減らず、アダプタ断片も比較的取り除かれている)

データやパラメータにもよると思いますが、cutadaptが比較的
バランスがとれているように思われます。
| hat | 次世代シーケンサー解析 | 14:58 | comments(0) | - |
bwa のバージョン検討 その1
 最も広く利用されているマッピングソフトの一つにbwaがあります。bwaは2011年にリリースされたバージョン0.6が広く使用されてきましたが、今年の2月末にバージョン0.7がリリースされました。0.7ではBWA-backtrack、BWA-SW、BWA-MEMの、三種類のアルゴリズムが選択できます。この中でBWA-backtrackはilluminaのショートリード用、残り二つがショートリードから1Mbまでのロングリードに対応しているようです。また、bwaのサイトによれば、BWA-SWとMEMでは、MEMの方がより新しく、精度も高いとのことです。
 こうなると今までの0.6と、0.7[MEM]を比較したくなりますね。というわけでやってみました。条件は以下の通りです。

【サンプル】
WholeExomeSequenceデータ(SRR077486、paired-end)

【データ解析】
1. QC(FastQC、QCleaner)
2. マッピング(bwa0.6.1、bwa0.7.4[MEM])
3. リアライメント/リキャリブレーション(GATK1.6)
4. 変異のコール(GATK1.6)

上記の解析の流れは、弊社製Reseq解析パイプラインとほぼ同じです。

 では、コールされた変異について簡単に比較してみたいと思います。




 bwa(バージョン0.6と0.7[MEM])を用いて得られた変異がどの程度重なっているかを、ベン図にしてみました。SNVとIndelで分けています。このように、かなりの変異が重なっていますが、どちらか一方でしか検出されなかった変異も数%存在しています。

次回はこの数%の変異に焦点を当ててみたいと思います。
| deda | 次世代シーケンサー解析 | 16:57 | comments(0) | - |
変異の絞り込み 【3】 候補の絞り込み方
変異の絞り込み 【1】論文紹介
変異の絞り込み 【2】変異検出

前回は、LCA患者とその両親のトリオのfastqデータをダウンロードし、変異検出まで行いました。
今回は、主題である疾患関連変異候補の絞り込みについてお話します。


絞り込みの流れは、以下の通りです。

[1]既知の変異の除外
[2]アレル頻度が高い変異の除外
[3]クオリティが低い変異の除外
[4]タンパク質に与える影響が低い変異の除外
[5]遺伝形式による絞り込み


既知の変異(dbSNP132)は、SNPのデータベースに登録されている変異です。データベースは、重大な疾患のない人々から収録されているため、登録されている変異は重篤な症状をもたらす変異ではないと考えられます。疾患関連変異の探索では除外します。
アレル頻度が高い変異は、1000人ゲノムプロジェクトにおけるアレル頻度が5%より高い変異としました。これも、既知の変異と同様の理由で除外します。
クオリティが低い変異とは、マッピングの結果から変異を検出する際に、カバレージやクオリティなど、経験的に正しい変異が多く含まれている条件に該当しない変異のことです。GATKを用いてその情報を変異に付与して、信頼性の低い変異を除外しました。
タンパク質に与える影響が低い変異の除外には、snpEffを使用しました。snpEffはよく使われるSNV/SNPのアノテーションプログラムです。SNV/SNPがタンパク質に与える影響の大きさを予測し、“High”, “Moderate”, “Low”, “Modifier”のいずれかをアノテーションします。このうち、“Low”, “Modifier”とアノテーションされた変異は重篤な症状をもたらしにくいとして除外しました。
最後に、遺伝形式による絞り込みを行っています。罹患者の両親がいずれもLCAを発症していないことから、この家系のLCAの原因変異は劣性遺伝していると考えられます。そこで、両親においてはヘテロ、罹患者においてはホモの変異を抜き出しました。



ステップの紹介に終始してしまいましたが、今回はここまでにします。
次回は上記の[1][2]にあたる、公開データベースを用いた絞り込みの、実際の結果をお見せします。

変異の絞り込み 【4】公開データベースを用いた候補の絞り込み
変異の絞り込み 【4.5】訂正とお詫び
変異の絞り込み 【5】変異のクオリティとインパクト
変異の絞り込み 【6】遺伝型による絞り込み
| kubo | 次世代シーケンサー解析 | 18:16 | comments(0) | - |
変異の絞り込み 【2】 変異検出
変異の絞り込み 【1】論文紹介

前回は、exome sequencingによりある家系においてレーバー先天黒内障(LCA)を引き起こしている原因変異を特定した論文についてご紹介しました。
今回は、その論文で使用されたデータから、変異検出を行った結果についてご紹介いたします。


まず、論文で用いられたデータは公開されているので、DRAsearch に登録されているexome sequencingのfastqデータをダウンロードします。データはこちらです。
この家系は少し複雑な婚姻を繰り返しており、その最初の世代をI、次の世代をIIとし、罹患者の親世代が世代III、罹患者が世代IVにあたります。他のIV世代(罹患者の兄弟や従兄弟、非罹患者とLCA罹患者を含む)のデータも登録されていますが、今回の解析では罹患者とその父母のトリオのデータのみ使用することにします。トリオのうち、罹患者はIV-1(Accession numberはSRS344411)、母親はIII-4(同SRS344410)、父親はIII-5(同SRS344412)です。
Paired-end readなので、データは3サンプルについてそれぞれforward readとreverse readの2つずつ、合計6つダウンロードしました。


まず、各サンプルの変異を検出します。

【変異検出解析の流れ】
1.リードのクリーニング: QCleaner(弊社製ツール)
2.マッピング・カバレッジ集計: bwa、samtools、picard
3.リアライメント・SNV/Indel検出、クオリティフィルタの付与: GATK
4.snpEff: SNV情報の付与


リードクリーニングの結果はこちらです。
Pakistan-clean
III-4のforwardの塩基減少率が高く、あまりきれいなデータではないようです。
実際、III-4のforwardのFastQC結果を確認すると、こんな感じでした。
fastqc
(縦軸がリードのクオリティで、赤い部分はクオリティ20未満の塩基です)
今回の解析結果にはあまり影響ありませんでした。


マッピング結果、SNV/Indel検出結果はこちらです。
pakistan-map

ここから疾患関連遺伝子を絞り込むのは大変そうですね。


それでは、次回はいよいよ本題である、疾患関連変異候補の絞り込み手法についてご紹介します。

変異の絞り込み 【3】候補の絞り込み方
変異の絞り込み 【4】公開データベースを用いた候補の絞り込み
変異の絞り込み 【4.5】訂正とお詫び
変異の絞り込み 【5】変異のクオリティとインパクト
変異の絞り込み 【6】遺伝型による絞り込み
| kubo | 次世代シーケンサー解析 | 14:28 | comments(0) | - |
変異の絞り込み 【1】
今回は、次世代シーケンサーを用いて検出した変異の絞り込み解析の具体例についてご紹介したいと思います。


次世代シーケンサーを用いた解析の難点のひとつに、多数の変異が検出されるため、目的の疾患に関連した遺伝子の探索が困難であることが挙げられます。そのため、目的変異を絞り込む手法が重要になってきます。
今回は、公開データベース及び推測される遺伝形式を用いた絞り込みについてお話します。

こちらの論文は、次世代シーケンサーを用いてある家系のexome sequencingを行い、その家系における遺伝性疾患の関連変異を特定した、ひいては、他の家系でも、同じ遺伝子上に変異がある罹患者が存在していること、その遺伝子が疾患に関わっていることを突き止めたという論文です。

NMNAT1 mutations cause Leber congenital amaurosis
Nat Genet. 2012 September; 44(9): 1040–1045

簡単に内容をご紹介しますと、

・レーバー先天黒内障(LCA、幼児期に視力が減退する病気)罹患者が多く生まれる家系があった
・LCAの罹患者やその血縁者でexome sequencingを行った
・絞り込み:
 1. 罹患者ではホモだが、視力に異常がない両親や罹患者の兄弟ではホモではない変異
 2. 非同義置換が起きている変異(コドンが変化しアミノ酸に置換が生じている)
 3. 新規の変異(dbSNP132、1000人ゲノムプロジェクト、NHlBI ESPに未登録)
 4. マウスの網膜RNA-seqで発現量が高いことが知られている遺伝子上の変異
 5. SIFT、Polyphen他、アミノ酸置換がタンパク質に与える影響を調べるプログラムで、タンパク質の機能を損なうとされた変異
→その家系のLCA原因変異を特定

特定された変異は
・NMNAT1遺伝子上にある
・25番目のヌクレオチドがGからAに置換(アミノ酸一次配列の9番目のValがMetに置換)

この研究の公開データを用いて、変異検出および絞り込みを行い、論文で報告された疾患の原因変異がを探し出せるか試してみました。

次回から、その変異検出と、検出後の絞り込みを行う過程をご紹介しようと思います。

変異の絞り込み 【2】変異検出
変異の絞り込み 【3】候補の絞り込み方
変異の絞り込み 【4】公開データベースを用いた候補の絞り込み
変異の絞り込み 【4.5】訂正とお詫び
変異の絞り込み 【5】変異のクオリティとインパクト
変異の絞り込み 【6】遺伝型による絞り込み
| kubo | 次世代シーケンサー解析 | 14:58 | comments(0) | - |
Exome家系解析論文
Plos Oneに出ていたExome家系解析の論文をご紹介します。

Wu CC, Lin YH, Lu YC, Chen PJ, Yang WS, Hsu CJ, Chen PL.
Application of massively parallel sequencing to genetic diagnosis in multiplex families with idiopathic sensorineural hearing impairment.
PLoS One. 2013;8(2):e57369. doi:10.1371/journal.pone.0057369.

解析手順は以下の通りです。
1.難聴12家系について、既知の難聴遺伝子(80)およびその周囲をエキソンキャプチャ
2.HiSeq2000 100bp Paired-Endでシーケンシング
3.BWAでヒトゲノムhg19にマッピング
4.GATKでSNVとIndelを検出
5.VariantToolsで変異を絞り込み
6.ANNOVARでアノテーション
7.絞り込み
・1000Genomes ProjectやESP5400におけるアレル頻度
・SIFTやPolyPhen2のスコア
・変異箇所のアミノ酸の進化的保存度

最終的に、既知の変異2つと、新規の変異2つが検出できたそうです。

それら以外の変異も検出されたのですが、
「疾患群以外でもアレル頻度が高い」
「アミノ酸置換による影響が少ない」
「症状の重さやタイプ(低音域が聞こえづらい、など)が異なる」
などの理由で「難聴には無関係の変異」と判断しています。

難聴のように原因遺伝子が多数あり、症状も一律でない疾患の場合は、単に「患者に共通する変異」でなく「同じ症状の患者に共通する変異」を見る必要があるのだということが、発見でした。

当たり前のことなのですが、その疾患について知っているほうがより適切なフィルタリングができると思います。
シーケンシングデータだけを見ているとこの辺を見落としがちなので、お客様から伺えた範囲の疾患やサンプルの情報を活用し、有効な解析結果をご提供できるようにしていきたいです。
| hat | 次世代シーケンサー解析 | 16:53 | comments(0) | - |
ソフトウェア結果比較【BAMソート編・2】
ソフトウェア結果比較【BAMソート編・1】のつづきです。

前回、同じBAMをsamtools sortとSortSam.jar(picard)でソートしたところ、結果のBAMが異なりました。

異なる箇所を調べるため、まず結果BAMをsamtools view -hでSAMに変換しました。

【samtools sortでソート】
・結果(SAMに変換)
sort_samt.sam(666,626 bytes)

【SortSam.jarでソート】
・結果(SAMに変換)
sort_pica.sam(666,651 bytes)

これらのSAMをdiffコマンドで比較したところ、2点違いがありました。

違い1:ヘッダー
SortSam.jarでソートしたBAMには、1行目に
@HD VN:1.0 SO:coordinate
という行がありましたが、samtools sortでソートしたBAMにはありませんでした。
@HDタグは、SAMの仕様によると「ヘッダー行の先頭につけるもの」とあります。
BWAやsamtoolsはデフォルトでは@HDタグをつけませんが、picardを使うと自動的に@HDタグがつくようです。

違い2:同一ポジションにおける並び順
同じポジションに複数のリードがマッピングされている場合、それらのリードの順番が異なる場合がありました。
一例を挙げます。

今回のデータでは、
chr8 75275246
にSRR077861.62とSRR077861.104の2リードがマッピングされていて、sort_samt.samではこれらの並びが
SRR077861.62
SRR077861.104
の順になっていましたが、sort_pica.samでは逆になっていました。
他にこのようなところが2か所ありました。


今回のデータに関しては、これ以外の違いはありませんでした。
この程度の違いなら、どちらを使っても問題なさそうです。

その他
今回のBAMのソートのようにメモリを大量に使う処理は、デフォルトで/tmpに中間ファイルを大量に出力する場合が多いようです。

サイズの大きいBAM(46GB)でもソートを試してみたところ、SortSam.jarの場合は「-Djava.io.tmpdir=mytmpdir」をつけて中間ファイル出力先を変更しないとメモリ不足で落ちましたが、samtools sortだとデフォルトでも実行できました。

テストデータやその時のマシンの使用状況にもよると思いますが、samtools sortのほうがもしかすると若干「大きいファイルに強い」のかもしれません。

今後、他の処理についてもソフトウェアによる違いを検証してみたいと思います。
| hat | 次世代シーケンサー解析 | 14:00 | comments(0) | - |
「Reseq解析GUIマニュアル」を公開しました
本日は、弊社が slide share に新しく公開した「次世代シーケンス解析サーバー Reseq解析GUIマニュアル」をご紹介いたします。




この資料は、弊社で販売しております次世代シーケンスデータ解析サーバーのマニュアルの一部です。
Linux上であってもWindowsと同じような操作感覚で、グラフィカルに次世代シーケンスデータの解析が出来ることが示されています。

是非一度ご覧下さい。
| deda | 次世代シーケンサー解析 | 10:14 | comments(0) | trackbacks(0) |
snpEffデータベースの作成方法
snpEffとは、SNVやIndelなどの変異にアノテーション付けをしてくれるソフトの事です。

既に広く解析されている生物種の場合は、snpEff内にアノテーションのデータベースが揃っていますが、多少マニアックなものだとデータベースが無い場合もあります。その際には、新しいデータベースを自作することができ、その方法はsnpEffのHPにも説明があります。ですが、少々不親切なところもあり、少し苦労をします。

今日は日本語で、その方法について簡単に紹介いたします。

1.ゲノムのリファレンス配列(fasta)を用意します。ここでは仮にhogenomeという生物種のゲノムを追加する事にします。

2.hogenomeのGFF、もしくはGTFファイルを用意します。

3. snpEffのインストールディレクトリ(/path_to_snpEff/)にあるsnpEff.configに、ゲノムの情報を追加します。詳細はこちら。ちなみにリンク先に記載されているsnpEffect.configとは、snpEff.configのことです。

4./path_to_snpEff/data/ ディレクトリに移動し、新しいデータベース名のディレクトリを作ります。ここでの名前は、上記3番で追加したhogenome.genomeの.より前(ここではhogenome)と同じものにした方が良いようです。

5.hogenome/の中に移動し、GFFもしくはGTFファイルを置きます。そして、ファイル名をgenes.gff(GTFの時は .gtf)と変換してください。

6.hogenome/から出て、/path_to_snpEff/data/ 以下にgenomes/というディレクトリを作成します。その中に、hogenomeのfastaを置いてください。ファイル名はhogenome.fa としてください。

snpEffのHPではgffとfaに関して、gz圧縮された例しか載っていませんが、圧縮されていなくても動きます。

7./path_to_snpEff/に移動して、以下のコマンドを打ちます。

$ java -jar snpEff.jar build -gff3 -v hogenome

これで、エラーが出なければ、新しいsnpEffのデータベースの完成です。あとは楽しくアノテーションをしてみてください。
| deda | 次世代シーケンサー解析 | 18:50 | comments(0) | trackbacks(0) |
miRNA-Seq解析の論文
最近のPLoS Oneに載っていたmiRNA-Seqの論文をご紹介します。

Bansal A, et al.
Discovery and Validation of Barrett's Esophagus
MicroRNA Transcriptome by Next Generation Sequencing.
PLoS One. 2013;8(1):e54240.
Epub 2013 Jan 23. PubMed PMID: 23372692.


バレット食道(Barrett's esophagus)は
食道下部の粘膜が胃と同じ円柱上皮に置き換わった状態で、
食道がんの前駆病変の一つと考えられており、
胃食道逆流症(GERD)により引き起こされるそうです。

この論文では、バレット食道患者と胃食道逆流症患者の間で
以下の手順でmiRNAの発現を比較しています。

1. バレット食道患者 6名と胃食道逆流症患者 5名の
Total RNAのうち70bp未満のものをSOLiDでシーケンシング
(リード長35bp)

2. アダプタ除去後、15bp以上残ったものだけを、
miRBase最新版にbowtieでマッピング

3. precursor miRNA配列にマッピング(〜1ミスマッチ)、
mature miRNA配列にマッピング(オーバーラップ7bp〜)

4. miRNAにマップできなかったものを
fRNAdbのNon-coding RNA配列、hg19ゲノム配列、
大腸菌ゲノム配列にマッピング(〜3ミスマッチ)

5. miRNAとアノテートできたものについて、Bioconductorの
DESeqパッケージでリード数を正規化、発現レベルを比較

6. 解析結果をRT-PCRで検証

7. miRNAのターゲットを複数プログラム(*1)で予測し、
それぞれ上位5%だけに絞り、2つ以上のプログラムで予測された
遺伝子をターゲット候補として残す
*1:microT, miRanda, miRTarget, PicTar, PITA, RNA22, TargetScan

8. ターゲット候補遺伝子についてEGANを用いてKEGG Pathway解析

同じような発現傾向を示すいくつかのmiRNAが
同じ遺伝子をターゲットにしていることがわかったそうです。

miRNA-Seq解析の論文を定期的にチェックしているのですが、
Exome解析やRNA-Seq解析と異なり、
miRNA-Seq解析はまだ定番の仕組みができておらず、
論文ごとに手順やソフトウェアが異なるようです。

miRNA-Seq解析のご要望にお応えできるよう、
引き続き論文をチェックしていきたいと思います。

個人的な話ですが、アノテーションで使われているfRNAdb
昔3年間ほどお仕事をさせていただいていたことがあります。
自分が関わったデータが最新の論文でも活用されていることが
とても嬉しかったです。
| hat | 次世代シーケンサー解析 | 17:49 | comments(0) | trackbacks(0) |
次世代シーケンスデータ解析サーバーのHPをリニューアルしました
弊社では研究目的に合わせてお選びいただける、下記の2つの解析サーバーをご用意しております。

☆スタンダード
☆エンタープライズ

2つのサーバーはCPU、メモリ、ストレージ、冷却システムなどの性能に違いがありますが、最も特徴的な差は、サーバーによって、お選びいただける解析パイプラインの数が異なる点です。

下記の4つの解析パイプラインからスタンダードでは1つエンタープライズでは2つ、お選びいただくことができます。

☆Resequence/Exome解析パイプライン → 解析フローはこちら
☆RNA-seq解析パイプライン
☆ChIP-seq解析パイプライン
☆microRNA解析パイプライン

詳細は弊社ホームページでもご覧いただけます。

お客様の研究目的に合わせて解析パイプライン、サーバー構成をカスタマイズいたします。
次世代シーケンスデータの解析を行っている方は、ぜひ一度ご検討下さい。
| akb | 次世代シーケンサー解析 | 16:20 | comments(0) | trackbacks(0) |
NGSの「顔」
チャーノフの顔グラフは、多次元のデータを「顔」のパーツの大きさや配置にあてはめて視覚的に示す手法です。

群馬大・青木繁伸先生が公開されている、Rによるチャーノフの顔グラフ描画プログラムを使って、次世代シーケンサーの「顔」を描いてみました。

この論文のTable1から、Illumina MiSeq、Ion Torrent PGM、PacBio RS、Illumina GAIIx、Illumina HiSeq2000について以下の8つの値を入力しました。

・Instrument Cost(K$)
・Sequence yield per run(Gb)
・Sequencing cost per Gb($)
・Run Time(hour)
・Observed Raw Error Rate(%)
・Read length(bp)
・Paired reads(0:No,1:Yes)
・Insert size(bp)

顔グラフには変数が18個必要なのですが、8個しか用意できなかったため、残りの10変数はすべて0としました。

描画結果です!(クリックで大きくなります)


作為的なことはしていないのですが、さすがライバル、IonPGMとPacBioがにらみ合う構図になりました。脇でMiSeqが心配しています。

GAIIがニコニコしてかわいいです。Illumina3兄弟は共通点が多いと思っていたのですが、意外と似なかったです。

変数を増やしたら全然違う結果になると思います。ご興味のある方はぜひお試しください。
| hat | 次世代シーケンサー解析 | 12:20 | comments(0) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
 皆様、こんにちは。detです。今回は、前回に引き続きまして、日本人全ゲノムシーケンスデータの解析についてご紹介いたします。

 解析には、弊社製の Reseq パイプラインを用いました。前回は、データのクリーニング結果について、ご紹介しました。今回は、マッピングとSNV/Indel検出の結果を、簡単にご紹介いたします。
 まず、ヒトゲノム(UCSC hg19 + Scaffold )およびデコイ配列(version5)に対して、BWAを用いてアライメントを行いました。その結果、92.77%のリードがマッピングされました。
 次に、CCDS(consensus coding sequence)について集計を行った結果を以下の表に示します。



表より、CCDS の 97.80% を平均カバレージ28でシーケンスできていることが分かります。ゲノム全体としては、89.15% (2,765,146,564 base)を、10以上のカバレージでシーケンスすることができました。前回のクリーニング結果でもそうでしたが、この日本人全ゲノムのシーケンシングは高い精度で行われたことが分かります。
 また、GATKを用いてリアライメント・SNV/Indel検出を行った結果を下の図に示します。



 SNV/Indel検出の結果、3,413,445個のSNVと311,413個のIndelを検出しました。そのうち、1000人ゲノムプロジェクト及びdbSNP135に登録されていない多型は375,656個ありました。さらに、BreakDancerを用いて、ペアリードの情報から構造多型を検出した結果、361個の逆位、800個の染色体内転座、215個の染色体間転座などの候補が検出されました。

 次回は、1000人ゲノムプロジェクトJPTサンプルの解析結果と、日本人全ゲノムシーケンスデータを簡単に比較してみたいと思います。

-----関連記事-----
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
| deda | 次世代シーケンサー解析 | 18:17 | comments(2) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
 皆様、こんにちは。detです。今回は、前回までの1000人ゲノムプロジェクトJPTデータの解析からは少し離れまして、日本人の全ゲノムシーケンスデータに関する解析についてご紹介いたします。

 2012年7月31日に慶応大の富田先生の全ゲノムが 日本DNAデータバンク(DDBJ) で公開されました。シーケンスは、Beijing Genomics Institute(BGI) によって Illumina HiSeq 2000を用いて行われ、合計1,079,459,974リードが得られています。今回は、この日本人全ゲノムデータを弊社製の Reseq パイプラインで解析した結果をご紹介いたします。解析の流れは、本連載の前の記事「1000 人ゲノムプロジェクトJPT データの活用◆でご紹介したものと同じです。それでは、データのクリーニング(QC)結果を見ていきたいと思います。



 上にQC結果を纏めました。このように、幾つかのステップを経て精度の低いリードを削除することで、以降の解析精度を保証することができます。データによっては、この段階でかなりの数のリードが削除されてしまうこともありますが、今回は 99.97 % のリード(1,079,147,434 リード)を残すことができました。この結果は実験の精度が高く、綺麗なリードであったことを示しています。このQC処理には弊社製のツールであるQCleanerを利用しています。
 次回は、マッピングや多型検出の結果についてご紹介いたします。

-----関連記事-----
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
| deda | 次世代シーケンサー解析 | 18:00 | comments(0) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
 皆様、こんにちは。detです。前回に引き続きまして、1000人ゲノムプロジェクトJPTデータの解析に関する記事を書かせていただきます。

 今回は、1000人ゲノムJPTサンプルのデータ解析結果から得られたアレル頻度分布と、他の人種におけるアレル頻度分布を比較した結果をご紹介いたします。比較には、1000 genome project が公開している下記のデータを利用しました。
1092サンプルの多型情報をまとめたVCFファイル
 このデータには人種毎のアレル頻度も記載されているため、今回の比較において有用であると判断しました。それでは、比較した結果を下の図に示します。
 
 濃い青色の線が、弊社で解析した1000人ゲノムJPTサンプルのアレル頻度分布であり、それ以外の線が公開VCFファイルの人種毎のアレル頻度分布になります。JPTとEast Asians(※)の分布が、他の人種に比べて比較的似ているのが分かります。つまり、欧米人とアジア人ではリファレンスゲノム上での多型の分布や、頻度が異なる可能性が高いことが分かります。このことは、日本人ゲノムのデータ解析をする際に、日本人のリファレンスゲノムに基づいたデータ解析が必要であることを示しているのかもしれません。

※East Asiansに含まれる人種:CHB(Han Chinese in Bejing), JPT(Japanese in Tokyo), CHS(Southern Han Chinese)

 次回からは、1000人ゲノムから少し離れて、日本人全ゲノムシーケンスデータ(慶応大冨田先生のデータ)の解析についてご紹介いたします。

-----関連記事-----
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
| deda | 次世代シーケンサー解析 | 13:40 | comments(0) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
 皆様、こんにちは。detです。前回に引き続きまして、1000人ゲノムプロジェクトJPTデータの解析に関する記事を書かせていただきます。

 前回は1000人ゲノムJPTサンプルのデータ解析から得られた、88人の多型情報をまとめたVCFファイルについてご紹介いたしました。本日の記事では、今回の解析で得られた多型情報と、IlluminaのSNPジェノタイピングチップ(Omni 2.5 chip)で読まれたJPTサンプルの多型情報を比較することで、SNVの検出精度を明らかにする方法についてご紹介します。この Omni 2.5 chip は1000人ゲノムプロジェクトで検出された多型情報を基に設計されているため、今回の評価に用いるchipとして有用であると考えられます。
 まず最初に、1000人ゲノムJPTサンプルと、Omni 2.5 chip の両方に含まれるサンプルを二つ選びます。次に、Omni 2.5 chipによるSNV検出結果と、弊社での reseqパイプラインを用いた1000人ゲノムJPTのSNV検出結果が一致する割合を解析します。この結果から解析の信頼性が議論できると考えられます。現在、さらなる精度の向上に向けて、ブラッシュアップを行っております。

 次回は人種間で解析結果を比較した結果をご紹介いたします。

-----関連記事-----
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
| deda | 次世代シーケンサー解析 | 18:44 | comments(0) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
 皆様、こんにちは。detです。前回に引き続きまして、1000人ゲノムプロジェクトJPTデータの解析に関する記事を書かせていただきます。

 前回は88人の1000ゲノムJPTサンプルに対して、Reseq パイプラインでデータ解析した結果から、リードのクリーニング結果とマッピング結果をご紹介しました。今回は、得られた88サンプルの多型情報を一つにまとめたVCFファイルについてご紹介いたします。(VCFとは多型情報のフォーマットの一つで、超高速シーケンサーのデータ解析に用いられています。詳細は過去のこちらの記事を参照ください。)

 VCFファイルは横長のため、部分毎に説明します。下の図は、VCFファイルの先頭部分になります。染色体(Chr)と位置情報(Pos)から始まり、ID、リファレンスの塩基、ALTアリル(ALT)の塩基、クオリティになります。



ALTが複数の場合は、カンマで区切って表示されます。またクオリティ情報はサンプル間の平均値になっています。次は、下の図に示すようにフィルタリング結果と、インフォメーションが続きます。



どのフィルターに引っ掛かったか、PASSしたかが分かります。インフォメーション行に関する説明は過去の記事をご参照ください。今回の解析ではインフォメーション行の最後に以下の図に示す、FreqALTallelsとNumbersの項目を追加しています。



FreqALTallelsは全サンプル中のALTの頻度、NumbersはALTをもっていたサンプル数を示します。最後に以下の図に示すように、各サンプル毎の多型情報が記されます。



フォーマットに続いて、サンプル毎に情報が記されます。GTで表現されるのがジェノタイプで、0/0がリファレンスのホモ、0/1がALTのヘテロ、1/1がALTのホモを表します。ALTが複数の場合は、1/2や、0/2などの表記になります。
 以上、今回の解析で得られたVCFファイルについてご紹介しました。次回は、結果の信頼性などを見ていきたいと思います。


-----関連記事-----
1000 人ゲノムプロジェクトJPT データの活用
1000 人ゲノムプロジェクトJPT データの活用
| deda | 次世代シーケンサー解析 | 11:40 | comments(0) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
皆様、こんにちは。detです。
akbさんに引き続きまして、1000人ゲノムプロジェクトJPTデータの解析に関する記事を書かせていただきます。
前回の記事の目的に従いまして、今回は、1000人ゲノムのデータベースに含まれる100程度の日本人サンプルから88サンプルを選び、解析対象としました。解析の流れは、以下のようになっています。

1.リードのクリーニング: QCleaner(弊社製ツール)
2.マッピング・カバレッジ集計: bwa、samtools、picard
3.リアライメント・SNV/Indel検出: GATK
4.アノテーション付け: QuickAnnotator(弊社製ツール)
5.snpEff: SNV情報の付与


QCleanerの詳細につきましては、こちらのSlideShareをぜひご覧ください。またQuickAnnotatorに関しましては、弊社HPに説明がございます。
またこれらの一連のデータ解析は、弊社製のReseq パイプラインを用いる事で簡単に実行することができます。

では、解析結果を紹介していきたいと思います。まずは、リードクリーニングの結果です。



最初の2行ではクリーニング前のファイルサイズとリード数が記してあります。サイズにはかなり幅があります。またそれ以降は、クリーニング後の結果です。かなり綺麗なリードから、クオリティの低いリードまで、様々なリードがJPTサンプルに含まれていることが分かります。次にマッピング結果について以下の表に示します。



SNV数の最大値は Whole Genome Sequence のデータ、最小値は Target reseq のデータの結果です。
次回は、VCFファイルの詳細についてみていきたいと思います。
| deda | 次世代シーケンサー解析 | 18:35 | comments(0) | trackbacks(0) |
1000 人ゲノムプロジェクトJPT データの活用
akbです。
「日本人類遺伝学会 第57回大会」に出展いたしました。
弊社のポスター発表、ならびにブースまで足を運んでくださった皆さま方に、この場をかりて厚く御礼申し上げます。

本日から「日本人類遺伝学会 第57回大会」で発表させていただいたポスター内容の連載を開始いたします。
タイトルは『超高速シーケンサーを用いた疾患関連遺伝子探索のデータ解析と今後の展望-1000 人ゲノムプロジェクトJPT データの活用-』です。
日本人のゲノムデータを重点的に収集して解析することは、個別医療の実現や臨床研究の促進に有用であると考えられます。
そこで弊社では、1000ゲノムプロジェクトのデータから全てのJPT(Japanese in Tokyo)サンプルを抽出し、Illuminaのシーケンサーを用いてペアエンドでシーケンスされた88サンプルにおけるアリル頻度の算出を行いました。
<JPT88サンプルの内訳>
Whole Genome Sequence(WGS): 29
Whole exome Sequence(WXS): 35
Target reseq: 24
さらに、慶應義塾大学 環境情報学部の冨田勝 教授の全ゲノムシーケンスデータとの比較も行いました。
今回の解析には、弊社で開発した超高速シーケンスデータ解析サーバを使用しております。

解析の詳細を本日より複数回に渡って連載していきますので、連載終了までご覧いただけましたら幸いです。
| akb | 次世代シーケンサー解析 | 13:06 | comments(0) | trackbacks(0) |
BWAでRNA-Seq解析
最近読んだRNA-Seqの論文をご紹介します。

Brooks MJ, Rajasimha HK, Roger JE, Swaroop A.
Next-generation sequencing facilitates quantitative analysis
of wild-type and Nrl(-/-) retinal transcriptomes.
Mol Vis. 2011;17:3034-54. Epub 2011 Nov 23.



マウス網膜のRNA-Seq解析を以下の2パターンで行い、結果を比較して、PCRで検証しています。

(1)
-マッピング:BWA
-アノテーション:ANOVA
-アノテーションデータ:UCSC GenomeBrowser mm9 refFlat.txt

(2)
-マッピング:Tophat
-アノテーション:Cufflinks
-アノテーションデータ:Ensembl

マッピング対象はどちらもゲノム mm9です。

おもしろいと思ったのが、BWAでもゲノムにマッピングしているところです。
リードデータを既知遺伝子配列にマッピングする場合にはBWAも使いますが、ゲノムにマッピングする場合はTophatを使うのがあたりまえのようになっています。
これは、BWAでは長いギャップを許容しないため、リードの途中にスプライシング箇所があるようなマッピングができないためです。

結果、当然ですが、(1)で検出された転写物数は(2)で検出された転写物数の半分くらいでした。
ただ、発現レベルの精度は(1)のほうが良かったそうです。
これは用いたアノテーションデータによるところもありそうで、(2)のEnsemblデータはアイソフォーム数が(1)のrefFlatの3倍もあるために、精度がぼやけてしまっているのではないかという考察がなされていました。

特定の遺伝子の発現レベルを正確に見るなら、(1)の選択肢もありうるわけです。
# その場合、ゲノム配列ではなく遺伝子配列にマッピングしたほうがもっといいような気はしますが...

「RNA-SeqならTophat」と固定概念にとらわれてしまいがちですが、場合に応じてその都度最適な解析手順を考える必要があると思いました。
用いるアノテーションデータも精査していかないといけないですね。
| hat | 次世代シーケンサー解析 | 18:44 | comments(0) | trackbacks(0) |
NCBI SRA Toolkitの使い方
NCBI SRAで公開されているデータをテストデータとして使うには、sraファイルをダウンロードした後、fastqファイルに変換する必要があります(注:fastqで公開されているものもあります)。

sra→fastq変換には、NCBIが提供しているSRA Toolkitに含まれるfastq-dumpを用います。

・実行コマンド
$ fastq-dump -A SRR290592 --split-files SRR290592.lite.sra
注:ペアエンドデータの場合、--split-filesオプションを付けることによりfastqファイル2ファイルにわけて出力します。

・出力fastq(Reverse側ファイルの先頭4行)
@SRR290592.1 FCB062TABXX:4:1101:1435:1945 length=90
CGTTCACGCATCAGCTTCACGGAGCCAGAGG ・・・(略)
+SRR290592.1 FCB062TABXX:4:1101:1435:1945 length=90
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH ・・・(略)


ここで問題があり、fastqの各リード1行目(ID行)の形式がIllumina CASAVA形式になっていないと動かないソフトウェアがたまにあります。
CASAVA形式のID表記は以下のようになります。

・CASAVAのID行形式(※CASAVA1.8の場合)
@ <instrument-name>:<run ID>:<flowcell ID>:<lane-number>:<tile-number>:<x-pos>: <y-pos> <read number>:<is filtered>:<control number>:<barcode sequence>


こんな時、fastq-dumpの--defline-seqオプションを使うことで、ID行をCASAVA形式風にしてfastqを出力することができます。

・fastq-dumpの--defline-seqで指定できる項目

$ac: accession【SRR290592】
$si: spot id【1】
$sn: spot name【FCB062TABXX:4:1101:1435:1945】
$sg: spot group【NAGNNNNN】
$sl: spot length in bases【180】
$ri: read number【2】
$rn: read name【(今回は空)】
$rl: read length in bases【90】  (【】内は例)


大体上記の項目でCASAVAとの対応が取れるのですが、該当する情報がないものもあります。
これらについては便宜的に適当な値を入れてしまいます。
CASAVAの<is filtered>は「フィルターをパスしていればN、していなければY」という項目なので'N'に、<control number>は「control bitsが1つも立っていない時は0、それ以外は偶数」という項目なので'18'にします。

また、--defline-qualオプションで各リード3行目の表示を変えることができます。
fastqの仕様では、3行目は「'+'だけ」でも「'+'とID」でもいいことになっていますので、'+'だけにすることによってファイルサイズを減らすことができます。
3行目にID行がないと動かないソフトウェアには私は今まで出会ったことがありませんが、もしかしたらあるかもしれませんので、一応ご注意ください。

--defline-seqオプションと--defline-qualオプションを付けてfastq-dumpを再実行します。

・実行コマンド
$ fastq-dump -A SRR290592 --split-files
--defline-seq '@$ac:$si:$sn $ri:N:18:$sg'
--defline-qual '+' SRR290592.lite.sra


・出力fastq(Reverse側ファイルの先頭4行)
@SRR290592:1:FCB062TABXX:4:1101:1435:1945 2:N:18:NAGNNNNN
CGTTCACGCATCAGCTTCACGGAGCCAGAGG ・・・(略)
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH ・・・(略)


IDがCASAVA風、3行目が'+'だけになりました。
あくまでCASAVA「風」で、意味的に全く同じではありませんが、fastqのID行の形式チェックにひっかかって動かないソフトウェアがある場合などに試していただければと思います。
| hat | 次世代シーケンサー解析 | 14:47 | comments(0) | trackbacks(0) |
BWA v0.5.xとv0.6.xの比較(2)
hatです。前回から随分開いてしまいましたが、BWA v0.5.xとv0.6.xの比較(1)のつづきです。

bwa v0.5.x系とv0.6.x系のマッピング結果の違いを、公開データを使って比較してみました。
それぞれhg19に対してマッピングして、samtools flagstatでマッピング結果を出しました。

■使用したデータ
ERR034601(日本人女性Exome)

■bwa v0.5.9のマッピング結果
118495810 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
117074076 + 0 mapped (98.80%:-nan%)
118495810 + 0 paired in sequencing
59247905 + 0 read1
59247905 + 0 read2
116292908 + 0 properly paired (98.14%:-nan%)
116571052 + 0 with itself and mate mapped
503024 + 0 singletons (0.42%:-nan%)
234772 + 0 with mate mapped to a different chr
217456 + 0 with mate mapped to a different chr (mapQ>=5)

■bwa v0.6.1のマッピング結果(減:青字、増:赤字
118495810 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
117071517 + 0 mapped (98.80%:-nan%)
118495810 + 0 paired in sequencing
59247905 + 0 read1
59247905 + 0 read2
116287794 + 0 properly paired (98.14%:-nan%)
116565920 + 0 with itself and mate mapped
505597 + 0 singletons (0.43%:-nan%)
234930 + 0 with mate mapped to a different chr
217494 + 0 with mate mapped to a different chr (mapQ>=5)

bwa v0.6.1では、若干マッピング率が下がりました。ペアエンドのマッピング数が減ってシングルトンのマッピング数が増えています。

しかし、個々のマッピング箇所を並べて見ていると、下図のような例もありました。

bwa v0.5.9ではペアエンドが別々の染色体にマッピングされていたものが、bwa v0.6.1ではペアエンドでマッピングできるようになった例です。ミスマッチは入ってしまっていますが、これくらいのミスマッチであればbwa v0.6.1の結果のほうを採用してもいい気がします。

bwa v0.5.x系とv0.6.x系の差はWebで調べてもまだあまり情報がないようです。微妙な違いなのでしょうか。なにか見つけましたら、またブログに書かせていただきます。
| hat | 次世代シーケンサー解析 | 14:33 | comments(0) | trackbacks(0) |
意外と知らない?FastQCの便利オプション
hatです。

次世代シーケンシングデータのQCツールとして、一番使われているのがFastQCだと思います。

FastQCは非常に高機能なツールなのですが、「ポジションごとのクオリティスコア分布図」で、ポジション10bp以降の領域がグループ化されてしまうのだけが残念だと思っていました。
100bpくらいのリードだと、リードの後ろのほうでどのくらいクオリティが低下するかディテールが見えづらくなってしまうのです。

ところが、今日FastQCのUsageを見ていて、FastQCに--nogroup オプションを付けて実行すると、グループ化を無効にできることを知りました。


このオプションは実は結構前(Version 0.9.1)から実装されていたようです。

今日の教訓:
使い慣れていると思っているツールこそ、Usageをよく読もう。


もしかしたら世間では常識かもしれないのですが、社内で最初にこれを発見したのが嬉しくて、ドヤ顔しています。
| hat | 次世代シーケンサー解析 | 17:55 | comments(0) | trackbacks(0) |
QCの道 その4
こんにちは。detです。
今日は前回のQCの道 その3の続きです。

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

・fastq_quality_trimmer
FASTQ ファイルの各リードの 3'側から、指定したクオリティ値(qv)未満の塩基を順番に削除していきます。指定したqvより高いqvの塩基が見つかったら、そこでストップします。そして削除の結果、リードの長さが指定した長さ未満になった時はそのリード自体が削除されます。5'側からの削除には対応していません。

【optionの説明】
-h: ヘルプを表示します。
-t N: qvの閾値を指定します。この閾値未満の塩基は削除されます。
-l N: リード長さの最小値を指定します。この最小値未満の長さのリードは削除されます。デフォルトでは0になっています。
-z: 出力をgzip形式で圧縮します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。
-v: 処理前後でのリード数などを出力してくれます。

【実行例1: 3'末端からqvが20未満の塩基削除】
$ fastq_quality_trimmer -t 20 -Q 33 -i test.fastq -o out.fastq

【実行例2: 3'末端からqvが20未満の塩基を削除し、長さが10塩基未満になったリードを削除】
$ fastq_quality_trimmer -t 20 -l 10 -Q 33 -i test.fastq -o out.fastq


・fastx_renamer
FASTA/FASTQ ファイルの、ID行とオプション行(1行目と3行目)を、リード毎に昇順にカウントされる数字かそのリードの塩基配列でリネームします。

【optionの説明】
-h: ヘルプを表示します。
-n TYPE: リネームタイプを指定します。COUNTにすると、昇順の数字に、SEQにすると、塩基配列でリネームされます。
-z: 出力をgzip形式で圧縮します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。

【実行例: 昇順の数字で置き換える場合】
$ fastx_renamer -n COUNT -i test.fastq -o out.fastq -Q33

【実行例: 塩基配列で置き換える場合】
$ fastx_renamer -n SEQ -i test.fastq -o out.fastq -Q33

それでは、また次回に。
| deda | 次世代シーケンサー解析 | 16:49 | comments(0) | trackbacks(0) |
QCの道 その3
こんにちは。detです。
今日は前回のQCの道 その2の続きです。

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

・fastx_trimmer
FASTA/FASTQ ファイルの各リードの先頭から、指定した塩基数を削除して出力することができます。全てのリードの塩基が同じ数だけ削除されます。

【optionの説明】
-h: ヘルプを表示します。
-f N: リードの先頭からN塩基の手前までを削除します。
-l N: リードの先頭からN塩基まで残して削除します。
-z: 出力をgzip形式で圧縮します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。

【実行例1: 先頭から4塩基削除】
$ fastx_trimmer -f 5 -i test.fastq -o out.fastq -Q33

【実行例2: 4から10塩基を抽出】
$ fastx_trimmer -f 4 -l 10 -i test.fastq -o out.fastq -Q33


・fastq_quality_filter
FASTA/FASTQ ファイルの各リードのクオリティをチェックし、一定のクオリティ以上の塩基が指定した割合未満であった時に、そのリードそのものを削除するツールです。

【optionの説明】
-h: ヘルプを表示します。
-q N: 最低でも維持しなければならないクオリティ値(QV)を指定します。
-p N: リードを占める指定したQVを維持する塩基の割合が、ここで指定した数値(%)未満のときはそのリードを削除します。
-z: 出力をgzip形式で圧縮します。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。
-v: 処理前後でのリード数などを出力してくれます。

【実行例: QVが20以上の塩基の割合が20%未満のリードを削除】
$ fastq_quality_filter -i test.fastq -o qual_fil_1.fastq -q 20 -p 20 -Q 33


それでは、また次回に。
| deda | 次世代シーケンサー解析 | 17:51 | comments(0) | trackbacks(0) |
QCの道 その2
こんにちは。detです。
今日は前回のQCの道 その1の続きです。

まずは、FASTX-Toolkitが持つ各機能について順番に説明して行きたいと考えています。

1.インストール方法
上記のFASTX_Toolkitのウェブサイトからダウンロードページに飛ぶと、一番上に、コンパイル済みのバイナリデータがあります。
あなたの環境が64bitのlinux環境であれば、Linux(64bit)を選択しダウンロードします。この記事執筆時点でのファイル名は以下のようになっています。

fastx_toolkit_0.0.13_binaries_Linux_2.6_amd64.tar.bz2

これをlinux上で以下のコマンドで解凍します。

$ tar jxvf fastx_toolkit_0.0.13_binaries_Linux_2.6_amd64.tar.bz2

これで同じ階層に bin/ というディレクトリができます。
あとはその bin/ にPATHを通せばどこでも利用できます。


2.使用方法
bin/の中を見て頂くと、たくさんのツールが存在しています。それら一つ一つがfastqファイルに対してさまざまな処理をする機能を持っています。それらについて一つずつ、紹介していきます。
注意事項:FASTX-Toolkitの各種ツールを実行するときは、オプションで -Q 33 を与えることを忘れないでください。

・fastq_to_fasta
fastqファイルに含まれる各リードをfastaに変換します。単純に考えれば、fastqファイルの3行目と4行目を削除したものになります。また、Nが一つでもあるリードは除去されてなくなってしまいます。その代り、どんなにクオリティ値が低いリードでも除去されません。

【optionの説明】
-h: ヘルプを表示します。
-r: fastqファイルの各リード1行目のシーケンスIDをただの番号に変えます。ファイルの先頭から各リードに1から順番に番号が付きます。
-n: デフォルトではNが一つでもあるリードは除去されます。このオプションを入れると、Nがあっても除去されません。
-v: 処理前後でのリード数などを出力してくれます。
-i: 入力ファイルを指定します。
-o: 出力ファイルを指定します。指定しない場合は標準出力に出力されます。

【実行例1: 普通に変換】
$ fastq_to_fasta -i test.fastq -o test.fasta -Q 33

【実行例2: 番号付出力、及びNがあっても除去しない】
$ fastq_to_fasta -i test.fastq -o test.fasta -r -n -Q 33


これ以外にもいろいろツールがあります。それはまた次回にご説明いたします。
| deda | 次世代シーケンサー解析 | 18:03 | comments(0) | trackbacks(0) |
BWA v0.5.xとv0.6.xの比較(1)
BWAの最新バージョンはv0.6.1ですが(2012年6月現在)、まだv0.5.x系をお使いの方も多いのではないでしょうか。
BWA v0.5.x系とv0.6.x系で何が違うかまとめてみました。

■BWA v0.6.xの変更点
・4GB以上のゲノムに対応
・ForwardとReverseのインデックスを統合することにより、BWA-shortで20%、BWA-SWで90%のスピードアップを実現
・BWA-SWがPaired-endに対応

■インデックスファイルのちがい
ヒトリファレンスゲノムhg19.faに対して、bwa v0.5.9とbwa v0.6.1でそれぞれインデックスを作成すると、生成されるファイルサイズは次のように異なります。



BWA v0.6.1ではReverse用のインデックスが生成されなくなっていて、bwtやsaといったインデックスファイルのサイズが大きくなっています。

■デフォルトパラメータのちがい
・bwa indexの-c(build color-space index)オプションが削除されました。
・bwa alnは-c(input sequences are in the color space)オプションが削除されて、-Y(filter Casava-filtered sequences)オプションが追加されました。
・bwa sampeは(Usageを見る限り)変更がないようです。

■互換性
・v0.5.x系はv0.5.x系のインデックスファイルを使って、v0.6.x系はv0.6.x系のインデックスファイルを使って、マッピングしないといけません。
・bwa v0.5.9で作成したインデックスファイルに対してbwa v0.6.1でalnを実行しても、bwa v0.6.1で作成したインデックスファイルに対してbwa v0.5.9でalnを実行しても、エラーになります。

■iGenomes
ちなみに最近のIllumina iGenomesには、bwa v0.5.x系と0.6.x系両方のインデックスファイルが含まれています。

次回は、マッピング結果のちがいについてまとめます。
| hat | 次世代シーケンサー解析 | 17:39 | comments(0) | trackbacks(0) |
DNAメチル化解析 その1
akbです。
本日は、DNAメチル化解析について説明したいと思います。

DNAメチル化解析の実験段階では、よくバイサルファイトと呼ばれる処理が用いられます。
バイサルファイト処理とは、メチル化されていないシトシン(非メチル化C)をウラシル(U)に変換する処理のことです。
この時、メチル化シトシンはシトシンとして残留しますので、バイサルファイト処理前とバイサルファイト処理後を比較解析することで、
シトシンのメチル化を1塩基レベルで判定することができます。

バイサルファイト後にPCRでシーケンスリードを増やすため、実際にはC→T(C→U→T)とreverse鎖のG→Aという置換がおきます。
つまりバイサルファイト後のシーケンスリードをマッピングする際には、このC→T、G→Aの置換を考慮したマッピングが必要になるわけです。

下記の論文では、バイサルファイト後リードのマッピングに特化したmappingソフトとして3のツールを比較しています。

1)BISMARK
2)BSMAP
3)RMAPBS

マッピング効率、時間、使用するCPUコア数などからBISMARKを推奨しています。

いくつかの連載にわけて、BISMARKをご紹介していきたいと思います。

参考文献
Comparison of alignment software for genome-wide bisulphite sequence data.
Nucleic Acids Res. 2012 May 1;40(10):e79. Epub 2012 Feb 16.
http://www.ncbi.nlm.nih.gov/pubmed/22344695
| akb | 次世代シーケンサー解析 | 18:39 | comments(0) | trackbacks(0) |
QC の道 その1
こんにちは。detです。

今日は、次世代シーケンサーから得られる生データのクオリティコントロールについて、お話したいと思います。

次世代シーケンサーから得られるデータ(例えば、Fastq形式のファイル)は、クオリティに問題があることが多いため、その後の解析にそのデータをそのまま用いることはあまり推奨されません。

そこで、得られたデータを処理する前に、そのリードのクオリティを確認し、必要に応じてトリミングやリードの除外などの処理(クオリティコントロール(QC))をする必要があります。

ひとつ前の記事の一番下の資料にも一部書いてありますが、いくつかのフリーツールを用いてFastqファイルのQCをすることができます。よく利用されるものを以下に挙げてみます。

FASTX-Toolkit
Fastq/Fasta形式に対応したプリプロセスツール群です。いくつかのツールがセットになっており、データの統計解析や、形式の変換、長さやクオリティなどに基づいたトリミング・フィルタリング等の豊富な機能を持ちます。

PRINSEQ
PRINSEQもFasta/Fastqに対応したQCツールです。豊富なQC項目だけでなく統計解析を簡単にグラフィカルに確認できます。web版とコマンドラインで動かすstandalone版があり、windowsユーザーでも気軽に利用できます。

Tagcleaner
シーケンス用ライブラリ調整時に結合するアダプター(タグ)配列が、Fastqファイルの中に紛れ込んでいることがあります。そのタグを検出して削除してくれるのがこのTagcleanerです。タグ配列が不明の場合はその予測機能も持ちます。Tagcleanerにもweb版とstandalone版があります。

cmpfastq
ペアエンドのデータをQCしていくと、ペアの片側が除去されて無くなってしまい、片側だけのデータができてしまいます。その片側だけのデータを除去して、ペアを揃えるのがこのcmpfastqです。perlのハッシュを用いて処理を高速化していますが、大容量データだとメモリ消費量が上昇するのが難点です。

以上、簡単にQCソフトを紹介しました。上記以外にもQCに関するツールは存在しており、それぞれさまざまな特徴があります。

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

-------
次の記事 QC の道 その2 はこちら。
| deda | 次世代シーケンサー解析 | 17:12 | comments(3) | trackbacks(0) |
SAM format
tokunagaです。

第一回目はSAM formatの概要について
第二回目はSAM formatのヘッダー部分についてお話ししました。
本日はSAM formatのアライメント部分についての説明をしたいと思います。

アライメント部分の形式は


というように
「QNAME」から「QUAL」までの11項目の必須フィールドと、その後ろに使ったツールや条件等で異なってくるオプショナルフィールドがあります。

・必須フィールドの項目




・オプショナルフィールド



オプショナルフィールドは、使用したツール等でタグの名前が決められていることがあるので詳しくはマッピングに使用したツールのwebページ等をご参照ください。

*参考URL
http://samtools.sourceforge.net/SAM1.pdf
| tokunaga | 次世代シーケンサー解析 | 13:27 | comments(0) | trackbacks(0) |
SAM format
tokunagaです。

前回はSAM formatの概要についてお話ししました。
本日はSAM formatのヘッダーについて説明します。

ヘッダーには、シーケンスやリードグループなどのタグ情報が記述されており、以下のような形式で書かれています。
@<TAG>   <TYPE>:<VALUE> <TYPE>:<VALUE>

(例)
@SQ SN:gi|49175990|ref|NC_000913.2| LN:4639675

@マークの次にアルファベット2文字でタグが記され、
その後ろに「タイプ(アルファベット2文字):バリュー」のセットが並んでいます。

そして、各タグおよびタイプの簡単な説明を下記の表にまとめました。



このタグとタイプの内容がバリューの部分に記されています。
これらの情報でヘッダー行以下のアライメント部分の定義をしています。

次はアライメント部分について説明したいと思います。

*参考URL
http://samtools.sourceforge.net/SAM1.pdf
| tokunaga | 次世代シーケンサー解析 | 17:47 | comments(0) | trackbacks(0) |
SAM format
tokunagaです。
本日からはSAM formatについてです。

SAMとはSequence Alignment/Mapの略で、
次世代シーケンサー解析に用いられるformatの一つで、マッピングの結果が書かれています。
BAM formatは、このSAMファイルをバイナリ化したものです。

SAM formatの中身は
@SQ SN:gi|49175990|ref|NC_000913.2| LN:4639675
@PG ID:bwa PN:bwa VN:0.6.1-r104
SRR022885.1 16 gi|49175990|ref|NC_000913.2| 2966086 25 36M * 0 0 ACATGAAGCAACTGGCGACGTTGAATAATTGGTACG /+&I0I53>I+I.IIII@IIIIIII>IIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:3C0C31
SRR022885.2 4 * 0 0 * * 0 0 AGATTTTTTCCTGTAACGCTGCCAGTTGGTGGGCTC IIIIIIIIIIIIIIIIIIIIIII SRR022885.3 0 gi|49175990|ref|NC_000913.2| 1766166 25 36M * 0 0 AGCGTACGCCAAGTACGTGATCTGACGTTTTGCCCT IIIIIIIIIIIIIIIII2IIIIII7IE,III1(%(4 XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:0C33A1
SRR022885.4 16 gi|49175990|ref|NC_000913.2| 3182924 37 36M * 0 0 AATCAGGGAGTTCGGGGAAGATGTGGAGAAAAAAAG I796II+IIII+IIIIIIIIIIIIIIIIIIIIIIII XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:1 XO:i:0 XG:i:0 MD:Z:35T0
・・・・・・・・・・(続く)

という形式になっています。
基本的にタブで区切られています。

@から始まっている行はヘッダー部分で、ここではヘッダー部分以降のファイルの中身を定義しています。
それ以外の行はアライメント部分で、各リードのマッピング情報について書かれています。

次回はヘッダー部分の形式について書きたいと思います。
| tokunaga | 次世代シーケンサー解析 | 09:30 | comments(0) | trackbacks(0) |
fastq format
tokunagaです。

少し1回目と間が空いてしまいましたが、
本日はfastq formatの2回目です。
各リードの4行目に記されているqualityについて説明したいと思います。

・Quality value
3行目に使われているquality valueは以下のように計算されています。




・ASCIIコード
各塩基ごとのquality valueを一文字で表すためにASCIIコードが使われています。




・比較
S - Sanger Phred+33, raw reads typically (0, 40)
X - Solexa Solexa+64, raw reads typically (-5, 40)
I - Illumina 1.3+ Phred+64, raw reads typically (0, 40)
J - Illumina 1.5+ Phred+64, raw reads typically (3, 40)
L - Illumina 1.8+ Phred+33, raw reads typically (0, 41)

()の中はquality valueの最大値と最小値です。
現行のIllumina formatでは、実際に使われているquality valueの最大値は41です。

参考URL
http://en.wikipedia.org/wiki/FASTQ_format

| tokunaga | 次世代シーケンサー解析 | 09:32 | comments(0) | trackbacks(0) |
VCF format
どうもakbです。
VCF formatの連載三日目です。
今日は、FORMATの"PL"について説明しようと思います。

まず、前回の復習ですが、FORMATの"GT"はGenotypeの略で
0:reference 1:alternate(sample1)のように数字に意味を持たせており、
0/0,0/1,1/1はそれぞれREFのホモ、REFとALTのヘテロ、ALTのホモを示します。

まず図をご覧ください。
(例)

--------------------------------------------
REF = A ALT = G GT:PL 0/1:191,0,255
--------------------------------------------
"PL"に対応する値は"191,0,255"のようにカンマ区切りで示されています。
なぜ三つあるかというと"REF = A ALT = G"から考えられるGenotypeごとに phred-scaledの尤度計算をしているからです。

つまり、REF = A ALT = G の場合

AA AG GG
191 0 255

この数字は小さいほど信頼性が高くなりますので、
AG = 0が採用され、0/1、REFとALTのヘテロとなるわけです。

では、REF=T  ALT=A,GのようにALTの候補が2つある場合です。

--------------------------------------------
REF=T  ALT=A,G   GT:PL 1/1:178,82,76,104,0,98
--------------------------------------------
この場合、"REF=T  ALT=A,G"から考えられるGenotypeは

TT TA AA TG AG GG (←この順番になります)
178 82 76 104 0 98

AG = 0が採用され、REFとALTのヘテロ、つまりA/Gのように見えます。
しかし、ここは注意が必要です。
なぜならGTは1/1を示しており、ALTのホモを表しているからです。
この場合、ALT同士の組み合わせになっているため
"1/1"という表記になっていますが、ALTのヘテロです。

(2014年11月14日訂正)


三回にわたってVCF formatについて連載してきました。
今回、ご紹介できなかったID等は、下記URLをご参照下さい。
http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/
vcf-variant-call-format-version-41

---------------------------------------------------
(2014年11月14日追記)
【訂正】
上記のGTの説明につきまして、誤りがあるとご指摘をいただきました。
ありがとうございます。
当ブログをご覧のみなさまには、ご迷惑をおかけしまして申し訳ありません。

GTの数字の意味は、以下の通りになります。
0: REFアリル
1: ALTの1つめのアリル
2: ALTの2つめのアリル
:
(以下同)

そのため、GTが1/1ということは、ここではA/Aのホモということになります。

GATKのガイドによりますと、PLは、GTで示したジェノタイプの確率を1.0(PLでは0と示す)とし、それをもとに他のジェノタイプの尤度を相対的に算出しているため、上記の例のようなケース(GTとPLでジェノタイプが食い違うケース)は実際にはないと考えられます。

一例ですが、下記のような組み合わせがあると考えられます。

【ジェノタイプがAAの場合】
--------------------------------------------------------
REF=T  ALT=A,G   GT:PL 1/1:178,82,0,104,76,98
--------------------------------------------------------

【ジェノタイプがAGの場合】
--------------------------------------------------------
REF=T  ALT=A,G   GT:PL 1/2:178,82,76,104,0,98
--------------------------------------------------------
| akb | 次世代シーケンサー解析 | 12:49 | comments(0) | trackbacks(0) |
VCF format

どうもakbです。
VCF formatの連載二日目です。

今日は、INFOカラムのIDとFORMATのIDの説明をします。

まずは、INFOカラムの説明です。

(例)





(例)


FORMATの"GT"はGenotypeの略で
0:reference 1:alternate(sample1)のように数字に意味を持たせており、
0/0,0/1,1/1はそれぞれREFのホモ、REFとALTのヘテロ、ALTのホモを示します。
REF = A ALT = G  GT = 0/1

これは、GenotypeがAGのヘテロであることを表しています。

次回は、FORMATの"PL"について説明しようと思います。
| akb | 次世代シーケンサー解析 | 18:07 | comments(0) | trackbacks(0) |
VCF format

どうもakbです。
今日からVCF formatの連載が始まります。

さて、VCFとは (Variant Call Format)の略で、次世代シーケンサー解析に用いられるformatの一つです。
samtoolsを用いて抽出した多型情報が、VCFファイルに格納されて出力されます。

まず、VCFファイルの中身を見てみます。



ファイルのメタデータは、"##"の文字列の後に記載されています。
例えば、最初の行には、VCFのバージョンが次のように示されます。
(例)
##fileformat = VCFv4.1

次にVCFファイルの本体で使用されるフィルタ(INFO)とFORMATが、
メタデータセクションに含まれます。
(例)
##INFO = <ID=DP,Number=1,Type=Integer,Description="Total Depth">

"##"から始まる行(メタデータ)の後の構造は、"#"から始まるヘッダー行と、そのヘッダー行のカラムに対応したデータ行から成ります。


ヘッダ行には9個のカラム+SAMPLEカラムがあります。これらのカラムは次のとおりです。





次回のブログでは、各カラムに対応するデータ行の説明をしたいと思います。
| akb | 次世代シーケンサー解析 | 09:16 | comments(0) | trackbacks(0) |
融合遺伝子検出ソフトウェア defuseを使ってみる【第三回】
融合遺伝子検出ソフトウェア defuseを使ってみる【第二回】のつづきです。

第一回ではデータの準備とソフトウェアのインストール、第二回では設定ファイルの編集とデータセットの作成を行いました。

今回はいよいよ融合遺伝子の検出を行います。

【第一回】
1. 解析データの準備
2. アノテーションデータの準備
3. ソフトウェアのインストール

【第二回】
4. 設定ファイルの編集
5. データセットの作成

【第三回】
6. 融合遺伝子予測



6. 融合遺伝子予測

defuseのscriptsディレクトリにあるdefuse.plを実行します。
この際、以下のオプションを指定します。
-c config.txtファイルのパス
-d 解析データ(Fastqファイル)があるディレクトリのパス
-o 結果を出力するディレクトリのパス
-p 実行プロセス数(指定しなければ1)

◆実行例
$ perl /usr/local/defuse-0.5.0/scripts/defuse.pl ¥
-c /home/hat/test_defuse/config.txt ¥
-d /home/hat/test_defuse/read ¥
-o /home/hat/test_defuse/out -p 2

-oで指定したディレクトリに結果が出力されます。
たくさんのファイルが出力されますが、重要なのは以下の3ファイルになります。

・result.tsv
 予測された融合遺伝子です。1行が1融合を示します。

・result.classify.tsv
 result.tsvにprobability列を追加したものです。

・result.filtered.tsv
 config.txtで指定したprobability_thresholdで、result.tsvをフィルタリングしたものです。

result.filtered.tsvの例を見てみましょう。
defuse-results.filtered.tsv
各遺伝子の染色体番号・染色体上の座標・向き、融合の位置・ロケーション(coding, UTR, intron, 上流下流など)、probabilityなどが記載されています。
本当に融合遺伝子らしいものがどれか、各情報で絞り込みを行ってください。

出力結果の詳しい説明はdeFuse manualをご参照ください。

これは第一回でご紹介した、融合遺伝子 TMPRSS2-ERGが入っているデータに対してdeFuseを実行した結果なのですが、ちゃんとTMPRSS2遺伝子-ERG遺伝子間の融合が予測できています。


これでdeFuseの記事は終わりです。お疲れ様でした!
無事に実行できましたでしょうか?

融合遺伝子解析が癌研究のさらなる発展につながることを期待しています。

ご注意
・本記事ではヒトの解析例を示します。同等の情報を用意できれば他の生物種でも実行可能と思われます。
・データを置く場所・ソフトウェアのインストール先などは、お使いの環境にあわせて読み替えてください。
・defuseは現時点の最新版(0.5.0)を使用しました。
| hat | 次世代シーケンサー解析 | 10:21 | comments(0) | trackbacks(0) |
融合遺伝子検出ソフトウェア defuseを使ってみる【第二回】
融合遺伝子検出ソフトウェア・defuse【第一回】のつづきです。

前回は、データの準備とソフトウェアのインストールを行いました。
今回は、設定ファイルの編集とデータセットの作成を行います。

【第一回】
1. 解析データの準備
2. アノテーションデータの準備
3. ソフトウェアのインストール

【第二回】
4. 設定ファイルの編集
5. データセットの作成


【第三回】
6. 融合遺伝子予測


4. 設定ファイルの編集

defuseのscriptsディレクトリ以下にあるconfig.txtを適当な場所にコピーし、以下の項目を編集します。

・source_directory にdeFuseのソースコードのディレクトリを記載
source_directory = /usr/local/defuse-0.5.0/



・dataset_directory にデータセットを置くディレクトリを記載
前回ダウンロードしたアノテーションデータを元に、解析に必要なデータを作成しておく必要があります。これをデータセットと呼びます。
dataset_directory = /home/hat/test_defuse/dataset


・gene_models に遺伝子モデルGTFファイルのパスを、genome_fasta にリファレンスゲノムFastaファイルのパスを記載します。
gene_models = /home/hat/test_defuse/annotation/Homo_sapiens.GRCh37.62.gtf
genome_fasta = /home/hat/test_defuse/annotation/Homo_sapiens.GRCh37.62.dna.chromosome.fa


・repeats_filename にリピート情報ファイルのパスを記載します。
repeats_filename = /home/hat/test_defuse/annotation/rmsk_hg19


・est_fasta にESTのFastaファイルのパスを、est_alignments にESTのIntron情報ファイルのパスを記載します。
est_fasta = /home/hat/test_defuse/annotation/est.fa
est_alignments = /home/hat/test_defuse/annotation/intronEst.txt


・unigene_fasta にUniGene情報ファイルのパスを記載します。
unigene_fasta = /home/hat/test_defuse/annotation/Hs.seq.uniq


・各ソフトウェアのパスを記載します。
bowtie_bin = /usr/local/bin/bowtie
bowtie_build_bin = /usr/local/bin/bowtie-build
blat_bin = /usr/local/bin/blat
fatotwobit_bin = /usr/local/bin/faToTwoBit
r_bin = /usr/local/bin/R
rscript_bin = /usr/local/bin/Rscript


5. データセットの作成

defuseのscriptsディレクトリにあるcreate_reference_dataset.plを実行します。この際、-cオプションでconfig.txtファイルのパスを指定します。

config.txtのdataset_directoryで指定したディレクトリ内にデータセットが作成されます。この処理は少し時間がかかります。

これで準備は整いました。
次回はいよいよ融合遺伝子の検出を行います。

ご注意
・本記事ではヒトの解析例を示します。同等の情報を用意できれば他の生物種でも実行可能と思われます。
・データを置く場所・ソフトウェアのインストール先などは、お使いの環境にあわせて読み替えてください。
・defuseは現時点の最新版(0.5.0)を使用しました。

| hat | 次世代シーケンサー解析 | 18:28 | comments(0) | trackbacks(0) |
融合遺伝子検出ソフトウェア defuseを使ってみる【第一回】
前回ご紹介した融合遺伝子検出ソフトウェア・defuseを使ってみます。

一回目の今回は、データの準備とソフトウェアのインストールを行います。

【第一回】
1. 解析データの準備
2. アノテーションデータの準備
3. ソフトウェアのインストール


【第二回】
4. 設定ファイルの編集
5. データセットの作成

【第三回】
6. 融合遺伝子予測

1. 解析データの準備

融合遺伝子を調べたいサンプルのリードデータ(Fastqファイル)を用意します。

とりあえず動かしてみるだけなら、deFuseの論文[1]で使われているテストデータで試してみましょう。
これは元々融合遺伝子検出ソフトウェア FusionSeq[2]で使われたヒト前立腺癌のRNA-Seqデータで、融合遺伝子 TMPRSS2-ERG が含まれています。

2. アノテーションデータの準備

以下のデータをダウンロードして解凍します。

ヒトリファレンスゲノムhg19のFastaファイル(染色体がchr1〜21、MT、X、Y以外のFastaファイルは不要なので削除し、ヘッダーのdescriptionを削除して1ファイルにマージしておきます)

ヒト遺伝子モデルのGTFファイル

ヒトESTのFastaファイル

ヒトESTのIntron情報

ヒトUniGene情報

ヒトのリピート情報(UCSC Genome BrowserのRepeatMaskerトラックをダウンロード)


3. ソフトウェアのインストール

以下のソフトウェアをインストールします。

deFuse

bowtie ※0.12.5以上が必要

blat バイナリ版

faToTwoBit バイナリ版

R


これで、必要なデータとソフトウェアが準備できました。
次回は「設定ファイルの編集」「データセットの作成」を行います。

ご注意
・本記事ではヒトの解析例を示します。同等の情報を用意できれば他の生物種でも実行可能と思われます。
・データを置く場所・ソフトウェアのインストール先などは、お使いの環境にあわせて読み替えてください。
・defuseは現時点の最新版(0.5.0)を使用しました。


[1] McPherson A, et.al
____deFuse: an algorithm for gene fusion discovery
____in tumor RNA-Seq data.
____PLoS Comput Biol. 2011 May;7(5):e1001138.
____http://www.ncbi.nlm.nih.gov/pubmed/21625565

[2] Sboner A, et.al
____FusionSeq: a modular framework for finding gene fusions
____by analyzing paired-end RNA-sequencing data.
____Genome Biol. 2010;11(10):R104.
____http://www.ncbi.nlm.nih.gov/pubmed/20964841
| hat | 次世代シーケンサー解析 | 09:19 | comments(0) | trackbacks(0) |
融合遺伝子検出ソフトウェア・defuse
最近この業界でホットなキーワードの一つが、「融合遺伝子(fusion gene)」ではないでしょうか?

融合遺伝子は、体細胞ゲノムで生じたリアレンジメント部位から転写された異常な遺伝子のことで、癌との関連が注目されています。
同じ染色体上で生じる場合も、異なる染色体間で生じることもあります。

融合遺伝子が発現するということは、ページがばらばらに混ざってしまったレシピ本を使って料理を作るようなことではないかと思います。
ハンバーグを作っていたのに途中からレアチーズケーキの手順になってしまったら、とんでもないものができてしまいますよね。

融合遺伝子検出ソフトウェアはいくつかありますが、今のところ最も使われているものの一つが、今回ご紹介するdefuseです。
defuseは、カナダのブリティッシュ・コロンビア州がん研究所(BC Cancer Agency)とサイモンフレーザー大学(SIMON FRASER UNIVERSITY)により開発・提供されているオープンソースのソフトウェアです。

明日から3回にわたって、defuseの使用手順をご紹介したいと思います。
明日は、データの準備と、ソフトウェアのインストールを行います。

【第一回】
1. 解析データの準備
2. アノテーションデータの準備
3. ソフトウェアのインストール

【第二回】
4. 設定ファイルの編集
5. データセットの作成

【第三回】
6. 融合遺伝子予測

ご注意
・本記事ではヒトの解析例を示します。同等の情報を用意できれば他の生物種でも実行可能と思われます。
・データを置く場所・ソフトウェアのインストール先などは、お使いの環境にあわせて読み替えてください。
・defuseは現時点の最新版(0.5.0)を使用しました。
| hat | 次世代シーケンサー解析 | 09:20 | comments(0) | trackbacks(0) |
Webinar
イルミナさんのWebinarに注目しています。
Webinarとは、Webサイト上での登録会員制オンラインセミナーです。
すでに、5回のセッションが開催されています。
2月23日にはTruSeqカスタムアンプリコンについてのセッションがあります。

見逃した方も、過去のセッションのPDFを閲覧することができますので、
下記のサイトを覗いてみてはいかがでしょうか。

http://www.illuminakk.co.jp/community/webinar.shtml?ID=BTJ120207-2
| きむ | 次世代シーケンサー解析 | 14:49 | comments(0) | trackbacks(0) |
vcfファイル
samtoolsは、次世代シーケンサデータ解析において、多型を抽出するツールとして広く使用されています。

このsamtoolsで抽出された多型情報は、“vcfファイル”のファイル形式で算出されます。

WETな先生方からよくお聞きする声の中に、

「このvcfファイルの中身がわかりずらい」

という意見がよくあります。

下図がvcfファイルの中身です


多くの情報が入っているため、一見ごちゃごちゃしているように見えますね。

ですので今週は、vcfファイルの中身について記載していこうかと思っています。
| | 次世代シーケンサー解析 | 19:26 | comments(0) | trackbacks(0) |
生命情報工学研究センターワークショップ 2011
1月24日から27日まで、生命情報工学研究センター(CBRC)ではBiWO2011 (Bioinformatics Week in Odaiba 2011) が開催されています。

今日は、「次世代シークエンサーのデータ解析」にフォーカスしたHPCIワークショップがあったので参加しました。
HPCI(High Performance Computing Infrastructure)戦略プログラムは、文部科学省の補助事業です。CBRCは、平成23年度から5年間、「HPCI戦略プログラムにおける人材養成プログラム」として、セミナーやワークショップを実施するそうです。

3月8日(木)〜9日(金)には、「次世代シークエンサー データ解析入門」として、1人1台のPCを使用する実習つきのチュートリアルセミナーが開催されます!!
URL: https://hpci.cbrc.jp/modules/tutorial/index.html

この機会を活かして、多くのノウハウを取り入れたいと思います。
| きむ | 次世代シーケンサー解析 | 15:17 | comments(0) | trackbacks(0) |
Ion Proton Sequencer
2012年1月10日(米国現地時間)、Life Technologies社がIon Protonシーケンサの先行受注を開始しました。
「1日と1000ドルでヒトゲノム解読を実現するシーケンサ」の登場です。
PGM(Ion Personal Genome Machine)の半導体シーケンシング技術を基盤としています。
すでに、ベイラー医科大学、イエール大学医学部、およびBroad Instituteの導入が決定しているそうですよ。

2種類のチップが、紹介されていました。
1. Ion Proton I Chip
 ヒトエクソーム2人分、10Gbase/ラン、1000ドル/ラン
2. Ion Proton II Chip
 ヒトゲノム1人分、100Gbase/ラン、1000ドル/ラン
スループットと費用は予定の段階だそうですが、短期間に低価格でスケールアップを実現するポテンシャルの高さは、半導体シーケンシング技術ならではということになるのでしょうか。

気になる価格ですが、下記のとおりです。
シーケンサ本体のみ: 149,000ドル
シーケンサ本体+サーバー+その他システム: 244,000ドル

プレスリリースはこちらです。
http://www.lifetechnologies.com/us/en/home/about-us/news-gallery/press-releases/2012/life-techologies-itroduces-the-bechtop-io-proto.html
日経バイオテクに、日本語の記事があります。
https://special.nikkeibp.co.jp/ts/article/ab0a/115203/

日本では、今年の秋にリリースされるそうですよ!
シーケンサ本体+サーバー+テンプレート調製自動システムIon OneTouch システム :3360万
詳細はこちらです。
http://www.appliedbiosystems.jp/website/jp/product/modelpage.jsp?BUCD=138106&PLCD=138105&MODELCD=142136
| きむ | 次世代シーケンサー解析 | 09:15 | comments(0) | trackbacks(0) |
アライメント続編
先日は、bwaのアライメントアルゴリズムに関して記述しました。
本日は、その続編を記述致します。
よろしくお願い致します。

bwaでアライメントを行う際は、必ずReferenceゲノムのindex化が必要です。つまり、Suffix Arrayを作成しておく必要があります。

そもそもアライメント(Alignment)とは、“位置合わせ”、“一列に揃える”という意味です。つまり、Referenceゲノムの配列と、Sampleの配列で位置合わせを行う事がアライメントです。
すごく簡単に言ってしまえば、Sample配列の文字列検索をするようなものだと考えます。

それでは、なぜ、ReferenceゲノムのSAを作る必要があったのでしょうか?
ということで、ここからSAの補足を行います。


ReferenceゲノムのSAを作っておくことで、検索が容易になり、かつ同時に位置情報まで判明します!
このようなメカニズムでアライメントを行うから、短時間で必要情報を得る事ができるのですね。
| | 次世代シーケンサー解析 | 18:55 | comments(0) | trackbacks(0) |
アライメントのアルゴリズム
次世代シーケンサの解析を行う際、アライメントという作業を行います。

本日は、アライメントツールとしてbwaを例にとり、アライメントのアルゴリズムに関して記述していきます。
よろしくお願いします。


bwaのInformationでは、Burrows-Wheeler Transform (BWT)をベースとしたアルゴリズムが組まれています(コチラを参照)。

さらに、bwaでアライメントを行うと、“calculate SA coordinate”と表示されます。

BWTって・・・?
SA(Suffix Array)とは・・・?

ということで下記のHPを参考にし、ゲノム配列に置き換えて考えてみました。

wiki ブロックソート
Burrows Wheeler TransformとSuffix Array
Suffix Arrayの概念

大量のデータを用いて検索を行う際に、用いられる方法のようです。
DNA断片の配列を繰り返し並列で並べ、印を付ける(index化)。
その後、ソートを行うことで、SAを作成することができるようです。



今日はここまでです!
| | 次世代シーケンサー解析 | 19:05 | comments(0) | trackbacks(0) |
未来に向けて。次世代シーケンサーと共に。
本日の日経産業新聞に“ITが導く医の進化論”という見出しがあり、次世代シーケンサーによる医学への貢献について特集されています。

記事には、主に下記の3つについて記載されています。

・ここ8年で次世代シーケンスに要するコスト、時間の著しい変動をグラフ化  当初は、ヒトゲノムを全て解読するのに8000日かかっていた。
  現在では1日。

・新薬開発への貢献
  突然変異による薬剤耐性が発し、その原因遺伝子、変異を解明した。

・国際プロジェクト!がん変異のリスト化!
  日本では、肝がんに関連する変異箇所を全てリスト化することが目標。
  ゲノム解析により、その変異箇所を特定する。

特に印象に残った言葉を下記に記します。
 「同じがん細胞でも種類によって様々である事が次世代シーケンサーによりわかってきた。今のペースでシーケンサーの能力が向上すれば、近い将来、患者の細胞のゲノムを検査し、検査結果に応じて治療法を決めるという時代が来る。」



上記の様な原因遺伝子の特定は、ゲノム解析が必須となります。次世代シーケンサー自体の性能ももちろんですが、ゲノム解析の向上も留まることなく行うべきです。
弊社では、学会や日々のお客様との関わりの中で、常に新しいゲノム解析へのヒントを頂き、解析の向上を図っています。上記の様な時代のサポートに貢献していきます。
| | 次世代シーケンサー解析 | 18:33 | comments(0) | trackbacks(0) |
NGSデータ解析Linuxサーバー販売開始しました
次世代シーケンサー(NGS)のデータ解析で利用するソフトには、専用ソフトやオープンソースソフトウェア、さらには自作のプログラムなど、いくつかの選択肢がありますが、自分の研究テーマにぴったりのソフトはなかなか見つからないのが現状かとおもいます。

社内ではNGSデータ解析受託サービスに利用している解析パイプラインがありますので、そのパイプライン自体を購入したいというお問い合わせを多く頂いていました。

そこで、自社のパイプラインを利用しやすいようにインターフェースを整えて、Linuxサーバーに全てを組み込んだ上で販売することといたしました。

→ ホームページはこちら

市販のソフトでは対応が難しいデータ解析パイプラインの出口の部分を研究に合わせてカスタマイズすることで、かゆいところに手が届くサービス・製品をご提供したいと考え、カスタマイズ(ハード・ソフトともに)にも対応させていただいております。

トレーニングなどの運用面でも順次サービスを整備していますので、遠慮無くお問い合わせいただけましたら幸いです。

NGSデータ解析Linuxサーバー
| 社長 | 次世代シーケンサー解析 | 18:41 | comments(0) | trackbacks(0) |
次世代シーケンシング結果を用いたデータ解析
本日は次世代シーケンシングのデータ解析を例として、アノテーションの方法を記述していきます。

次世代シーケンシングのデータは、まずリファレンスゲノムと比較する必要があります。

以下に大まかなデータ解析方法を記載していきます。
公開されているリファレンスゲノムをダウンロードする(例:hg19.fa)。

リファレンスゲノムのデータと次世代シーケンシングのデータを比較(マッピング)します。
- バイオインフォマティクスツールを使用します(例:bwa)。Linux上でコマンドを用いる必要があります。

多型の抽出を行います(例:samtools)。
- 多型情報を得ます。
- 多型情報の結果は、通常膨大なデータとして出力されます。

多型情報に位置や変異など目的に合わせた注釈を入れていきます(アノテーション)。
- comon SNPを除外します。
- exone intron, intergenicなどのアノテーションを付与します。
- さらにexoneに含まれるSNPに対して、synonymous、nonsynonymousなどの詳細なアノテーションを付与します。


上記のプロトコールやアノテーションの内容などは、実験の目的や手法によって大きく変わります。実験によっては、適当なバイオインフォマティクスツール(オープンソースツール)がなく、自前でプログラムを作成する必要がある場合もあります。
そのため次世代シーケンシングによるデータ解析には、お客様が求めている事をしっかりと把握し、その目的に合わせて解析手法を組み立てる必要があり、弊社もそれに対応できるよう心がけています。
| | 次世代シーケンサー解析 | 17:25 | comments(0) | trackbacks(0) |
アノテーション付け(入門編)について
本日は、データ解析を行っているとよく耳にする“アノテーション付け”について記していきたいと思います。よろしくお願いします。

アノテーションという言葉を調べると、“塩基配列データに遺伝子構造や遺伝子機能の情報、また文献情報などを注釈付けする事を指す”とありました。注釈というのは補足的な説明を加えるという意味です。具体的にはどういうことを指すのかを疾患遺伝子解析を例として図に表わします。

次世代シーケンサーから得られたデータを解析する際に困難なことは、以下のことだと思います。
最適なリファレンスゲノムを選択し、次世代シーケンサーから得られるデータと比較。データの量が膨大で容易にできない。
多型情報を得るためにバイオインフォマティクスツールを駆使する必要がある。初めてだとどのようなツールをどうやって使うのかがわからない。
△覇世蕕譴紳新疹霾鵑鮃覆蟾むため、アノテーション付けをして多型の分類をする。

本日はここまで!
明日は具体的なアノテーション付けについて記載します。
| | 次世代シーケンサー解析 | 19:17 | comments(0) | trackbacks(0) |
次世代シーケンサーの参考動画
こんにちは。タノです。
本日は、私が次世代シーケンサーの勉強資料として参考にしている動画を紹介します。

次世代シーケンサを開発していらっしゃる各企業様のHPには、シーケンシング原理の紹介動画が掲載されています。
たとえば
illumina社
http://www.illumina.com/technology/sequencing_technology.ilmn(各機種のサイト内に動画があります)

Applied Biosystems社
http://www.appliedbiosystems.com/absite/us/en/home/applications-technologies/solid-next-generation-sequencing/videos-webinars.html

私にとって動画は、大まかな原理を学習するための必須ツールです。
私のお気に入り動画を下記に記します。

1. 基礎中の基礎から各シーケンサーの原理まで。
http://www.wellcome.ac.uk/Education-resources/Teaching-and-education/Animations/DNA/index.htm

2. Illumina社 Solexaのシーケンシング方法。
http://www.youtube.com/watch?v=77r5p8IBwJk&feature=related

3. Applied Biosystems社 SOLIDのシーケンシング方法。
http://www.youtube.com/watch?v=nlvyF8bFDwM&NR=1

4. DNA microarray法。
http://www.youtube.com/watch?v=VNsThMNjKhM&feature=related

他にもありますが、どれも原理を可視化することができ、大変わかりやすいと思います。
ぜひ次世代シーケンサー初心者という方は参考にしてみてはいかがでしょうか?
| | 次世代シーケンサー解析 | 17:50 | comments(0) | trackbacks(0) |
次世代シーケンサーのデータ解析受託サービス 2
次世代シーケンサーのデータ解析受託サービスを開始して、1年が経とうとしています。同様のサービスを行う会社さまもどんどん増えています。NGS Surfer's Wikiの次世代シーケンサー受託データ解析企業一覧だけでも、12社もの名前が挙がっていました。
そこで今日は、弊社サービスの特徴をまとめてみます。
データ解析受託サービスを企業間で比較する際の、一助となれましたら幸いです。


オーダーメイド
研究に合った解析をご提案します。大量の結果のフィルタリング、データ間比較、GO解析などを行います。


レポート作成
アノテーションを付与した解析データをご提供します。
解析手順のレポートを作成します。


教育・システム
解析手順やビューワの使用方法をレクチャーします。
解析パイプラインの納入およびトレーニングも承ります。

上記の3つ以外の特徴は、創業当初から「バイオインフォマティクス勉強会」を開催している点です。今年の夏は、NGSによるデータ解析をテーマに勉強会を数回行う予定です。ご興味がございましたら、ぜひご参加ください。

| きむ | 次世代シーケンサー解析 | 14:39 | comments(0) | trackbacks(0) |
NGS勉強会@熱海に参加します!
今週末の土日(5/28〜29)に開かれるNGS現場の会 第一回研究会に参加します。
会社からは私を含め2名での参加です。

参加者は、当初予定の60人からなんと100人余りになっているそうです!
しかも全員がポスター発表をしないといけないというルールまであります。

一度に100通りのノウハウを吸収できるなんて、とてもありがたいイベントです。しかも1泊2日の温泉付きです!!!

夜には懇親会もあるので、飲み過ぎないように今から注意しておきます♪
| 社長 | 次世代シーケンサー解析 | 19:09 | comments(0) | trackbacks(0) |
次世代シーケンサーのデータ解析受託サービス 1
アメリエフ株式会社ホームページの
「次世代シーケンサーのデータ解析受託サービス」ページ
より詳細なサービス内容をご覧いただけるようになりました。

弊社では、データ解析の受託に加え、「パイプライン構築/解析ツール開発」も承ります。研究室内に処理系を構築することにより、素早いデータ解析と結果の検証、その後の再解析を行うことが可能になります。

詳細は、下記URLよりお気軽にお問い合わせください。
URL:http://amelieff.jp/enquiry/index.html
| きむ | 次世代シーケンサー解析 | 18:01 | comments(0) | trackbacks(0) |
BEDTools 2
以前、BEDToolsについて簡単にご紹介しましたが、今日は「BEDファイル」についてお話します。
UCSC Genome Browserではアノテーションが図のように表示されますね。
オリジナルのアノテーションを追加したいとき、BEDファイルを作成してカスタムトラックとして追加すると、緑の枠で囲ったアノテーションと同じように表示することが出来ます。

ファイルに記入すべき必須の情報は、3つです。
1.染色体番号
2.開始位置
3.終了位置
他に、9つのオプションを記入することもできます。
必須の情報を見ても分かる通り、BEDファイルは高さの情報をもたないファイル形式です。線の長さと太さで表現しているんですね。

詳しくはこちらのサイトをご覧ください。
URL:http://genome.ucsc.edu/FAQ/FAQformat

【参考】
・BEDTools
URL:http://code.google.com/p/bedtools/
| きむ | 次世代シーケンサー解析 | 17:59 | comments(0) | trackbacks(0) |
BEDTools 1
今日は、BEDToolsを簡単にご紹介いたします。
URL:http://code.google.com/p/bedtools/
BEDToolsは、NGSのデータ解析で役に立つツールです。

まず、bamToBedを使って、BAMファイルをBEDファイルに変換する作業をしてみましょう。
igvで表示すると下記のようになります。
SRX033211.sorted.bedという灰色で表示されているファイルが、
作成したBEDファイルです。

折をみて、詳しい説明をしていきたいと思います。
| きむ | 次世代シーケンサー解析 | 14:57 | comments(0) | trackbacks(0) |
NGSデータ解析勉強会
明日は、第5回バイオインフォマティクス勉強会の日です。
今回の内容は、次世代シーケンサーのデータ解析入門です。

これまでの勉強会の内容を思い出してみると‥
第1回 統計解析ソフトR入門
第2回 SNP解析の実践
第3回 ゲノムマップのviewer
第4回 バイオ研究のためのLinux入門
などなど、内容は様々です。

この場を通して、多くの皆様と出会えたことを大変幸せに思います。
勉強会後に開催される懇親会で、学ぶことも多く、とても楽しく有意義な時間を過ごしています。懇親会からの参加もお待ちしておりますラブラブてれちゃうラッキー
| きむ | 次世代シーケンサー解析 | 18:33 | comments(0) | trackbacks(0) |
DDBJ Sequence Read Archive (DRA) 3
DRASeqarch
http://trace.ddbj.nig.ac.jp/DRASearch/

で、オブジェクトごとにAccession番号のついたリストを表示できると紹介しました。
ちなみに、オブジェクトごとにAccession番号のプレフィックスが異なります。
DRA : Submission(登録者など)
DRP : Study(研究のタイトルやタイプ、文献へのリンクなど)
DRX : Experiment(用いたライブラリーやシーケンサーなど)
DRS : Sample(Taxonomy ID、取得方法など)
DRR : Run(ランが行われた日付、シーケンシングを行った組織名など)
プレフィックスが、DRXやERXだと「実験の情報だな」とすぐに判断できます。

詳しく知りたい方はこちらをどうぞ
http://trace.ddbj.nig.ac.jp/dra/documentation.shtml

| きむ | 次世代シーケンサー解析 | 18:16 | comments(0) | trackbacks(0) |
NGSデータ解析勉強会の準備中
今週末23日(土)に行われる「バイオインフォマティクス勉強会 次世代シーケンサーのデータ解析入門」の準備をしています。

おかげさまで増席後も満席となっておりますので、期待の大きさに少々プレッシャーを感じつつも、自分たちが行っているNGSのデータ解析のフローを簡単にではありますがご紹介したいとおもいます。

また、その際の簡易的なスクリプトもお持ち帰りいただいて、ご利用いただけたらと考えております。


私が学生時代にゲル板でシーケンシング(といえるかわかりませんが)をしていた頃は、3日間の実験で100bpとかそのオーダーでした。

それが今や、NGS1台を1週間ランさせるだけ(ラン途中の大半の時間は放置ですし)で、100Gbpのオーダーで配列が読めるわけですから。時間や精度を買っていると考えられますね。


歳をとると昔話が多くなってしまいます。(笑)
いかん・いかん。

身も心も若くいたいとおもっていますが、皇居周辺を9kmほどランニングしただけで、2日間筋肉痛です。(-_-;)

脳みそも肉体も鍛錬が必要ですね。
| 社長 | 次世代シーケンサー解析 | 18:09 | comments(0) | trackbacks(0) |
DDBJ Sequence Read Archive (DRA) 2
DRASeqarch
http://trace.ddbj.nig.ac.jp/DRASearch/

で欲しいデータを検索しましょう!
Statisticsを眺めると、どのくらいの数のどんなデータがどのように管理されているか、概要をつかむ事が出来ます。

「Released Entries」の各Typeをクリックすると、Accession番号のついたリストが表示されます。
各Typeの大まかな内容は下記の通りです。
Submission : 登録者など
Study : 研究のタイトルやタイプ、文献へのリンクなど
Experiment : 用いたライブラリーやシーケンサーなど
Sample : Taxonomy ID、取得方法など
Run : ランが行われた日付、シーケンシングを行った組織名など

詳しく知りたい方はこちらをどうぞdowndowndown
http://trace.ddbj.nig.ac.jp/dra/documentation.shtml
| きむ | 次世代シーケンサー解析 | 18:59 | comments(0) | trackbacks(0) |
DDBJ Sequence Read Archive (DRA) 1
「試しに、次世代シーケンサーで得られたデータの解析をしたいなぁ」
という時、皆さまはどうしていますか??

日本にいる皆さまは、ココdowndowndown
DDBJ Sequence Read Archive (DRA)
http://trace.ddbj.nig.ac.jp/dra/index.shtml

から、データをダウンロードしている方が多いのではないでしょうか。

「Search」タブをクリックすると、新しいウィンドウが開きます。

説明が必要ないくらい、見やすくて使いやすいサイトですが、
少しだけ使い方を説明します。
| きむ | 次世代シーケンサー解析 | 18:51 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第四世代シーケンサー5」
Ion Personal Genome Machine (PGM™) sequencerについて、この場を借りて、簡単にメモしたいと思います。

・原理について
半導体チップ上の各ウェルでシーケンシングが行われる。
ヌクレオチドがポリメラーゼによってDNA鎖に取り込まれる時に、副産物として水素イオンが放出される。PHの変化は、イオンセンサーに感知され電圧に変換される。電圧の変化は、半導体に伝えられ出力される。
わかりやすい図や説明はここで見られます。
downdowndowndown
http://www.iontorrent.com/the-simplest-sequencing-chemistry/

Ion 314 Sequencing Chipについて
スケール:130万ウェル
データ量:10Mb以上
リード長:100base

クローニングはemulsion PCRで行うようです。

【参考】
・Ion Torrent Systems
http://www.iontorrent.com/
| きむ | 次世代シーケンサー解析 | 19:40 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第四世代シーケンサー4」
Life Technologies社が、12月14日にIon Personal Genome Machine (PGM™) sequencerの発売を開始しましたね。

製品紹介には、3つのポイントが大きく書いてありました。
・scalability
massively scalable semiconductor technology
・simplicity
natural biochemistry, no optical components
・speed
runs in about two hours, not days or weeks

私が注目しているのは、これまで公表されていなかったlibrary and target selection protocolsです。
折を見て、詳細を見ていきたいと思います。

【参考】
・invitrogen by Life Technologies(プレスリリース)
http://www.lifetechnologies.com/news-gallery/press-releases/2010/life-techologies-lauches-io-pgm-sequecer.html
・iontorrent(製品紹介)
http://www.iontorrent.com/lib/images/PDFs/ion_prod_a_121410.pdf
| きむ | 次世代シーケンサー解析 | 20:54 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第三世代シーケンサー4」
Next Generation Sequencing Platforms「第三世代シーケンサー」として
PACBIO RSについて少し書きました。
Next Generation Sequencing Platforms「第三世代シーケンサー1」
Next Generation Sequencing Platforms「第三世代シーケンサー2」
Next Generation Sequencing Platforms「第三世代シーケンサー3」

PACBIOのSMRT(Single Molecule Real Time)シーケンシングと同時期に話題になった、Life TechnologiesのQuantum Dot FRETシーケンシングについて触れたいと思います。
Quantum Dotは半導体物質からなるユニークな蛍光物質です。2006年から様々な製品として提供されています。
このInvitrogenのQuantum Dot技術と、VisiGen BiotechnologiesのFRET技術を組み合わせることで、高い精度の1分子シーケンシングを目指しています。
確かに、1分子シーケンシングでネックになる蛍光検出の問題は、CCDカメラの性能向上を目指す方法もありますが、明るく安定した蛍光を用いる方法も取り入れることで、解決の道が開けそうですね。

近いうちに、Quantum Dot FRETシーケンシングシステムの発表があることを期待して、来週おさらいしたいと思います。

【参考】
・invitrogen by Life Technologies
http://www.invitrogen.jp/qdot/products.shtml
| きむ | 次世代シーケンサー解析 | 22:11 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第四世代シーケンサー3」
今日は、Ion Torrent SystemsのPostLight sequencing technologyについて紹介したいと思います。

チップは、CMOS(complementary metal-oxide semiconductor:相補型金属酸化膜半導体)という半導体のチップを使用します。
CMOS(シーモスと呼ぶらしい)構造は、メモリやマイクロプロセッサなどにも広く利用されているそうです。電圧駆動型のトランジスタで構成されているので、ヌクレオチドがDNAポリメラーゼに取り込まれる時に放出される水素イオンの情報を電圧に変換して、CMOSに渡すと、信号に変換してくれるということになります。

Iontorrentのホームページでは、
半導体技術とシーケンシング技術の出会いを
Watson meets Moore
と表現していました。
半導体技術の導入は、シーケンシングチップの低コスト化やスループットの増大に大きな恩恵をもたらしてくれるのではないでしょか。

【参考】
・Ion Torrent Systems
http://www.iontorrent.com/
| きむ | 次世代シーケンサー解析 | 20:00 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第四世代シーケンサー2」
Ion Torrent Systemsのシーケンサーについて書こうと思います。

ケミカルデータを直接デジタルデータに変換するため、第三世代シーケンサーでは重要な技術であったカメラ/レーザー/スキャンが要らないのです。
そのため第四世代シーケンサーで使われている技術は「PostLight sequencing technology」と呼ばれています。

技術革新により、全ての研究室やクリニックが導入できる価格やサイズを目指すことが出来るようになったそうです。
ちなみに、サイズはデスクトッププリンターと同じくらいで、重さは50poundsおよそ22.7キロだとか。

PostLight sequencing technologyを支えているのは、半導体技術です。
来週、勉強していきたいと思います!

それでは、良い週末を〜

【参考】
・Ion Torrent Systems
http://www.iontorrent.com/
| きむ | 次世代シーケンサー解析 | 20:29 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第四世代シーケンサー1」
キムです。皆さまこんにちは。
今日は、オフィス家具の搬入やPCの設置などで、
経理のスタッフさんと大わらわでした。

とりあえずは、床に転がっている物はなくなったので、
分子生物学会に出張している山口さんに顔向けできそうです(笑)

昨日まで、第三世代のSingle Molecule Real Time sequencing technologyについて書きました。
明日からは、第四世代のPostLight sequencing technologyについて書きたいと思います。
| きむ | 次世代シーケンサー解析 | 18:13 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第三世代シーケンサー3」
先週、PACBIO RSのSMRT(Single Molecule Real Time)シーケンシングを実現するための3つの技術革新シーケンスオプションについて少し書きました。

彼らは、
・flushing, scanning, washingステップが不要
・PCR増幅が不要
・長いシーケンス長(平均1000base)
・時間短縮(sample preparation:one day, run time:30min)
・使いやすさ(SMRTbell sample preparation protocolはシンプルで早い、タッチスクリーンのインターフェース、標準と互換性の高いデータフォーマット)
・柔軟性と精度(多様なプロトコル、コストとスループットを調節可能、SMRT Cellの数を調節)
・kineticsの測定と評価が可能
を優れている点として挙げています。

カメラの性能向上が課題だと思いますが、技術の進歩は目覚ましいですね。
外気や振動に左右されず安定して性能を発揮できるのか、
またコストはどのくらいまで下がりそうなのかが気になるところです。

【参考】
・Pacific Biosciences
http://www.pacificbiosciences.com/
| きむ | 次世代シーケンサー解析 | 19:05 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第三世代シーケンサー2」
昨日から、PACBIO RSについて書いています。
テンプレートとなるSMRTbell libraryは、DNAサンプルの両端でループ状になるアダプターをライゲーションすることで作成されます。シーケンスオプションは下記の通り。

・Standard sequencing
標準的なシーケンシング
長い一本のリード

・Circular consensus sequencing
ぐるぐる循環しながら同じ配列を何度も読むことで、精度がアップ。
複数のリード(forwardとreverse含む)

・Strobe sequencing
レーザーを断続的に照射することで、レーザーによるポリメラーゼへのダメージを抑えることが出来る。長くシーケンシングできので、構造多型の解析に有効。
paired readやmate pairsの様に、リード間の距離情報を推定できるリード

・Combining sequencing modes
複数のテクニックを使用して配列データを作り出せる。

まず、Strobe sequencingでscaffold作成
次に、Standard sequencingで長く読む
最後に、Circular consensus sequencingで変異を確認
という流れの解析に一台で対応できるというわけです。

【参考】
・Pacific Biosciences
http://www.pacificbiosciences.com/
・Strobe sequencing
http://bioinformatics.oxfordjournals.org/content/26/10/1291.abstract
| きむ | 次世代シーケンサー解析 | 09:10 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第三世代シーケンサー1」
PACBIO RSは、DNAポリメラーゼを利用することで、リアルタイム一分子シーケンシングを可能にしています。
ホームページに、SMRT(Single Molecule Real Time)シーケンシングを実現するための3つの技術革新について書かれていました。
簡単にまとめてみます。

・The SMRT Cell
シーケンシングは、SMRT Cell上で行われる。
150,000個のZMW(zero-mode waveguide)がSMRT Cellに並んでいる。(同時にモニターできるのは75,000個)
ZMWは、ガラス基板上の金属フィルムに空けられた直径数十ナノメートルの穴であり、底にはDNAポリメラーゼが固定されている。
レーザー(600nm)を照射すると、ガラスを通過して、ZMWに入る。
ZMW内で指数関数的に減衰するため、底面の30nmのみが照らし出される。

・Phospholinked nucleotides
ヌクレオチドではなくリン酸基に蛍光物質を付加している。
DNAポリメラーゼによって、ヌクレオチドがDNA鎖に取り込まれると、リン酸基に付加されていた蛍光物質は解き放たれる。

・The PacBio RS
大きな開口数の対物レンズと単一光子が検出可能な4つの高感度カメラを備えている。蛍光を検出して、リアルタイムで塩基配列データに変換する。

こうした技術の組み合わせにより、ヌクレオチドが取り込まれる瞬間の蛍光を観測することに成功したんですね。
| きむ | 次世代シーケンサー解析 | 20:21 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第三・四世代シーケンサー」
今週は、第三世代シーケンサーについて少し調べて書いていこうと思います。

3つのテクノロジー革新でSingle Molecule Real Time sequencingを実現したPacific BiosciencesのPACBIO RSや、
PostLight sequencing技術で期待を集めているIon Torrent Systemsに焦点を当てていきます。

【参考】
・Pacific Biosciences
http://www.pacificbiosciences.com/
・Ion Torrent Systems
http://www.iontorrent.com/
| きむ | 次世代シーケンサー解析 | 18:04 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第二世代シーケンサー4」
今週は、第二世代シーケンサーを比較する際のポイントについて、とりとめなく書いてきたので、そろそろまとめたいと思います。

・データ量(リード長、目標リード長、タグ数)
・コスト(サイズ、導入、維持、ランニング)
・時間(ハンズオンタイム、ランニングタイム)
・精度
・データ解析(ソフト、パイプライン)
・データ管理/保管
・テクニカルサポート/コミュニティー

特にデータ解析に関しては、かなり苦労しているというお話をよく伺います。
今、大きな変革期を迎えているこの業界でお仕事をさせていただいているのは、本当に幸せなことだなとつくづく思います。今週は、のんきなつぶやきで締めくくらせていただきます〜てれちゃうポッイヒヒ



| きむ | 次世代シーケンサー解析 | 18:32 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第二世代シーケンサー3」
今日も第二世代シーケンサーに関して書きたいと思います。

Complete Genomics社のサービスが、とても興味深いので触れておこうと思います。CGA™ Service (Complete Genomics Analysis Service) は、包括的な受託サービスを目指しています。
つまり、
・シーケンシング装置や試薬を販売
・解析サーバを案内
・解析ソフトを提供/紹介
・トレーニング/テクニカルサポート/コミュニティーの場を提供
などを行う「シーケンシング装置と試薬の開発・販売」を主体とした機器メーカーさんとは異なり、
・サンプル調整
・シーケンシング
・データ解析
・データ管理/保管
(amazon web services)
を組み合わせたまさに「次世代シーケンシンサー研究の包括的支援サービス」を目指しているのです。

対象はヒトゲノム、アプリケーションは下記の4つなので、
かなり特化(限定?)したメニューを提供しています。
・Cancer Research
・Mendelian Disease Research
・Rare Variant Disease Research
・Clinical Trial Optimization
HPを見ただけなので、他にもメニューがあるかもしれません。
ご興味がある方はチェックしてみてください。

【参考】
・Complete Genomics
http://www.completegenomics.com/
・amazon web services
http://aws.amazon.com/jp/
| きむ | 次世代シーケンサー解析 | 14:49 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第二世代シーケンサー2」
第二世代シーケンサーに関して書きたいと思います。
各社ともにリード長を延ばす努力を行っています。
今年の目標は
Roche 454 1000base
Illumina 150base
ABI SOLiD 100baseでした。
他にも、
GS Junior Bench Top System
SOLiD PI Systemなど、
ベンチトップサイズへの小型化も盛んですね。

以前、比較するポイントとして
・データ量(リード長、タグ数)
・コスト(導入、維持、ランニング)
・時間(ハンズオンタイム、ランニングタイム)
・精度
を挙げました。
これらのテーマはいわば「シーケンサー開発の必須課題」
なのかなと思います。

一方、
・目標リード長
・ベンチトップ型の有無
は、研究用途やその市場まで視野に入れた問題であり、
最適化されていく中で定着したり淘汰されるのだと思います。

さて、私個人的にとても楽しみなのはPolonatorです。
かずさDNA研究所が導入したことで話題になりましたが、プロトコルとソフトウェアがオープンソース(!)、試薬も低価格で販売するそうです。
しかも、ユーザーが自社製品以外の試薬を使うのも大歓迎。低価格でより良い試薬があるのならば、むしろ情報を共有しましょうというのです。
Open evolution is the accelerant that will drive this platform forward.Please join us!
とHPにも力強い言葉がありますから、ユーザーコミュニティーが育てば爆発的に普及しそうな予感がします。
| きむ | 次世代シーケンサー解析 | 19:04 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms「第二世代シーケンサー1」
第二世代シーケンサーと呼ばれる
Roche 454
Illumina
ABI SOLiD
を比較する際のポイントを挙げてみます。

・1ランあたりのデータ量(リード長、タグ数)
・コスト(導入、維持、ランニング)
・時間(ハンズオンタイム、ランニングタイム)
・精度
でしょうか。

| きむ | 次世代シーケンサー解析 | 15:06 | comments(0) | trackbacks(0) |
Next Generation Sequencing Platforms
社長のリクエストにお応えして、まずはPlatformごとに書いていこうと思います。

Roche 454
Illumina HiSeq 2000
ABI SOLiD
Helicos
Complete Genomics
Polonator
Pacific Biosciences SMRT
Ion Torrent

ざっと書いてみました。
次世代シーケンサーも多様になりましたね。

明日から、少しずつ見ていきたいと思います。
| きむ | 次世代シーケンサー解析 | 20:17 | comments(0) | trackbacks(0) |
次世代シーケンサーのデータ解析
私の探し方が悪いのかもしれませんが、次世代シーケンサーのデータ解析について、日本語で網羅的に記述されたものが無いように感じます。

弊社でも次世代シーケンサーのデータ解析や、二次解析・三次解析用のパイプラインを構築しているので、社内でまとめた文章をブログ記事にまとめていきたいとおもいます。
| 社長 | 次世代シーケンサー解析 | 23:09 | comments(0) | trackbacks(0) |
次世代シーケンサーSOLiDのデータ解析
次世代シーケンサーの一つ、ライフテクノロジーズのSOLiDのデータ解析をしています。

オープンソースの解析ソフトは多数公開されていますので、そのソフトをつなぎ合わせて目的の解析をするパイプラインを構築するのですが、そのパイプラインの途中に組み込む解析ソフトやファイル書式のコンバーターなどが存在しない場合は(存在しない場合が多数ですが・・・)、ソフトウェアを自前で作成しています。

今日は作成したソフトウェアをお客様にプレゼンしました。数千万のデータから目的の情報を抽出することは通常のパソコンでは困難なので、専用ソフトを作成しているわけですが、お客様は抽出された結果(統計処理した結果)を初めて見て、予想外の結果に驚いていました。

研究上は新たな課題の出現ですが、ソフトを作成した弊社としては研究に貢献できた喜びを味わっています。(笑)
| 社長 | 次世代シーケンサー解析 | 23:26 | comments(0) | trackbacks(0) |
   1234
567891011
12131415161718
19202122232425
262728    
<< February 2017 >>

このページの先頭へ