checks the reference genome and create indices if necessary

Author: priesgo

.gitignore

@@ -6,4 +6,6 @@ report.html*
 timeline.html*
 trace.txt*
 dag.dot*
-*.swp
+*.swp
+/test_data/ucsc.hg19.minimal.without_indices.dict
+/test_data/ucsc.hg19.minimal.without_indices.fasta.fai

Makefile

@@ -17,3 +17,4 @@ test:
 	bash tests/test_08.sh
 	bash tests/test_09.sh
 	bash tests/test_10.sh
+	bash tests/test_11.sh

main.nf

@@ -7,6 +7,7 @@ include { MARK_DUPLICATES; SPLIT_CIGAR_N_READS } from './modules/02_mark_duplica
 include { METRICS; HS_METRICS; COVERAGE_ANALYSIS; FLAGSTAT } from './modules/03_metrics'
 include { REALIGNMENT_AROUND_INDELS } from './modules/04_realignment_around_indels'
 include { BQSR; CREATE_OUTPUT } from './modules/05_bqsr'
+include { CREATE_FAIDX; CREATE_DICT } from './modules/00_reference_indices'
 
 params.help= false
 params.input_files = false
@@ -82,9 +83,32 @@ else if (params.input_files) {
     .set { input_files }
 }
 
+workflow CHECK_REFERENCE {
+    take:
+        reference
+
+    main:
+        // checks the reference and its indexes, if the indexes are not there creates them
+        reference_file = file(reference)
+        if (reference_file.isEmpty()) {
+            log.error "--reference points to a non existing file"
+            exit 1
+        }
+        faidx = file("${reference}.fai")
+        if (faidx.isEmpty()) {
+            CREATE_FAIDX(reference)
+        }
+        dict =  file("${reference_file.getParent() }/${reference_file.baseName }*.dict")
+        if (dict.isEmpty()) {
+            CREATE_DICT(reference)
+        }
+}
+
 
 workflow {
 
+    CHECK_REFERENCE(params.reference)
+
     PREPARE_BAM(input_files, params.reference)
 
     if (!params.skip_deduplication) {

modules/00_reference_indices.nf

@@ -0,0 +1,30 @@
+
+process CREATE_FAIDX {
+    cpus "1"
+    memory "4g"
+    tag "${name}"
+
+    conda (params.enable_conda ? "bioconda::samtools=1.12" : null)
+
+    input:
+    val(reference)
+
+    """
+    samtools faidx ${reference}
+    """
+}
+
+process CREATE_DICT {
+    cpus "1"
+    memory "4g"
+    tag "${name}"
+
+    conda (params.enable_conda ? "bioconda::gatk4=4.2.5.0" : null)
+
+    input:
+    val(reference)
+
+    """
+    gatk CreateSequenceDictionary --REFERENCE ${reference}
+    """
+}

nextflow.config

@@ -44,7 +44,7 @@ process.shell = ['/bin/bash', '-euo', 'pipefail']
 
 cleanup = true
 
-VERSION = '2.0.1'
+VERSION = '2.1.0'
 DOI = 'https://zenodo.org/badge/latestdoi/358400957'
 
 manifest {

test_data/ucsc.hg19.minimal.without_indices.fasta

@@ -0,0 +1,9 @@
+#!/bin/bash
+
+
+source tests/assert.sh
+output=output/test11
+nextflow main.nf -profile test,conda --output $output --reference `pwd`/test_data/ucsc.hg19.minimal.without_indices.fasta
+
+test -s `pwd`/test_data/ucsc.hg19.minimal.without_indices.fasta.fai || { echo "Missing output FAI index!"; exit 1; }
+test -s `pwd`/test_data/ucsc.hg19.minimal.without_indices.dict || { echo "Missing output dict index!"; exit 1; }