Edited at

cromwell を利用してWorkflows for germline short variant discovery with GATK4を実行する


Cromwellを利用してgatk4-germline-snps-indelsを実行する

https://github.com/gatk-workflows/gatk4-germline-snps-indels

を実行する

このワークフローを実行している環境は


  • OS: macOS Mojave Version 10.14.2

  • CPU: Core i5 (4 Core)

  • Memory: 32GB

下記の記述内容の続き

cromwell でGATKのワークフローgatk4-data-processing-1.0.1/processing-for-variant-discovery-gatk4.wdlを実行する


手順概要


  • gatk4-germline-snps-indelsをダウンロード

  • 必要なファイルをgoogle storageから手元にダウンロード


    • 下記の2つのファイルからユーティリティーdownload_gs.pyを利用してファイルを取得する

    • haplotypecaller-gvcf-gatk4.hg38.wgs.inputs.json に記載のgs://から始まるファイル

    • gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txt に記載のファイル




  • 入力ファイルを修正


    • gs://部分を加工

    • haplotypecaller-gvcf-gatk4.hg38.wgs.inputs.json


      • 手元で実行するため、小さい方のbamファイルを選択する

      • gvcfを作成する






  • ワークフローの修正


    • 実行時にreferenceファイルのサイズから動的にディスク容量を計算する部分をコメントアウト




gatk4-germline-snps-indelsをダウンロード

https://github.com/gatk-workflows/gatk4-germline-snps-indels/releases

Source Code (zip)をダブルクリックしてダウンロード

gatk4-germline-snps-indels-1.0.4

$ mv ~/Downloads/gatk4-germline-snps-indels-1.0.4 .


必要なファイルをgoogle storageから手元にダウンロード

$ python3 download_gs.py gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs.json

ダウンロードしたgatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txtには解析領域を示す一連のファイルが記載されているので、このファイルに記載されている一連のファイルもダウンロードする

$ python3 download_gs.py gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txt --non_json


入力ファイルを修正

gs://を次のファイルから取り除いて

gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs.json

gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs_local.json

を作成する

$ cat gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs.json | sed 's/gs:\/\///g' > gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs_local.json

gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txtの内部に記載されているファイルについても同様

元のファイルをgatk-test-data/intervals/hg38_wgs_scattered_calling_intervals_bak.txtにリネームして、ファイルの記述からgs://を除いたファイルをgatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txtとして作成する

$ mv gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txt  gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals_bak.txt 

$ cat gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals_bak.txt | sed 's/gs:\/\///g' > gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals.txt


入力用jsonファイルの変更

gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs_local.jsonファイルのサンプルbamはお試し用には大きすぎるので

  "##_COMMENT1": "INPUT BAM",

"#HaplotypeCallerGvcf_GATK4.input_bam": "gatk-test-data/wgs_bam/NA12878_24RG_hg38/NA12878_24RG_small.hg38.bam",
"#HaplotypeCallerGvcf_GATK4.input_bam_index": "gatk-test-data/wgs_bam/NA12878_24RG_hg38/NA12878_24RG_small.hg38.bai",
"HaplotypeCallerGvcf_GATK4.input_bam": "broad-public-datasets/NA12878/NA12878.cram",
"HaplotypeCallerGvcf_GATK4.input_bam_index": "broad-public-datasets/NA12878/NA12878.cram.crai",

小さい方のサンプルを利用するために、前2行のコメントを外し、後2行をコメントアウトする

  "##_COMMENT1": "INPUT BAM",

"HaplotypeCallerGvcf_GATK4.input_bam": "gatk-test-data/wgs_bam/NA12878_24RG_hg38/NA12878_24RG_small.hg38.bam",
"HaplotypeCallerGvcf_GATK4.input_bam_index": "gatk-test-data/wgs_bam/NA12878_24RG_hg38/NA12878_24RG_small.hg38.bai",
"#HaplotypeCallerGvcf_GATK4.input_bam": "broad-public-datasets/NA12878/NA12878.cram",
"#HaplotypeCallerGvcf_GATK4.input_bam_index": "broad-public-datasets/NA12878/NA12878.cram.crai",

次の#HaplotypeCallerGvcf_GATK4.HaplotypeCaller.make_gvcfのコメントを外す

  "##_COMMENT4": "MISCELLANEOUS PARAMETERS",

"HaplotypeCallerGvcf_GATK4.HaplotypeCaller.make_gvcf": "True",

java -Dbackend.providers.Local.config.concurrent-job-limit=2 -jar cromwell.jar run  gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.wdl -i gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs_local.json

しかし、このままでは動作しない

 [error] WorkflowManagerActor Workflow 30e9fce5-4ef2-4974-abf7-fa7496cb7022 failed (during ExecutingWorkflowState): cromwell.engine.workflow.lifecycle.execution.job.preparation.JobPreparationActor$$anonfun$1$$anon$1: Call input and runtime attributes evaluation failed for HaplotypeCaller:

Failed to evaluate input 'disk_size' (reason 1 of 1): [Attempted 1 time(s)] - NoSuchFileException: /Users/xxxx/Documents/cromwell/cromwell-executions/HaplotypeCallerGvcf_GATK4/30e9fce5-4ef2-4974-abf7-fa7496cb7022/call-HaplotypeCaller/shard-14/execution/gatk-test-data/wgs_bam/NA12878_24RG_hg38/NA12878_24RG_small.hg38.bam
Failed to evaluate input 'ref_size' (reason 1 of 1): [Attempted 1 time(s)] - NoSuchFileException: /Users/xxxx/Documents/cromwell/cromwell-executions/HaplotypeCallerGvcf_GATK4/30e9fce5-4ef2-4974-abf7-fa7496cb7022/call-HaplotypeCaller/shard-14/execution/broad-references/hg38/v0/Homo_sapiens_assembly38.fasta
at cromwell.engine.workflow.lifecycle.execution.job.preparation.JobPreparationActor$$anonfun$1.applyOrElse(JobPreparationActor.scala:66)
at cromwell.engine.workflow.lifecycle.execution.job.preparation.JobPreparationActor$$anonfun$1.applyOrElse(JobPreparationActor.scala:62)
at scala.runtime.AbstractPartialFunction.apply(AbstractPartialFunction.scala:34)
at akka.actor.FSM.processEvent(FSM.scala:684)
at akka.actor.FSM.processEvent$(FSM.scala:681)
at cromwell.engine.workflow.lifecycle.execution.job.preparation.JobPreparationActor.processEvent(JobPreparationActor.scala:39)
at akka.actor.FSM.akka$actor$FSM$$processMsg(FSM.scala:678)
at akka.actor.FSM$$anonfun$receive$1.applyOrElse(FSM.scala:672)
at akka.actor.Actor.aroundReceive(Actor.scala:517)
at akka.actor.Actor.aroundReceive$(Actor.scala:515)
at cromwell.engine.workflow.lifecycle.execution.job.preparation.JobPreparationActor.aroundReceive(JobPreparationActor.scala:39)
at akka.actor.ActorCell.receiveMessage(ActorCell.scala:588)
at akka.actor.ActorCell.invoke(ActorCell.scala:557)
at akka.dispatch.Mailbox.processMailbox(Mailbox.scala:258)
at akka.dispatch.Mailbox.run(Mailbox.scala:225)
at akka.dispatch.Mailbox.exec(Mailbox.scala:235)
at akka.dispatch.forkjoin.ForkJoinTask.doExec(ForkJoinTask.java:260)
at akka.dispatch.forkjoin.ForkJoinPool$WorkQueue.runTask(ForkJoinPool.java:1339)
at akka.dispatch.forkjoin.ForkJoinPool.runWorker(ForkJoinPool.java:1979)
at akka.dispatch.forkjoin.ForkJoinWorkerThread.run(ForkJoinWorkerThread.java:107)


ワークフローの修正

wdl内部で必要なディスク容量を計算するためにリファレンスファイルのサイズを求めているが、inputの指定の段階では該当のファイルがコピーされていないためにNoSuchFileException:が発生する

結局wdlファイルに手を入れない限り、このワークフローを動作させられない

そこで、該当領域をコメントアウトしたwdlを作成する


変更箇所について

gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.wdlを修正してgatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4_local.wdlを作成する

103行目の

# TASK DEFINITIONS 以降の一部をコメントアウトしている

オリジナルのwdlの該当行以降を下記で置き換えてgatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4_local.wdlを作成


gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4_local.wdlの一部

# TASK DEFINITIONS

task CramToBamTask {
# Command parameters
File ref_fasta
File ref_fasta_index
File ref_dict
File input_cram
String sample_name

# Runtime parameters
String docker
Int? machine_mem_gb
Int? disk_space_gb
Boolean use_ssd = false
Int? preemptible_attempts

#Float output_bam_size = size(input_cram, "GB") / 0.60
#Float ref_size = size(ref_fasta, "GB") + size(ref_fasta_index, "GB") + size(ref_dict, "GB")
#Int disk_size = ceil(size(input_cram, "GB") + output_bam_size + ref_size) + 20

command {
set -e
set -o pipefail

samtools view -h -T ${ref_fasta} ${input_cram} |
samtools view -b -o ${sample_name}.bam -
samtools index -b ${sample_name}.bam
mv ${sample_name}.bam.bai ${sample_name}.bai
}
runtime {
docker: docker
#memory: select_first([machine_mem_gb, 15]) + " GB"
#disks: "local-disk " + select_first([disk_space_gb, disk_size]) + if use_ssd then " SSD" else " HDD"
preemptibe: preemptible_attempts
}
output {
File output_bam = "${sample_name}.bam"
File output_bai = "${sample_name}.bai"
}
}

# HaplotypeCaller per-sample in GVCF mode
task HaplotypeCaller {
File input_bam
File input_bam_index
File interval_list
String output_filename
File ref_dict
File ref_fasta
File ref_fasta_index
Float? contamination
Boolean make_gvcf

String gatk_path
String? java_options
String java_opt = select_first([java_options, "-XX:GCTimeLimit=50 -XX:GCHeapFreeLimit=10"])

# Runtime parameters
String docker
Int? mem_gb
Int? disk_space_gb
Boolean use_ssd = false
Int? preemptible_attempts

Int machine_mem_gb = select_first([mem_gb, 7])
Int command_mem_gb = machine_mem_gb - 1

#Float ref_size = size(ref_fasta, "GB") + size(ref_fasta_index, "GB") + size(ref_dict, "GB")
#Int disk_size = ceil(size(input_bam, "GB") + ref_size) + 20

command <<<
set -e

${gatk_path} --java-options "-Xmx${command_mem_gb}G ${java_opt}" \
HaplotypeCaller \
-R ${ref_fasta} \
-I ${input_bam} \
-L ${interval_list} \
-O ${output_filename} \
-contamination ${default=0 contamination} ${true="-ERC GVCF" false="" make_gvcf}
>>>

runtime {
docker: docker
#memory: machine_mem_gb + " GB"
#disks: "local-disk " + select_first([disk_space_gb, disk_size]) + if use_ssd then " SSD" else " HDD"
preemptible: select_first([preemptible_attempts, 3])
}

output {
File output_vcf = "${output_filename}"
File output_vcf_index = "${output_filename}.tbi"
}
}
# Merge GVCFs generated per-interval for the same sample
task MergeGVCFs {
Array[File] input_vcfs
Array[File] input_vcfs_indexes
String output_filename

String gatk_path

# Runtime parameters
String docker
Int? mem_gb
Int? disk_space_gb
Boolean use_ssd = false
Int? preemptible_attempts

Int machine_mem_gb = select_first([mem_gb, 3])
Int command_mem_gb = machine_mem_gb - 1

command <<<
set -e

${gatk_path} --java-options "-Xmx${command_mem_gb}G" \
MergeVcfs \
--INPUT ${sep=' --INPUT ' input_vcfs} \
--OUTPUT ${output_filename}
>>>

runtime {
docker: docker
memory: machine_mem_gb + " GB"
#disks: "local-disk " + select_first([disk_space_gb, 100]) + if use_ssd then " SSD" else " HDD"
preemptible: select_first([preemptible_attempts, 3])
}

output {
File output_vcf = "${output_filename}"
File output_vcf_index = "${output_filename}.tbi"
}
}



ワークフローの実行

$ java -Dbackend.providers.Local.config.concurrent-job-limit=2 -jar cromwell.jar run  gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4_local.wdl -i gatk4-germline-snps-indels-1.0.4/haplotypecaller-gvcf-gatk4.hg38.wgs.inputs_local.json

出力ファイル

{

"outputs": {
"HaplotypeCallerGvcf_GATK4.output_vcf": "/Users/xxxx/Documents/cromwell/cromwell-executions/HaplotypeCallerGvcf_GATK4/c87278d5-f092-4179-86d9-3bf5a6b04a11/call-MergeGVCFs/execution/NA12878_24RG_small.hg38.g.vcf.gz",
"HaplotypeCallerGvcf_GATK4.output_vcf_index": "/Users/xxxx/Documents/cromwell/cromwell-executions/HaplotypeCallerGvcf_GATK4/c87278d5-f092-4179-86d9-3bf5a6b04a11/call-MergeGVCFs/execution/NA12878_24RG_small.hg38.g.vcf.gz.tbi"
},
"id": "c87278d5-f092-4179-86d9-3bf5a6b04a11"
}

今回はここまで:smile: