htseq-countを用いた遺伝子発現量のカウント【Python】

2020-11-192020-11-29

RNA-Seqのシークエンスデータをマッピングし、それぞれの遺伝子領域にいくつのリードがマッピングされたかをカウントすることで遺伝子発現量を定量化することができます。マッピングファイルから遺伝子発現量をカウントするプログラムはいくつかありますが、ここではPythonを用いて作成されたプログラムであるhtseq-countについて説明していきます。

動作確認環境

OS

Windows10 + WSL (WindowsユーザーのためのPythonを用いたゲノム解析環境)

Python

Python 3.7.8

モジュール

BCBio-GFF 0.6.6
HTSeq 0.12.4

[toc]

HTSeqモジュールとhtseq-countコマンド

NGS関連の様々なデータ形式を扱うためのパーサーやクラスを提供するモジュールであるHTSeqでは、RNA-Seqの遺伝子発現量カウントに特化したスクリプトとしてhtseq-countコマンドが用意されています。(公式ドキュメント)

htseq-countコマンド自体のソースコードは以下の通りで、

https://github.com/htseq/htseq/blob/master/HTSeq/scripts/count.py

これを＜こちら＞から実行しています。ソースコードを見てわかる通り、HTSeqのモジュールを用いてRNA-Seqの遺伝子発現量カウント処理を実装したスクリプトがhtseq-countコマンドだと言えます。なお、Pythonから実行するAPIは用意されていないようなので、コマンドラインから実行する必要があります。

ここでは、このhtseq-countコマンドを用いて、RNA-Seqのマッピングデータの発現量をカウントする方法を見ていきましょう。

htseq-countコマンドを実行するための準備

リファレンスとなる遺伝子のリスト(アノテーションファイル)

htseq-countではGTF形式で遺伝子のアノテーションファイルを指定することが想定されていますが、GFF形式でアノテーションファイルを指定することも可能です。

ただし、GTF形式ではfeature領域を表すIDが必ず「gene_id」というタグで与えられていますが、GFF形式ではfeature領域を表すIDは必須項目ではないので、feature領域を表すIDがない場合は自分でGFFファイルにfeature領域を表すIDを予め書き込んでおく必要があります。GFF形式のfeature領域で親子関係のあるところは同じ「gene_id」となるようにしたプログラムの例は以下のようになります。なお、このプログラムは「先進ゲノム支援 2018年度情報解析中級者講習会」のプログラム例を参考にしています。

from BCBio import GFF

# 読み込むGFFファイル名と出力するファイル名を指定する
file_name = 'sample.gff'
out_file_name = 'sample_mod.gff'

# 処理したレコードを一時的に格納する
records = []

# 指定されたGFFファイルを開く
with open(file_name) as handle:
    gene_count = 1

    # GFFファイルのレコードごとに処理をする
    for record in GFF.parse(handle):
        for f in record.features:
            if f.type == "gene" or f.type == "pseudogene":
                gene_id = "gene_" + str(gene_count).zfill(4)
                f.qualifiers["gene_id"] = [gene_id]
                for sf in f.sub_features:
                    sf.qualifiers["gene_id"] = [gene_id]
                    for ssf in sf.sub_features:
                        ssf.qualifiers["gene_id"] = [gene_id]
                gene_count += 1
        records.append(record)

# 指定したGFFファイルを作成し、書き込む
with open(out_file_name, "w") as fh:
    GFF.write(records, fh)

GFF形式のファイルの扱いについては以下の記事もご覧ください。

また、feature領域を表すIDが「gene_id」以外のタグ名で与えられている場合は、コマンドの-iオプションに指定します。エクソン領域以外の発現解析を行うときは、-tオプションにfeature typeを指定します。

発現量をカウントするゲノムデータ(マッピングファイル)

発現量をカウントしたいゲノムデータをマッピング済みのデータ形式(SAM形式／BAM形式、またはCRAM形式)で準備します。paired-end リードの場合はBAM形式であらかじめソートしておく必要があります。

マッピング済みのデータ形式については以下の記事もご覧ください。

htseq-countコマンドを実行する

htseq-countコマンドの書式は以下のようになります。

htseq-count [options] <alignment_files> <gff_file>

コマンドオプション

htseq-countコマンドの主なオプションは以下の通りになります。

オプション	説明
-r	どのような順番でアラインメントがソートされているかを示します。(デフォルト：name) 次の値を指定可能です：name (名前順にソートした場合) / pos (ゲノム位置順にソートした場合) ※ single-end リードでは不要ですが、paired-end リードの場合はアラインメントがソートされている必要があります。
-s	ストランド特異的なRNA-Seqかどうかを指定します。(デフォルト：yes)
-a	マッピングクオリティスコアが与えられた閾値以下の場合はスキップします。(デフォルト：10)
-t	発現量をカウントするfeature typeを指定します。(デフォルト：exon)
-i	feature IDのタグ名を指定します。(デフォルト：gene_id)
-m	複数のfeatureにまたがっているリードの扱いを指定します。(デフォルト：union) 次のモードから選択します：union / intersection-strict / intersection-nonempty
-n	使用するプロセッサ数。(デフォルト：1) 1つのマッピングファイルの発現量カウントには1つのCPUしか使用されませんが、複数のマッピングファイルを指定したときに使用するプロセッサが多いと高速化できます。

なお、以前は-fオプションで入力ファイル形式がSAM形式かBAM形式かなどを指定する必要があったようですが、現在では自動で認識するようになっているので不要となっています。

single-end リードの場合

single-end リードの場合はSAMファイルやBAMファイルをそのまま指定するだけで構いません。例えば、annotation.gffをアノテーションファイルとしてsample.bamの発現量をカウントするためには次のようなコマンドになります。

htseq-count sample.bam annotation.gff > result.txt

カウント結果をresult.txtにリダイレクションで保存しています。

Python上から実行する場合は、subprocessモジュールのrun関数の引数にコマンドをそのまま書いて渡しましょう。

import subprocess
subprocess.run('htseq-count sample.bam annotation.gff > result.txt', shell=True)

出力形式はタブ区切りのテキスト形式で、以下のようにgene_idごとにカウント結果が示されたものになります。

paired-end リードの場合

paired-end リードの場合はマッピングファイルをソートしておく必要があります。そのため、SAMファイルの場合はあらかじめBAMファイルに変換してソートを行います。samtoolsのsortコマンドでソートする際に、引数に何もつけなければゲノム位置でソートされ、-n引数をつけると名前順にソートされます。htseq-countコマンドでもソートの順番を指定する必要があり、デフォルトでは名前順となっているので、ゲノム位置順でソートした場合は-r引数にposを指定します。

htseq-count -r pos sample_sorted.bam annotation.gff > result.txt

Python上から実行する場合は、subprocessモジュールのrun関数の引数にコマンドをそのまま書いて渡しましょう。

import subprocess
subprocess.run('htseq-count -r pos sample_sorted.bam annotation.gff > result.txt', shell=True)

なお、インデックスは作成していなくても実行できますが、あらかじめインデックスファイル(「***.bai」)を作成して同一フォルダ内に配置しておいた方が若干早いようです。

複数のマッピングファイルの発現を解析する場合

single-end リードでもpaired-end リードでも同様ですが、複数のマッピングファイル(sample1_sorted.bam, sample2_sorted.bam, sample3_sorted.bam)の遺伝子発現のカウントをする場合は次のようにマッピングファイルを連続で指定するだけで構いません。この時、-nオプションでプロセッサの数を指定すると複数のマッピングファイルを同時に処理できるので高速化できます。

htseq-count -r pos -n 16 sample1_sorted.bam sample2_sorted.bam sample3_sorted.bam annotation.gff > all_count_data.txt

Python上から実行する場合は、subprocessモジュールのrun関数の引数にコマンドをそのまま書いて渡しましょう。

import subprocess
subprocess.run('htseq-count -r pos -n 16 sample1_sorted.bam sample2_sorted.bam sample3_sorted.bam annotation.gff > all_count_data.txt', shell=True)

出力形式はタブ区切りのテキスト形式で、以下のようにコマンドに指定したマッピングファイルの順番でgene_idごとにカウント結果が示されたものになります。