実験室より一言申し上げます

初心者からのNGS解析 ~悪戦苦闘日記６~

さて、今回はTopHat ~ Cufflinks準備の部分での修正Ver.となります。

前回実際に解析を行ったのですが、少々トラブルがありました。
今回はそのトラブルシューティングと実際にやったことについてです。

まず何があったのかというと、
TopHatで染色体21本を選んだにも関わらず、５本しかデータが参照されなかった
ということです。

Cuffdiffまで実際にやってみて、結果を見ながらおかしな所を探し、一個ずつ遡ってデータを見ていったところ、最初のTopHatでやらかしていたという涙目な結果だったのです。

この解析では
１．DDBJ pipelineのTopHatでGenBank IDを入力し、各染色体のデータセットをインポート
２．全部インポートしたらまとめてひとつのデータセットとして保存
３．これをリファレンスに指定してTopHatを実行
としたわけです。

が、どうもインポートしたゲノムデータセットがちゃんとしていなくて、参照された５本以外がおかしくなっていたようです。
DDBJの担当者によると、染色体をまとめて指定するとこういうことが起きるそうです。
以前に一括インポートを試みたことがあって、ダメなことはわかっていたので今回は全部個別にインポートしたのですが、それでもこういうことがあるんですね。。。
現在はこの件を踏まえ、デフォルトのMajor genome setにラットの最新のアセンブル (Rn6.0) が標準で選択できるようにしてくれたので、以後はこれを用いると安全でしょう。

そして、染色体名の書き換えです。

こちらは悪戦苦闘日記３でも書きましたが、上記の方法でTopHatを行った際も必要です。
(作業行程はこっちを参照してください)

標準ゲノムセットでTopHatを実行し、生成されたsamファイルをみてみると、例えば１番染色体が
gi|6619202784|gb|CM000072.5|
となっています。

一方、UCSCからインポートしたRefFlatファイルはchr1と書かれています。

なので、まずはUSCSからインポートした各染色体のRefFlatファイルを[Concatenate]で結合し、ひとまとめにします。→ これをファイル１とします。
(ゲノム全体をまとめてとってきてもY染色体の分ははいっていないので、Y染色体を別途インポートし、concatenateで結合させましょう)

このファイルをPCにダウンロードし、excel等を用いて文字列 [chr1] を [gi|6619202784|gb|CM000072.5|] に置換していきましょう。

置換したらtxt形式で保存し、再度Galaxyにアップロードします。→ ファイル２
で、ファイル２にファイル１を[Paste]します。
すると、ファイル１のデータの横にファイル２のデータがつながったファイルができあがります。→ ファイル３
上記の一連の操作の際に余計なタブ区切り認識で9列目がおかしくなっているので、[Cut]を用いてファイル３の[C1-C8,C18]を選択すると、まともな形式のGTFファイルになっているはずです。

あとは前回の様にCufflinksを実行してください。

初心者からのNGS解析 ~悪戦苦闘日記５~

さて、いよいよCuffdiffです。
Cuffdiffではサンプル間の発現量比較を行います。

そこで、またまた準備です。

・リファレンスデータを作成する
何回もやったじゃねーか！！と思いたくなるのをぐっと我慢です。
今回用意するリファレンスデータは、「Cufflinksで各サンプル毎に出てきたtranscriptデータを、比較できるように揃える」作業です。

Cufflinksはサンプル毎に個別に実行しているので、出力された結果では各転写産物がサンプル毎に名前が付けられています。
発現量を比較する際には、同じ転写産物同士で比較する必要があるので、名前を揃えなきゃいけないわけです。

そこで、「Cuffcompare」あるいは「Cuffmerge」を使って、統合したリファレンスデータを作成するわけです。

[Cuffcompare]と[Cuffmerge]の違いについては原著論文の筆頭著者であるCole Trapnellによる解説があり、彼はCuffmergeを推奨しているようです。

"Cole Trapnellによる説明"
=======================================================================
both are designed to gracefully merge full-length and partial transcript assemblies without ever merging transfrags that disagree on splicing structure. Consider two transfrags, A and B, each with a couple exons. If A and B overlap, and they don't disagree on splicing structure, we can (and according to Cufflinks' assembly philosophy, we should) merge them. The difference between Cuffcompare and Cuffmerge is that Cuffcompare will only merge them if A is "contained" in B, or vice versa. That is, only if one of the transfrags is essentially redundant. Otherwise, they both get included. Cuffmerge on the other hand, will merge them if they overlap, and agree on splicing, and are in the same orientiation. As turnersd noted, this is done by converting the transfrags into SAM alignments and running Cufflinks on them.

The other thing that distinguishes these two options is how they deal with a reference annotation. You can read on our website how the Cufflinks Reference Annotation Based Transcript assembler (RABT) works. Cuffcompare doesn't do any RABT assembly, it just includes the reference annotation in the combined.gtf and discards partial transfrags that are contained and compatible with the reference. Cuffmerge actually runs RABT when you provide a reference, and this happens during the step where transfrags are converted into SAM alignments and assembled. We do this to improve quantification accuracy and reduce errors downstream. I should also say that Cuffmerge runs cuffcompare in order annotate the merged assembly with certain helpful features for use later on.

So we recommend #3 for a number of reasons, because it is the closest in spirit to #1 while still being reasonably fast. For reasons that I don't want to get into here (pretty arcane details about the Cufflinks assembler) I also feel that option #3 is actually the most accurate in most experimental settings.
=======================================================================

実際にはどちらかを使えば良いわけですが、これを読む限りCuffmergeの方が良さそうですね。
(そもそもCufflinksが良い手法なのかどうかについては議論があるようですが...)
というわけでCuffmergeを実行します。

Cufflinksで出力されたファイルを選び、Cufflinksで用いたGTFファイルをリファレンスにして実行すればOK.
すぐに[merged transcripts]というファイルが生成されます。

・Cuffdiffを実行する
いよいよCuffdiffです。
左のタブからCuffdiffを選択してやるだけです。

サンプル名を入力し、sorted samファイルを選びましょう。
replicateとなっているものは[add replicate]から追加しましょう。

最後にパラメーターの設定です。
normalizationには[Geometric]、[Quartile]、[classic-fpkm]の３種類があります。
Geometricは中央値に対して、Quartileは上位25%をはじいたものの平均値に対して、classic-fpkmは生のfpkmを出力します。
GeometricもQuartileもデータをRobustnessをあげるためなので、どちらを使っても良いでしょうが、生は怖いですね。
私はGeometric (中央値) を使っています。

Dispersion estimation methodはreplicateのデータの合わせ方です。
調べたのですが細かいことは分からなかったです。
[pooled]は平均分散モデルというものを用いたもので、これで良さそうです。
もしreplicateがない場合は[blind]を選びましょう。

あとは[multi-read correction]と[Bias correction]です。
どちらもやった方が定量性が増すようですが、私はやっていません。
*multi-read correctionをやると結果が思わしくなかったです (因果はまだ不明)。現在詳細について解析中。

そしてExecuteをポチれば終わり！！！

しばらく (1日くらい) 待てば結果が出てくるはずです。

出てきたデータは次回、さばきます。

初心者からのNGS解析 ~悪戦苦闘日記４~

第４回。いよいよCufflinksです。

前回までで準備万端整いました。

まずは左のタブ[NGS: RNA Analysis]から[Cufflinks]を選びましょう。

そして、事前に準備しておいた[Sorted sam file]を選択します。

あとはオプションを選んでいくだけです。
以下は私の場合です。デフォルトからの変更点だけ記載します。

イントロンについてはデフォルトで200kbとなっています。
論文を調べていると、どうやら最大で50kbくらいのものはありそうですが、それ以上はないようなので私は[50,000]としています。
あまり設定が長いと変な物までひっつけてしまうかもしれませんからね。

次にquartile Normalizationを[yes]にしてます。
Cufflinksでは発現量をfpkmで表しており、簡単に言うと発現量を [100万リードあたりいくら分？]という形で表すわけです。
これはサンプル間で比較するに当たってとてもわかりやすいのですが、欠点として「発現量の多い遺伝子が増減すると、全体に影響を与える」というものです。
極端な例を挙げると、興味のある遺伝子Aならびに殆どの遺伝子の発現量がサンプル間で変わっていないとする。しかし、元々発現量が100万リード中10万リード分だった遺伝子Bが５倍の50万リード分になった。すると、Bの変化のせいでAを含む他の遺伝子の発現量が相対的に小さくなってしまう。
というように、高発現遺伝子の変動が低発現遺伝子の変動に影響を与えてしまう、という問題があります。
これを回避するために、総リード数を求める際に[上位25%分はなしにしよう]というのがこの正規化法です。
私的に、この正規化はやらないよりやった方がいいだろう、という判断で実行しています。

最後にExecuteをポチっとやってやればCufflinksが動きます。
しばらく待っていると
・gene expression
・transcript expression
・assembled transcript
という3つのファイルが生成されます。

あとはサンプル分、この作業を繰り返します。

次回は発現データをみるためのCuffdiffと、その準備のためのCuffcompare, Cuffmergeです。

初心者からのNGS解析 ~悪戦苦闘日記３~

第３回目となります。

前回まででCufflinksに用いるSAMファイルの準備までやりました。
次はCufflinksです。

実際にCufflinksを実行する前にもう少し準備が必要になります。
それは「リファレンスゲノムデータ」の準備です。

あれ？それはTopHatの時にやったのでは？と思われるかも知れませんが、残念ながらCufflinksではもう一度準備が必要です。

ぶっちゃけ面倒です。
なぜかというと、今回もデータをインポートして使うのですが、そのデータがそのままでは使えないからです。

以下、手順を追って説明します。

1. USCSからゲノムデータを取得する
まずはページ左のタブの[Get Data]より[USCS Main table browser]を選んで下さい。
(ページが飛ばないときは右クリックから新規タブで開くとOKです)

するとUSCSのページが表示されます。

そしたら上記の例のように[生物種]、[Assemblyの種類]を選び、Groupを[Gene and Gene Predictions]、trackを[RefSeqGenes]、tableを[refFlat]とします。
特に大事なのがtableの種類で、refFlatを用いないと悲しいことになります。
用いるデータテーブルの種類を決定したら、output formatで[GTF]を選択して下さい。
あとは[get output]を押すことでGalaxyにimportされます。

2. GTFファイルの染色体名を変更する

DDBJ pipelineに用意されているゲノムセットを用いた場合はこちら

このままCufflinks！といきたいところですが、やってしまうと悲しいことになります。
Cufflinksでは[TopHat]のデータ (SAMファイル) とリファレンスゲノムデータの「染色体名」が一致していないと別のものとして判定されてしまい、データが染色体にひも付けされないままになってしまいます。

TopHatで用いた染色体名 (悪戦苦闘日記２を参照) は、例えば１番染色体に関しては[CM000072.5]となっているのに対し、USCSからインポートしたものは[chr1]となっています。
これじゃいけません。

というわけで、書き換えます。
やり方は自由ですが、私はExcelを使いました。
インポートしたUSCSのrefFlatファイルを自分のPCにダウンロードし、Excelで開きます。
次に、[chr1]を[CM000072|CM000072.5]に置換します。
(検索文字列は完全一致を用いて下さい。chr11とかも置換されてしまいますので)
これを染色体の数分繰り返せばOK。
2~3分で終わると思います。

しかし、油断してはいけません。
私がTopHatに用いたリファレンスゲノムはchromosome 1 ~ 20, X, Y, mitochondriaなのですが、このUSCSのデータは 1 ~ 20とXしかありません。

なので、先ほどと同様にUSCSからY染色体とmitochondriaのデータを取得してきて、同じ作業をします。
Y染色体はpositionを[chrY]、ミトコンドリアは[chrM]として[lookup]ボタンを押せばOKです。

しかし、さらなるトラップがあります。
YについてはrefFlatをそのままとってくればいいですが、ミトコンについてはtrackを[other RefSeq]、tableを[xenoRefFlat]を選んでください。
そのままrefFlatをとってくると、中身がからっぽになっています。
(YについてもxenoRefSeqを使うのもありかも知れません)
*xenoRefFlat: RefFlatは対象生物の遺伝子名だけですが、xenoRefFlatは他生物のオルソログも参照しています。

全染色体の染色体名を書き換えたらファイルをひとまとめにします。
最初の編集したRefFlatの行下にY, ミトコンのものをコピペで貼り付けます。
そしたらこのexcelファイルを[テキスト(txt)]形式で保存します。

次にこのファイルをGalaxyの[Get Data]の[Upload file from your computer]でアップロードします。
ファイルを選んで、形式を[GTF]にしてください。
自動でtxtファイルがGTF形式に変換されます。

3. GTFファイルのフォーマットを整える
そして、最後のトラップが待ち受けています。
変換されたファイルをみてみると、一番右 (9列目) のAttributionが変なタブ認識されて、過剰に["]が追加されています。
まずいです。

タブの設定を色々やったのですが、うまくいかなかったのでごり押しで直します。

まず、UCSCからGalaxyにインポートした元のrefFlatファイル (1~20 +X, Y, Mの３つ) がありますよね。
これを左タブの[Text manipulation]から[Concatenate datasets tail-to-head]を選んで、excelでやったのと同様に全部ひとまとめにします。
このファイルはAttibutionは正常なのですが、染色体名が書き換える前のものになっています。

次に、先ほどアップロードしたファイル (染色体名を書き換えたもの) と上記のファイルを[Paste two files side by side]で合体させてください。すると９列＋９列 = 18列のファイルになります。
そして、[Cut columns from a table]でC1-C8, C18を選びます。
つまり、「狂った９列目のAttibutionを元のRefFlatファイルの9列目に置き換える」わけです。
最後に、この編集したファイルは[Tabular形式]になっているので、ファイルの鉛筆印 (Edit Attributes) から[Data type]を選んで、GTFファイルに変換すれば完成です。

かなりぐだったやり方ですが、要は以下のポイントができていれば他のやり方でもOKです。
 1. Y, ミトコンドリアを含むTopHatに用いた全ての染色体のRefFlatファイルを用いる
 2. 染色体名がTopHatの時のものと同じ
 3. GTFファイルの表示がいじくる前と同じ


これで準備が整ったので次回はCufflinksです。

初心者からのNGS解析 ~悪戦苦闘日記２~

素人でもできるNGS解析！というわけで続編です。(１年以上空いてるけど、しれっと書きます)

前回のところではFASTQファイルのプレプロセシングまでやりました。
今回はここからTopHatによるマッピングです。

私が利用しているのは「DDBJ read annotation pipeline (以下DDBJ pipeline)」で、ここには各種演算のツールがおいてあります。

なので、ページ左側のタブから「Step.1 Mapping」をぽちっとして、方法で「TopHat」を選びます。
その後、事前にプレプロセシングしたデータを選択して、あとはページを順番に進めていくだけです。

ここで大事なのが「リファレンスゲノム」のデータです。
DDBJ pipelineではすでに「主要生物のゲノムデータ」がサーバーにおいてあり、選択肢からポチッとするだけでTopHatを動かすことができます。
しかし、私が解析しているラットに関しては昨年までありませんでした。
(今はおいてありますが、アノテーションが１バージョン前のものです)

というわけで、データを用意する必要があります。
TopHatでリファレンスゲノムを選ぶページにくると、「Download or upload reference」という項目があるので、ここをポチります。
そして、「使いたいゲノムデータのGenBank ID」を入力してください。
するとゲノムデータがとれます。
私の場合はラットの最新のもの (Rnor_6.0) が使いたかったので、CM000072.5からCM000092.5 (chr1 ~ chrX) とCM002824.1 (chrY)、AY172581.1 (chrM) をダウンロードしました。
(NCBIのAssemblyのデータベースにまとめてあるので、そこを参照)
染色体5本分くらいはまとめてダウンロードできますが、エラーがでることがあるので一つずつやった方がよいでしょう。
ダウンロードした各染色体ゲノムのデータは[Create Data set]でひとまとめにして、わかりやすい名前をつけておきましょう。

あとはそのままTopHatを実行するだけです。
するとSAMファイルが生成されます。(SAMファイルは次のCufflinksで使います)
ちなみに、FASTQファイルが分割されている場合は、まとめてファイルを選んでMappingすることでひとつのSAMファイルにすることができます。
ファイルの大きさやサンプル数にもよりますが、数時間 ~ 1日もあれば全部終わるでしょう。


次は「Cufflinks」になるわけですが、まずはページの左側[step-2]の[Workflow]をポチってください。
するとGalaxy/P-GALAXYというページに飛びます。

今後の作業はここでやることになるので、先ほど生成されたSAMファイルをこちらに移します。
左側タブの[Workflow]の中に[Import samfile  from DDBJ Pipeline]というのがあるので、こちらをポチってください。
するとSAMファイルがGalaxyの方に移せます。

最後に、このSAMファイルを「sorting」します。
Cufflinksはアルゴリズムの都合上、sortingしたSAMファイルでないと扱えないようなので、忘れずにやりましょう。

やり方は簡単です。
同じく[Workflow]のタブの中に[Cufflinks preprocessing]というのがあるので、これを選んで先ほどimportしたSAMファイルを選んで実行してください。
すると、[Sorted SAM file] というのができあがるはずです。

ここまででSAMファイルの準備はおしまいです。

次回 (連続して投稿しますが) はいよいよCufflinksです。