- 著者: Simon Anders, Paul Theodor Pyl, Wolfgang Huber
- Corresponding author: Simon Anders (European Molecular Biology Laboratory, Heidelberg, Germany)
- 雑誌: Bioinformatics
- 発行年: 2015
- Epub日: 2014-09-25
- Article種別: Original Article (Methods / Tool paper)
- PMID: 25260700
背景
HTS (High-throughput sequencing、高速並列シーケンス) の急速な普及により、RNA-seq (RNA sequencing)・ChIP-seq (chromatin immunoprecipitation sequencing)・variant calling 等の多様な assay が出現し、それぞれ専用解析パイプラインを要求するようになった。アラインメント・アセンブリ等の “big tasks” については BWA (Burrows-Wheeler Aligner)・Bowtie・Tophat 等の標準ツールが利用可能であり、開発者は適切な選択肢からツールを組合せて pipeline を構築できる。しかし、ChIP-seq の TSS (transcription start site、転写開始点) 周辺 coverage profile 集計、RNA-seq read counting における ambiguous overlap 処理、複数 file format 間のデータ変換等、標準パイプラインから外れる “specialized tasks” にはユーザーが自前 script を書く必要がある。具体的に扱う必要のある format には FASTA (FAST-All、ヌクレオチド reference 配列形式)、FASTQ (FASTA + base-call quality)、SAM (Sequence Alignment/Map)、BAM (Binary Alignment/Map)、GFF (General Feature Format)、GTF (Gene Transfer Format)、VCF (Variant Call Format)、BED (Browser Extensible Data)、Wiggle (ゲノム track score 形式)、CIGAR (Compact Idiosyncratic Gapped Alignment Report) 文字列等がある。
先行研究の文脈では、Biopython (Cock et al. Bioinformatics 2008) は古典的 sequence analysis を Python で提供し PMID 19304878 で広く参照されたが HTS file format 対応は限定的、PySam (htslib ラッパー) は SAM/BAM 専用、pybedtools (Dale et al. Bioinformatics 2011) は BED interval 操作専用と、いずれも単一機能であり多種 format を統合した汎用 base layer は存在しなかった。R 側では Bioconductor の Rsamtools・GenomicRanges (Lawrence et al. PLoS Comput Biol 2013、PMID 23950696) が同等機能を提供していたが、Python ユーザーには等価ライブラリが不足しており、custom script を書く際の base layer を再発明する負担が大きかった点が当時の knowledge gap であった。
すなわち先行研究では「Python で HTS の多種 file format を統一的に扱い、genomic interval ベースの効率的クエリを行える統合フレームワーク」が決定的に不足しており、ユーザーは format ごとに分散したライブラリを組合せ座標系変換を自前実装するという生産性低下を強いられていた。HTSeq の初期版は 2010 年に公開され htseq-count は RNA-seq differential expression 解析の前処理として広く使われるようになったが、フレームワーク全体の系統的記述と最近の改良点 (paired-end 処理の sort-free 化等) を公式論文として発表する必要性があった。
目的
HTS データに対する custom script 作成の生産性を高める Python ライブラリ HTSeq を提示し、その設計思想・主要構成要素・代表的応用例 (htseq-qa・htseq-count) を記述する。具体的には、(1) 主要 HTS file format に対する統一された parser と record クラス、(2) genomic position・interval にデータを紐付ける効率的な container クラス GenomicArray、(3) 複雑な alignment (gapped・spliced・paired-end) のハンドリング、(4) 中規模 Python 経験者でも custom script を書けるドキュメント・API 設計、(5) RNA-seq read counting の de facto standard tool である htseq-count の現行仕様、を 1 本の論文として公開することを目指す。
結果
6 系統の file format を統合する parser layer: HTSeq は FASTA/FASTQ・SAM/BAM (legacy formats 含む)・GFF/GTF・VCF・BED・Wiggle の 6 系統以上に対する parser を統合し (Fig 1a)、record class (SequenceWithQualities・GenomicPosition・GenomicInterval・SAM_Alignment・CigarOperation 等) を parser 間で共通利用する設計により、parser 間でのデータ受け渡しが seamless となる。zero-based half-open interval convention を Python 標準に揃えつつ、入力 format ごとの座標系差異 (BED は 0-based、GFF は 1-based 等) を自動変換するため、ユーザーは座標系の差を意識せずに interval 演算を記述できる。SAM_Alignment object 内の cigar field は CigarOperation のリストとして展開され、5 種類の operation type (M・N・I・D・S 等) を size・query 範囲・reference 範囲の 3 属性で構造化保持する (Fig 1b)。
GenomicArray の二重バックエンド設計: step-based + tree-backed 実装 (C++ red-black tree map) と NumPy array backend の透過切替により、histone mark coverage profile (sparse、step-based 有利)・gene model annotation (sparse、object-valued)・read count distribution (dense、NumPy 有利) 等、データ密度や型 (scalar / object) が異なる用途に同一 API で対応可能となる (Fig 2)。memmap を活用すれば genome-wide single-base resolution データ (ヒトゲノム約 30 億 bp 規模) も RAM 制約下で扱える。Fig 2 の GenomicArrayOfSets の例では、interval [1000-1020] に gene A、[1020-1320] に intron、[1320-1350] に gene B 等が割当てられ、query interval [1000-1358] を投げると 5 step に分割された feature set が返却される構造が示される。
htseq-count の処理性能ベンチマーク (n=2 構成での測定、n=1 マシンを 2 入力 sort 順序で評価): 標準ノート PC 上で htseq-count v0.6.1 は 約 1.2 million reads/min (paired-end は 約 0.6M read pairs/min = 半減) を、ヒト遺伝子アノテーション (約 60,000 genes) を RAM 約 250 MB で保持しつつ処理可能である (本文 Sec 5 で記述、n=1 マシン構成での実測)。SAM ファイルが read name 順でなく position 順 (mate が隣接しない) の場合でも、buffer 機構により処理時間は 約 2 倍以内 (約 2-fold 増) に抑制され、26 GB SAM file でも追加メモリは 450 MB 未満 (= 元データの約 1.7%) で済む。Fonseca et al. PLoS ONE 2014 (PMID 25268973) のベンチマーク (n=5 counting tool 比較) で htseq-count は accuracy 面で良好と評価された一方、後発の Bioinformatics et al. Basic 2014 (featureCounts) は C 実装ゆえ runtime で htseq-count を 概ね 5-10 倍 (= 5-10-fold) 上回り、以降の RNA-seq pipeline では速度で featureCounts、柔軟性で HTSeq という棲み分けが生じた。HTSeq の総合性能は n=3 評価指標 (read throughput / memory footprint / buffer overhead) で実証されている。
Loop-based query paradigm の独自性: HTSeq の paradigm は “一方の interval list を GenomicArrayOfSets に格納し、もう一方を loop で 1 件ずつ query する” 方式で、Bioinformatics et al. Basic 2010 (BEDtools) や IRanges/GenomicRanges の “two-list overlap” paradigm と異なる。この explicit looping は、split read・gapped alignment・ambiguous mapping 等、inner loop 内での branching を 多数の特殊ケース (典型的に 5 種以上の条件分岐) が必要となる複雑処理で開発者の思考に直感的に沿うため、custom script 開発の生産性向上に寄与する。
応用範囲と RNA-seq pipeline 統合: 同梱の htseq-qa は read 位置別 nucleotide composition と base-call quality (Phred score 分布) をプロットし初期 QC に有用である。htseq-count は GenomeBiol et al. Basic 2014 (DESeq2) や Bioinformatics et al. Basic 2010 (edgeR) の入力 count matrix 生成として de facto standard 化し、Bioinformatics et al. Basic 2014 (Trimmomatic) と組合せて典型的 RNA-seq pipeline (Trim → align → htseq-count → DESeq2/edgeR) を完結できる。HTSeq は htseq-count 以外にも、ChIP-seq の peak 周辺 coverage 集計、bulk RNA-seq の novel transcript detection 用 read partitioning、VCF アノテーション集計等の custom script 基盤として活用される。
考察/結論
① 先行研究との違い: HTSeq は performant だが monolithic な専用ツール (SAMtools・BEDtools の C 実装等) と完全自前 script の中間ニッチを埋めるという、これまでの Python bioinformatics ライブラリと異なる位置づけを取る。Biopython (Cock et al. 2009) が古典的 sequence analysis に強い一方 HTS 対応が手薄であった点と対照的に、HTSeq は HTS の各種 file format と genomic coordinate 演算に集中特化する。R 側の Rsamtools/GenomicRanges (Lawrence et al. 2013) と機能的に並ぶが、Python ベース解析環境という別エコシステムで等価機能を提供した点で相違が明確である。BEDtools の two-list overlap paradigm とは異なり、HTSeq は一方を GenomicArrayOfSets に格納して loop query する設計を取る点で根本的に paradigm が違う。
② 新規性: 本研究で初めて、HTS データ操作に必要な 6 系統以上の file format parser と GenomicArray・GenomicArrayOfSets という genomic coordinate データ構造を Python で統合的に提供する新規なフレームワークが体系的に記述された。これまで報告されていない設計として、step-based tree + NumPy array の透過切替バックエンド、CigarOperation による gapped alignment の構造化表現、buffer-based paired-end pairing の 3 点が挙げられる。htseq-count についても、ambiguous overlap を 3 mode (union/intersection-strict/intersection-nonempty) で明示的に区別する設計は当時の counting tool では新規であった。
③ 臨床応用: 本ライブラリの臨床応用としては、bench-to-bedside の橋渡し的役割として、希少疾患 whole-genome sequencing の variant filtering、cancer panel sequencing の custom annotation aggregation、clinical RNA-seq の fusion transcript detection、long-read sequencing データの preliminary パーシング等、標準商用パイプラインに収まらないシナリオで継続的に有用である。臨床的意義としては、患者検体に特化した custom QC pipeline 開発の生産性向上を通じて、診断精度向上に間接的に寄与する。translational research としては、cancer biomarker discovery において tumor-normal pair の read counting や allele-specific expression 解析の prototyping ツールとして活用される。
④ 残された課題: 今後の検討課題として、(1) interpreted Python のオーバヘッドにより featureCounts のような C 実装には速度で劣るため、large cohort 解析では runtime が bottleneck となる、(2) single-cell RNA-seq・spatial transcriptomics・multi-omics 等 2015 年以降の新パラダイムに対する直接サポートは限定的で、scanpy・anndata 等の後発エコシステムが台頭した、(3) ambiguous overlap・multi-mapping read 処理の選択肢は htseq-count では union / intersection-strict / intersection-nonempty の 3 mode に限定され、より柔軟な EM-based 推定 (RSEM・Salmon) には対応しない、という limitation が挙げられる。今後の研究方向としては、long-read alignment 用 CigarOperation 拡張、Cython 以外の高速化 (Numba・PyPy 対応) の検討、scverse エコシステムとの interop 強化が望まれる。
方法
実装は Python (Python 2.6+/3.x 対応) + Cython (型情報を付与し C 変換することで性能改善、Behnel et al. Comput Sci Eng 2011) で構成される。BAM_Reader / BAM_Writer は PySam (htslib ラッパー) を依存ライブラリとして利用するため、BAM 読み書きには pysam インストールが必須となる一方、テキスト系 (FASTA/FASTQ/SAM/GTF/VCF/BED/Wiggle) はネイティブ実装で外部依存なし。Parser はすべて iterator generator として実装され、for record in HTSeq.SAM_Reader(filename): の形で record object を順次産出する設計を取る。
GenomicArray クラス は piecewise constant な genomic position 依存データを効率的に保持するため C++ の map template (red-black tree、Josuttis 1999) を SWIG (Beazley et al. 1996) 経由で Python から利用可能にし、step ベースの圧縮表現を採用する。dense data 用に NumPy array バックエンド (van der Walt et al. 2011、memmap 対応で巨大データを disk swap 可能) も切替可能で、storage 選択は API レベルで透過的に行える。
GenomicArrayOfSets は overlap 可能な annotation (exon・gene model 等) を Python set で各 step に保持し、interval query で重なる feature 全体を union として返すデータ構造である (Fig 2)。SAM_Alignment クラス は CigarOperation のリストを通じて gapped alignment (cigar string 20M300N30M2I8M のような split read 等) を構造化表現し、pair_SAM_alignments_with_buffer 関数で SAM の paired record を read name 順でなくとも buffer ベースで pairing する。
htseq-count スクリプト は SAM/BAM + GTF/GFF を入力に、各遺伝子 (meta-feature) の exon に重なる read を集計し、ambiguous overlap (複数遺伝子重なり) や multi-mapping read は除外する設計を取り、union / intersection-strict / intersection-nonempty の 3 mode を提供する。Identifier として GitHub repo (http://www-huber.embl.de/HTSeq) と PyPI package index (https://pypi.python.org/pypi/HTSeq) で配布、GNU General Public Licence の OSS としてリリースされる。最近のバージョン (v0.6.1+) では read name sort なしで paired-end 処理が可能となり、buffer 機構で mate 検索を吸収する設計が追加された。