マッピング結果を処理する Python モジュール

HTSeq

HTSeq (Anders et al., 2014) は BAM あるいは SAM ファイルからカウントデータを取得する際に利用されるプログラムである。その機能に特化したコマンドとしては htseq-count がある。ここで、このコマンドの使い方を説明せずに、Python のモジュールの一つとして使い方を示す。

HTSeq は Python3 に対応しておらず、Python 2.7 の環境で実行する必要がある。以下に示す Python スクリプトはすべて Python 2.7 で動作するスクリプトである。

インストール

HTSeq は Python3 に対応していない。そのため、Python3 が存在する環境におて、 pip コマンドで HTSeq をインストールしようとするとエラーが発生する場合がある。

pip install htseq
## Collecting htseq
##   Using cached HTSeq-0.7.1.tar.gz
##     Complete output from command python setup.py egg_info:
##     Error in setup script for HTSeq:
##     See the python3 branch on github for Python 3 support.
##     Please use Python 2.7.
## 
##     ----------------------------------------
## Command "python setup.py egg_info" failed with error code 1 in /private/var/folders/8s/m3wznvcx3zgc7d9vs_1l5rb80000gn/T/pip-build-xry2ffyg/htseq/

このようなエラーが起きたとき、pip コマンドの代わりに pip2 を用いればよい。

pip install htseq

GFF3 ファイルの処理

HTSeq の GFF_Reader メソッドを利用することで、遺伝子のアノテーションデータが記載されている GTF あるいは GFF ファイルを処理することができる。例えば、次の Python スクリプトでは、GFF_Reader メソッドを利用して、A. thaliana の GFF ファイルを読み込み、mRNA のデータを画面上に出力する例である。

import HTSeq
gff_file = HTSeq.GFF_Reader('Arabidopsis_thaliana.TAIR10.34.gff3.gz')

for feature in gff_file:
    if feature.type == "mRNA":
        print("--------------------------------")
        print(feature.name)
        print(feature.attr)
        print(feature.iv)
## --------------------------------
## transcript:AT1G01010.1
## {'transcript_id': 'AT1G01010.1', 'ID': 'transcript:AT1G01010.1', 'Parent': 'gene:AT1G01010', 'biotype': 'protein_coding'}
## 1:[3630,5899)/+
## --------------------------------
## transcript:AT1G01020.6
## {'transcript_id': 'AT1G01020.6', 'ID': 'transcript:AT1G01020.6', 'Parent': 'gene:AT1G01020', 'biotype': 'protein_coding'}
## 1:[6787,8737)/-
## --------------------------------
## transcript:AT1G01020.2
## {'transcript_id': 'AT1G01020.2', 'ID': 'transcript:AT1G01020.2', 'Parent': 'gene:AT1G01020', 'biotype': 'protein_coding'}
## 1:[6787,8737)/-

マッピング結果からカウントデータを取得

GFF ファイル内に書かれている feature として mRNA のほかに、gene、exon、CDS や miRNA などがある。GFF ファイルから特定の feature のデータだけを取り出して 1 つのオブジェクトに保存する場合は、次のようにする。

import HTSeq

gtf_file = HTSeq.GFF_Reader('Cama.genome.gene.gff')
exons = HTSeq.GenomicArrayOfSets( "auto", stranded=True )

for feature in gtf_file:
    if feature.type == "exon":
       exons[ feature.iv ] += feature.attr["gene_id"]


gff_file = HTSeq.GFF_Reader('Arabidopsis_thaliana.TAIR10.34.gff3.gz')

for feature in gff_file:
    if feature.type == "mRNA":
        print("--------------------------------")
        print(feature.name)
        print(feature.attr)
        print(feature.iv)

GFF ファイル中に多くのアノテーションが記載されている。ここで、exon のアノテーションのみを取り出して gff_exon オブジェクトに保存する。その際に、GenomicArrayOfSets を利用してひな形を作成する。RNA-seq が stranded-specific の場合、stranded = True も合わせて指定する。

gff_exon = HTSeq.GenomicArrayOfSets("auto", stranded = True)

for feature in gff_file:
    if feature.type == "exon":
        gff_exon[feature.iv] += feature.name

list(gff_exon.steps())[0:2]
## [(<GenomicInterval object 'ctg22_00009', [0,35), strand '+'>, set([])), (<GenomicInterval object 'ctg22_00009', [35,389), strand '+'>, set(['exon:EKP97577-1']))]

さらに gff_exon 中の指定範囲のアノテーションデータを取り出す際には、以下のように GenomicInterval メソッドを利用して範囲を作成して取得する。

iv = HTSeq.GenomicInterval("ctg28_00010", 10000, 20000, "+")

list(gff_exon[iv].steps())
## [(<GenomicInterval object 'ctg28_00010', [10000,18115), strand '+'>, set([])), (<GenomicInterval object 'ctg28_00010', [18115,19096), strand '+'>, set(['exon:EKP96087-1'])), (<GenomicInterval object 'ctg28_00010', [19096,19501), strand '+'>, set([])), (<GenomicInterval object 'ctg28_00010', [19501,20000), strand '+'>, set(['exon:EKP96088-1']))]

マッピング結果からカウントデータを取得

次にマップされたリードをカウントする例を示す。カウントする際にディクショナリを利用するので、予めカウントをすべて 0 にリセットする。次に、リードを 1 つ読む度に、GFF のアノテーションデータにフィッティングする(どの exon に由来するかを調べる) intersection_update(step_set)。フィッティングした結果、1 箇所にしかフィッティングしなかったリードに対して、そのカウントを 1 増やす len(iset) == 1

counts = {}
for feature in gff_file:
    if feature.type == "exon":
        counts[feature.name] = 0

bam_file = HTSeq.BAM_Reader("./SRR616268.sort.bam")

for read in bam_file:
    if read.aligned:
        iset = None
        for step_iv, step_set in gff_exon[read.iv].steps():
            if iset is None:
                iset = step_set.copy()
            else:
                iset.intersection_update(step_set)
        if len(iset) == 1:
            counts[list(iset)[0]] += 1

for name in sorted(counts.keys()):
    print name, counts[name]
## exon:EKP96074-1 448
## exon:EKP96075-1 198
## exon:EKP96076-1 156
## exon:EKP96077-1 5
## exon:EKP96078-1 1041
## exon:EKP96079-1 151

References

  • Anders S, Pyl PT, Huber W. HTSeq - a Python framework to work with high-throughput sequencing data. Bioinformatics. 2014, 31(2):166-9. DOI: 10.1093/bioinformatics/btu638 PMID: 25260700