もしかしたら、常識かもしれないので、 情報お持ちのかたがおられたら、お願いします。 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
選択肢2
一応、BioStarやSeqAnswersで上記難点の解決策も探しては見たのですが、 自分では見つけられませんでした。 coverageBEDで、イントロンも含めてカウントした後に、 全イントロンのタグ数を引くとか、 小賢しい手も考えられますが、 境界のリードなどの扱いで、厳密では無くなる点に難があります。 もし、うまい手やら、既存のスクリプトでこんなものがある等、 関係する情報お知らせいただけますと幸いです。 直接の解答でなくとも、ご助言ありましたら、お願いします。 |
自己レスですが、現在までに分かったことを記載します。
もし、間違っていたり、関連情報ありましたらコメント等引き続きお願いします。 |