[DRAFT] Add NGSCheckMate in as part of a cram sampleQC subworkflow

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - [#1113](https://github.com/nf-core/sarek/pull/1113) - Adding CNVkit genemetrics module
 - [#1193](https://github.com/nf-core/sarek/pull/1193) - Adding support for Sentieon's DnaScope for germline variant-calling including joint-germline
+- [#1252](https://github.com/nf-core/sarek/pull/1252) - Added NGSCheckMate tool for checking that samples come from the same individual
 - [#1271](https://github.com/nf-core/sarek/pull/1271) - Back to dev
 - [#1290](https://github.com/nf-core/sarek/pull/1290) - Add nf-test for whole pipeline.
 

conf/igenomes.config

@@ -35,6 +35,7 @@ params {
             known_indels_tbi      = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi"
             known_indels_vqsr     = '--resource:1000G,known=false,training=true,truth=true,prior=10.0 1000G_phase1.indels.b37.vcf.gz --resource:mills,known=false,training=true,truth=true,prior=10.0 Mills_and_1000G_gold_standard.indels.b37.vcf.gz'
             mappability           = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/Control-FREEC/out100m2_hg19.gem"
+            ngscheckmate_bed      = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/NGSCheckMate/SNP_GRCh37_hg19_wChr.bed"
             snpeff_db             = 87
             snpeff_genome         = 'GRCh37'
             vep_cache_version     = 110
@@ -68,6 +69,7 @@ params {
             known_indels_tbi      = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/{Mills_and_1000G_gold_standard.indels.hg38,beta/Homo_sapiens_assembly38.known_indels}.vcf.gz.tbi"
             known_indels_vqsr     = '--resource:gatk,known=false,training=true,truth=true,prior=10.0 Homo_sapiens_assembly38.known_indels.vcf.gz --resource:mills,known=false,training=true,truth=true,prior=10.0 Mills_and_1000G_gold_standard.indels.hg38.vcf.gz'
             mappability           = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/Control-FREEC/out100m2_hg38.gem"
+            ngscheckmate_bed      = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/NGSCheckMate/SNP_GRCh38_hg38_wChr.bed"
             pon                   = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/1000g_pon.hg38.vcf.gz"
             pon_tbi               = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/1000g_pon.hg38.vcf.gz.tbi"
             snpeff_db             = 105
@@ -79,6 +81,7 @@ params {
         'Ensembl.GRCh37' {
             bwa                   = "${params.igenomes_base}/Homo_sapiens/Ensembl/GRCh37/Sequence/BWAIndex/version0.6.0/"
             fasta                 = "${params.igenomes_base}/Homo_sapiens/Ensembl/GRCh37/Sequence/WholeGenomeFasta/genome.fa"
+            ngscheckmate_bed      = "${params.igenomes_base}/Homo_sapiens/Ensembl/GRCh37/Annotation/NGSCheckMate/SNP_GRCh37_hg19_woChr.bed"
             readme                = "${params.igenomes_base}/Homo_sapiens/Ensembl/GRCh37/Annotation/README.txt"
             snpeff_db             = 87
             snpeff_genome         = 'GRCh37'
@@ -89,6 +92,7 @@ params {
         'NCBI.GRCh38' {
             bwa                   = "${params.igenomes_base}/Homo_sapiens/NCBI/GRCh38/Sequence/BWAIndex/version0.6.0/"
             fasta                 = "${params.igenomes_base}/Homo_sapiens/NCBI/GRCh38/Sequence/WholeGenomeFasta/genome.fa"
+            ngscheckmate_bed      ="${params.igenomes_base}/Homo_sapiens/NCBI/GRCh38/Annotation/NGSCheckMate/SNP_GRCh38_hg38_wChr.bed"
             snpeff_db             = 105
             snpeff_genome         = 'GRCh38'
             vep_cache_version     = 110

conf/modules/ngscheckmate.config

@@ -0,0 +1,26 @@
+process {
+    withName: ".*BAM_NGSCHECKMATE:BCFTOOLS_MPILEUP" {
+
+        publishDir = [
+            mode: params.publish_dir_mode,
+            path: { "${params.outdir}/reports/ngscheckmate/vcfs" },
+            pattern: "*{vcf.gz}"
+        ]
+        ext.prefix = { "${meta.id}.ngscheckmate" }
+        ext.when   = { params.tools && params.tools.split(',').contains('ngscheckmate') }
+        ext.args2 = '--no-version --ploidy 1 -c'
+        ext.args3 = '--no-version'
+    }
+
+    withName: ".*BAM_NGSCHECKMATE:NGSCHECKMATE_NCM" {
+        ext.args = '-V'
+
+        publishDir = [
+            mode: params.publish_dir_mode,
+            path: { "${params.outdir}/reports/ngscheckmate/" },
+            pattern: "*"
+        ]
+
+    }
+
+}

conf/test.config

@@ -40,6 +40,7 @@ params {
     bcftools_annotations_index  = "https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/sarscov2/illumina/vcf/test2.vcf.gz.tbi"
     bcftools_header_lines = "${projectDir}/tests/config/bcfann_test_header.txt"
 
+ngscheckmate_bed = "https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/homo_sapiens/genome/chr21/germlineresources/SNP_GRCh38_hg38_wChr.bed"
     // default params
     split_fastq       = 0         // no FASTQ splitting
     tools             = 'strelka' // Variant calling with Strelka

conf/test/cache.config

@@ -48,6 +48,7 @@ params {
     vep_cache_version = 110
     vep_genome        = 'WBcel235'
     vep_species       = 'caenorhabditis_elegans'
+    ngscheckmate_bed  = params.test_data['homo_sapiens']['genome']['ngscheckmate_bed']
 
     // default params
     split_fastq       = 0         // no FASTQ splitting

docs/output.md

@@ -61,6 +61,7 @@ The pipeline is built using [Nextflow](https://www.nextflow.io/) and processes d
     - [FastQC](#fastqc)
     - [FastP](#fastp)
     - [Mosdepth](#mosdepth)
+    - [NGSCheckMate](#ngscheckmate)
     - [GATK MarkDuplicates reports](#gatk-markduplicates-reports)
     - [Sentieon Dedup reports](#sentieon-dedup-reports)
     - [samtools stats](#samtools-stats)
@@ -985,6 +986,25 @@ Plots will show:
   - CSI index for per-base depth for targeted data, per-window (500bp) depth of WGS
   </details>
 
+#### NGSCheckMate
+
+[NGSCheckMate](https://github.com/parklab/NGSCheckMate) is a tool for determining whether samples come from the same genetic individual, using a set of commonly heterozygous SNPs. This enables for the detecting of sample mislabelling events. The output includes a text file indicating whether samples have matched or not according to the algorithm, as well as a dendrogram visualising these results.
+
+<details markdown="1">
+<summary>Output files for all samples</summary>
+
+**Output directory: `{outdir}/reports/ngscheckmate/`**
+
+- `ngscheckmate_all.txt`
+  - Tab delimited text file listing all the comparisons made, whether they were considered as a match, with the correlation and a normalised depth.
+- `ngscheckmate_matched.txt`
+  - Tab delimited text file listing only the comparison that were considered to match, with the correlation and a normalised depth.
+- `ngscheckmate_output_corr_matrix.txt`
+  - Tab delimited text file containing a matrix of all correlations for all comparisons made.
+- `vcfs/<sample>.vcf.gz`
+  - Set of vcf files for each sample. Contains calls for the set of SNP positions used to calculate sample relatedness.
+  </details>
+
 #### GATK MarkDuplicates reports
 
 More information in the [GATK MarkDuplicates section](#gatk-markduplicates)

modules.json

@@ -24,7 +24,7 @@
                     "bcftools/mpileup": {
                         "branch": "master",
                         "git_sha": "8fc1d24c710ebe1d5de0f2447ec9439fd3d9d66a",
-                        "installed_by": ["modules"]
+                        "installed_by": ["bam_ngscheckmate", "modules"]
                     },
                     "bcftools/sort": {
                         "branch": "master",
@@ -334,6 +334,11 @@
                         "git_sha": "8fc1d24c710ebe1d5de0f2447ec9439fd3d9d66a",
                         "installed_by": ["modules"]
                     },
+                    "ngscheckmate/ncm": {
+                        "branch": "master",
+                        "git_sha": "32d6725f584ebf460de39b7c1c53a29d5384d697",
+                        "installed_by": ["bam_ngscheckmate"]
+                    },
                     "samblaster": {
                         "branch": "master",
                         "git_sha": "8fc1d24c710ebe1d5de0f2447ec9439fd3d9d66a",
@@ -483,6 +488,11 @@
             },
             "subworkflows": {
                 "nf-core": {
+                    "bam_ngscheckmate": {
+                        "branch": "master",
+                        "git_sha": "32d6725f584ebf460de39b7c1c53a29d5384d697",
+                        "installed_by": ["subworkflows"]
+                    },
                     "vcf_annotate_ensemblvep": {
                         "branch": "master",
                         "git_sha": "cfd937a668919d948f6fcbf4218e79de50c2f36f",

nextflow.config

@@ -375,6 +375,9 @@ includeConfig 'conf/modules/recalibrate.config'
 includeConfig 'conf/modules/trimming.config'
 includeConfig 'conf/modules/umi.config'
 
+//ngscheckmate
+includeConfig 'conf/modules/ngscheckmate.config'
+
 // variant calling
 includeConfig 'conf/modules/ascat.config'
 includeConfig 'conf/modules/cnvkit.config'
@@ -393,7 +396,6 @@ includeConfig 'conf/modules/sentieon_haplotyper.config'
 includeConfig 'conf/modules/sentieon_haplotyper_joint_germline.config'
 includeConfig 'conf/modules/strelka.config'
 includeConfig 'conf/modules/tiddit.config'
-
 includeConfig 'conf/modules/post_variant_calling.config'
 
 //annotate

nextflow_schema.json

@@ -100,7 +100,8 @@
                     "fa_icon": "fas fa-toolbox",
                     "description": "Tools to use for duplicate marking, variant calling and/or for annotation.",
                     "help_text": "Multiple tools separated with commas.\n\n**Variant Calling:**\n\nGermline variant calling can currently be performed with the following variant callers:\n- SNPs/Indels: DeepVariant, FreeBayes, GATK HaplotypeCaller, mpileup, Sentieon Haplotyper, Strelka\n- Structural Variants: Manta, TIDDIT\n- Copy-number: CNVKit\n\nTumor-only somatic variant calling can currently be performed with the following variant callers:\n- SNPs/Indels: FreeBayes, mpileup, Mutect2, Strelka\n- Structural Variants: Manta, TIDDIT\n- Copy-number: CNVKit, ControlFREEC\n\nSomatic variant calling can currently only be performed with the following variant callers:\n- SNPs/Indels: FreeBayes, Mutect2, Strelka2\n- Structural variants: Manta, TIDDIT\n- Copy-Number: ASCAT, CNVKit, Control-FREEC\n- Microsatellite Instability: MSIsensorpro\n\n> **NB** Mutect2 for somatic variant calling cannot be combined with `--no_intervals`\n\n**Annotation:**\n \n- snpEff, VEP, merge (both consecutively).\n\n> **NB** As Sarek will use bgzip and tabix to compress and index VCF files annotated, it expects VCF files to be sorted when starting from `--step annotate`.",
-                    "pattern": "^((ascat|bcfann|cnvkit|controlfreec|deepvariant|freebayes|haplotypecaller|sentieon_dnascope|sentieon_haplotyper|manta|merge|mpileup|msisensorpro|mutect2|sentieon_dedup|snpeff|strelka|tiddit|vep)?,?)*(?<!,)$"
+
+                    "pattern": "^((ascat|bcfann|cnvkit|controlfreec|deepvariant|freebayes|haplotypecaller|sentieon_dnascope|sentieon_haplotyper|manta|merge|mpileup|msisensorpro|mutect2|ngscheckmate|sentieon_dedup|snpeff|strelka|tiddit|vep)?,?)*(?<!,)$"
                 },
                 "skip_tools": {
                     "type": "string",
@@ -752,6 +753,12 @@
                     "hidden": true,
                     "help_text": "If you use AWS iGenomes, this has already been set for you appropriately."
                 },
+                "ngscheckmate_bed": {
+                    "type": "string",
+                    "fa_icon": "fas fa-file",
+                    "description": "Path to SNP bed file for sample checking with NGSCheckMate",
+                    "help_text": "If you use AWS iGenomes, this has already been set for you appropriately."
+                },
                 "snpeff_db": {
                     "type": "string",
                     "fa_icon": "fas fa-database",

subworkflows/local/cram_sampleqc/main.nf

@@ -0,0 +1,30 @@
+include { BAM_NGSCHECKMATE } from '../../../subworkflows/nf-core/bam_ngscheckmate/main'
+
+workflow CRAM_SAMPLEQC {
+
+    take:
+    ch_cram             // channel: [ val(meta), cram, crai ]
+    ngscheckmate_bed    // channel: [ ngscheckmate_bed ]
+    fasta               // channel: [ fasta ]
+
+    main:
+
+    ch_versions = Channel.empty()
+
+    ch_ngscheckmate_bed = ngscheckmate_bed.map{bed -> [[id: "ngscheckmate"], bed]}
+
+    ch_fasta            = fasta.map{fasta -> [[id: "genome"], fasta]}
+
+    BAM_NGSCHECKMATE ( ch_cram.map{meta, cram, crai -> [meta, cram]}, ch_ngscheckmate_bed, ch_fasta)
+    ch_versions = ch_versions.mix(BAM_NGSCHECKMATE.out.versions.first())
+
+    emit:
+    corr_matrix  = BAM_NGSCHECKMATE.out.corr_matrix  // channel: [ meta, corr_matrix ]
+    matched      = BAM_NGSCHECKMATE.out.matched      // channel: [ meta, matched ]
+    all          = BAM_NGSCHECKMATE.out.all          // channel: [ meta, all ]
+    vcf          = BAM_NGSCHECKMATE.out.vcf          // channel: [ meta, vcf ]
+    pdf          = BAM_NGSCHECKMATE.out.pdf          // channel: [ meta, pdf ]
+
+    versions = ch_versions                     // channel: [ versions.yml ]
+}
+