implement rtgtools for somatic small bench

bin/merge_reports.py

@@ -260,13 +260,16 @@ def main(args=None):
 			summ_table2.reset_index(drop=True, inplace=True)
 			summ_table2.to_csv(args.output + ".regions.csv", index=False)
 
+		elif args.bench == "rtgtools":
+			summ_table = get_rtgtools_resuls(args.inputs)
+
 		elif args.bench == "truvari":
 			summ_table = get_truvari_resuls(args.inputs)
 
 		elif args.bench == "svbenchmark":
 			summ_table = get_svbenchmark_resuls(args.inputs)
 		else:
-			raise ValueError('Only truvari/svbenchmark/sompy results can be merged for somatic analysis!!')
+			raise ValueError('Only truvari/svbenchmark/sompy/rtgtools results can be merged for somatic analysis!!')
 
 		summ_table.reset_index(drop=True, inplace=True)
 		summ_table.to_csv(args.output + ".summary.csv", index=False)

conf/modules.config

@@ -276,7 +276,17 @@ process {
 
     withName: "RTGTOOLS_VCFEVAL" {
         ext.prefix = {"${meta.id}.${params.truth_id}.${meta.caller}"}
-        //ext.args   = {"--ref-overlap --all-record --output-mode ga4gh"}
+        ext.args   = {"--all-record"}
+        publishDir = [
+            path: {"${params.outdir}/${params.variant_type}/${meta.id}/benchmarks/rtgtools"},
+            pattern: "*{.vcf.gz,vcf.gz.tbi,tsv.gz,txt}",
+            mode: params.publish_dir_mode
+        ]
+    }
+    // squash-ploidy is necessary to be able to match het-hom changes
+    withName: "RTGTOOLS_VCFEVAL_SOMATIC" {
+        ext.prefix = {"${meta.id}.${params.truth_id}.${meta.caller}"}
+        ext.args   = {"--all-record --squash-ploidy"}
         publishDir = [
             path: {"${params.outdir}/${params.variant_type}/${meta.id}/benchmarks/rtgtools"},
             pattern: "*{.vcf.gz,vcf.gz.tbi,tsv.gz,txt}",

conf/tests/somatic_indel.config

@@ -28,11 +28,11 @@ params {
     genome          = 'GRCh38'
     analysis        = 'somatic'
     truth_id        = "SEQC2"
-    method          = 'sompy'
+    method          = 'rtgtools'
     preprocess      = "filter_contigs"
     variant_type    = "indel"
 
     truth_vcf       = "${params.test_data_base}/testdata/hg38/truth/somatic/sINDEL_truth_set_v1.0.chr21.vcf.gz"
-    regions_bed     = "${params.test_data_base}/testdata/hg38/truth/somatic/high-confidence_sINDEL_in_HC_regions_v1.2.chr21.vcf.gz"
+    regions_bed     = ""
     rename_chr      = "${projectDir}/assets/rename_contigs/grch37_grch38.txt"
 }

subworkflows/local/small_somatic_benchmark/main.nf

@@ -3,17 +3,25 @@
 //
 
 include { HAPPY_SOMPY           } from '../../../modules/nf-core/happy/sompy/main'
+include { RTGTOOLS_FORMAT  } from '../../../modules/nf-core/rtgtools/format/main'
+include { RTGTOOLS_VCFEVAL  as RTGTOOLS_VCFEVAL_SOMATIC  } from '../../../modules/nf-core/rtgtools/vcfeval/main'
+include { BCFTOOLS_REHEADER as BCFTOOLS_REHEADER_1       } from '../../../modules/nf-core/bcftools/reheader'
+include { BCFTOOLS_REHEADER as BCFTOOLS_REHEADER_2       } from '../../../modules/nf-core/bcftools/reheader'
+include { BCFTOOLS_REHEADER as BCFTOOLS_REHEADER_3       } from '../../../modules/nf-core/bcftools/reheader'
+include { BCFTOOLS_REHEADER as BCFTOOLS_REHEADER_4       } from '../../../modules/nf-core/bcftools/reheader'
 
 workflow SMALL_SOMATIC_BENCHMARK {
     take:
     input_ch    // channel: [val(meta), test_vcf, test_index, truth_vcf, truth_index, bed]
     fasta       // reference channel [val(meta), ref.fa]
     fai         // reference channel [val(meta), ref.fa.fai]
+    sdf       // reference channel [val(meta), sdf]
 
     main:
 
     versions        = Channel.empty()
     summary_reports = Channel.empty()
+    tagged_variants = Channel.empty()
 
     if (params.method.contains('sompy')){
         // apply sompy for small somatic variant benchmarking
@@ -35,7 +43,100 @@ workflow SMALL_SOMATIC_BENCHMARK {
             .set{ report }
         summary_reports = summary_reports.mix(report)
     }
+
+    if (params.method.contains('rtgtools')){
+
+        if (!params.sdf){
+
+            // Use rtgtools format to generate sdf file if necessary
+            RTGTOOLS_FORMAT(
+                fasta.map { meta, file -> [ meta, file, [], [] ] }
+            )
+            versions = versions.mix(RTGTOOLS_FORMAT.out.versions)
+            sdf = RTGTOOLS_FORMAT.out.sdf
+
+        }
+
+        // apply rtgtools eval method
+        RTGTOOLS_VCFEVAL_SOMATIC(
+            input_ch.map { meta, vcf, tbi, truth_vcf, truth_tbi, bed ->
+                [ meta, vcf, tbi, truth_vcf, truth_tbi, bed, [] ]
+            },
+            sdf
+        )
+        versions = versions.mix(RTGTOOLS_VCFEVAL_SOMATIC.out.versions.first())
+
+        // collect summary reports
+        RTGTOOLS_VCFEVAL_SOMATIC.out.summary
+            .map { _meta, file -> tuple([vartype: params.variant_type] + [benchmark_tool: "rtgtools"], file) }
+            .groupTuple()
+            .set{ report }
+
+        summary_reports = summary_reports.mix(report)
+
+        // reheader benchmarking results properly and tag meta
+        BCFTOOLS_REHEADER_1(
+            RTGTOOLS_VCFEVAL_SOMATIC.out.fn_vcf.map{ meta, vcf ->
+            [ meta, vcf, [], [] ]
+            },
+            fai
+        )
+        versions = versions.mix(BCFTOOLS_REHEADER_1.out.versions.first())
+
+        BCFTOOLS_REHEADER_1.out.vcf
+            .join(BCFTOOLS_REHEADER_1.out.index)
+            .map { _meta, file, index -> tuple([vartype: params.variant_type] + [tag: "FN"] + [id: "rtgtools"], file, index) }
+            .set { vcf_fn }
+
+        BCFTOOLS_REHEADER_2(
+            RTGTOOLS_VCFEVAL_SOMATIC.out.fp_vcf.map{ meta, vcf ->
+            [ meta, vcf, [], [] ]
+            },
+            fai
+        )
+        versions = versions.mix(BCFTOOLS_REHEADER_2.out.versions)
+
+        BCFTOOLS_REHEADER_2.out.vcf
+            .join(BCFTOOLS_REHEADER_2.out.index)
+            .map { _meta, file, index -> tuple([vartype: params.variant_type] + [tag: "FP"] + [id: "rtgtools"], file, index) }
+            .set { vcf_fp }
+
+        BCFTOOLS_REHEADER_3(
+            RTGTOOLS_VCFEVAL_SOMATIC.out.baseline_vcf.map{ meta, vcf ->
+            [ meta, vcf, [], [] ]
+            },
+            fai
+        )
+        versions = versions.mix(BCFTOOLS_REHEADER_3.out.versions)
+
+        BCFTOOLS_REHEADER_3.out.vcf
+            .join(BCFTOOLS_REHEADER_3.out.index)
+            .map { _meta, file, index -> tuple([vartype: params.variant_type] + [tag: "TP_base"] + [id: "rtgtools"], file, index) }
+            .set { vcf_tp_base }
+
+        BCFTOOLS_REHEADER_4(
+            RTGTOOLS_VCFEVAL_SOMATIC.out.tp_vcf.map{ meta, vcf ->
+            [ meta, vcf, [], [] ]
+            },
+            fai
+        )
+        versions = versions.mix(BCFTOOLS_REHEADER_4.out.versions)
+
+        BCFTOOLS_REHEADER_4.out.vcf
+            .join(BCFTOOLS_REHEADER_4.out.index)
+            .map { _meta, file, index -> tuple([vartype: params.variant_type] + [tag: "TP_comp"] + [id: "rtgtools"], file, index) }
+            .set { vcf_tp_comp }
+
+        tagged_variants = tagged_variants.mix(
+            vcf_fn,
+            vcf_fp,
+            vcf_tp_base,
+            vcf_tp_comp
+        )
+    }
+
     emit:
     summary_reports // channel: [val(meta), reports]
+    tagged_variants // channel: [val(meta), vcfs]
     versions        // channel: [versions.yml]
 }

tests/germline_small.nf.test.snap

@@ -449,7 +449,7 @@
                 "test1.HG002.strelka.roc.all.csv.gz:md5,cbf4db06f9cbdd0b3b3a29ee4bfc362b",
                 "test1.HG002.strelka.summary.csv:md5,8a5fff1b336e93cee21b3a81f49ab290",
                 "test1.HG002.strelka.phasing.txt:md5,8ed2c215bc14dc461a10e923d0127347",
-                "test1.HG002.strelka.summary.txt:md5,537e7469d43f94647018e173928f088e",
+                "test1.HG002.strelka.summary.txt:md5,98aa938559029495ee3ac7e99b4e6a67",
                 "test1.strelka.bcftools_stats.txt:md5,32586b4f0851fb4f386534a29fe8658f",
                 "test2.HG002.bcftools.extended.csv:md5,9ecb6d4bd55bc7e1ac18bd6cddaa15d9",
                 "test2.HG002.bcftools.roc.Locations.INDEL.PASS.csv.gz:md5,8200b1ab6ab8d4665967ce24b28df103",
@@ -467,6 +467,6 @@
             "nf-test": "0.9.2",
             "nextflow": "24.10.4"
         },
-        "timestamp": "2025-01-29T16:47:42.84439134"
+        "timestamp": "2025-03-04T09:28:13.60011649"
     }
 }

tests/somatic_indel.nf.test

@@ -7,7 +7,7 @@ nextflow_pipeline {
     tag "indel"
     config "../conf/tests/somatic_indel.config"
 
-    test("Params: --analysis 'somatic' --variant_type 'indel' --method 'sompy'") {
+    test("Params: --analysis 'somatic' --variant_type 'indel' --method 'rtgtools'") {
 
         when {
             params {