htseq-count 使い方 gene単位・transcript単位の数え方

バイオインフォ道場、くまぞうです。

RNA-Seq解析では、ゲノムに張りついたリードの数を数えます。gene単位・transcript単位の数え方は、オプションで切り替えることができます。

スポンサーリンク



関連記事 HTSeq htseq-count 使い方 インストール v0.7.2
関連記事 HTSeq htseq-count 使い方 よく使うオプション

htseq-count 使い方 gene単位・transcript単位のカウント方法

htseq-countのデフォルト動作はgene単位のカウントです。exon領域に対してgene_id毎に集計します。オプションの-i,--idattrを明示的に指定してGTFファイルの参照情報を変更することで、カウント単位をgene単位やtranscript単位に変更することができます。

GTFファイルの9列目(attribute)の情報で切り替える

GTFファイルは、ゲノムや遺伝子の「位置・遺伝子名・エクソン番号など」の情報を1行毎に記述したファイルです。1行には9つの列があり、タブで区切られています。GTFの9列目には、type/valueで対になった情報(必須:gene_id, transcript_id)を記載します。

htseq-countの-i,--idattrオプションで、GTFの9列目のgene_idまたはtranscript_idを指定します(定義があれば他のIDも指定可能です)。このオプションで指定した単位について、3列目のfeature領域のリードをカウントします。これらは連動するので、どちらかが正しく指定されていない(または定義がない)とエラーが発生します。

1.seqname   : 例 Chr1, 1 ...
2.source    : 例 hg19(生成プログラム名など)
3.feature   : 例 gene, CDS, exon, ...
4.start     : 例 1, 1000 (開始位置。先頭は1)
5.end       : 例 100,999 (終了位置)
6.score     : 例 '.' (任意のスコア) 
7.strand    : 例 +, -, '.'
8.frame     : 例 0, 1, 2, '.' (翻訳開始塩基位置。0はコドン1番目)
g.attribute : 例 遺伝子名, gene_id="xxx", transcript_id="xxx"

htseq-count 使い方 gene単位・transcript単位のオプション

機能 説明
-i, –idattr リード数をまとめる単位を変更する場合に指定。デフォルトはgene単位。gtfファイルの9列目で定義されているIDを指定。このオプションを変更することで、カウントの単位をgene単位やtranscript単位で切り替える。IDが定義されていれば他のIDも指定可能。
-t, –type どの領域(exon, mRNA, gene…)で数えるかを指定。デフォルトはexon。gtfファイルの3列目で定義されている情報を指定。カウント領域を変更したい場合に変更。idattrと連動し、idattrで参照する行に指定したの定義がないとエラーになる。

gene単位・exon領域の指定

デフォルトなので、通常は特に指定しなくてもgene単位・exon領域になります。(bam・samtoolsでソート・non-specific-strandの場合)

htseq-count -f bam -r pos -t exon -i gene_id accepted_hits.bam ref.gtf > count_gene.txt
htseq-count -f bam -r pos accepted_hits.bam ref.gtf > count_gene.txt #同じ設定(デフォルト)

transcript単位・exon領域の指定

transcript単位・exon領域をカウントします。(bam・samtoolsでソート・non-specific-strandの場合)

htseq-count -f bam -r pos -t exon -i transcript_id accepted_hits.bam ref.gtf > count_gene.txt
スポンサーリンク