ショートリードの憂鬱 - 次世代シーケンサー

第2世代高速シーケンサー（Illumina, SOLiD, 454）とそのアプリケーションについて。できるだけ最新に近いことと、トレンド、面白いと思った論文やアプリケーションを書いています。シーケンサーに興味のあるひとは是非覗いていってください。コメントも大歓迎です。　一分子シーケンサー　PacBioについては、姉妹ブログ「パックマンの挑戦　http://pacbiobrothers.blogspot.com/　」　を、10X Genomicsについては「くろみうんの冒険」を覗いてみてください！

2011年7月29日金曜日

SeattleSeq Annotation : WebベースのHuman SNP機能予測システム

SeattleSeq　というサイトをご存じでしょうか？
ワシントン大学が運営する、ヒトSNV機能予測データベース＆システムです。

Web上でデータをやり取りします。
ユーザはSNVのリストファイルをアップロードして、欲しいアノテーションにチェックを入れ、メールアドレスを入力して、ボタンを押して、後は待つだけ。
インターネット上にデータを投げるわけですから、機密性が低い、公共データやデモデータ、自分で責任の負えるデータ、でやるのが無難です。

ファイルのフォーマットは割合と自由です。ここに詳しく
公共データのReadファイルから、マッピング、SNV検出してきた結果ファイルVCFフォーマットを入れてみましょうか。
VCFフォーマットの場合はちょっとコツがいりまして、SNVのみの結果にしなければいけません。
もしSNVとInDelが混ざって出力されているVCFの場合、InDelのデータは除きます。
VCFからSNVとInDelを分けて別々のファイルに保存するには、
こんなawkファイルを作って、

----------------Separate_SNV_and_InDel.awk     ここから------------------------------
/^#/    {
    print $0 > "snv_only.vcf";
    print $0 > "indel_only.vcf";
    next;
    }

/^[^\t]+\t[0-9]+\t[^\t]*\t[atgcATGC]\t[a-zA-Z]\t/   {
    print $0 > "snv_only.vcf";
    next;
    }

    {
    print $0 > "indel_only.vcf";
    next;
    }
------------------------ここまで---------------------------------

awk -f separate_snp_indel.awk [あなたのファイル.vcf]

なんてコマンドを打つと、2つのファイルができて、snv_only.vcf ファイルにはSNVのみのデータが、
indel_only.vcf　にはInDelのみのデータが作られます。
このコマンドはあるBio-info仲間が教えてくれたものですが、簡単で重宝しています。

さて、もうひとつ、VCF 4.0 ファイルであることを示すため、できた snv_only.vcf ファイルの先頭行に
##fileformat=VCFv4.0
という一行を入れてあげます。

sed -e "1i ##fileformat=VCFv4.0" [あなたのファイル.vcf] > [あなたのファイル2.vcf]
これで準備OK。

SeattleSeqの画面を見てみましょう。　http://snp.gs.washington.edu/SeattleSeqAnnotation131/
ヒトゲノムのバージョンが、hg18 と　hg19 の２つあるので、hg19を選択。

メールアドレスを入力し、ファイルをUpします。
あとは欲しいアノテーションにチェックを入れて、Submitボタンを押せばOK！
画面が切り替わり、ジョブの進行状況が確認できます。
しばらくすると、メールアドレス宛てに、メールで結果ファイルのURLが届きます。

そこをクリックして、ファイルをダウンロードします。

これがアップロードしたVCFファイル

そしてこれが、SeattleSeqの結果

New　で示したところが、SeattleSeqで新しくつけられたアノテーションですね。
VCFフォーマットを保ったまま、INFO　のところにrs番号をはじめ、SNVの分類、アミノ酸置換の種類、配列、GranthamScore、PolyPhen分類、などが新たについています。

以下の例で言うとこれらは、rs番号はrs3748597、SNVの分類(FD)はNon-synonymousでmissense、アミノ酸置換（AC）はイソロイシンからバリン（ILEからVAL）、GranthamScore（GS）は29、PolyPhen分類（PH）はpossibly-damaging、というようなことが書かれています。

1 888659 rs3748597 T C 149 . DP=11;AF1=1;CI95=1,1;DP4=0,0,3,8;MQ=60;
FQ=-60;DN;DT;DA=C/T;GM=NM_015658;GL=NOC2L;FG=missense;FD=missense;

AC=ILE/VAL;PP=300/750;GS=29;PH=possibly-damaging;CP=0.5990;CG=2.300;AA=C;
CN=2294,3274,30362;HA=12.1;HE=6.7;HC=6.2;DG;DV=by-frequency,by-cluster;
PS=MAAAGSRKRRLAELTVDEFLASGFDSESESESENSPQAETREAREAARSPDK
PGGSPSASRRKGRASEHKDQLSRLKDRDPEFYKFLQENDQSLLNFSDSDSSEEEE
GPFHSLPDVLEEASEEEDGAEEGEDGDRVPRGLKGKKNSVPVTVAMVERWKQA
AKQRLTPKLFHEVVQAFRAAVATTRGDQESAEANKFQVTDSAAFNALVTFCIRD
LIGCLQKLLFGKVAKDSSRMLQPSSSPLWGKLRVDIKAYLGSAIQLVSCLSETTV
LAAVLRHISVLVPCFLTFPKQCRMLLKRMVVVWSTGEESLRVLAFLVLSRVCRH
KKDTFLGPVLKQMYITYVRNCKFTSPGALPFISFMQWTLTELLALEPGVAYQHA
FLYIRQLAIHLRNAMTTRKKETYQSVYNWQYVHCLFLWCRVLSTAGPSEALQPL
VYPLAQVIIGCIKLIPTARFYPLRMHCIRALTLLSGSSGAFIPVLPFILEMFQQVDFN

RKPGRMSSKPINFSVILKLSNVNLQEKAYRDGLVEQLYDLTLEYLHSQAHCIGFP

ELVLPVVLQLKSFLRECKVANYCRQVQQLLGKVQENSAYICSRRQRVSFGVSEQ

QAVEAWEKLTREEGTPLTLYYSHWRKLRDREIQLEISGKERLEDLNFPEIKRRKM

ADRKDEDRKQFKDLFDLNSSEEDDTEGFSERGILRPLSTRHGVEDDEEDEEEGEE

DSSNSEDGDPDAEAGLAPGELQQLAQGPEDELEDLQLSEDD* GT:PL:GQ 1/1:182,33,0:63

ここで登場した、granthamScoreとpolyPhen、ほかにもscorePhastCons（CP）、consScoreGERP（CG）などは、塩基の変異によりアミノ酸が置換されるとき、これがタンパク質にどれくらい影響を及ぼすかの予測された指標です。

Grantham Scoreは、アミノ酸が変わるときに隣近所のアミノ酸との間にどれだけ化学的な影響を及ぼすか、分子量や極性などをもとに数値化しています。

conservative (0-50)、moderately conservative (51-100)、moderately radical (101-150)、radical (>151) というような分類もあるようです（Li et al. J Mol evol. 21, 58-71(1984).）

PolyPhen は、他の生物種とのホモログの配列情報から、変異があった箇所の配列保存性、化学的な特徴、タンパク質構造やドメイン情報をもとに、この変異がタンパク質の機能に及ぼす影響を数値化します。

probably damaging (>2.00)、possibly damaging (1.50-1.99)、potentially damaging (1.25-1.49)、borderline (1.00-1.24)、benign (0.00-0.99) という分類があるようで、SeattleSeqではこの分類表記になっていますね。

ちなみに、PolyPhenは現在、PolyPhen-2というのもありますが、こちらはアカデミック・非営利機関のみ利用できます。

SNVを探してから、missense のSNVのみを抽出し、そのアミノ酸変異がタンパク質にどれくらい影響を及ぼすのか、という研究は昔からありました。

現在はNGSのおかげでどんどん新しいSNVが見つかっています。

その塩基の変異が、タンパク質レベルでどのように働くのか、というテーマは今後も続いていくでしょう。

まずはPolyPhenなどの予測アルゴリズムを用いて、分子レベルの機能変化をスコア化し、数個に絞った後で、（今はここで終わっている論文がほとんどですが）今後はタンパク質の3次元構造モデリングなどを行ったり、アッセイに持ち込んで本当にタンパク質の活性が変化するのかを確かめたり、する研究も増えていくかもしれません。

2011年7月27日水曜日

Exome 解析　non-synonymous SNVを見つけた後は・・・

Agilent社のSureSelect、Illumina社のTruSeq、NimbleGen社のSeqCap といえば、最近盛んに宣伝されているExomeの実験用キットです。
Exomeというのは、ゲノム全体ではなく、その中のExon、実際はExon＋α　の配列のみを選択してシーケンスする解析方法です。
原理はまず、Exonの部分配列DNA（プローブ）が用意されていて、断片化されたゲノムと、溶液中またはアレイ上でハイブリし、Exonプローブにハイブリされた（キャプチャーされた）ゲノム断片だけを抽出します。
そのExon部分配列のみをシーケンスするというわけです。
実際Exonの上を何塩基ごとにプローブが設定されているのか、という情報は非公開ですが、Exon全てをキャプチャーできるようには設計されているそうです。

このキットが普及されるにつれ、Exome、特にヒトExomeの研究が盛んになってきました。
それまでは1000GenomeやCancer Genome Instituteのようなゲノムセンターでの研究が一般的でしたが、昨年末ごろから一般の大学の研究室レベルでも、ヒトExome実験が行われるようになってきたようです。

実験プロトコールもいろいろ大変なのでしょうが、私はデータ解析屋なので、データが出てからの料理に興味があります。

世の中に、フリーツールはたくさんありますが、Exome解析は、それらをいくつも組み合わせていく、一見面倒くさいけれどもエレガントなワークフローです。

ちょっと検索すると、
NGS Surfer's Wiki のリシークエンス
BioStar
などにパイプライン（ワークフロー）が出てきます。

まとめると、

リードをゲノムにマッピングし　（BWA）
冗長性のあるリード・Duplicateを除去し　（SamtoolsやPicard）
キャプチャーしたExon領域だけを取り出し　（Bedtools）
SamtoolsでSNVを抽出し
SNVにアノテーションをつける non-synonymous SNV, missense, frame-shift など

というところまでを行っています。
SNVというのはSNPとほぼ同意語で、変異か多型かの違いです。厳密には異なりますが、ここでは同じ意味とします。

さて、Exomeの論文はたくさんありますが、その中でいいな、と思ったものを3つ

Vissers et al. A de novo paradigm for mental retardation. Nature Genetics 42, 1109-12 (2010).
Timmermann et al. Somatic mutation profiles of MSI and MSS colorectal cancer identified by whole exome next generation sequencing and bioinformatics analysis. PLoS One 22, e15661 (2010).
Ng et al. Exome sequencing identifies MLL2 mutations as a cause of Kabuki syndrome. Nature Genetics 42, 790-3 (2010).

3のKabukiシンドロームは、10人の患者データセットがあるのですが、10人全部に共通する新規レアSNVを見つけようとしたらうまくいかなかったのですね。　そこで10人全員ではなく9人に共通、8人に共通、7人に共通というふうにレベルを下げていったのです。
そうするとNonsense置換またはフレームシフトを起こしていたMLL2遺伝子上の変異が、10人中7人に見つかったのです。
後はちゃんとCGHアレイとサンガーシーケンスで確認しています。

2のMSI/MSSは、直腸癌の種類（microsatellite instable / stable）が違う6人の患者サンプルで、MSIとMSSとで見つかった変異にどう違いがあるか、を見ています。
それぞれのサンプルで55,000程のSNVをリストした後、遺伝子の中にあるか、タンパクコード領域にあるか、dbSNP/1000Genomeに登録されていないか、カバレージは十分にあるか、Somaticな変異かどうか、とフィルタリングしていき、最終的に「タンパク質の機能を変化させる変異であるか」というところまで行き、数十程のSNVまで絞り込んでいます。
せっかくタンパク質の3D立体構造まで表示しているので、もう少しここから先が欲しかったです。

1、2、3、共通するのが de novo SNVを見つけた後にその変異が引き起こすであろうアミノ酸置換またはフレームシフトが、タンパク質の機能にどのような影響を及ぼすか、をスコア化していることです。

1はPhyloPとGrantham、2はPolyPhen（Polymorphism Phenotyping）とMutation Taster、3はGERP　というプログラムを使ってスコアを計算しています。
このような潜在的なnon-synonymous SNVの機能予測プログラムは他にも、SIFT（Sorting Intolerant from Tolerant）、PolyPhen‐2などがあります。
私も全部を知っているわけではないので（知っていなくてはいけないのでしょうが）、これらのうちの
いくつか、多分PolyPhenとGranthamを、SeattleSeq というWebツールと一緒に紹介します。

2011年7月22日金曜日

Trinity: de novo Transcriptome の救世主？

前回の最後にチョイ出ししたTrinityというソフト、使ってみました。

http://trinityrnaseq.sourceforge.net/
Broad Instituteで開発されたこのソフトは、Perlで動かして使い、Linux上でテストされているそうです。
先のサイトからダウンロードしてインストールすると、

こんなディレクトリができあがり、/trinityrnaseq/sample_data/test_Trinity_Assembly の下に、サンプル用のペアエンドデータが作られているはずです。

とりあえずどんな結果が出てくるのかが待ち遠しいひとは、まずこのサンプルデータをランしてみましょう。　runMe.sh を打つようにマニュアルでは言っていますが、あえて（わざわざ）このシェルの中のコマンドを打ってみました。　（以下、一行です）

Trinity.pl --seqType fq --left reads.left.fq --right reads.right.fq --SS_lib_type RF --paired_fragment_length 280 --min_contig_length 305 --run_butterfly --CPU 2 --bfly_opts "-V 10 --stderr"

--seqType fq で、このリードがFastqだということを示し、
--left / --right で、ペアエンドの配列ファイルを指定し、
--SS_lib_type RF で、このリードがストランドスペシフィック（日本語で何と訳すんだろう）でリバース・フォワードの向きだと宣言し、
--paired_fragment_length 280 はフラグメントの長さが280ベースだと言い、
--min_contig_length 305　で、コンティグは最低305塩基必要とし、
--run_butterfly --CPU 2 で、一緒にbutterfly もランして、その時はCPUを2個使えと言い、
--bfly_opts "-V 10 --stderr"　これは何かよくわかりません

ちなみにシングルエンドリードのときは、--left / --right　を　--single にすればOKです。
さて、このコマンドを流すと、trinity_out_dir　というディレクトリが作られて結果が保存されますが、--output で出力先を決めることもできます。

結果で大事なのは、この中でTrinity.fastaです。　これが最終的なコンティグです。

これがそのファイルを開いたところです。　名前（ID）と、Contigの長さ、リードをマッピングした時のFPKM、がヘッダーに書かれているMulti-Fasta形式ですね。

Trinityは、Inchworm（シャクトリ虫）、Chrysalis（サナギ）、Butterfly（チョウ）という名前の3つに分かれています。

Inchwormでは、k-mer ベースでの de Bruijn graph でアセンブります。　Alternative Variantがある時でも、そのうちの代表的な転写産物をつくるようにアセンブるので、完全なTranscriptome、Isoformを反映したTranscriptomeができるわけではありません。　

デフォルトは k=25 です。　k-1 の重なりでContigを伸長していき、これ以上伸長しなくなったら止まるそうです。

Inchworm.K25.L48.fa という名前のファイルができます。　

メモリが足りないとこのファイルすらできません。

次のChrysalisは、名前の通り、動いているのか止まっているのかわからない、ユーザーを不安にさせるコマンドです。

冗談です。

Inchwormで作られたInchworm.K25.L48.fa のコンティグたちから似ている配列のContigを集め、Alternative Splicingや、Paralogっぽい配列をまとめる、極めて重要な作業をしている、そうです。

k-1 の重なり具合をもとに、Inchwormで出来上がったコンティグで、de Bruijn graph を作っていくつものパスを作るらしいです。

最後のButterflyは、いかにもチョウらしく、何百行ものコマンドラインが現れては流れていきます。　ChrysalisでまとめられたIsoform/Paralogのコンティグをもとに、リード配列を使って、あいまいさの残る部分を修正し、全体を再構築しているそうです。　

この過程でFPKMも算出されます。

文献によると、ベンチマークでは、ABySS、TransABySS、SOAP denovo、と比べると、完全長の遺伝子、Isoformの数ともにTrinityは優れているそうです。

しかーし、

このツールは恐らくショートリード向けでしょう。　k-mer が25と短いのがその理由です。

また、デモデータではなく、本番データをアセンブるには、結構な量のメモリが必要です。

76塩基程度のイルミナペアリード100万本あたり、1GB のメモリが必要だそうです。

とうことは、1,000万本で10GB、1億本で100GB !!!

彼らのお勧めは、256GB の、22コアCPUのサーバーとのこと。

とてもじゃないですが、常識的なお値段で買える代物ではありません。

と、いうことで、私は多分、リード数が2,000万本を超えるような de novo Transcriptome でこのツールを使うことはしないでしょう。　

分散化などをうまく施せば別ですが。

2011年7月18日月曜日

de novo Transcriptome; 454用のベストなアセンブラーはどれだ !?　【論文紹介】

Roche 454 と言えば、パイロシーケンスでロングリードを読める代表格です。
数百塩基も読めるロングリードは、未知の配列決定にも良く使われます。
de novo のシーケンスですね。
以前、de novo Transcriptome を話題にしましたが、リファレンスが未知の生物で、転写産物配列を決定したい場合、ロングリードの454が、真っ先に使うシーケンサーの候補に挙げられるでしょう。
もちろんショートリードのペアエンドで読む、ということも可能です。
しかし、長く読めるということは、未知の配列決定において非常に大きなアドバンテージですので、454を選ぶ方は多いのです。

ゲノムではなく、転写産物を読む場合、Isoformの存在が気になります。
スプライシングのバリアントを見たい場合、やっぱりロングリードが有利でしょうか。

実は、リファレンス未知の生物のTranscriptomeを行っている研究は、結構あります。
論文になっているものだけでも、昨年紹介した、ほかにも、

Coral larval (サンゴ): Meyer et al. BMC Genomics 10, 219 (2009).
Eucalyptus grandis (グランディスユーカリ): Novaes et al. BMC Genomics 9, 312 (2008).
Sarcophaga crassipalpis (ニクバエ) : Hahn et al. BMC Genomics 10, 234 (2009).
Populus trichocarpa　(ブラックコットンウツド): Geraldes et al. Mol.Ecol.Resour 11, 81 (2011).

なんかがあるようです。　（全部は読んでいませんが）

454のリードは一般的には、シーケンサーの機械に付属しているアセンブラー、Newblerでアセンブルすることが多いと思います。
でも、Newbler以外にも、ロングリードをアセンブルできるアセンブラーはいくつかあるんですよ。
そこで、454のリードをアセンブルできるアセンブラーを比較した論文を紹介します。

Kumar et al. Comparing de novo assemblers for 454 transcriptome data. BMC Genomics 11, 571 (2010).

この論文は2010年に出されていますから、今年は各ソフトもバージョンアップして、論文投稿時とは若干状況が違っていると思いますが、参考になると思います。

比較しているアセンブラーは、

Newbler 2.3
Newbler 2.5
CAP3
CLC Assembly Cell 3.0
MIRA 3.0
SeqMan NGen 2.1

の6つです。　CLCとSeqManは商用ソフト、Newblerはアカデミックはフリー、CAP3は非営利ならフリー、MIRAはどこでもフリー、だそうです。　
結論からいうと、Newbler 2.5が最も良いとのこと。
「454のシーケンスメーカーが作っているんだから当たり前だろ！」
私も思わず突っ込みましたが、ほかにも面白いことが書いてあるかと。

アルゴリズムの違いとして重要なところは、
CLCはde Bruijn graphを、その他の5つはOverlap-Layout-Consensus (OLC) を使っていること。

de Bruijn graphはVelvetやABySSといったアセンブラーでも使われています。
de Bruijn graphの特徴は、リードをk-merという決められた長さの塩基ブロックに切り、この部分の重なりをもとに、アセンブルをします。
k-mer = 31なら31塩基の重なりを見て、リード同士をつなげていきます。
k-merの範囲でのみ重なり具合を見るんですね。　本当はショートリードの、たくさんカバレージがあるデータに向いています。
一方OLCはクラシカルな方法？で、ペアワイズアライメントを元にしています。
より慎重に伸長していくのですね。
もちろん、6つのアセンブラーはそれぞれ、計算アルゴリズムを工夫しているので、結果はちがうのですが。

この論文では、彼らは線虫の一種をサンプルに、Roche 454 FLX を使って de novo Transcriptome をやっていました。
アダプタートリムした　741,387本のリード、約2億塩基のデータを、先の6つのアセンブラーでアセンブルしています。

最も速くアセンブルが終わったのは、CLCで4分、次がNewbler 2.5の45分。　一番遅いのがMIRAの3日。
1kb以上の長さのContigが一番多くできたのは、Newbler 2.5で7,661本、次がNewbler 2.3の6,320本。　一番少ないのはCLCの4,174本。
全Contigの合計塩基数は、一番多いのがSeqManの2,136万塩基、次がMIRAの2,134万塩基、3位はNewbler2.5の2,007万塩基。　最も少ないのはNewbler 2.5の1,446万塩基。
（詳細は同論文のTable4を参照）

Contigの合計塩基数でNewblerのバージョン2.5と2.3で成績が全然違うのにはびっくりします。
2,007万塩基と1,446万塩基ですからねえ。

6つのアセンブラーでできたContigをそれぞれBLATして、相似性・非相似性を調べ、新規配列がどれだけ作られたか、を彼らは次に見ています。
結果、MIRAとNewbler2.5は、より多くの（長い）ContigをCLCやCAP3より作ったが、その余分の配列は特に新規の配列ではなかった。　曰く、MIRAとNewbler2.5のContigには冗長性があるのではないか、と。
CLCのcontigは最も冗長性が低かったらしい。　ということは余分な配列は作らなかった。　
「待てよ、それって良いことなのか？　TranscriptomeではIsoformはつきもの。」

繰り返しですが、
de Burijn graph は高カバレージの大量リードアセンブルに向いています。　なぜならアセンブルがめっちゃ速いから。
OLCは低カバレージのアセンブルに向いています。　なぜなら正確だけど、ペアワイズアライメントなので計算量がかかるから。

私の経験でも、CLCの de Burijn graph アセンブラーは他と比べて速いです。
でも、速いがゆえに、犠牲にしていることがあります。

アセンブルに使われたリードの情報（どのリードがどのContig形成に用いられたか等の情報）は持たない。　他のアセンブラーではACEファイルというものを持つようです。
k-mer は最大でも31です。　これはCLCの設定です。　31塩基は、ショートリードには向いているかもしれませんが、400塩基を超すようなロングリードでは短いですし、結果、Contigが断片化してしまうことがあります。
Isoformの認識は失われる。　CLCのde Burijn graphの宿命です。　Contigが途中から枝分かれするようなとき、「枝」は別のContigとされます。　これを修正しているアルゴリズムは別途紹介します。

しかし、4分という速さは驚くべきことで、「とりあえずやってみよう」的な使い方には向いているでしょう。

現在、MiSeqやPGMに代表されるように、より多くの、ロングリードが出てくるようになってきています。　OLCアルゴリズムでは今後、計算速度が大きな壁になるでしょう。
ちなみにこの論文では、64bit Linuxの、普通のワークステーションレベルのPCを使っていました。

さて、同じ de Bruijn graph でも、Isoformを見つけよう！　パラログを認識しよう！というコンセプトのアセンブラーが最近出てきました。
454のようなロングリードではなく、ショートリード用ですが、Trinityというものです。

http://trinityrnaseq.sourceforge.net/

最近使い始めたばかりですので、このツールについてはまたの機会に。

2011年7月6日水曜日

MiSeqのデータセット

イルミナ社も、いよいよデスクトップ型シーケンサーをリリースするそうです。
それに先立ち、早くもデモデータセットが公開されています。　
http://www.illumina.com/systems/miseq/ecoli.ilmn

このページの、右上あたりにある、Download the data　から、FastqファイルとBAMファイルが落とせます。　E.coliのデータです。
サイズが大きく、ダウンロードに結構時間がかかるので、wget で取得すると良いと思います。

（wget http://.........ファイルのURL............）

ついでに、というかこっちの方が面白いかもしれませんが、BAMファイルの下の、プレゼンテーションファイルも落として見てみましょう。
びっくりしますが、スライドの後半部分に、MiSeqデータとIon Torrent PGMデータの比較、があるのです！

ここから先は実際にダウンロードしてからのお楽しみ？
当たり前ですが、MiSeqの方がいい！　と結論しています。

カバレージで言うと、1ランで1.7Gb読めるMiSeqは、393x
一方PGMは6~8ラン読んでも11~24Mbで、2~5x

クオリティはどうでしょう？
プレゼンにも書いてあるのですが、せっかくReadファイルを落とせるので、fastqファイルを取得し、FastQCにかけてみました。
ここから先は、イルミナさんのプレゼンではなく、私が行った内容です。
ちなみにPaired-Endです。以下、R1とR2の順に示しています。

Quality Score

GC Content per Sequence

その他、リード長は151、平均クオリティはQ35/Q34、GCはどちらも50%でした。

私が気になったのは、クオリティのばらつきです。
150塩基まで読めるとはいっても、リードによっては、後ろの方は結構Qが低い塩基も含まれるようですね。
とは言っても100塩基あたりまで、平均してQ30あったのは驚きました。

実際、このデータでアセンブルをすると、どれくらいつながるのでしょうか。
先のプレゼンの中に、じつは答えはあるのですが。

2011年7月2日土曜日

Ion Torrent データセット(2)　クオリティチェック

先日、Ion Torrent PGM のデータセットが、Ion Torrent 社のホームページまたはBGIのFTPからダウンロードして誰でも使えるということを紹介しました。　Ion Torrent データセット
そのデータとは、E.coli O104-H4 の de novo sequence です。　
このことを書いた後、ある方と、クオリティの話をしました。

ということで、本日はクオリティチェックの結果を書きます。
使ったファイルは、Ion Torrent 社のサイトから取得した、LB226692株のfastq ファイルの、64.fastq ファイルです。　本当は64.fastq - 71.fastq まで、8ファイル全部やりましたが、ここに示すのは64.fastq ファイルのみということで。

使ったツールは、FastQC と、PRINSEQ　です。
両方とも同じようなツールです。
FastQCはPCの上で動くのに対し、PRINSEQはクラウド上で動くのが違いです。
ちなみにPRINSEQは、Internet Explore 9 ではまだうまく動かないみたいです。
FireFox、Google ChromeではOKでした。

FastQCの方が有名？