DSL2: Qualimap

CITATIONS.md

@@ -85,8 +85,13 @@
   > Jun, G., Wing, M. K., Abecasis, G. R., & Kang, H. M. (2015). An efficient and scalable analysis framework for variant extraction and refinement from population-scale DNA sequence data. Genome Research, 25(6), 918–925. doi: [10.1101/gr.176552.114](https://doi.org/10.1101/gr.176552.114)
 
 - [DamageProfiler](https://doi.org/10.1093/bioinformatics/btab190)
+
   > DamageProfiler Neukamm, J., Peltzer, A., & Nieselt, K. (2020). DamageProfiler: Fast damage pattern calculation for ancient DNA. In Bioinformatics (btab190). doi: [10.1093/bioinformatics/btab190](https://doi.org/10.1093/bioinformatics/btab190). Download: https://github.com/Integrative-Transcriptomics/DamageProfiler
 
+- [QualiMap](https://doi.org/10.1093/bioinformatics/btv566)
+
+  > QualiMap Okonechnikov, K., Conesa, A., & García-Alcalde, F. (2016). Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics , 32(2), 292–294. Download: http://qualimap.bioinfo.cipf.es/
+
 ## Software packaging/containerisation tools
 
 - [Anaconda](https://anaconda.com)

conf/modules.config

@@ -726,4 +726,12 @@ process {
             enabled: true
         ]
     }
+
+    withName: "QUALIMAP_BAMQC" {
+        publishDir = [
+            path: { "${params.outdir}/qualimap/" },
+            mode: params.publish_dir_mode,
+            enabled: true
+        ]
+    }
 }

docs/output.md

@@ -380,18 +380,28 @@ Within nf-core/eager, when BAM trimming is activated alongside PMD-filtering, tr
 [DamageProfiler](https://github.com/Integrative-Transcriptomics/DamageProfiler)
 is a tool which calculates a variety of standard 'aDNA' metrics from a BAM file. The primary plots here are the misincorporation and length distribution plots. Ancient DNA undergoes depurination and hydrolysis, causing fragmentation of molecules into gradually shorter fragments, and cytosine to thymine deamination damage, that occur on the subsequent single-stranded overhangs at the ends of molecules.
 
-### Contamination estimation for Nuclear DNA
-
-#### ANGSD nuclear contamination
+#### QualiMap
 
 <details markdown="1">
 <summary>Output files</summary>
 
-- `contamination_estimation/angsd/`
+- `qualimap/`: this contains a sub-directory for every sample, which includes a qualimap report and associated raw statistic files. You can open the `.html` file in your internet browser to see the in-depth report (this will be more detailed than in MultiQC). This includes stuff like percent coverage, depth coverage, GC content and so on of your mapped reads.
+
+</details>
+
+[QualiMap](http://qualimap.bioinfo.cipf.es/)
+is a tool which provides statistics on the quality of the mapping of your reads to your reference genome. It allows you to assess how well covered your reference genome is by your data, both in 'fold' depth (average number of times a given base on the reference is covered by a read) and 'percentage' (the percentage of all bases on the reference genome that is covered at a given fold depth). These outputs allow you to make decision if you have enough quality data for downstream applications like genotyping, and how to adjust the parameters for those tools accordingly.
 
-  - `*.txt`: Text file containing the results of nuclear contamination estimation with ANGSD for each library.
-  - `nuclear_contamination.txt`: Text file containing a summary table of contamination estimates for all libraries.
-  - `nuclear_contamination_mqc.json`: JSON file containing a summary table of contamination estimates for all libraries.
+> NB: Neither fold coverage nor percent coverage on there own is sufficient to assess whether you have a high quality mapping. Abnormally high fold coverages of a smaller region such as highly conserved genes or un-removed-adapter-containing reference genomes can artificially inflate the mean coverage, yet a high percent coverage is not useful if all bases of the genome are covered at just 1x coverage.
+
+**Note that many of the statistics from this module are displayed in the General Stats table, as they represent single values that are not plottable.**
+
+You will receive output for each sample. This means you will statistics of deduplicated values of all types of libraries combined in a single value (i.e. non-UDG treated, full-UDG, paired-end, single-end all together).
+warning: If your library has no reads mapping to the reference, this will result in an empty BAM file. Qualimap will therefore not produce any output even if a BAM exists!
+
+### Contamination estimation for Nuclear DNA
+
+#### ANGSD nuclear contamination
 
 </details>
 

modules.json

@@ -150,6 +150,11 @@
                         "git_sha": "911696ea0b62df80e900ef244d7867d177971f73",
                         "installed_by": ["modules"]
                     },
+                    "qualimap/bamqc": {
+                        "branch": "master",
+                        "git_sha": "911696ea0b62df80e900ef244d7867d177971f73",
+                        "installed_by": ["modules"]
+                    },
                     "samtools/faidx": {
                         "branch": "master",
                         "git_sha": "911696ea0b62df80e900ef244d7867d177971f73",

nextflow.config

@@ -164,6 +164,10 @@ params {
     damageprofiler_threshold     = 15
     damageprofiler_yaxis         = 0.3
 
+    // Qualimap options
+    skip_qualimap                = false
+    snpcapture_bed               = null
+
     // Damage Manipulation options
     run_mapdamage_rescaling                                          = false
     damage_manipulation_rescale_seqlength                            = 12

nextflow_schema.json

@@ -966,19 +966,20 @@
                     "description": "Specifies the length filter for DamageProfiler.",
                     "help_text": "> Modifies DamageProfile parameter: -l "
                 },
-                "damageprofiler_threshold": {
-                    "type": "integer",
-                    "default": 15,
-                    "description": "Specify number of bases of each read to consider for DamageProfiler calculations.",
-                    "fa_icon": "fas fa-ruler-horizontal",
-                    "help_text": "Specifies the length of the read start and end to be considered for profile generation in DamageProfiler.\n\n> Modifies DamageProfile parameter: -t\n"
-                },
+
                 "damageprofiler_yaxis": {
                     "type": "number",
                     "default": 0.3,
                     "fa_icon": "fas fa-ruler-vertical",
                     "description": "Specify the maximum misincorporation frequency that should be displayed on damage plot. Set to 0 to 'autoscale'.",
-                    "help_text": "Specifies what the maximum misincorporation frequency should be displayed as, in the DamageProfiler damage plot. This is set to 0.30 (i.e. 30%) by default as this matches the popular [mapDamage2.0 program](https://ginolhac.github.io/mapDamage/). However, the default behavior of the DamageProfiler is to 'autoscale' the y-axis maximum to zoom in on any possible damage that may occur (e.g. if the damage is about 10%, the highest value on the y-axis would be set to 0.12). This 'autoscale' behavior can be turned on by specifying the number to 0.\n\n> Modifies DamageProfile parameter: -yaxis_dp_max"
+                    "help_text": "Specifies what the maximum misincorporation frequency should be displayed as, in the DamageProfiler damage plot. This is set to 0.30 (i.e. 30%) by default as this matches the popular [mapDamage2.0 program](https://ginolhac.github.io/mapDamage/). However, the default behavior of the DamageProfiler is to 'autoscale' the y-axis maximum to zoom in on any possible damage that may occur (e.g. if the damage is about 10%, the highest value on the y-axis would be set to 0.12). This 'autoscale' behavior can be turned on by specifying the number to 0. Default: 0.30.\n\n> Modifies DamageProfile parameter: -yaxis_dp_max"
+                },
+                "damageprofiler_threshold": {
+                    "type": "integer",
+                    "default": 15,
+                    "description": "Specify number of bases of each read to consider for DamageProfiler calculations.",
+                    "fa_icon": "fas fa-ruler-horizontal",
+                    "help_text": "Specifies the length of the read start and end to be considered for profile generation in DamageProfiler. By default set to 15 bases.\n\n> Modifies DamageProfile parameter: -t"
                 }
             }
         },
@@ -1061,6 +1062,26 @@
                     "fa_icon": "fas fa-map"
                 }
             }
+        },
+        "statistics": {
+            "title": "Statistics",
+            "type": "object",
+            "description": "Additional options.",
+            "default": "",
+            "properties": {
+                "skip_qualimap": {
+                    "type": "boolean",
+                    "fa_icon": "fas fa-forward",
+                    "description": "Turns off QualiMap and thus does not compute coverage and other mapping metrics."
+                },
+                "snpcapture_bed": {
+                    "type": "string",
+                    "default": null,
+                    "help_text": "Can be used to set a path to a BED file (3/6 column format) of SNP positions of a reference genome, to calculate SNP captured libraries on-target efficiency. This should be used for array or in-solution SNP capture protocols such as 390K, 1240K, etc. If supplied, some on-target metrics are automatically generated for you by qualimap in the 'Globals inside' section of the 'genome_results.txt' file in the qualimap results directory. These statistics are currently NOT displayed in MultiQC!",
+                    "fa_icon": "fas fa-magnet"
+                }
+            },
+            "fa_icon": "fas fa-calculator"
         }
     },
     "allOf": [
@@ -1112,6 +1133,9 @@
         {
             "$ref": "#/definitions/adna_damage_analysis"
         },
+        {
+            "$ref": "#/definitions/statistics"
+        },
         {
             "$ref": "#/definitions/contamination_estimation"
         },

workflows/eager.nf

@@ -91,6 +91,7 @@ include { PRESEQ_CCURVE               } from '../modules/nf-core/preseq/ccurve/m
 include { PRESEQ_LCEXTRAP             } from '../modules/nf-core/preseq/lcextrap/main'
 include { FALCO                       } from '../modules/nf-core/falco/main'
 include { MTNUCRATIO                  } from '../modules/nf-core/mtnucratio/main'
+include { QUALIMAP_BAMQC              } from '../modules/nf-core/qualimap/bamqc/main' 
 include { HOST_REMOVAL                } from '../modules/local/host_removal'
 
 /*
@@ -128,6 +129,13 @@ workflow EAGER {
         if ( params.preprocessing_tool == 'fastp' && !adapterlist.extension.matches(".*(fa|fasta|fna|fas)") ) error "[nf-core/eager] ERROR: fastp adapter list requires a `.fasta` format and extension (or fa, fas, fna). Check input: --preprocessing_adapterlist ${params.preprocessing_adapterlist}"
     }
 
+    // QualiMap
+   if ( params.snpcapture_bed ) {
+    ch_snpcapture_bed = Channel.fromPath(params.snpcapture_bed, checkIfExists: true).collect()
+    } else {
+    ch_snpcapture_bed = Channel.empty()
+    }
+
     // Contamination estimation
     hapmap_file = file(params.contamination_estimation_angsd_hapmap, checkIfExists:true)
 
@@ -352,6 +360,23 @@ workflow EAGER {
     }
 
     //
+    // MODULE: QUALIMAP_BAMQC
+    //
+
+    if ( !params.skip_qualimap ) {
+
+        qualimap_input = ch_dedupped_bams
+        .map {
+            meta, bam, bai ->
+            [ meta, bam ]
+        }
+
+        QUALIMAP_BAMQC( qualimap_input, ch_snpcapture_bed )
+        ch_versions = ch_versions.mix( QUALIMAP_BAMQC.out.versions )
+        ch_multiqc_files = ch_multiqc_files.mix(QUALIMAP_BAMQC.out.results.collect{it[1]}.ifEmpty([]))
+
+    }
+
     // SUBWORKFLOW: Contamination estimation
     //