Platform for Drug Discovery


HTSeq count


概要


  • pythonで記載されたタグカウントプログラム。
  • モードを切り替えることによって、tophat等のsplicing-junctionをまたいだリード適切にカウントすることができる。
  • 複数のアノテーションに重複してアノテートされるリードは、カウントしない。
  • アノテーション情報はGFF(GTFも可)。リードはsamを利用する。ペアリードの場合は、名前順ソートが必要。
  • デフォルトではGTFのgene_idをもとにグルーピングする。isoform単位(transcript_id単位)のタグ数は計上しない。

使用例


実行前に


  • pythonのパスを通しておきます。
    export PATH=/toolDir/HTSeq-0.5.3p9/bin:@CIP_HOME@/tool/samtools-0.1.18:${PATH};
    
    export PYTHONPATH=/toolDir/HTSeq-0.5.3p9/lib/python2.6/site-packages/:$PYTHONPATH;
    • 他のプログラムとpythonのバージョンが競合しなければ毎回実行する必要はないと思います。

sam入力の場合


  • htseq-count input.sam annotation.gtf -m intersection-strict
  • オプション指定した例
    htseq-count input.sam annotation.gtf -m intersection-nonempty -s yes
    • directional RNA-Seqの場合は"-s yes"。
    • "-m intersection-nonempty"はアノテーション情報の境界が不正確で、境界領域のリードを救済したい場合に付加します。

bam(染色体座標順ソート済み)入力の場合


  • TopHat bamの場合
    samtools view -b input.bam | samtools sort -on - tmporaryFilePrefix | samtools view - | htseq-count -q - annotation.gtf -m intersection-strict
    • samtools sort -onはペアエンドのbamでなければ必ずしも必要ありません。
    • もちろんsamtoolsだけ先に実行しておいても良いです。
    • "-s yes"と"-m intersection-strict"は、入力がbam
  • BWA paired-end bamの場合
    #Propter pairに限定するフィルタを実行する例
    
    samtools view -b input.bam | samtools sort -on - tmporaryFilePrefix | samtools view - -f 0x2 | htseq-count -q - annotation.gtf -m intersection-strict
    
    #ペアの片鎖あるいは両鎖がunmappedなものを除外する例
    
    samtools view -b input.bam | samtools sort -on - tmporaryFilePrefix | samtools view - -F 0xC | htseq-count -q - annotation.gtf -m intersection-strict
    • BWA由来のBAMの場合など、ペアのどちらかがマップされないものを除外するフィルタをかまさないと、エラーが発生するようです。
      • proper pairに絞ると、適切なinsertサイズ、適切な方向に両鎖がマップされたペアのみが計上されます。
      • ペアの片鎖、あるいは両鎖がunmappedなものを除外すると、"not_aligned"のタグ数が0になります。
    • 関連記事SeqAnswers(HTSeq sort sam)(external link)

実行結果の例と結果の検討


実行例と出力例


  • 実行コマンドの例
    $ /toolDir/HTSeq-0.5.3p9/bin/htseq-count randRna20M.name.sam hg19.refFlat.20130205.gtf > randRna20M.tc.refFkat.txt 2> randRna20M.tc.refFkat.err
  • 標準エラー出力の、進捗ログ出力
    $ less randRna20M.tc.refFkat.err
    
    100000 GFF lines processed.
    
    200000 GFF lines processed.
    
    :
    
    (中略)
    
    :
    
    800000 GFF lines processed.
    
    868195 GFF lines processed.
    
    Warning: Malformed SAM line: MRNM != '*' although flag bit &0x0008 set
    
    Warning: Malformed SAM line: RNAME != '*' although flag bit &0x0004 set
    
    200000 sam line pairs processed.
    
    400000 sam line pairs processed.
    
    :
    
    (中略)
    
    :
    
    19900000 sam line pairs processed.
    
    20000000 sam line pairs processed.
    
    20000000 sam line pairs processed.
    • コマンドで実行した場合、コマンドの進捗状況に応じて、GTFのエントリー読み込み数、SAMの処理件数が逐次出力され、進捗状況が確認できます。
  • 標準出力の、タグカントファイル
    $ less randRna20M.tc.refFkat.txt
    
    A1BG    5
    
    A1BG-AS1        4
    
    A1CF    2
    
    A2M     2
    
    A2M-AS1 0
    
    :
    
    ZYG11A  0
    
    ZYG11B  1
    
    ZYX     0
    
    ZZEF1   1
    
    ZZZ3    0
    
    no_feature      17035385
    
    ambiguous       1942
    
    too_low_aQual   0
    
    not_aligned     2753965
    
    alignment_not_unique    0
    • 末尾の5行はアノテーション付けされなかったリード数のカウントです。

検討1


  • 入力のアノテーションGTF
    $ head -n 5 hg19.refFlat.20130205.gtf
    
    chr1    hg19_refFlat    exon    11874   12227   0.000000        +       .       gene_id "DDX11L1"; transcript_id "DDX11L1";
    
    chr1    hg19_refFlat    exon    12613   12721   0.000000        +       .       gene_id "DDX11L1"; transcript_id "DDX11L1";
    
    chr1    hg19_refFlat    exon    13221   14409   0.000000        +       .       gene_id "DDX11L1"; transcript_id "DDX11L1";
    
    chr1    hg19_refFlat    exon    14362   14829   0.000000        -       .       gene_id "WASH7P"; transcript_id "WASH7P";
    
    chr1    hg19_refFlat    exon    14970   15038   0.000000        -       .       gene_id "WASH7P"; transcript_id "WASH7P";
    • GTFファイルは、人間には分かりにくい書式ですが、基本的にexonごとに1行あり、第9列のIDで、行間のグループ化ができるようになっています。このままでは、遺伝子がいくつリストされているのか分かりにくいので、以下のようにIDだけ抜き出してみます。
  • アノテーションGTFのgene_id抜き出し
    $ cut -f 9 hg19.refFlat.20130205.gtf | sed 's/gene_id "\([^"]\+\)".*$/\1/' | sort | uniq > tmpGtf.ids
    
    $ head -n 5 tmpHtseq.ids
    
    A1BG
    
    A1BG-AS1
    
    A1CF
    
    A2M
    
    A2M-AS1
    
    $ tail tmpHtseq.ids
    
    ZYG11A
    
    ZYG11B
    
    ZYX
    
    ZZEF1
    
    ZZZ3
    
    no_feature
    
    ambiguous
    
    too_low_aQual
    
    not_aligned
    
    alignment_not_unique
    • GTFファイル中のgene_idが格納されている第9列から、sedでIDを抜き出し、ABC順に並び替え(sort)して、重複除去(uniq)しています。
  • 上記のオリジナルのアノテーションのIDと、HTSeq countの結果出力のIDを以下のように比較
    $ cut -f 1 randRna20M.tc.refFkat.txt > tmpHtseq.ids
    
    $ head -n 5 tmpHtseq.ids
    
    A1BG
    
    A1BG-AS1
    
    A1CF
    
    A2M
    
    A2M-AS1
    
    $ tail -n 5 tmpHtseq.ids
    
    no_feature
    
    ambiguous
    
    too_low_aQual
    
    not_aligned
    
    alignment_not_unique
    
    $ diff tmpGtf.ids tmpHtseq.ids
    
    23706a23707,23711
    
    > no_feature
    
    > ambiguous
    
    > too_low_aQual
    
    > not_aligned
    
    > alignment_not_unique
    • diffでの比較の結果、アノテートされなかったものに関する最後の5行以外は完全一致しました。
    • このことは以下の特徴を示しています。
      #結果の特徴
      1HTSeqの結果は、gene_id単位でタグ数を計上する。
      2HTSeqの結果はgene_idのABC順で出力され、入力に対する過不足がない。(0カウントでも必ず出力行がある)
    • 出力行がずれないことは、結果をマージしてカウントテーブルを作る際に、重宝するので重要です。

検討2


  • HTSeq countでは、基本的に1リードが複数アノテーションに重複して計上されることはないため、アノテートされなかった最後の5行も含め2列目の数値を合計すると、入力のリード数の総合計に一致することになります。これを証明するため、以下の検討を行いました。
  • 入力の行数を数えます。
    $ grep -v "^@" randRna20M.name.sam | wc -l
    
    40000000
  • ヘッダー行をgrep -vで除外した後、行数を数えています。
  • ペアのそれぞれが1行で出力されるため、40,000,000/2=20,000,000ペアあることが分かります。
  • 2列名のカウントの合計値を求めます。
    $ awk 'BEGIN{a=0}{a+=$2}END{print a}' randRna20M.tc.refFkat.txt
    
    20000000
  • awkコマンドで、順次2列目を変数aに加算していって総合計を出力しています。
  • 以上の結果から
    #結果の特徴
    3ペアエンドの場合は、ペアで1カウントされる。
    4全出力行の合計値は、入力行の合計値に等しい。
    4HTSeqの出力では、1リードが重複してアノテートされることはない。
    • この性質はモード(-mオプション)が'union','intersection-strict','intersection-nonempty'のいずれでも同じ。
  • 検討3


    • HTSeq countは"-i transcript_id"オプションを付加することで、Isoform(またはSplicing variant)単位のタグカウントも行うことができます。
    • 一方HTSeqの計算方法は、複数のfeatureで重複する領域のタグは計上されないため、下図のような共通するexonを持つIsoformのタグは重複部分だけ計上されないと予想されます。この挙動の確認を行いました。
    • テストデータには"SRR307920"(ペアエンド、ディレクショナルRNA-seq)からランダムに200万リード抜き出したものを用いています。このデータではディレクショナルのリードは逆になるため、"-s reverse"オプションを付加します。マッピングはBWAで行ったものから、片鎖のみマップしたペアをアンマップに変換して使用ししています。
    • 通常の遺伝子単位の計上
      $ cat SRR307920.bam | htseq-count -q - ucsc_genes_20120309032441.gtf -m intersection-strict -s reverse   > SRR307920.bam.exp
      
      $ tail SRR307920.bam.exp
      
      ZYG11B  965
      
      ZYX     2127
      
      ZZEF1   1091
      
      ZZZ3    723
      
      tAKR    0
      
      no_feature      12943599
      
      ambiguous       143052
      
      too_low_aQual   0
      
      not_aligned     10517490
      
      alignment_not_unique    0
    • Isoform単位の計上
      $ cat SRR307920.bam | htseq-count -q - ucsc_genes_20120309032441.gtf -m intersection-strict -s reverse -i transcript_id  > SRR307920.bam.exp;
      
      $ tail SRR307920.bam.exp
      
      NR_046375       0
      
      NR_046376       0
      
      NR_046377       0
      
      NR_046378       0
      
      NR_046379       0
      
      no_feature      12956989
      
      ambiguous       7038734
      
      too_low_aQual   0
      
      not_aligned     10517490
      
      alignment_not_unique    0
      • 集計の単位がtranscript_idになったのが確認できます。
    • transcript_idで計測した場合のambiguousの割合が高くなっています。
      • -遺伝子(gene_id)単位Isoform(transcript_id)単位
        (A)ambiguous1430527038734
        (B)any feature1782586910916797
        (A)+(B)1796892117955531
        (A)/((A)+(B))0.8%39.2%
    • 適当な遺伝子"ZSCAN2"に着目し、遺伝子単位と、Isoform単位のタグ数を比較しました。
      • Image
      • 上図の吹き出し内は、下表の"略号":"領域内タグ数"("内境界タグ数")を表します。
      • 下図は手作業で計測したセグメントごとのタグ数の内訳です。上図の吹き出しのタグ数と一致します。
      • Image
      • 遺伝子単位では、"-m intersection-strict"で除外される境界領域のリードを除いた全てのIsoformの和集合が計上されています。
      • Isoform単位では領域が重複しないセグメントのみ計上されるため、Isoform1ではJ2+E4, Isoform2はJ3+E4, Isoform3はJ4+E5のみ計上されています。
    • 境界領域の例Image
      • 今回用いたBWAではsplicing junctionを考慮しない分tophatなどより高速ですが、境界をまたぐリードはsoft clippingという非アライメント領域として扱われます。図中グレーの三角(三角内の数値は非アライメント領域の長さ)。
      • 稀に、境界領域をまたぐリードが含まれるため、"-m intersection-strict"オプション付加時はそのようなリードは計上されません。
      • exonの境界以降がsoft clippingされている場合は、境界リードとみなされず、適切に計上されるようです。

    インストール方法


    参考



    Contact us
    Copyright © 2009-2017 National Institute of Genetics  [Site Policy] [Privacy Policy]