- 著者: Yang Liao, Gordon K. Smyth, Wei Shi
- Corresponding author: Wei Shi (Walter and Eliza Hall Institute of Medical Research, Parkville, Australia)
- 雑誌: Bioinformatics
- 発行年: 2014
- Epub日: 2013-11-13
- Article種別: Original Article
- PMID: 24227677
背景
次世代シーケンス (NGS) 技術は、DNAシーケンスを前例のない速度で可能にし、生物学に革命をもたらした。数百万の短いシーケンスリードが生成され、これらは通常、リファレンスゲノムにアラインされる。多くのアプリケーションにおいて、下流解析に必要とされる主要な情報は、各ゲノム特徴（例えば、各エキソンや各遺伝子）にマッピングされるリードの数である。このリードを数えるプロセスは「リードサマライゼーション」と呼ばれ、差次的発現解析や差次的結合解析のためのカウントベース統計手法（例えば、Robinson et al. Bioinformatics 2010、DESeq2、limmaなど）の入力として不可欠である。リードサマライゼーションは多様なゲノム解析に必要とされるが、これまで文献において比較的少ない注目しか受けてこなかった。この問題は表面上は単純に見えるかもしれないが、実際には多くの微妙な点が存在する。リードカウントプログラムは、DNAシーケンスとRNAシーケンスの両方、およびシングルエンドリードとペアエンドリードの両方に対応する必要がある。カウントされるリードまたはペアエンドフラグメントは、リファレンスゲノムに対して挿入、欠失、または融合を含む可能性があり、これらの複雑さは各リードまたはフラグメントの位置を各可能なターゲットゲノム特徴と比較する際に考慮される必要がある。特徴の数が多い場合、リードカウンティングの計算コストはリードのアライメントステップに匹敵する可能性がある。
ChIP-seqのようなDNAシーケンスリードは、転写因子結合部位の同定（Zhang et al. GenomeBiol 2008）、ヒストン修飾の解析（Park, 2009）、DNAメチル化の検出（Harris et al., 2010）など、さまざまな技術から生じる。DNAリードの対象となるゲノム特徴は、通常、単純なゲノム区間として指定できる。例えば、Pal et al. (2013)は、遺伝子プロモーター領域および遺伝子全体におけるヒストン修飾に関連するリードをカウントした。RNA-seqリードのカウントは、エキソンスプライシングに対応する必要があるため、やや複雑である。一つの方法は、各アノテーションされたエキソンにオーバーラップするリードをカウントすることであり、これは実験条件間の選択的スプライシングをテストするために使用できる（Anders et al., 2012; Reyes et al., 2013）。もう一つの一般的なアプローチは、各遺伝子の任意のエキソンにオーバーラップするすべてのリードをカウントすることにより、遺伝子レベルでカウントを要約することである（Anders et al., 2013; Bhattacharyya et al., 2013; Man et al., 2013）。
当時、主要なリードサマライゼーションツールとしては、BioconductorのGenomicRangesパッケージに含まれるsummarizeOverlaps関数やIRanges (Infrastructure for Manipulating Intervals on Sequences) パッケージのcountOverlaps関数がR言語で実装され、HTSeq (High-throughput Sequencing) フレームワークのhtseq-countスクリプトがPython言語で実装されていた。これらのツールはインタープリタ言語ベースであり、速度とメモリ効率に最適化されていなかったため、大規模なデータセットの処理において計算コストがボトルネックとなることが指摘されていた。一方、Quinlan et al. Bioinformatics 2010のBEDToolsはC++で実装され高速であったが、RNA-seq解析に不可欠なエキソンジャンクションへの対応が不足していた。アライナーの高速化が進む一方で（Li et al. Bioinformatics 2009、Langmead et al. GenomeBiol 2009、Liao et al., 2013）、後段のリードカウンティングがボトルネックとなりうる状況であり、より効率的なリードサマライゼーションツールの開発が喫緊の課題として認識されていた。特に、大規模なデータセットや不完全なゲノムアセンブリに対する効率的な処理能力が不足しており、この点が未解明な課題として残されていた。
目的
本研究の目的は、次世代シーケンスデータ（RNA-seq、ChIP-seq、およびその他のDNA-seq由来リード）を、エキソン、遺伝子、プロモーター領域、遺伝子ボディなどの任意のゲノム特徴に対して、効率的かつ正確にリードサマライゼーションを行うための汎用プログラムfeatureCountsを開発することである。具体的には、以下の要件を満たすことを目指した。(1) 既存のリードカウンティングツール（summarizeOverlaps、countOverlaps、htseq-countなど）と同等以上の高い精度を達成すること。(2) 既存ツールと比較して、処理速度において桁違いの高速化を実現し、かつ大幅に少ないコンピュータメモリで動作すること。(3) シングルエンドリードとペアエンドリードの両方のシーケンスデータに対応し、多様なシーケンスアプリケーションに適したオプションを提供すること。(4) リファレンススカフォールド数が大規模な不完全アセンブリゲノムや、トランスクリプトームアライメントのように参照配列数が多い状況でも安定して効率的に動作すること。これらの目的を達成することで、リードサマライゼーション工程における計算上のボトルネックを解消し、下流の統計解析パイプラインへのシームレスな統合を可能にする、実用性の高いツールを提供することを目指した。
結果
RNA-seq遺伝子レベルカウント（シングルエンド）における効率性: SEQC (SEquencing Quality Control) データを用いたシングルエンドリードの遺伝子レベルカウントにおいて、featureCountsとsummarizeOverlapsは全25,702遺伝子で完全に一致するカウント数（4,385,354リード）を生成した。htseq-countはわずかに少ないリード数（4,385,207リード）をカウントしたが、これはGFF (General Feature Format) アノテーションの右端座標の解釈の違いに起因することが判明した。GFF仕様では開始位置と終了位置は両方とも含まれるとされているため、featureCountsとsummarizeOverlapsの解釈が正しいと考えられた。この差異を補正した後、3つのツールは完全に一致した。featureCountsは1.0分・16 MBという極めて高速かつ低メモリで処理を完了した。これは、summarizeOverlaps（全ゲノム一括処理で12.1分・3,400 MB）、summarizeOverlaps（染色体ごとの処理で41.7分・661 MB）、htseq-count（22.7分・101 MB）と比較して、10倍以上高速であり、メモリ使用量は桁違いに少なかった (Table 1)。この評価には n=6.8Mペアのリードが使用された。
RNA-seqペアエンドフラグメントカウントにおける精度と効率性: ペアエンドリードを用いたフラグメントカウントにおいて、featureCountsは4,796,948フラグメントを1.0分・16 MBで集計した。summarizeOverlapsは両端がマップされたフラグメントのみを要求するため、カウント数が3,942,439フラグメントと少なく、片端のみがマップされたフラグメントも回収できるfeatureCountsの方が下流解析に有用であることが示唆された。多くの主要アライナー（Subread、Subjunc、Langmead et al. GenomeBiol 2009、TopHat）は片端のみがマップされたフラグメントも報告しており、これらのフラグメントをカウントすることは下流解析に利益をもたらすと考えられる。htseq-count（4,769,913フラグメント）との差分は、featureCountsがオーバーラップの大きさ（リード数ベース）で曖昧なフラグメントを解決する戦略によるものであった。例えば、あるフラグメントが2つの遺伝子にオーバーラップする場合、featureCountsは両方のリードがオーバーラップする遺伝子にフラグメントを割り当てるが、htseq-countはこれを曖昧として割り当てない。Venn解析により、featureCountsが追加でカウントしたフラグメントの86%以上は、少なくとも100個の曖昧さのないフラグメントを持つ発現遺伝子に割り当てられ、誤割り当ては0.2%未満であることが確認された (Figure 2b)。この結果は、featureCountsがより多くの関連フラグメントを正確に回収できることを示している。
ChIP-seq特徴レベルカウントにおける性能比較: H3K27me3 ChIP-seqデータ（広範な遺伝子ボディ+3 kb上流領域）を用いた特徴レベルカウントにおいて、featureCountsとcountOverlaps (IRanges) は完全に一致するカウント数（5,392,155フラグメント）を生成した。このアプリケーションでは、フラグメントが複数の遺伝子にオーバーラップする場合、複数回カウントされるべきである。featureCountsは0.9分・4 MBで処理を完了し、countOverlaps（全ゲノム一括処理で24.4分・7,000 MB）、countOverlaps（染色体ごとの処理で36.6分・783 MB）と比較して約5倍高速であり、メモリ使用量は100倍以上低減された (Table 2)。htseq-count（4,978,050フラグメント、unionモード）とcoverageBED（5,366,902フラグメント）はわずかに少ないカウント数を示し、これは各ツールのマルチオーバーラップ処理やペアエンド処理戦略の違いに起因する。coverageBEDは、フラグメント全体を特徴に割り当てるために、各フラグメントの最初のリードのみを使用していたため、カウント数が少なかった。htseq-countは、マルチオーバーラップフラグメントをカウントしないため、7〜8%少ないフラグメントをカウントした。この評価には n=15Mペアのリードが使用された。
大規模リファレンススカフォールド環境下での性能維持: Budgerigarゲノム（2,850スカフォールド、16,204遺伝子、153,724エキソン）のシミュレーションデータを用いた評価でも、featureCountsは0.6分・15 MBで7,924,065リードをsummarizeOverlapsと完全に一致させつつ集計した。この環境下でも、summarizeOverlaps（全ゲノム一括処理で12.6分・2,400 MB）、summarizeOverlaps（スカフォールドごとの処理で53.3分・262 MB）、htseq-count（12.1分・78 MB）と比較して、featureCountsの優位性は維持された (Table 3)。これは、featureCountsのChromosome hashing技術が大規模なリファレンスシーケンス数に対して効果的に機能することを示している。この評価には n=800万個のシミュレートされたシングルエンドリードが使用された。
理論計算量と実測性能の相関: featureCountsの時間計算量はO(f log f + r·sqrt(k1))であり、リード数rに対してサブリニアに近い特徴検索コストを実現している。これは、実測された桁違いの高速化の根拠となっている。空間計算量はO(f + b1)であり、特徴数fとビン数b1に対して線形に増加する (Table 4)。featureCountsは、比較対象のアルゴリズムの中で最も低い時間計算量と空間計算量を持つことが理論的に裏付けられた。htseq-countの赤黒木検索アルゴリズムは、featureCountsが使用するビンあたりの特徴数の平方根よりもlog fが大きいため、より高い複雑性を持つ。countOverlapsとsummarizeOverlapsはリードをマッピング位置でソートするため、r log rの項が導入され、特にコストが高い。この理論的優位性は、C言語による高効率な実装と直接的なメモリ操作によってさらに強化されている。
考察/結論
先行研究との違い: featureCountsは、リードカウンティングという一見単純な工程に対し、chromosome hashingとfeature blockingの組み合わせというアルゴリズム的革新を導入した。これにより、summarizeOverlapsやcountOverlapsと同等の高い精度を維持しつつ、処理速度とメモリ使用量を1〜2桁改善した点で、これまでの先行研究と異なり、明確な優位性を持つ。特に、SEQC (SEquencing Quality Control) のような大規模なRNA-seqデータセットでも、16 MBという極めて低いメモリフットプリントで1分以内に処理を完了できることは、ノートブックPCレベルでも全ゲノム解析を可能にする実用的な価値が高い。また、マルチスレッド処理をサポートする唯一の既存リードカウントツールである点も、大規模シーケンス研究におけるデータ要約に特に有用である。
新規性: 本研究で初めて、染色体ハッシングと特徴ブロッキングという二つの効率的な技術を組み合わせることで、次世代シーケンスデータのリードサマライゼーションにおける計算効率を劇的に向上させることに成功した。これは、これまで報告されていないアプローチであり、リードカウンティングのボトルネックを解消する新規な解決策を提供する。さらに、本プログラムは、マッピング品質スコア、フラグメントマッパビリティ、フラグメント長、ストランド特異性、キメラ性など、さまざまな指標を用いてリードをフィルタリングするオプションや、マルチマッピングリードのカウントを制御するオプションなど、他のプログラムではサポートされていない多くの有用な機能を提供する。
臨床応用: featureCountsは、SubreadパッケージのUnixコマンドとして、またRsubreadパッケージとしてRから直接呼び出せるため、Robinson et al. Bioinformatics 2010、limma、voomなどのBioconductor統計解析パイプラインへシームレスに接続できる。これにより、ベンチサイエンスの研究者にとって解析の障壁が低い。本ツールは、大規模コホートのRNA-seqデータ（TCGA、GTExなど）の再解析、シングルセル擬似バルクプロファイリングの前処理、ChIP-seqによるエピゲノム定量など、現代の癌ゲノミクスにおけるカウントマトリックス生成のデファクトスタンダードツールとして広く採用されている。これらの応用は、臨床現場での診断や治療法開発に繋がるトランスレーショナルな研究に貢献する。
残された課題: 今後の検討課題としては、(1) トランスクリプトアイソフォームレベルの定量（Li et al. BMCBioinformatics 2011、Salmon、kallistoなどの領域）は本ツールの主要な対象外であること、(2) UMI (Unique Molecular Identifier) 集計やシングルセル専用の最適化は未実装であること、(3) NanoporeやPacBio Iso-Seqなどのロングリードデータへの直接的な対応は限定的であることが挙げられる。これらの課題は、今後の研究で解決されるべき点である。本論文は被引用数1万を超えるバイオインフォマティクス分野のサイテーションクラシックとなり、Subread/Rsubreadエコシステムの基盤を成している。
方法
featureCountsはC言語で実装され、SubreadパッケージのUnixコマンドとして提供されている。また、R言語のラッパー（Rsubreadパッケージ、Bioconductor）も提供されており、R環境からのアクセスも可能である。入力ファイルは、アラインされたリードデータを含むSAM (Sequence Alignment/Map) またはBAM (Binary Alignment/Map) 形式のファイル（自動判別）と、ゲノム特徴のアノテーション情報を含むGFF (General Feature Format) またはSAF (Simplified Annotation Format) 形式のファイルである。SAF形式は、feature ID、chromosome、start、end、strandの5列からなる簡素な形式であり、リード定量に必要最小限の情報を提供する。
featureCountsのアルゴリズムは、主に以下の2つの革新的な技術に基づいている。
Chromosome hashing: 最初に、リファレンスシーケンス名のハッシュテーブルを生成する。これにより、SAMファイルとGFFアノテーションファイル間でリファレンス名（染色体名やコンティグ名）のマッチングを高速化する。この技術は、特にスカフォールド数が膨大な不完全アセンブリゲノムや、トランスクリプトームアライメントのように参照配列数が多い場合に有効である。マッチング後、各リファレンスシーケンスごとに独立して解析を進めることができる。
Feature blocking (二階層構造): 各リファレンスシーケンス内の特徴は、まず開始位置でソートされる。次に、リファレンスシーケンスは重複しない128 kbのビンに分割され、特徴はその開始位置に基づいてビンに割り当てられる。各ビン内では、連続する特徴がsqrt(ビン内の特徴数)個ずつのブロックにグループ化される（図1）。この階層的なデータ構造（ビン内のブロック内の特徴）は、クエリリードとオーバーラップする可能性のあるゲノム領域を迅速に絞り込むことで、リードの割り当てを高速化する。クエリリードは、まずゲノムビンと比較され、次にオーバーラップするビン内の特徴ブロックと比較され、最後にオーバーラップするブロック内の特徴と比較される。この階層的検索により、時間計算量はO(f log f + r·sqrt(k1))となる（f=特徴数、r=リード数、k1=ビン内の特徴数）。
マルチオーバーラップリードの処理: 複数の特徴にオーバーラップするリードやフラグメントの処理は、実験タイプに応じてオプションで制御できる。RNA-seq実験では、単一のフラグメントは単一の遺伝子に由来するため、複数の遺伝子にオーバーラップするリードはカウントから除外することが推奨される。一方、ChIP-seq実験では、エピジェネティック修飾が複数のオーバーラップする遺伝子の生物学的機能を調節する可能性があるため、マルチオーバーラップリードを各特徴にカウントすることが推奨される。また、マルチスレッド処理にも対応しており、大規模データセットでのさらなる速度向上が可能である。
検証は以下のデータセットで実施された。
- RNA-seqデータ: SEQC (SEquencing Quality Control) プロジェクト（MAQC-III）のRNA-seqデータ（n=6.8Mペア × 101 bpリード、Illumina HiSeq 2000で生成）。ヒトゲノムGRCh37 RefSeq build 37.2（25,702遺伝子/225,071エキソン）にSubjuncアライナーでアラインされた。
- ChIP-seqデータ: Pal et al. (2013)のH3K27me3 ChIP-seqデータ（n=15Mペア × 35 bpリード、Illumina Genome Analyzer IIxで生成）。マウスゲノムmm9にSubreadアライナーでアラインされた。
- 大規模リファレンススカフォールドデータ: Assemblathon 2プロジェクトで生成されたBudgerigarゲノムアセンブリ（16,204遺伝子/153,724エキソン/2,850スカフォールド）からシミュレートされたn=800万個の100 bpシングルエンドリード。
比較対象ツールは、BioconductorのsummarizeOverlaps（GenomicRangesパッケージ）とcountOverlaps（IRangesパッケージ）、HTSeqのhtseq-count、およびBEDToolsのcoverageBEDである。すべてのプログラムは、64 AMD Opteron 2.3GHz CPUと512 GBメモリを搭載したHP Bladeスーパーコンピュータ上で、シングルCPU（マルチスレッドなし）で実行された。使用されたソフトウェアパッケージは、Subread 1.4.2、Rsubread 1.12.2、GenomicRanges 1.12.5、IRanges 1.18.4、htseq-count 0.5.4p3、およびBEDTools 2.17.0である。統計解析には、リードカウントデータに対して差次的発現解析に一般的に用いられるRobinson et al. Bioinformatics 2010パッケージが使用されることを想定している。本研究では特定の細胞株やマウス系統は使用せず、既存の公開データセットを用いた計算機科学的評価に重点を置いた。統計的有意性の評価には、各ツールのカウント結果の比較にFisher’s exact testが用いられた。