Merge pull request #220 from ENCODE-DCC/dev

v1.8.0

Author: leepc12

chip.wdl

@@ -1,16 +1,16 @@
 version 1.0
 
 workflow chip {
-    String pipeline_ver = 'v1.7.1'
+    String pipeline_ver = 'v1.8.0'
 
     meta {
-        version: 'v1.7.1'
+        version: 'v1.8.0'
         author: 'Jin wook Lee (leepc12@gmail.com) at ENCODE-DCC'
         description: 'ENCODE TF/Histone ChIP-Seq pipeline'
         specification_document: 'https://docs.google.com/document/d/1lG_Rd7fnYgRpSIqrIfuVlAz2dW1VaSQThzk836Db99c/edit?usp=sharing'
 
-        caper_docker: 'encodedcc/chip-seq-pipeline:v1.7.1'
-        caper_singularity: 'docker://encodedcc/chip-seq-pipeline:v1.7.1'
+        caper_docker: 'encodedcc/chip-seq-pipeline:v1.8.0'
+        caper_singularity: 'docker://encodedcc/chip-seq-pipeline:v1.8.0'
         croo_out_def: 'https://storage.googleapis.com/encode-pipeline-output-definition/chip.croo.v5.json'
 
         parameter_group: {
@@ -146,6 +146,8 @@ workflow chip {
         String aligner = 'bowtie2'
         File? custom_align_py
         Boolean use_bwa_mem_for_pe = false
+        Int bwa_mem_read_len_limit = 70
+        Boolean use_bowtie2_local_mode = false
         Int crop_length = 0
         Int crop_length_tol = 2
         String trimmomatic_phred_score_format = 'auto'
@@ -176,7 +178,7 @@ workflow chip {
         # group: resource_parameter
         Int align_cpu = 6
         Float align_bowtie2_mem_factor = 0.15
-        Float align_bwa_mem_factor = 0.3
+        Float align_bwa_mem_factor = 1.0
         Int align_time_hr = 48
         Float align_bowtie2_disk_factor = 8.0
         Float align_bwa_disk_factor = 8.0
@@ -639,9 +641,19 @@ workflow chip {
             help: 'There is a template included in the documentation for inputs. Defining this parameter will automatically change "chip.aligner" to "custom". You should also define "chip.custom_aligner_idx_tar".'
         }
         use_bwa_mem_for_pe: {
-            description: 'For paired end dataset with read length >= 70bp, use bwa mem instead of bwa aln.',
+            description: 'For paired end dataset with read length >= chip.bwa_mem_read_len_limit (default 70) bp, use bwa mem instead of bwa aln.',
             group: 'alignment',
-            help: 'Use it only for paired end reads >= 70bp.'
+            help: 'Use it only for paired end reads >= chip.bwa_mem_read_len_limit (default 70) bp. Otherwise keep using bwa aln.'
+        }
+        bwa_mem_read_len_limit: {
+            description: 'Read length limit for bwa mem (for PE FASTQs only).',
+            group: 'alignment',
+            help: 'If chip.use_bwa_mem_for_pe is activated and reads are shorter than this limit, then bwa aln will be used instead of bwa mem.'
+        }
+        use_bowtie2_local_mode: {
+            description: 'Use bowtie2\'s local mode (soft-clipping).',
+            group: 'alignment',
+            help: 'This will add --local to bowtie2 command line so that it will replace the default end-to-end mode.'
         }
         crop_length: {
             description: 'Crop FASTQs\' reads longer than this length.',
@@ -1133,6 +1145,11 @@ workflow chip {
             msg = 'To use chip.use_bwa_mem_for_pe, choose bwa for chip.aligner.'
         }
     }
+    if ( aligner_ != 'bowtie2' && use_bowtie2_local_mode ) {
+        call raise_exception as error_use_bowtie2_local_mode_for_non_bowtie2 { input:
+            msg = 'To use chip.use_bowtie2_local_mode, choose bowtie2 for chip.aligner.'
+        }
+    }
     if ( aligner_ == 'custom' && ( !defined(custom_align_py) || !defined(custom_aligner_idx_tar) ) ) {
         call raise_exception as error_custom_aligner { input:
             msg = 'To use a custom aligner, define chip.custom_align_py and chip.custom_aligner_idx_tar.'
@@ -1185,6 +1202,8 @@ workflow chip {
                     else custom_aligner_idx_tar,
                 paired_end = paired_end_,
                 use_bwa_mem_for_pe = use_bwa_mem_for_pe,
+                bwa_mem_read_len_limit = bwa_mem_read_len_limit,
+                use_bowtie2_local_mode = use_bowtie2_local_mode,
                 ref_fa = ref_fa_,
 
                 trimmomatic_java_heap = align_trimmomatic_java_heap,
@@ -1282,7 +1301,9 @@ workflow chip {
                     else if aligner=='bowtie2' then bowtie2_idx_tar_
                     else custom_aligner_idx_tar,
                 paired_end = false,
-                use_bwa_mem_for_pe = use_bwa_mem_for_pe,
+                use_bwa_mem_for_pe = false,
+                bwa_mem_read_len_limit = 0,
+                use_bowtie2_local_mode = use_bowtie2_local_mode,
                 ref_fa = ref_fa_,
 
                 cpu = align_cpu,
@@ -1413,6 +1434,8 @@ workflow chip {
                     else custom_aligner_idx_tar,
                 paired_end = ctl_paired_end_,
                 use_bwa_mem_for_pe = use_bwa_mem_for_pe,
+                bwa_mem_read_len_limit = bwa_mem_read_len_limit,
+                use_bowtie2_local_mode = use_bowtie2_local_mode,
                 ref_fa = ref_fa_,
 
                 trimmomatic_java_heap = align_trimmomatic_java_heap,
@@ -2028,6 +2051,8 @@ task align {
         File? idx_tar            # reference index tar
         Boolean paired_end
         Boolean use_bwa_mem_for_pe
+        Int bwa_mem_read_len_limit
+        Boolean use_bowtie2_local_mode
 
         String? trimmomatic_java_heap
         Int cpu
@@ -2036,7 +2061,7 @@ task align {
         Float disk_factor
     }
     Float input_file_size_gb = size(fastqs_R1, "G") + size(fastqs_R2, "G")
-    Float mem_gb = 5.0 + mem_factor * input_file_size_gb
+    Float mem_gb = 5.0 + size(idx_tar, "G") + mem_factor * input_file_size_gb
     Float samtools_mem_gb = 0.8 * mem_gb
     Int disk_gb = round(40.0 + disk_factor * input_file_size_gb)
 
@@ -2102,6 +2127,7 @@ task align {
                 ${if paired_end then 'R2$SUFFIX/*.fastq.gz' else ''} \
                 ${if paired_end then '--paired-end' else ''} \
                 ${if use_bwa_mem_for_pe then '--use-bwa-mem-for-pe' else ''} \
+                ${'--bwa-mem-read-len-limit ' + bwa_mem_read_len_limit} \
                 ${'--mem-gb ' + samtools_mem_gb} \
                 ${'--nth ' + cpu}
 
@@ -2112,6 +2138,7 @@ task align {
                 ${if paired_end then 'R2$SUFFIX/*.fastq.gz' else ''} \
                 ${'--multimapping ' + multimapping} \
                 ${if paired_end then '--paired-end' else ''} \
+                ${if use_bowtie2_local_mode then '--local' else ''} \
                 ${'--mem-gb ' + samtools_mem_gb} \
                 ${'--nth ' + cpu}
         else

dev/test/test_task/test_bowtie2.json

@@ -12,5 +12,9 @@
         "chip-seq-pipeline-test-data/input/se/fastqs/rep1/rep1.subsampled.25.fastq.gz"
     ],
     "test_bowtie2.ref_pe_flagstat" : "chip-seq-pipeline-test-data/ref_output/test_bowtie2/pe/rep1-R1.subsampled.67.samstats.qc",
-    "test_bowtie2.ref_se_flagstat" : "chip-seq-pipeline-test-data/ref_output/test_bowtie2/se/rep1.subsampled.25.samstats.qc"
+    "test_bowtie2.ref_se_flagstat" : "chip-seq-pipeline-test-data/ref_output/test_bowtie2/se/rep1.subsampled.25.samstats.qc",
+
+    "test_bowtie2.ref_pe_local_flagstat" : "chip-seq-pipeline-test-data/ref_output/test_bowtie2/pe/local/rep1-R1.subsampled.67.srt.samstats.qc",
+    "test_bowtie2.ref_se_local_flagstat" : "chip-seq-pipeline-test-data/ref_output/test_bowtie2/se/local/rep1.subsampled.25.srt.samstats.qc"
+   
 }

dev/test/test_task/test_bowtie2.wdl

@@ -11,7 +11,9 @@ workflow test_bowtie2 {
         # we don't compare BAM because BAM's header includes date
         # hence md5sums don't match all the time
         String ref_pe_flagstat
-        String ref_se_flagstat
+        String ref_se_flagstat    
+        String ref_pe_local_flagstat
+        String ref_se_local_flagstat
 
         String pe_bowtie2_idx_tar
         String se_bowtie2_idx_tar
@@ -30,6 +32,8 @@ workflow test_bowtie2 {
         fastqs_R2 = pe_fastqs_R2,
         paired_end = true,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = 0,
         crop_length_tol = 0,
 
@@ -46,6 +50,45 @@ workflow test_bowtie2 {
         fastqs_R2 = [],
         paired_end = false,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
+        crop_length = 0,
+        crop_length_tol = 0,
+
+        cpu = bowtie2_cpu,
+        mem_factor = bowtie2_mem_factor,
+        time_hr = bowtie2_time_hr,
+        disk_factor = bowtie2_disk_factor,
+    }
+
+    call chip.align as pe_bowtie2_local { input :
+        aligner = 'bowtie2',
+        idx_tar = pe_bowtie2_idx_tar,
+        mito_chr_name = 'chrM',
+        fastqs_R1 = pe_fastqs_R1,
+        fastqs_R2 = pe_fastqs_R2,
+        paired_end = true,
+        use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = true,
+        crop_length = 0,
+        crop_length_tol = 0,
+
+        cpu = bowtie2_cpu,
+        mem_factor = bowtie2_mem_factor,
+        time_hr = bowtie2_time_hr,
+        disk_factor = bowtie2_disk_factor,
+    }
+    call chip.align as se_bowtie2_local { input :
+        aligner = 'bowtie2',
+        idx_tar = se_bowtie2_idx_tar,
+        mito_chr_name = 'chrM',
+        fastqs_R1 = se_fastqs,
+        fastqs_R2 = [],
+        paired_end = false,
+        use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = true,
         crop_length = 0,
         crop_length_tol = 0,
 
@@ -59,14 +102,20 @@ workflow test_bowtie2 {
         labels = [
             'pe_bowtie2',
             'se_bowtie2',
+            'pe_local_bowtie2',
+            'se_local_bowtie2',
         ],
         files = [
             pe_bowtie2.samstat_qc,
             se_bowtie2.samstat_qc,
+            pe_bowtie2_local.samstat_qc,
+            se_bowtie2_local.samstat_qc,
         ],
         ref_files = [
             ref_pe_flagstat,
             ref_se_flagstat,
+            ref_pe_local_flagstat,
+            ref_se_local_flagstat,
         ],
     }
 }

dev/test/test_task/test_bwa.wdl

@@ -30,6 +30,8 @@ workflow test_bwa {
         fastqs_R2 = pe_fastqs_R2,
         paired_end = true,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = 0,
         crop_length_tol = 0,
 
@@ -46,6 +48,8 @@ workflow test_bwa {
         fastqs_R2 = [],
         paired_end = false,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = 0,
         crop_length_tol = 0,
 

dev/test/test_task/test_trimmomatic.wdl

@@ -36,6 +36,8 @@ workflow test_trimmomatic {
         fastqs_R2 = pe_fastqs_R2,
         paired_end = true,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = pe_crop_length,
         crop_length_tol = pe_crop_length_tol,
 
@@ -52,6 +54,8 @@ workflow test_trimmomatic {
         fastqs_R2 = [],
         paired_end = false,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = se_crop_length,
         crop_length_tol = se_crop_length_tol,
 
@@ -69,6 +73,8 @@ workflow test_trimmomatic {
         fastqs_R2 = pe_fastqs_R2,
         paired_end = true,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = pe_crop_length,
         crop_length_tol = pe_crop_length_tol,
         trimmomatic_phred_score_format = 'phred33',
@@ -86,6 +92,8 @@ workflow test_trimmomatic {
         fastqs_R2 = [],
         paired_end = false,
         use_bwa_mem_for_pe = false,
+        bwa_mem_read_len_limit = 70,
+        use_bowtie2_local_mode = false,
         crop_length = se_crop_length,
         crop_length_tol = se_crop_length_tol,
         trimmomatic_phred_score_format = 'phred64',

docs/input.md

@@ -153,7 +153,9 @@ Parameter|Type|Default|Description
 `chip.crop_length` | Int | 0 | Crop FASTQs with Trimmomatic (using parameters `CROP`). 0: cropping disabled.
 `chip.crop_length_tol` | Int | 2 | Trimmomatic's `MINLEN` will be set as `chip.crop_length` - `abs(chip.crop_length_tol)` where reads shorter than `MINLEN` will be removed, hence not included in output BAM files and all downstream analyses.
 `chip.trimmomatic_phred_score_format` | String | auto | Base encoding (format) for phred score in FASTQs. Choices: `auto`, `phred33` or `phred64` (without hyphen). This is used for Trimmomatic only. It is `auto` by default, which means that Trimmomatic automatically detect it from FASTQs. Otherwise `-phred33` or `-phred64` will be passed to the Trimmomatic command line. Use this parameter if you get an error `Error: Unable to detect quality encoding` in Trimmomatic.
-`chip.use_bwa_mem_for_pe` | Boolean | false | For PE dataset, uise bwa mem instead of bwa aln.
+`chip.use_bwa_mem_for_pe` | Boolean | false | For `chip.aligner` as `bwa` and PE datasets only, use `bwa mem` instead of `bwa aln`. If read length of R1 FASTQ is shorter than `chip.bwa_mem_read_len_limit` (70 by default) then `bwa aln` will be used instead.
+`chip.bwa_mem_read_len_limit` | Int | 70 | For `chip.aligner` as `bwa` and PE dataset only, R1 FASTQ's read length limit for `bwa mem`.
+`chip.use_bowtie2_local_mode` | Boolean | false | Use bowtie2's local mode (adding `--local` to bowtie2 command line). If not defined, the default mode (end-to-end) will be used.
 `chip.custom_align_py` | File | | Python script for your custom aligner. See details about [how to use a custom aligner](#how-to-use-a-custom-aligner)
 
 
@@ -245,8 +247,8 @@ Base memory/disk is 4GB/20GB for most tasks.
 Parameter|Default|Description
 ---------|-------|-----------
 `chip.align_cpu` | 6 |
-`chip.align_bowtie2_mem_factor` | 0.15 | Multiplied to size of FASTQs to determine required memory
-`chip.align_bwa_mem_factor` | 0.30 | Multiplied to size of FASTQs to determine required memory
+`chip.align_bowtie2_mem_factor` | 0.15 | Multiplied to size of FASTQs to determine required memory. 5.0 + bowtie2_index_file_size + sum(all_fastqs) GB.
+`chip.align_bwa_mem_factor` | 1.0 | Multiplied to size of FASTQs to determine required memory. 5.0 + bwa_index_file_size + sum(all_fastqs) GB.
 `chip.align_time_hr` | 48 | Walltime (HPCs only)
 `chip.align_bowtie2_disk_factor` | 8.0 | Multiplied to size of FASTQs to determine required disk
 `chip.align_bwa_disk_factor` | 8.0 | Multiplied to size of FASTQs to determine required disk

example_input_json/template.full.json

@@ -22,6 +22,10 @@
     "chip.ctl_fastqs_rep2_R1" : [ "ctl2_R1.fastq.gz" ],
     "chip.ctl_fastqs_rep2_R2" : [ "ctl2_R2.fastq.gz" ],
 
+    "chip.use_bwa_mem_for_pe" : false,
+    "chip.bwa_mem_read_len_limit" : 70,
+    "chip.use_bowtie2_local_mode" : false,
+
     "chip.crop_length" : 0,
 
     "chip.mapq_thresh" : 30,
@@ -53,23 +57,23 @@
 
     "chip.align_cpu" : 6,
     "chip.align_bowtie2_mem_factor" : 0.15,
-    "chip.align_bwa_mem_factor" : 0.15,
+    "chip.align_bwa_mem_factor" : 1.0,
     "chip.align_time_hr" : 48,
     "chip.align_bowtie2_disk_factor" : 8.0,
     "chip.align_bwa_disk_factor" : 8.0,
 
     "chip.filter_cpu" : 4,
     "chip.filter_mem_factor" : 0.4,
     "chip.filter_time_hr" : 24,
-    "chip.filter_disk_factor" : 6.0,
+    "chip.filter_disk_factor" : 8.0,
 
     "chip.bam2ta_cpu" : 2,
     "chip.bam2ta_mem_factor" : 0.35,
     "chip.bam2ta_time_hr" : 6,
     "chip.bam2ta_disk_factor" : 4.0,
 
-    "chip.spr_mem_factor" : 4.5,
-    "chip.spr_disk_factor" : 6.0,
+    "chip.spr_mem_factor" : 13.5,
+    "chip.spr_disk_factor" : 18.0,
 
     "chip.jsd_cpu" : 4,
     "chip.jsd_mem_factor" : 0.1,
@@ -81,17 +85,17 @@
     "chip.xcor_time_hr" : 24,
     "chip.xcor_disk_factor" : 4.5,
 
-    "chip.subsample_ctl_mem_factor" : 7.0,
-    "chip.subsample_ctl_disk_factor" : 7.5,
+    "chip.subsample_ctl_mem_factor" : 14.0,
+    "chip.subsample_ctl_disk_factor" : 15.0,
 
     "chip.call_peak_cpu" : 6,
     "chip.call_peak_spp_mem_factor" : 5.0,
-    "chip.call_peak_macs2_mem_factor" : 2.5,
+    "chip.call_peak_macs2_mem_factor" : 5.0,
     "chip.call_peak_time_hr" : 72,
     "chip.call_peak_spp_disk_factor" : 5.0,
-    "chip.call_peak_macs2_disk_factor" : 15.0,
+    "chip.call_peak_macs2_disk_factor" : 30.0,
 
-    "chip.macs2_signal_track_mem_factor" : 6.0,
+    "chip.macs2_signal_track_mem_factor" : 12.0,
     "chip.macs2_signal_track_time_hr" : 24,
-    "chip.macs2_signal_track_disk_factor" : 40.0
+    "chip.macs2_signal_track_disk_factor" : 80.0
 }

scripts/requirements.txt

@@ -48,6 +48,7 @@ ptools_bin
 
 jsondiff ==1.1.1
 libgcc
+tbb ==2021.1.1
 requests
 ncurses
 gnuplot