SAM/BAM ファイルからリードカウントを取得するソフトウェア

HTSeq

HTSeq は Python で書かれているプログラムで、主に BAM あるいは SAM 形式のファイルを処理する機能を提供している。Bowtie や TopHat などのマッピングプログラムを利用してマッピングを行った結果からカウントデータを取得する際に利用できる。特に、HTSeq をインストールすると htseq-count コマンドが利用できるようになり、コマンド 1 行分で BAM あるいは SAM ファイルからカウントデータを取得できる。

HTSeq のインストール

HTSeq は Python のモジュールとして提供されている。HTSeq は NumPy モジュールを利用しているので、NumPy を先にインストールしておく必要がある。また、Pysam モジュールは必須ではないが、HTSeq で BAM ファイルを取り扱うときに必要になってくるので、ここでインストールしておく後々が楽になる。

pip install numpy --user
pip install htseq --user
pip install pysam --user

htseq-count の基本的な使い方

htseq-count を利用してリファレンスゲノム上にマッピングされたリードを遺伝子ごとに数える際に必要な入力ファイルとして、マッピング結果である SAM ファイルと遺伝子の座標情報などが保存されているアノテーションファイル(.gff, .gtf)が必要である。

single-end

single-end リードの場合は以下のように実行すれば良い。bowtie_result.sam と annotation.gtf を入力として、数え上げたカウントデータは count_data.txt ファイルに保存される。

htseq-count bowtie_result.sam annotation.gtf > count_data.txt

これにオプションなどを加えことでより細かい調整が行うことができる。例えば入力ファイルが SAM ではなく、BAM である場合は以下のようにオプション -f を与える。

htseq-count -f bam bowtie_result.bam annotation.gtf > count_data.txt

paired-end

paired-end リードの場合、SAM あるいは BAM ファイルが予めソートさせておく必要がある。ソートされていない場合はカウントデータを正しく取得できない。また、paired-end の場合は、strand specific に注意してオプションを指定する必要がある。

マッピング結果が SAM ファイルである場合、ソートを行うために BAM ファイルへ変換を行う。マッピング結果が BAM ファイルである場合は、このステップを行う必要はない。

samtools view -bS align.sam > align.bam

次にソートを行う。samtools でオプションを指定しないでソートすると、リードが遺伝子座標順にソートされる。

samtools sort align.bam align_sorted

コマンドの実行が終了すると、この場合は align_sorted.bam が生成される。これを htseq-count の入力ファイルとして用いる。上のソートは遺伝子座標順であるから -r pos と指定しておく必要がある。

htseq-count -f bam -r pos align_sorted.bam annotation.gtf > count_data.txt

取得したカウントデータの合計値がマッピングプログラムから出力されたログ中の mapped read 数と大きく離れている場合、-s オプションなどを調整してみてください。

htseq-count -f bam -r pos -s reverse align_sorted.bam annotation.gtf > count_data.txt

htseq-count のオプション

htseq-count で利用できるオプションなどは以下のようになっている。詳しくは公式サイトに詳しく書かれている。

option default
-f sam 入力するマッピング結果はのフォーマットを指定する。sam または bam を指定できる。
-r name ペアエンドリードの場合、SAM あるいは BAM ファイルを名前順または座標順にソートする必要がある。入力した SAM あるいは BAM はどのようにソートされたかを指定する。name または pos を指定する。
samtools sort でソートした場合は -n pos とする必要がある。また、samtools sort -n でソートした場合は -r name と指定する。ただし省略した場合は -r name とみなす。
-s yes ストランド特異的な RNA-seq の場合に指定する。yes を指定した時、single-end ならばすべてのリードがセンス鎖にマッピングされる場合を考慮する。paired-end ならば最初の _1 のリードがすべて同じ鎖にマッピングされ、_2 のリードがすべてアンチセンス鎖にマッピングされる場合を考慮する。yes の代わりに reverse を指定すると、yes とは逆の設定でカウントされる。
-a 10 クオリティが与えられた値未満の場合に無視する。
-t exon feature type を指定する。GFF ファイルの 3 列目に書かれている feature type を指定する。exon の他に、geneCDSrRNA などがある。また、使用する GFF / GTF ファイルにより、feature type が異なる。RNA-seq による発現変動解析を行うとき、利用するアノテーションファイルが Ensembl GTF の場合は exon にする。
-i gene_id feature type の属性 ID を指定する。遺伝子ごとにカウントしたいか、トランスクリプトごとにカウントしたいかなどを指定する。
-m union リードをカウントする際に利用するアルゴリズムを指定する。unionintersection-strictintersection-nonempty を指定できる。(詳細

References

  1. Anders S, Theodor PP, Huber W. HTSeq — A Python framework to work with high-throughput sequencing data. Bioinformatics 2015, 31(2):166-9. PubMed Abstract