Answers to: RNA-SeqのrRNAの除去

Comment by nob_fj on nob_fj's 回答

nob_fj — Wed, 29 Dec 2010 14:13:18 +0900

私自身が確認したわけではないですが、cufflinksを使用したフローでもTopHatが出力するSAMを上記処理でフィルタリング可能という話を聞きましたので、質問は一旦閉じさせていただきます。もし間違っていたら、質問の再開をお願いします。

Comment by nob_fj on nob_fj's 回答

nob_fj — Wed, 08 Dec 2010 20:25:19 +0900

hacchyさんgaouさんからいただいたヒントを元に、 BAMファイル中のrRNA、tRNAを除去する方法が分かりました。以下のwikiにまとめてあります。

rRNA,tRNA除去forBAMフロー

rRNA,tRNAの座標情報をBEDフォーマットでUCSC Table Browserから取得する

Comment by nob_fj on gaou's 回答

nob_fj — Wed, 08 Dec 2010 16:44:30 +0900

教えていただいて恐縮ですが、現在cufflinksはまだ導入できておらず、 cufflinksベースの解析フローが出来たら再度gaouさんのフローを試してみる予定です。対応には暫くお時間いただきたいと思います。
現在のBioConductorのDEGseqを使用した既知遺伝子アノテーションをベースとしたスプライスジャンクションを考慮しない解析フローの場合は、アノテーション済みのUCSCから取得したリピートの座標情報を使用して BEDtoolsを使用するやり方の方が相性は良さそうです。

Comment by nob_fj on gaou's 回答

nob_fj — Wed, 08 Dec 2010 11:30:42 +0900

UCSCのアノテーションがない生物種でも、RepeatMaskerベースの方法論であれば、ゲノムが決定しており、マッピング出来さえすれば、その情報を使用して、
1.genomeからマップされた箇所の配列切り出し
2.切り出した配列へ、Repbaseでリピートマスク
3.ヒットしたリピートがrRNA,tRNAのもののみ、除外
というステップがとれるという理解でよいでしょうか。
今回の対象が、ヒトなので、まずはUCSCアノテーションに依存した配列ベースの情報を介さない座標ベースのやり方でまずは試してみたいと考えております。

Comment by gaou on gaou's 回答

gaou — Wed, 08 Dec 2010 02:33:34 +0900

基本的にはそうだと思います。汎用的に使うためにはRepeatMasker自前でかけた方が楽かとは思いますが。

Comment by nob_fj on gaou's 回答

nob_fj — Tue, 07 Dec 2010 22:55:54 +0900

$ awk '{print $12}' allFieldsFromSelectedTable.txt | sort | uniq -c  
 461751 DNA  
   1881 DNA?  
1498690 LINE  
     51 LINE?  
 717656 LTR  
    122 LTR?  
 371543 Low_complexity  
   3733 Other  
   2236 RC  
    729 RNA  
1793723 SINE  
    425 SINE?  
   9566 Satellite  
 417913 Simple_repeat  
   7036 Unknown  
     97 Unknown?  
   1769 rRNA  
      1 repClass  
   1340 scRNA  
   4386 snRNA  
   1481 srpRNA  
   2002 tRNA

Comment by nob_fj on gaou's 回答

nob_fj — Tue, 07 Dec 2010 22:55:38 +0900

gaouさんもしご存知でしたら Repbaseから直接取得したわけではないのですが、 UCSCのTable browserで "Variation and Repeats"のカテゴリに"RepeatMasker"というトラックがあり、このカラムの "repClass"のカラムが以下のような構成になっており、このうちtRNA,rRNAのもののみ抜き出せば、 gaouさんのフィルタリング対象と同じになるのではないかと思いますが、違いますでしょうか。

Comment by gaou on gaou's 回答

gaou — Fri, 03 Dec 2010 17:50:27 +0900

クエリはcufflinksをかけたあと、gtfの座標をもとに簡単なスクリプトでtranscriptのfastaファイルを作って使っています。Repeatmaskerのリファレンスは普通にRepbaseで、いろんなものがヒットしますが、RNAに該当するものだけを除去する（こちらも簡単なスクリプト、というかgrep）ようにしています。

Comment by nob_fj on gaou's 回答

nob_fj — Fri, 03 Dec 2010 17:36:22 +0900

素人質問で申し訳ないですが、RepeatMaskerにかける際の、クエリとリファレンスはどうなるのでしょうか。
クエリ->ヒットのあった遺伝子の配列?
リファレンス->rRNA,tRNAの配列リスト?
でしょうか。

Comment by nob_fj on hacchy's 回答

nob_fj — Fri, 03 Dec 2010 17:31:48 +0900

申し訳ないですが、Ensembl APIですが、サーバ環境が直接外に出られないので、ちょっと敷居が高そうです。ローカルにデータをもって来てMySQL dbを構築する手もありますが、このためだけに行うにはちょっと作業負荷が高そうなので、他の方法を探したいと思います。

Comment by nob_fj on nob_fj's 回答

nob_fj — Fri, 03 Dec 2010 17:29:27 +0900

tRNAの座標に関しては、UCSCのTable browserでgroup="Genes and Gene Prediction Tracks"のtrack="tRNA Genes"からBEDフォーマットで入出できそうです。 rRNAがどこにあるかはまだ分かっていません。

Answer by gaou

gaou — Fri, 03 Dec 2010 12:17:56 +0900

私はお手軽にはとりあえずCufflinksかけてしまって、予測されたtranscriptに対してRepeatMaskerかけちゃったりしています。

Comment by nob_fj on nob_fj's 回答

nob_fj — Thu, 02 Dec 2010 22:41:22 +0900

どちらがやりやすいか検討してから対応を行ってみます。

Answer by nob_fj

nob_fj — Thu, 02 Dec 2010 22:40:31 +0900

良く調べたら、同じ質問がSEQAnswersにもありました。こちらの回答では、BEDtoolsを使用するようですが、あらかじめBED形式のrRNA,tRNA座標情報を用意しておく必要があるようです。

Comment by nob_fj on hacchy's 回答

nob_fj — Thu, 02 Dec 2010 22:32:41 +0900

生物種はヒトなので、Ensembl APIを試してみます。できたらやり方をwikiにアップしようと思います。

Answer by hacchy

hacchy — Thu, 02 Dec 2010 22:19:14 +0900

生物種はマウスでしょうか？一般的かどうかは分かりませんが，私は Ensembl API を使ってアノテーションの取得をしています。 Ensembl Genome Browser にある全てのアノテーションを柔軟に取得することができるので，tRNA の位置情報を取得したり，その配列を取得したり，mRNA-seq データから遺伝子ごとの RPKM を計算するときに必須なエキソンの位置を取得したり，という処理を短いコードで書けるので，愛用してます。

ヒト/マウスであればアノテーションはある程度しっかりしていると思いますが，Ensembl 生物種の中にも tRNA のアノテーションが足りなかったりする場合もあるようです。そのような場合は，tRNA-SE などを自前で動かして，ncRNA 領域の予測を行うと良いそうです。また，Ensembl にないような生物種は，GenBank ファイルか GFF ファイルがある場合が多く，BioPerl ベースのパーサを作って情報を取得しています。