Merge pull request #93 from NBISweden/mask_swegen

Mask swegen feature

Author: viklund

README.md

@@ -65,7 +65,10 @@ Options:
     --steps         Specify what steps to run, comma separated: (default: manta, vep)
                 Callers: manta, fermikit
                 Annotation: vep, snpeff
-                Extra: normalize (with vt)
+                Extra: normalize (with vt),
+                       filter (with bed files in masks_filters/, by default swegen is used)
+    --sg_mask_ovlp  Fractional overlap for use with the filter option
+    --no_sg_reciprocal  Don't use a reciprocal overlap for the filter option
     --outdir        Directory where resultfiles are stored (default: results)
     --prefix        Prefix for result filenames (default: no prefix)
 ```
@@ -101,14 +104,23 @@ usage message for information on them.
 The next section has the reference assembly to use, both as fasta and assembly
 name.
 
-You may want to use different versions of the modules used by this workflow, and
-the `modules` section contains all of them and their
-versions. Customize the modules here, NOT in the `main.nf` file. 
-
-Finally the `runtime` section has the different runtimes for the different
-parts of the workflow. `fermikit` has it's own timespec since that is a very
-long running program, otherwise the workflow differentiates between `callers`
-and other supporting `simple` single-core jobs.
+You may want to use different versions of the modules used by this workflow,
+currently you will have to edit the profiles to do that. On uppmax we have the
+milou profile which specifies all the modules and versions, see the
+`config/milou.config`.
+
+The runtimes of the different programs is set in the `config/standard.config`
+file. That file also specifies how to deal with errors and the interaction
+with the Slurm scheduler, you probably don't want to change those unless you
+know what you are doing.
+
+The two folders `masks_artifacts` and `masks_filters` contain bed files to
+filter the vcf-files from the callers. The artifact directory contains files
+that should mask out problematic regions, it removes everything that has an
+overlaps at least 25% with a region in the artifact mask. The filter one is
+for more stringent filtering of already known variants, and here the default
+filter threshold is instead a reciprocal overlap of 95%. It can be customized
+with the two options `sg_mask_ovlp` (default 0.95) and `no_sg_reciprocal`.
 
 
 ## Support

config/devcore.config

@@ -1,14 +1,13 @@
 process {
+    clusterOptions = {
+        "-A $params.project"
+    }
 
-   clusterOptions = {
-       "-A $params.project"
-   }
-   errorStrategy = { task.exitStatus == 143 ? 'retry' : 'terminate' }
-   queue = 'devcore'
-   executor = 'slurm'
-
-   time = '1h'
+    errorStrategy = { task.exitStatus == 143 ? 'retry' : 'terminate' }
+    queue = 'devcore'
+    executor = 'slurm'
 
+    time = '1h'
 }
 
 executor {

config/milou.config

@@ -1,43 +1,46 @@
 process {
+    $index_bamfile {
+        module = ['bioinfo-tools', 'samtools/1.3']
+    }
 
-  $index_bamfile {
-    module = ['bioinfo-tools', 'samtools/1.3']
-  }
+    $manta {
+        maxRetries = 3
+        cpus = 16
+        module = ['bioinfo-tools', 'manta/1.0.3']
+    }
 
-  $manta {
-    maxRetries = 3
-    cpus = 16
-    module = ['bioinfo-tools', 'manta/1.0.3']
-  }
+    $create_fastq {
+        module = ['bioinfo-tools', 'samtools/1.3']
+    }
 
-  $create_fastq {
-    module = ['bioinfo-tools', 'samtools/1.3']
+    $fermikit {
+        maxRetries = 3
+        cpus = 16
+        module = ['bioinfo-tools', 'samtools/1.3', 'fermikit/r178', "vcftools/0.1.14", "tabix/0.2.6"]
     }
 
-  $fermikit {
-    maxRetries = 3
-    cpus = 16
-    module = ['bioinfo-tools', 'samtools/1.3', 'fermikit/r178', "vcftools/0.1.14", "tabix/0.2.6"]
+    $artifact_mask_vcfs {
+        module = ['bioinfo-tools', 'BEDTools/2.25.0']
     }
 
-  $mask_vcfs {
-    module = ['bioinfo-tools', 'BEDTools/2.25.0']
+    $swegen_mask_vcfs {
+        module = ['bioinfo-tools', 'BEDTools/2.25.0']
     }
 
-  $intersect_files {
-    module = ['bioinfo-tools', 'BEDTools/2.25.0']
+    $intersect_files {
+        module = ['bioinfo-tools', 'BEDTools/2.25.0']
     }
 
-  $normalize_vcf {
-    module = ['bioinfo-tools', 'vt/0.5772']
+    $normalize_vcf {
+        module = ['bioinfo-tools', 'vt/0.5772']
     }
 
-  $variant_effect_predictor {
-    cpus = 4
-    module = ['bioinfo-tools', 'vep/84']
+    $variant_effect_predictor {
+        cpus = 4
+        module = ['bioinfo-tools', 'vep/84']
     }
 
-  $snpEff {
-    module = ['bioinfo-tools', 'snpEff/4.2', 'vt/0.5772']
+    $snpEff {
+        module = ['bioinfo-tools', 'snpEff/4.2', 'vt/0.5772']
     }
 }

config/standard.config

@@ -1,5 +1,4 @@
 process {
-
     clusterOptions = {
         "-A $params.project"
     }
@@ -20,5 +19,5 @@ process {
 
     $fermikit {
         time = params.runtime.fermikit * 2**( task.attempt -1 )
-
+    }
 }

main.nf

@@ -156,36 +156,75 @@ process fermikit {
 
 // 3. Create summary files
 
-mask_dir = file("$baseDir/masks")
 ch_vcfs = ch_manta_vcf.mix( ch_fermi_vcf )
 
-process mask_vcfs {
+process artifact_mask_vcfs {
     input:
-        each mask_dir from mask_dir
         set file(svfile), val(uuid), val(dir) from ch_vcfs
     output:
-        set file('*_masked.vcf'), val(uuid), val(dir) into ch_masked_vcfs
+        set file(svfile), val(uuid), val(dir) into ch_artifact_masked_vcfs
 
     tag "$uuid $svfile"
 
     executor choose_executor()
 
     """
     BNAME=\$( echo $svfile | cut -d. -f1 )
-    MASK_FILE=\${BNAME}_masked.vcf
-    MASK_DIR=$mask_dir
+    MASK_DIR=$params.mask_dirs.masks_artifacts
 
+    # We don't want to change the filename in this process so we copy the
+    # infile and remove the symbolic link. And then recreate the file at the
+    # end.
     cp $svfile workfile
+    rm $svfile # It's a link, should be ok :)
     for mask in \$MASK_DIR/*; do
+        if [ ! -f "\$mask" ]; then
+            break
+        fi
         cat workfile \
             | bedtools intersect -header -v -a stdin -b \$mask -f 0.25 \
             > tempfile
         mv tempfile workfile
     done
+    mv workfile $svfile
+    """
+}
+
+
+ch_noswegen_mask = ch_artifact_masked_vcfs.tap { ch_swegen_mask_in }
+reciprocal = params.no_sg_reciprocal ? '': '-r'
+
+process swegen_mask_vcfs {
+    input:
+        set file(svfile), val(uuid), val(dir) from ch_swegen_mask_in
+    output:
+        set file('*_swegen_masked.vcf'), val(uuid), val(dir) into ch_swegen_masked_vcfs
+
+    tag "$uuid $svfile"
+
+    executor choose_executor()
+    when 'filter' in workflowSteps
+
+    """
+    BNAME=\$( echo $svfile | cut -d. -f1 )
+    MASK_FILE=\${BNAME}_swegen_masked.vcf
+    MASK_DIR=$params.mask_dirs.masks_filters
+
+    cp $svfile workfile
+    for mask in \$MASK_DIR/*; do
+        if [ ! -f "\$mask" ]; then
+            break
+        fi
+        cat workfile \
+            | bedtools intersect -header -v -a stdin -b \$mask -f $params.sg_mask_ovlp $reciprocal \
+            > tempfile
+        mv tempfile workfile
+    done
     mv workfile \$MASK_FILE
     """
 }
 
+ch_masked_vcfs = ch_noswegen_mask.mix(ch_swegen_masked_vcfs)
 
 // To make intersect files we need to combine them into one channel with
 // toList() and then sort in the map so that fermi is before manta in the
@@ -199,7 +238,7 @@ process intersect_files {
     input:
         set file(vcfs), val(uuid), val(dir) from ch_intersect_input
     output:
-        set file('combined_masked.vcf'), val(uuid), val("${dir[0]}") into ch_intersections
+        set file('combined_*.vcf'), val(uuid), val("${dir[0]}") into ch_intersections
 
     tag "$uuid"
 
@@ -217,21 +256,22 @@ process intersect_files {
         manta_vcf=${vcfs[0]}
     fi
 
+    OUTNAME=`basename \$fermi_vcf|sed 's/fermikit_//;s/.vcf//'`
     ## Create intersected vcf files
     for WORD in DEL INS DUP; do
         intersectBed -a <( grep -w "^#.\\+\\|\$WORD" \$fermi_vcf) \
                      -b <( grep -w "^#.\\+\\|\$WORD" \$manta_vcf) \
             -f 0.5 -r \
-            | sort -k1,1V -k2,2n > combined_masked_\${WORD,,}.vcf
+            | sort -k1,1V -k2,2n > cmb_\${OUTNAME}_\${WORD,,}.vcf
     done
 
     cat <( grep -v -w '^#.\\+\\|DEL\\|INS\\|DUP' \$fermi_vcf ) \
         <( grep -v -w '^#.\\+\\|DEL\\|INS\\|DUP' \$manta_vcf ) \
         | cut -f 1-8 \
-        | sort -k1,1V -k2,2n > combined_masked_OTHER.vcf
+        | sort -k1,1V -k2,2n > cmb_\${OUTNAME}_OTHER.vcf
 
     ( grep '^#' \$fermi_vcf; \
-        sort -k1,1V -k2,2n combined_masked_*.vcf ) >> combined_masked.vcf
+        sort -k1,1V -k2,2n cmb_\${OUTNAME}_*.vcf ) >> combined_\${OUTNAME}.vcf
     """
 }
 
@@ -333,7 +373,7 @@ process variant_effect_predictor {
         --total_length             \
         --canonical                \
         --ccds                     \
-        --fork "\$SLURM_JOB_CPUS_PER_NODE" \
+        \$( test "\$SLURM_CPUS_ON_NODE" -gt 1 && echo "--fork \$SLURM_CPUS_ON_NODE" ) \
         --fields Consequence,Codons,Amino_acids,Gene,SYMBOL,Feature,EXON,PolyPhen,SIFT,Protein_position,BIOTYPE \
         --assembly "\$ASSEMBLY" \
         --offline
@@ -398,7 +438,10 @@ def usage_message() {
     log.info '    --steps         Specify what steps to run, comma separated: (default: manta, vep)'
     log.info '                Callers: manta, fermikit'
     log.info '                Annotation: vep, snpeff'
-    log.info '                Extra: normalize (with vt)'
+    log.info '                Extra: normalize (with vt),'
+    log.info '                       filter (with bed files in masks_filters/, by default swegen is used)'
+    log.info '    --sg_mask_ovlp  Fractional overlap for use with the filter option'
+    log.info '    --no_sg_reciprocal  Don't use a reciprocal overlap for the filter option'
     log.info '    --outdir        Directory where resultfiles are stored (default: results)'
     log.info '    --prefix        Prefix for result filenames (default: no prefix)'
     log.info ''

masks/LCR-hs37d5.bed.gz masks_artifacts/LCR-hs37d5.bed.gzmasks/LCR-hs37d5.bed.gz renamed to masks_artifacts/LCR-hs37d5.bed.gz

@@ -8,6 +8,10 @@ params {
     prefix = ''
     help = False
 
+    //optional swegen masking
+    sg_mask_ovlp = 0.95 // --sg_mask_ovlp in case other overlap value than default 0.95 for swegen
+    no_sg_reciprocal = False // --no_sg_reciprocal in case swegen masking should not include reciprocal overlap option
+
     // Reference assemblies
     ref_fasta = "/sw/data/uppnex/ToolBox/ReferenceAssemblies/hg38make/bundle/2.8/b37/human_g1k_v37.fasta"
     assembly = 'GRCh37'
@@ -18,6 +22,12 @@ params {
         fermikit = '48h' // Fermikit is the longest running of them all
         caller = '24h'  // The rest are a lot quicker
     }
+
+    // Dirs that contain BED files for masking
+    mask_dirs {
+        masks_artifacts = "$baseDir/masks_artifacts"
+        masks_filters = "$baseDir/masks_filters"
+    }
 }
 
 profiles {