このトピックでは、E-HPC クラスターを使用して BWA、GATK、Samtools を実行し、遺伝子シーケンシングを行う方法について説明します。
背景情報
生命科学における遺伝子シーケンシング技術の急速な発展により、遺伝子シーケンスデータの量は指数関数的に増加しています。この大量のデータを処理するには、ホモロジー検索、アライメント、変異検出などのタスクのために、大規模なデータ処理と並列計算が必要です。高性能コンピューティング (HPC) は、これらのタスクに必要な強力な計算能力を提供します。HPC は、さまざまなスケジューラを使用して並列性を向上させ、GPU を使用して計算を高速化します。
このトピックでは、典型的な全ゲノムシーケンシング (WGS) ワークフローを例として使用します。GATK ソフトウェアを使用したヒト WGS の一般的なプロセスを示します。実際のバイオインフォマティクス解析では、特定のニーズに基づいてプロセスを調整する必要があります。たとえば、変異検出のための品質管理とフィルタリングを追加する必要がある場合があります。
遺伝子シーケンシングでは、BWA を使用してインデックスを構築し、レコードをアライメントできます。次に、Samtools を使用してアライメントレコードをソートします。最後に、GATK を使用して重複シーケンスを除去し、塩基クオリティスコアを再キャリブレーションし、変異を検出します。
BWA (Burrows-Wheeler Alignment Tool) は、DNA シーケンスをヒトゲノムなどのリファレンスゲノムにマッピングするソフトウェアです。
GATK (The Genome Analysis Toolkit) は、第二世代の再シーケンシングデータを解析するためのソフトウェアツールキットです。主に、重複シーケンスの除去、塩基クオリティスコアの再キャリブレーション、変異検出に使用されます。
Samtools は、SAM および BAM フォーマットを処理するためのツールです。バイナリファイルの表示、ファイル形式の変換、ファイルのソートとマージが可能です。また、SAM ファイルのフラグとタグの情報を使用して、アライメント結果を要約します。
事前準備
E-HPC クラスターを作成します。詳細については、「Standard Edition クラスターの作成」をご参照ください。
次の表に、クラスター構成の例を示します。
設定項目
構成
シリーズ
Standard Edition
デプロイモード
パブリッククラウドクラスター
クラスタータイプ
SLURM
ノード構成
次の仕様のマスターノード 1 台、ログインノード 1 台、コンピュートノード 1 台が含まれます。
マスターノード:ecs.c7.large インスタンスタイプを使用し、2 vCPU と 4 GiB のメモリを備えています。
ログインノード:ecs.c7.large インスタンスタイプを使用し、2 vCPU と 4 GiB のメモリを備えています。
コンピュートノード:ecs.g7.16xlarge インスタンスタイプを使用し、64 vCPU と 256 GiB のメモリを備えています。
クラスターイメージ
CentOS 7.6 パブリックイメージ
クラスターユーザーを作成します。詳細については、「ユーザー管理」をご参照ください。
クラスターユーザーは、クラスターにログインしてソフトウェアのコンパイル、ジョブのサブミット、その他の操作を実行するために使用されます。次の例は、このトピックで作成されるユーザーを示しています。
ユーザー名:testuser
ユーザーグループ:sudo 権限グループ
ソフトウェアをインストールします。
次の表に、インストールするソフトウェアの詳細を示します。
ソフトウェア名
バージョン
ダウンロード URL
BWA
4.2.0.0
GATK
0.7.17
Samtools
1.14
手順1:E-HPC クラスターへの接続
E-HPC クラスターに接続します。詳細については、「クラスターへの接続」をご参照ください。
次のコマンドを実行して、ジョブ実行ディレクトリを作成します。このトピックでは、
/home/data/humanを例として使用します。説明データをダウンロードして処理するための十分なディスク領域を確保するために、クラスターにマウントされている
/home共有ストレージディレクトリにジョブ実行ディレクトリを作成してください。sudo mkdir -p /home/data/human
手順2:データのダウンロードと前処理
次のコマンドを実行して、
data_pre.slurmという名前のジョブスクリプトファイルを作成して開きます。sudo vim data_pre.slurmスクリプトは次のとおりです。
#!/bin/bash -l #SBATCH -J data_pre #SBATCH --partition comp #SBATCH -N 1 #SBATCH -n 1 #SBATCH --output=data_pre.log cd /home/data/human/ # ヒトのリファレンスゲノムと 2 つのヒト遺伝子シーケンシングサンプルをダウンロードして解凍します。 wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_human_g1k_v37.fasta wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_1000G_phase1.indels.b37.vcf wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_Mills_and_1000G_gold_standard.indels.b37.vcf wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_dbsnp_138.b37.vcf.zip unzip b37_dbsnp_138.b37.vcf.zip wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/gatk-examples_example1_NA20318_ERR250971_1.filt.fastq.gz gunzip gatk-examples_example1_NA20318_ERR250971_1.filt.fastq.gz wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/gatk-examples_example1_NA20318_ERR250971_2.filt.fastq.gz gunzip gatk-examples_example1_NA20318_ERR250971_2.filt.fastq.gz # bgzip を使用してファイルを圧縮します。 bgzip -f b37_1000G_phase1.indels.b37.vcf bgzip -f b37_Mills_and_1000G_gold_standard.indels.b37.vcf bgzip -f b37_dbsnp_138.b37.vcf # tabix を使用してインデックスを構築します。 tabix -p vcf b37_1000G_phase1.indels.b37.vcf.gz tabix -p vcf b37_Mills_and_1000G_gold_standard.indels.b37.vcf.gz tabix -p vcf b37_dbsnp_138.b37.vcf.gz次のコマンドを実行してジョブをサブミットします。
sbatch data_pre.slurm次の出力は、ジョブが正常にサブミットされたことを示しています。

次のコマンドを実行してジョブステータスを確認します。
squeueジョブが完了したら、次のコマンドを実行して
data_pre.logファイルで結果を表示します。cat data_pre.log
手順3:リファレンスゲノムインデックスの構築
次のコマンドを実行して、
bwa_index.slurmという名前のジョブスクリプトファイルを作成して開きます。sudo vim bwa_index.slurmスクリプトは次のとおりです。
#!/bin/bash -l #SBATCH -J bwa_index #SBATCH --partition comp #SBATCH -N 1 #SBATCH -n 1 #SBATCH --output=bwa_index.log cd /home/data/human/ # BWA アライメントのためにリファレンスゲノムのインデックスデータベースを構築します。 time bwa index b37_human_g1k_v37.fasta # FASTA シーケンスフォーマットのインデックスを作成します。 samtools faidx b37_human_g1k_v37.fasta次のコマンドを実行してジョブをサブミットします。
sbatch --dependency=afterok:jobid1 bwa_index.slurm次のコマンドを実行してジョブステータスを確認します。
squeueジョブが完了したら、次のコマンドを実行して
bwa_index.logファイルで結果を表示します。cat bwa_index.log
手順4:遺伝子シーケンシング解析の実行
次のコマンドを実行して、
wgs_main_process.slurmという名前のジョブスクリプトファイルを作成して開きます。sudo vim wgs_main_process.slurmスクリプトは次のとおりです。
重要6 行目の
#SBATCH -w compute000で、compute000をご利用のコンピュートノードの名前に置き換えてください。#!/bin/bash -l #SBATCH -J wgs_main_process #SBATCH --partition comp #SBATCH -N 1 #SBATCH -n 64 #SBATCH -w compute000 #SBATCH --output=wgs_main_process.log cd /home/data/human/ # サンプルシーケンシングデータのリードをヒトのリファレンスゲノムにアライメントし、出力ファイルを BAM フォーマットに変換してディスク領域を節約します。 time bwa mem -t 64 \ -R '@RG\tID:ehpc\tPL:illumina\tLB:library\tSM:b37' \ b37_human_g1k_v37.fasta \ gatk-examples_example1_NA20318_ERR250971_1.filt.fastq \ gatk-examples_example1_NA20318_ERR250971_2.filt.fastq | samtools view -S -b - > ERR250971.bam # アライメントレコードを昇順でソートします。 time samtools sort -@ 64 -O bam -o ERR250971.sorted.bam ERR250971.bam # DNA シーケンスの PCR 増幅によって生成された重複リードを除去します。 time gatk MarkDuplicates \ -I ERR250971.sorted.bam -O ERR250971.markdup.bam \ -M ERR250971.sorted.markdup_metrics.txt # markdup シーケンシング BAM ファイルのインデックスを構築します。 samtools index ERR250971.markdup.bam # リファレンスゲノムの dict ファイルを生成します。 time gatk CreateSequenceDictionary \ -R b37_human_g1k_v37.fasta \ -O b37_human_g1k_v37.dict # 既存の SNP および indel データベースを使用してモデルを構築し、再キャリブレーションテーブルを作成します。 time gatk BaseRecalibrator \ -R b37_human_g1k_v37.fasta \ -I ERR250971.markdup.bam \ --known-sites b37_1000G_phase1.indels.b37.vcf.gz \ --known-sites b37_Mills_and_1000G_gold_standard.indels.b37.vcf.gz \ --known-sites b37_dbsnp_138.b37.vcf.gz -O ERR250971.BQSR.table # 元の塩基クオリティスコアを調整します。 time gatk ApplyBQSR \ --bqsr-recal-file ERR250971.BQSR.table \ -R b37_human_g1k_v37.fasta \ -I ERR250971.markdup.bam \ -O ERR250971.BQSR.bam # BQSR シーケンシング BAM ファイルのインデックスを構築します。 time samtools index ERR250971.BQSR.bam # 変異検出結果を VCF ファイルに出力します。 time gatk HaplotypeCaller \ -R b37_human_g1k_v37.fasta \ -I ERR250971.BQSR.bam \ -O ERR250971.HC.vcf.gz次のコマンドを実行してジョブをサブミットします。
sbatch --dependency=afterok:jobid2 wgs_main_process.slurm次のコマンドを実行してジョブステータスを確認します。
squeueジョブが完了したら、次のコマンドを実行して
wgs_main_process.logファイルで結果を表示します。cat wgs_main_process.log
シーケンシング時間
このトピックで使用されるコンピュートノードは ecs.g7.16xlarge です。これはデモンストレーションのみを目的としています。全ゲノムシーケンシングプロセス全体には約 6~7 時間かかります。次の表に、参考として各手順のおおよその時間を示します。これらの時間は、最適なパフォーマンスのベンチマークではありません。
プロセス | 機能 | 実行時間 (分) |
bwa index | インデックスの構築 | 54 |
bwa mem | アライメント | 25 |
samtools sort | ソート | 2 |
gatk MarkDuplicates | 重複の除去 | 13 |
gatk CreateSequenceDictionary | 辞書の作成 | 0.15 |
gatk BaseRecalibrator | キャリブレーションテーブルの構築 | 40 |
gatk ApplyBQSR | 再キャリブレーション | 21 |
gatk HaplotypeCaller | 変異の検出 | 211 |