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

CHANGELOG.md

@@ -13,6 +13,7 @@ Ráhpajávvre is the Lule Sámi spelling of Rapaselet.
 
 - [#1246](https://github.com/nf-core/sarek/pull/1246) - Back to dev
 - [#1259](https://github.com/nf-core/sarek/pull/1259) - nf-prov plugin
+- [#1252](https://github.com/nf-core/sarek/pull/1252) - Added NGSCheckMate tool for checking that samples come from the same individual
 
 ### Changed
 

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/tools.config

@@ -16,6 +16,7 @@ params {
     germline_resource      = params.test_data['homo_sapiens']['genome']['gnomad_r2_1_1_21_vcf_gz']
     intervals              = params.test_data['homo_sapiens']['genome']['genome_21_multi_interval_bed']
     pon                    = params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_21_vcf_gz']
+    ngscheckmate_bed       = "https://raw.githubusercontent.com/parklab/NGSCheckMate/master/SNP/SNP_GRCh38_hg38_wChr.bed"
     nucleotides_per_second = 20
     step                   = 'variant_calling'
     tools                  = null

conf/test/tools_germline.config

@@ -18,6 +18,7 @@ params {
     known_indels_vqsr      = "--resource:1000G,known=false,training=true,truth=true,prior=10.0 mills_and_1000G.indels.hg38.vcf.gz"
     known_snps             = params.test_data['homo_sapiens']['genome']['hapmap_3_3_hg38_21_vcf_gz']
     known_snps_vqsr        = "--resource:hapmap,known=false,training=true,truth=true,prior=10.0 hapmap_3.3.hg38.vcf.gz"
+    ngscheckmate_bed       = "https://raw.githubusercontent.com/parklab/NGSCheckMate/master/SNP/SNP_GRCh38_hg38_wChr.bed"
     nucleotides_per_second = 20
     step                   = 'variant_calling'
     tools                  = null

docs/output.md

@@ -1054,6 +1054,25 @@ The plots will show:
   - Summary of the VEP run to be visualised with a web browser
   </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>
+
 ### Reporting
 
 #### MultiQC

modules.json

@@ -18,7 +18,7 @@
                     "bcftools/mpileup": {
                         "branch": "master",
                         "git_sha": "911696ea0b62df80e900ef244d7867d177971f73",
-                        "installed_by": ["modules"]
+                        "installed_by": ["bam_ngscheckmate", "modules"]
                     },
                     "bcftools/sort": {
                         "branch": "master",
@@ -322,6 +322,11 @@
                         "git_sha": "a6e11ac655e744f7ebc724be669dd568ffdc0e80",
                         "installed_by": ["modules"]
                     },
+                    "ngscheckmate/ncm": {
+                        "branch": "master",
+                        "git_sha": "32d6725f584ebf460de39b7c1c53a29d5384d697",
+                        "installed_by": ["bam_ngscheckmate"]
+                    },
                     "samblaster": {
                         "branch": "master",
                         "git_sha": "603ecbd9f45300c9788f197d2a15a005685b4220",
@@ -461,6 +466,11 @@
             },
             "subworkflows": {
                 "nf-core": {
+                    "bam_ngscheckmate": {
+                        "branch": "master",
+                        "git_sha": "32d6725f584ebf460de39b7c1c53a29d5384d697",
+                        "installed_by": ["subworkflows"]
+                    },
                     "vcf_annotate_ensemblvep": {
                         "branch": "master",
                         "git_sha": "dedc0e31087f3306101c38835d051bf49789445a",

nextflow.config

@@ -93,6 +93,9 @@ params {
     vep_spliceai              = null // spliceai plugin disabled within VEP
     vep_spliceregion          = null // spliceregion plugin disabled within VEP
 
+    // NGSCheckMate
+    ngscheckmate_bed          = null
+
     // MultiQC options
     multiqc_config             = null
     multiqc_title              = null
@@ -387,6 +390,9 @@ includeConfig 'conf/modules/post_variant_calling.config'
 //annotate
 includeConfig 'conf/modules/annotate.config'
 
+//ngscheckmate
+includeConfig 'conf/modules/ngscheckmate.config'
+
 // Function to ensure that resource requirements don't go beyond
 // a maximum limit
 def check_max(obj, type) {

nextflow_schema.json

@@ -100,7 +100,7 @@
                     "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|cnvkit|controlfreec|deepvariant|freebayes|haplotypecaller|sentieon_haplotyper|manta|merge|mpileup|msisensorpro|mutect2|sentieon_dedup|snpeff|strelka|tiddit|vep)?,?)*(?<!,)$"
+                    "pattern": "^((ascat|cnvkit|controlfreec|deepvariant|freebayes|haplotypecaller|sentieon_haplotyper|manta|merge|mpileup|msisensorpro|mutect2|ngscheckmate|sentieon_dedup|snpeff|strelka|tiddit|vep)?,?)*(?<!,)$"
                 },
                 "skip_tools": {
                     "type": "string",
@@ -713,6 +713,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 ]
+}
+