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ファイルの９列目（attribute）の情報で切り替える

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

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

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

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

機能 説明
-i, –idattr リード数をまとめる単位を変更する場合に指定。デフォルトはgene単位。gtfファイルの９列目で定義されているIDを指定。このオプションを変更することで、カウントの単位をgene単位やtranscript単位で切り替える。IDが定義されていれば他のIDも指定可能。
-t, –type どの領域（exon, mRNA, gene…）で数えるかを指定。デフォルトはexon。gtfファイルの３列目で定義されている情報を指定。カウント領域を変更したい場合に変更。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
スポンサーリンク