A step-by-step guide to DNA sequencing data analysis

Kolabtreeの経験豊富なバイオインフォマティクスコンサルタント、Javier Quilez Oliete博士が、データの読み込みに使うツールやソフトなど、DNAシーケンスデータ分析の総合ガイドを提供します。

はじめに

Deoxyribonucleic acid (DNA) は、生物の遺伝情報の大部分を伝える分子です。 (ウイルスの種類によっては、リボ核酸(RNA)によって遺伝情報が運ばれるものもあります)。 DNA分子の基本単位はヌクレオチド（従来はA、C、G、Tの文字で表される）である。 概念的には、DNAを構成するヌクレオチド（例えば、”GCAAACCAAT “は10塩基のDNA文字列）を読み取ることをDNA塩基配列解析というが、DNA塩基配列解析は、この10塩基のDNA文字列を読み取ることで、DNAの塩基配列を決定する。 現在のシーケンシング技術では、このようなDNAリードを数百万個、適度な時間と比較的低いコストで生成することができる。 参考までに、ヒトゲノム（ゲノムとは生物のDNA分子の完全な集合）の配列決定費用は100ドルの壁を越え、数日で完了する。 これは、10 年で完了し、約 27 億ドルのコストがかかったヒト ゲノム配列決定の最初の取り組みとは対照的です。

高スループットと低コストで DNA の配列を決定するこの能力により、配列決定に基づく手法やアプリケーションがますます多く開発されるようになってきています。 例えば、病気や健康な人のゲノム全体やそのタンパク質コード領域の配列決定（それぞれ全ゲノム配列決定とエクソーム配列決定として知られる2つのアプローチ）により、病気の原因となるDNAの変化が明らかになる。 また、DNAから転写されるRNAの配列決定（RNA-sequencing）により、遺伝子活性を定量化し、異なる条件下（例えば、未治療と治療）でそれがどのように変化するかを調べることができる。 一方、染色体コンフォメーションキャプチャーシーケンスでは、近傍のDNA分子間の相互作用を検出し、細胞内の染色体の空間分布を決定するのに役立つ。

これらおよび他のDNAシーケンスのアプリケーションに共通しているのは、数百万の読み取り配列からなるギガバイト単位のデータセットを作成することである。 したがって、ハイスループットシーケンス（HTS）実験を理解するためには、かなりのデータ解析能力が必要である。 幸いなことに、ほとんどの種類の HTS データについて、専用の計算・統計ツールや比較的標準的な解析ワークフローが存在する。 初期）解析ステップのいくつかはほとんどのシーケンサーデータタイプに共通ですが、より下流の解析はデータの種類や解析の最終目標に依存します。 以下では、HTSデータの解析における基本的なステップの入門編を提供し、一般的なツールも紹介します。

以下のセクションのいくつかは、ショートリードのシーケンス技術（主にイルミナ）から生成されたデータの解析に焦点を当てています。 しかし、より長いリードを生成する新しい技術（例：Oxford Nanopore Technologies、PacBio）が急速に台頭してきています。 ロングリードシーケンスにはいくつかの特殊性（例：高いエラー率）があるため、この種のデータを解析するための特別なツールが開発されている。

Quality control (QC) of raw reads

熱心な解析者は、FASTQファイルから解析を開始する。FASTQ形式は、長い間ショートリードのシーケンスデータを保存する標準となっている。 FASTQファイルには、何百万ものリードの塩基配列と塩基ごとのコール品質が含まれています。 ファイルサイズは実際のリード数に依存しますが、FASTQファイルは一般的に大きく（メガバイトやギガバイトのオーダー）、圧縮されています。 なお、FASTQ ファイルを入力とするほとんどのツールは圧縮された状態で扱うことができるので、ディスク容量を節約するためには、解凍しないことを推奨する。

FastQC は、生リードの QC を実行する最も一般的なツールです。 このツールは、ビジュアルインターフェースまたはプログラムによって実行することができます。 コマンドライン環境が苦手なユーザーには前者の方が便利かもしれませんが、後者は比類ないスケーラビリティと再現性を提供します（数十ファイルに対してツールを手動で実行するのがどれほど退屈でエラーが起こりやすいか考えてみてください）。 いずれにせよ、FastQCの主な出力は、与えられたサンプルからの生シーケンスリードの全体的な品質に関する主要な要約統計情報を報告するHTMLファイルです。 数十のFastQCレポートを1つずつ検査するのは面倒で、サンプル間の比較も複雑になります。 そのため、FastQC からの HTML レポート (およびアダプター トリミングやアライメントなど、ダウンストリームで使用される他のツールからのレポート) を 1 つのレポートに集約した MultiQC を使用することをお勧めします。 残念ながら、FastQC メトリクスに基づく、サンプルを良品質か不良品質かに分類するためのコンセンサス閾値は存在しません。 私が使用しているアプローチは以下の通りです。 同じ手順（DNA抽出、ライブラリー調製など）を経たすべてのサンプルは、品質統計が類似しており、「合格」フラグが大半を占めると予想します。 一部のサンプルが平均より低い品質であったとしても、このことを念頭に置いて下流の解析に使用することにしています。 一方、実験中のすべてのサンプルが複数のメトリクスで体系的に「警告」または「失敗」のフラグを立てた場合（この例を参照）、私は実験中に何か問題があったのではないかと考え（例えば、悪いDNA品質、ライブラリ準備など）、実験の繰り返しを推奨します。 そのためには、リードをトリミングしてテクニカルシーケンスや低品質の末端を取り除く必要があります。 このような配列がリードの真の配列に隣接している場合、アライメント（後述）によりリードがゲノム内の誤った位置にマップされたり、与えられたアライメントの信頼度が低下したりすることがあります。 技術的な配列の他に、生物学的由来の配列がリードの中に多く存在する場合、それらを除去したい場合があります。 例えば、最適でないDNA調製手順により、サンプル中にDNA変換されたリボソームRNA（rRNA）が高い割合で残っていることがあります。 このタイプの核酸がシーケンス実験のターゲットでない限り、rRNAに由来するリードを保持することは、下流のステップの計算負荷を増加させるだけであり、結果を混乱させる可能性があります。

ショートリードシーケンスでは、DNA配列は一度に1ヌクレオチドずつ（厳密にはシーケンスサイクルごとに1ヌクレオチドずつ）決定されます。 言い換えれば、シーケンシングサイクルの数がリードの長さを決定します。 HTSシーケンシング法の課題として、シーケンシングサイクルを重ねるごとに塩基配列の決定精度が低下していくことが知られています。 これは、特にリードの終点に向かうにつれて、塩基あたりのコール品質が全体的に低下することに反映される。 技術的な配列と同様に、低品質の末端を含むリードをアライメントしようとすると、ミスプレースメントやマッピングの質の低下につながることがあります。 このようなツールは、技術的な配列（内部で利用可能なもの、またはユーザーが提供するもの）を除去し、リードの長さを最大化しながら、品質に基づいてリードをトリミングします。 トリミング後に短くなりすぎたリードは破棄されます（例えば、<36ヌクレオチドなど過度に短いリードは、ゲノム内の複数の部位にマッピングされる可能性が高く、アライメントステップを複雑にしています）。 廃棄されたリードの割合が高い場合は、データの品質が低い可能性があるため、トリミングに耐えたリードの割合に注目するとよいでしょう。

最後に、トリミングしたリードに対してFastQCを再実行し、このステップが有効であること、QC指標が系統的に改善されたことを確認します。

アライメント

例外（例えばde novoアセンブリ）を除き、ほとんどのHTSデータタイプや用途ではアライメントは次のステップです（マッピングとも呼ばれる）。 リードのアライメントは、リードの配列が由来するゲノム上の位置（通常、染色体：始端として表現されます）を決定することからなります。 したがって、このステップでは、リードをアライメント/マッピングするための参照配列が必要となります。 1 つは、配列決定された DNA が由来する種です。 高品質の参照配列が利用できる種の数は増えていますが、あまり研究されていない生物では、まだそうでない場合があります。 このような場合、参照ゲノムが利用可能な進化的に近い生物種にリードをアラインメントすることが考えられます。 例えば、コヨーテのゲノムは参照配列がないため、近縁種のイヌのゲノムを用いてリードのアライメントを行うことができます。 同様に、より質の高い参照配列が存在する近縁種にリードをアライメントしたい場合もあります。 例えば、テナガザルのゲノムは公開されていますが、これは数千の断片に分かれており、数十本の染色体からなるゲノムの構成を完全に再現しているわけではありません。 重要なのは、あるアラインメントの座標がバージョン間で異なる可能性があることである。 例えば、UCSC Genome Browserでは、ヒトゲノムの複数のバージョンを見ることができます。 どのような生物種であっても、アセンブリが完全にリリースされたら、最新のバージョンに移行することを強く推奨します。 これは、既存の結果が古いバージョンと比較されるため、移行中に多少の迷惑をかけるかもしれませんが、長い目で見れば得策です。

さらに、シーケンスデータの種類も重要です。 DNA-seq、ChIP-seq、またはHi-Cプロトコルから生成されたリードは、ゲノム参照配列にアライメントされます。 一方、DNAから転写されたRNAはさらにmRNAに加工される（イントロンが除去される）ため、多くのRNA-seqリードはゲノム参照配列にアラインメントされない。 その代わりに、トランスクリプトーム参照配列にアライメントするか、ゲノム配列を参照にする場合は分割を考慮したアライナー（下記参照）を使用する必要があります。 これに関連して、参照配列のアノテーション、つまり、遺伝子、転写産物、セントロメアなどの座標を持つデータベースのソースの選択があります。

多くの短鎖配列アライメントツールが開発されています（短鎖配列アライメントのセクションを参照）。 これらのツールをレビューすることは、この記事の範囲を超えている（これらのツールの背後にあるアルゴリズムの詳細については、ここで見つけることができます）。 私の経験では、Bowtie2、BWA、HISAT2、Minimap2、STAR、TopHatが最もよく使われているようだ。 私のお勧めは、HTSデータの種類やアプリケーション、コミュニティでの受け入れ状況、ドキュメントの質、ユーザー数などの重要な要素を考慮してアライナーを選択することです。 例えば、RNA-seqをゲノムにマッピングする場合、エクソン-エクソン結合を意識したSTARやBowtie2のようなアライナーが必要です。

ほとんどのマッパーに共通するのは、実際のアライメントが行われる前に、参照として使用される配列にインデックスを付ける必要があることです。 このステップは時間がかかるかもしれませんが、各参照配列に対して 1 回だけ行う必要があります。 ほとんどのマッパーは、SAM/BAMファイル（BAMファイルはSAMファイルのバイナリバージョンです）にアライメントを保存します。 アラインメントはシーケンスデータの解析の中で最も計算と時間がかかるステップの一つであり、SAM/BAMファイルは重い（ギガバイトオーダー）です。 したがって、妥当な時間でアライメントを実行し、結果を保存するために必要なリソース（以下の最終セクションを参照）があるかどうかを確認することが重要です。 同様に、BAMファイルのサイズとバイナリ形式のため、テキストエディタで開くことは避け、UnixコマンドやSAMtoolsのような専用ツールを使用します。

アラインメントから

この時点では、各HTSデータの種類やアプリケーションが異なるため、アラインメント後の共通のステップはないといってよいでしょう。

DNA-seqデータの一般的なダウンストリーム解析は、バリアントコール、つまり、ゲノム参照と個人間で異なるゲノムの位置を特定することです。 このアプリケーションのための一般的な分析フレームワークは、一塩基多型（SNP）または小さな挿入/欠失（indel）用のGATKです（図2）。 DNAの大きな塊からなるバリアント（構造バリアントとも呼ばれる）には、専用のコールメソッドが必要です（包括的な比較はこちらの記事をご覧ください）。 アライナーと同様、バリアントの種類 (SNP、インデル、構造バリアント)、コミュニティによる受け入れ、ドキュメントの品質、ユーザー数などの重要な要因を考慮して、適切なツールを選択することをお勧めします。 歴史的には、リードは参照配列にアラインされる必要があり、その後、与えられた遺伝子または転写物にアラインされたリードの数が、その発現レベルを定量化するためのプロキシとして使用されました。 このアラインメントと定量化の手法は、Cufflinks、RSEM、featureCountsなどのツールで実行されています。 しかし、この方法は、KallistoやSalmonのようなソフトウェアに実装された新しい方法に取って代わられつつあります。 概念的には、このようなツールでは、リードの全塩基配列を参照配列にアライメントする必要はありません。 その代わり、あるリードがある転写産物に由来していると確信できる程度の塩基配列のアライメントが必要である。 簡単に言えば、アライメントと定量を1つのステップで済ませることができるのです。 このアプローチは擬似マッピングと呼ばれ、遺伝子発現の定量化の速度を大幅に向上させます。 一方、フルアライメントが必要なアプリケーション（例：RNA-seqデータからのバリアントコール）には、pseudo-mappingは適さないことに留意してください。

シーケンスベースのアプリケーションにおけるダウンストリーム解析ステップと必要なツールの違いのもう一つの例が、ChIP-seqです。 このような手法で生成されたリードは、ターゲットタンパク質が結合している場所を示すリードが著しく過剰な領域をゲノムから検出するピークコールに使用されます。 ピークコーリングツールはいくつか存在するが、本書ではそれらを調査した。 最後の例として、Hi-Cデータについて述べます。Hi-Cデータは、相互作用行列を決定するツールの入力として使用され、これらからゲノムの3次元的な特徴を決定します。 この記事の範囲を超えて、すべてのシーケンスベースのアッセイについてコメントします（比較的完全なリストは、この記事を参照してください）

Before you start…

この記事の残りの部分は、HTSデータの解析のステップとして厳密に考慮されていないかもしれない、ほとんど無視されている側面について触れています。 表1

考えてみよう	提案する行動
分析に必要なサンプルの情報はすべて揃っていますか？	実験のメタデータを系統的に収集する
サンプルを明確に識別できるか	各サンプルに固有の識別子を付与するシステムを構築する
データと結果はどこにあるか
		Data of Results in Where are you?	データの構造化、階層化
複数のサンプルをシームレスに処理できますか？	コードの拡張性、並列化、自動構成、モジュール化
あなたや他の誰かが結果を再現できますか？	Documentment your code and procedures!

前述のように、HTS 生データとその解析中に生成されるファイルのいくつかはギガバイトオーダーであり、数十サンプルを含むプロジェクトがテラバイト単位のストレージを必要とすることは例外的なことではありません。 その上、HTSデータの解析のいくつかのステップは計算集約的である（例えば、アライメント）。 しかし、HTS データの解析に必要なストレージやコンピューティングインフラは重要な検討事項であり、見過ごされたり、議論されなかったりしがちである。 例として、最近の分析の一部として、我々は、phenome-wide association analysis (PheWAS)を実行している数十の発表された論文をレビューした。 最新のPheWASでは、100-1,000の遺伝的変異と表現型の両方を解析するため、重要なデータストレージと計算能力が必要とされます。 しかし、我々がレビューした論文の中で、PheWASの解析に必要なインフラについてコメントしているものはほぼ皆無でした。 当然のことながら、私が推奨するのは、あなたが直面するストレージとコンピューティングの要件を前もって計画し、それをコミュニティと共有することです」

DNA 配列データの解析でお困りですか？ Kolabtreeでフリーランスのバイオインフォマティシャンやゲノミクスの専門家と連絡を取ってみてください。