HTSeq

Python 3 に HTSeq をインストールすることができなかったため、以下のサンプルコードはすべて Python 2.7 で実行した。

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

GFF ファイルの読み込み

GFF あるいは GTF ファイルは GFF_Reader メソッドを利用して読み込む。

import HTSeq

gff_file = HTSeq.GFF_Reader('./v22/Lactobacillus_casei_12a.GCA_000309565.1.22.gff3')

counts = {}
for feature in gff_file:
    if feature.type == "exon":
        print("--------------------------------")
        print(feature.name)
        print(feature.attr)
        print(feature.iv)
## --------------------------------
## exon:EKP96509-1
## {'Parent': 'transcript:EKP96509', 'rank': '1', 'ID': 'exon:EKP96509-1', 'constitutive': '1', 'ensembl_end_phase': '0', 'ensembl_phase': '0'}
## ctg28_00010:[465890,466664)/-
## --------------------------------
## exon:EKP96510-1
## {'Parent': 'transcript:EKP96510', 'rank': '1', 'ID': 'exon:EKP96510-1', 'constitutive': '1', 'ensembl_end_phase': '0', 'ensembl_phase': '0'}
## ctg28_00010:[466841,467447)/+
## --------------------------------
## exon:EKP96511-1
## {'Parent': 'transcript:EKP96511', 'rank': '1', 'ID': 'exon:EKP96511-1', 'constitutive': '1', 'ensembl_end_phase': '0', 'ensembl_phase': '0'}
## ctg28_00010:[467786,468011)/+

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