ログイン 概要 よくある質問

リード長の長いRNA-Seqのタグカウント

0

もしかしたら、常識かもしれないので、 情報お持ちのかたがおられたら、お願いします。

topHatやmapSpliceなど、splicingを考慮したマッパーの結果から、 遺伝子レベルのタグカウントを計算するうまいやり方など、 情報お持ちのかたがおられましたら、 お知らせいただけますと助かります。

以下背景

topHat-cufflinks-cuffdiffや、それに類似のフローが デファクトスタンダード的解析であると知りつつ、 リード長の長いRNA-SeqにedgeR等をかけるやりかたを考えています。

edgeRなどのDEG検出用のRパッケージ等は、 fpkm/rpkmは使わないように推奨しているかと思います。

EdgeRのedgeRUsersGuide.pdfより引用 http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf "There are good reasons why the NB model is appropriate for the raw count data, but transforming the data using RPKM (or FPKM or similar) renders our distributional assumptions invalid and we cannot guarantee that our methods will be reliable for such transformed data."

35塩基や、50塩基程度であれば、 BWAやBowtieなどで、splicingを考慮せずに、 マップしたものを、 DEGseqパッケージ中のgetGeneExpや、 BEDtoolsのcoverageBEDなどで、タグカウントすればよかったのですが、 最近リード長が伸び、150塩基のRNA-Seqなども よく見かけるようになったかと思います。 リード長が伸びるとexon長を超えてしまう場合も 多々あるので、junctionをまたぐタグの数を無視するには、 支障が出てきているように思います。

topHatが出力するaccepted_hits.bamに対して、

選択肢1

  • Step1 BEDtoolsのbamToBedに"-bed12"オプション付でかけBAMをBEDに変換
  • Step2 BEDを入力とするDEGseqのgetGeneExpでタグカウント
  • アノテーション情報は、refFlat.txt
  • => イントロン領域のタグは計上されない。
  • => splicing-junctionをまたぐタグがタグカウントに計上されないため、難あり。

選択肢2

  • Step1 coverageBEDに対して、"-abam"でBAM入力で"-splice"オプションを付けて実行。
  • アノテーション情報はrefFlat.bed
  • => splicing-junctionをまたぐタグはタグカウントに計上されるので良し。
  • => イントロン領域のタグ(第一イントロンの制御ノンコーディングRNAなど)もカウントされるので難あり。

一応、BioStarやSeqAnswersで上記難点の解決策も探しては見たのですが、 自分では見つけられませんでした。

coverageBEDで、イントロンも含めてカウントした後に、 全イントロンのタグ数を引くとか、

小賢しい手も考えられますが、 境界のリードなどの扱いで、厳密では無くなる点に難があります。

もし、うまい手やら、既存のスクリプトでこんなものがある等、 関係する情報お知らせいただけますと幸いです。 直接の解答でなくとも、ご助言ありましたら、お願いします。

質問日 Apr 03 '12 at 16:59

nob_fj's gravatar image

nob_fj ♦
507●8●16●28

edited Apr 04 '12 at 01:08

3 Answers:
回答順最新支持されている順
0

こんばんは。
僕は、nob_fjさんが箇条書きなされている方法の中で、５番の方法を使っています。
Cufflinksを使えば解決すると思うのですが、間違っていたらご指摘していただければと思います。
なお、下記は、cuffdiffをつかうため、比較データAとBを仮定します。
また、オプションは、pair-end directionalの例を書いていますが、気にしないでください。

  1. tophatをつかってgenome mappingとjunction mappingをする
    eval tophat --library-type fr-secondstrand -o ./dataA -r 800 -p 10 -G $mouse_ref $mouse_index1 dataA.fastq_1 dataA.fastq_2
    eval tophat --library-type fr-secondstrand -o ./dataB -r 800 -p 10 -G $mouse_ref $mouse_index1 dataB.fastq_1 dataB.fastq_2

  2. フォルダ内にある transcripts.gtfファイル(転写産物ごとのexon intron情報のファイル)をcuffmergeを用いて、統一したjunctionのreference用gtfにする。
    cuffmerge -s $mm9genome assembly.txt
    (assembly.txt内は以下です。)
    [suimye@Brahms testdata]$ cat assembly.txt
    ./dataA/transcripts.gtf
    ./dataB/transcripts.gtf

cuffmergeの結果、merged.gtfというA、Bのデータに基づく新しいjunction情報がえられますので、
こいつにcountingしてやれば、intronデータを過大評価するようにはならないと思います。

  1. merged.gtfを用いて、merged.gtfに対するカウンティングをしてやる。
    cuffdiff -p 10 --library-type fr-secondstrand -o ./resCufdiff ./merged.gtf ./dataA/dataA.bam ./dataB/dataB.bam

すると、cuffdiffでは、FPKM値を算出しますが、FPKM値は無視して
read数だけの行列をcutコマンドで作ってやって、最後は、cuffcompareを
merged.gtfにかけたやつでannotationしてedgeRやれば大丈夫だと思いますが、
いかがでしょうか。間違ってたらごめんなさい。

回答日 Jul 13 '12 at 20:24

suimye's gravatar image

suimye
296●1●8●15

edited Jul 13 '12 at 20:37

お返事遅れました。以前来ていたスパム的投稿の影響でメーラがLSQAの通知メールを迷惑行きにしてしまっておりました。 また詳細な情報大変ありがとうございます。お知らせいただいたやりかた、試してみます。 cuffdiffのタグカウントのjunction readの比率や割り振り方が気になるので、HTSeqをカスタムしたやり方と比べて評価しようと思っています。(大分先になりそうですが)

(Jul 19 '12 at 21:52) nob_fj ♦ nob_fj's gravatar image
0

同bamをcoverageBEDで"-splice"オプションを付けても、 イントロン領域のタグも含めてカウントしてしまいます。

coverageBedの -b オプションには遺伝子アノテーションのファイルを指定しているのでしょうか？ bam を bamToBed -split でBEDに変換してからcoverageBedをかけたら、うまくいきそうな気がします。

それでもダメなら、イントロン領域にあるタグ（の部分）を intersectBedで削るとうまく行くかもしれません。 - mergeBed -i annotation.gff > annotation.merge.bed - bamToBed -split -i accepted_hits.bam > accepted_hits.bed - intersectBed -a accepted_hits.bed -b annotation.merge.bed > accepted_hits.exon.bed - accepted_hits.exon.bed にcoverageBedをかける

試してないので、うまく動くか分かりませんが、

回答日 Apr 03 '12 at 20:27

yuifu's gravatar image

yuifu
51●1●1●4

試してみました。 以下の条件全て結果は同一で、intron領域もカウントされてしまいました。 結果のカウントは全て同じで、lengthのカラムのみ、-splitオプションが有効な、 2と6の条件のみ、intronを除外した値が表示されていました。

  1. リードの入力BAM -splitなしcoverageBED
  2. リードの入力BAM -splitありcoverageBED
  3. リードの入力BED(-bed12オプションなしbamToBed変換後) -splitなしcoverageBED
  4. リードの入力BED(-bed12オプションなしbamToBed変換後) -splitありcoverageBED
  5. リードの入力BED(-bed12オプションありbamToBed変換後) -splitなしcoverageBED
  6. リードの入力BED(-bed12オプションありbamToBed変換後) -splitありcoverageBED

intersectBedを使う件についてですが、intron領域に別の遺伝子が含まれる場合等、根こそぎ削ってしまったりしまいませんでしょうか。アノテーションにも依るような気がしますが、refFlatなんかはnonCodingなども含まれていたような気がするので、削られてしまいそうな気はします。それらは別にやるのも手ですが。

(Apr 03 '12 at 22:47) nob_fj ♦ nob_fj's gravatar image

補足ですが、-bには、BED形式でダウンロードしたrefFlatが入っています。上記の、入力BAMと書いてあるのは"-abam"で渡し、入力BEDと書いてあるのが"-a"で渡したものです。

(Apr 03 '12 at 23:55) nob_fj ♦ nob_fj's gravatar image
0

自己レスですが、現在までに分かったことを記載します。

  1. coverageBEDはRNA-Seqのタグカウントには、極めて不向きです。 理由としては、exonと非転写領域、exonとintronの境界のリードをカウントしてしまいます。intron領域のリードもカウントする他、splicing junctionのリードもカウントします。最初にexon+intronをcoverageBEDで数えたあとに、intron領域に張り付いたリードを引くことを試しましたが、intron領域のタグカウント時に、splicing junctionをまたぐ単一リードが複数回カウントされる可能性がある上、この影響は、junctionの数や、exon長などいろいろな要素に影響を受ける為、このカウントのしかたは、厳密性に問題がありすぎます。
  2. getGeneExpについて 遺伝子単位での計測には、ほぼ問題ないものの、ゲノム上で重複する遺伝子(例えばAPITD1とAPITD1-CORT)の重複領域のタグを"保守的"に数えないので、解析の目的によっては支障があるかもしれません。無理やり、splicing valiantに別IDをふってカウントすると、variantのある遺伝子の殆どのタグが0になってしまったりします。
  3. DESeqのマニュアルで、お勧めされているHTSeqは、splicing junctionをまたぐリードを計上してくれます。しかし、exonをskipするjunctionなども、領域に内に含まれていれば計上されるので、cufflinksのような、splicing junctionに特化した厳密性はないことを考慮する必要がありそうです。また、ゲノム中で重複する遺伝子のカウントは、getGeneExp同様、カウントしてくれません。splicing variant用のカウントには不向きかと。
  4. cufflinks2でタグカウント値を出してくれます。まだ詳しく把握できていませんが、ばらつきの情報などもあるようです。これを使うのがやはり最も素直でしょうか。ただ、小数のjunctionをまたぐリードの比率で、比率が大きく変わってしまう点が、検定する側のプログラム(edgeRやらDEseqやら)の前提を満たしているのかが気になります。もし、どこそこにその辺の議論が書いてあるなどの情報ありましたら、お知らせいただけると幸いです。
  5. cufflinksなどで、予測されたトランスクリプトや、既知トランスクリプトの領域情報から配列を切り出し、bowtie等でリードをトランスクリプト配列に再マップしてカウントする方法を知人に示唆されました。この方法では、splicing-junctionにマップされずらいなどのバイアスを考えなくても良いので、再マップの実行時間(それほど大したことはないですが)を抜きにすれば一番シンプルかもしれません。トランスクリプトの情報に依存する点と、マッピング時のマルチヒットの扱いには留意が必要だとは思います。

もし、間違っていたり、関連情報ありましたらコメント等引き続きお願いします。

回答日 Jul 11 '12 at 21:01

nob_fj's gravatar image

nob_fj ♦
507●8●16●28
あなたの回答
プレビューをトグルする

この質問をフォローする

By Email:

Once you sign in you will be able to subscribe for any updates here

By RSS:

回答

回答とコメント

タグ:

×47
×13
×6
×2
×1

質問日: Apr 03 '12 at 16:59

閲覧数: 10,054 回

最終更新日: Jul 19 '12 at 21:52

関係した質問

Tophatのbutterflyオプションについて

strand specific RNA-Seqデータ解析におけるtophatの利用について

Trinityのアセンブリ結果の見方

RNA-SeqのrRNAの除去

TopHatのreferenceデータについて

BMB2010 のNGS関連の発表

FindPeaksの使い方

C-value からゲノムサイズを推定する

Celera WGA Assembler(CABOG)のインストールと使い方

NGSライブラリのPCR増幅

概要 | よくある質問 | プライバシーポリシー | お問い合わせ

powered by OSQA