Merge pull request #231 from ENCODE-DCC/dev

v1.9.0

Author: leepc12

.circleci/config.yml

@@ -72,7 +72,7 @@ jobs:
           name: build image
           command: |
             source ${BASH_ENV}
-            export DOCKER_CACHE_TAG=v1.4.0.1
+            export DOCKER_CACHE_TAG=v1.8.1
             echo "pulling ${DOCKER_CACHE_TAG}!"
             docker pull encodedcc/chip-seq-pipeline:${DOCKER_CACHE_TAG}
             docker login -u=${DOCKERHUB_USER} -p=${DOCKERHUB_PASS}

chip.wdl

@@ -1,16 +1,16 @@
 version 1.0
 
 workflow chip {
-    String pipeline_ver = 'v1.8.1'
+    String pipeline_ver = 'v1.9.0'
 
     meta {
-        version: 'v1.8.1'
+        version: 'v1.9.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.8.1'
-        caper_singularity: 'docker://encodedcc/chip-seq-pipeline:v1.8.1'
+        caper_docker: 'encodedcc/chip-seq-pipeline:v1.9.0'
+        caper_singularity: 'docker://encodedcc/chip-seq-pipeline:v1.9.0'
         croo_out_def: 'https://storage.googleapis.com/encode-pipeline-output-definition/chip.croo.v5.json'
 
         parameter_group: {
@@ -162,6 +162,7 @@ workflow chip {
         Int xcor_subsample_reads = 15000000
         Int xcor_exclusion_range_min = -500
         Int? xcor_exclusion_range_max
+        Int pseudoreplication_random_seed = 0
 
         # group: peak_calling
         Int ctl_depth_limit = 200000000
@@ -728,7 +729,11 @@ workflow chip {
             group: 'alignment',
             help: 'For run_spp.R -s. If not defined default value of `max(read length + 10, 50)` for TF and `max(read_len + 10, 100)` for histone are used'
         }
-
+        pseudoreplication_random_seed: {
+            description: 'Random seed (positive integer) used for pseudo-replication (shuffling reads in TAG-ALIGN and then split it into two).',
+            group: 'alignment',
+            help: 'Pseudo-replication (task spr) is done by using GNU "shuf --random-source=sha256(random_seed)". If this parameter == 0, then pipeline uses input TAG-ALIGN file\'s size (in bytes) for the random_seed.'
+        }
         ctl_depth_limit: {
             description: 'Hard limit for chosen control\'s depth.',
             group: 'peak_calling',
@@ -1262,6 +1267,7 @@ workflow chip {
             call spr { input :
                 ta = ta_,
                 paired_end = paired_end_,
+                pseudoreplication_random_seed = pseudoreplication_random_seed,
                 mem_factor = spr_mem_factor,
                 disk_factor = spr_disk_factor,
             }
@@ -2280,6 +2286,7 @@ task spr {
     input {
         File? ta
         Boolean paired_end
+        Int pseudoreplication_random_seed
 
         Float mem_factor
         Float disk_factor
@@ -2292,6 +2299,7 @@ task spr {
         set -e
         python3 $(which encode_task_spr.py) \
             ${ta} \
+            ${'--pseudoreplication-random-seed ' + pseudoreplication_random_seed} \
             ${if paired_end then '--paired-end' else ''}
     }
     output {

dev/test/test_task/test_spr.json

@@ -5,5 +5,10 @@
     "test_spr.ref_pe_ta_pr1" : "chip-seq-pipeline-test-data/ref_output/test_spr/pe/rep1-R1.subsampled.67.merged.nodup.pr1.tagAlign.gz",
     "test_spr.ref_pe_ta_pr2" : "chip-seq-pipeline-test-data/ref_output/test_spr/pe/rep1-R1.subsampled.67.merged.nodup.pr2.tagAlign.gz",
     "test_spr.ref_se_ta_pr1" : "chip-seq-pipeline-test-data/ref_output/test_spr/se/rep1.subsampled.25.merged.nodup.pr1.tagAlign.gz",
-    "test_spr.ref_se_ta_pr2" : "chip-seq-pipeline-test-data/ref_output/test_spr/se/rep1.subsampled.25.merged.nodup.pr2.tagAlign.gz"
+    "test_spr.ref_se_ta_pr2" : "chip-seq-pipeline-test-data/ref_output/test_spr/se/rep1.subsampled.25.merged.nodup.pr2.tagAlign.gz",
+
+    "test_spr.ref_pe_seed_10_ta_pr1" : "chip-seq-pipeline-test-data/ref_output/test_spr/pe/pseudoreplication_random_seed_10/rep1-R1.subsampled.67.merged.nodup.pr1.tagAlign.gz",
+    "test_spr.ref_pe_seed_10_ta_pr2" : "chip-seq-pipeline-test-data/ref_output/test_spr/pe/pseudoreplication_random_seed_10/rep1-R1.subsampled.67.merged.nodup.pr2.tagAlign.gz",
+    "test_spr.ref_se_seed_10_ta_pr1" : "chip-seq-pipeline-test-data/ref_output/test_spr/se/pseudoreplication_random_seed_10/rep1.subsampled.25.merged.nodup.pr1.tagAlign.gz",
+    "test_spr.ref_se_seed_10_ta_pr2" : "chip-seq-pipeline-test-data/ref_output/test_spr/se/pseudoreplication_random_seed_10/rep1.subsampled.25.merged.nodup.pr2.tagAlign.gz"
 }

dev/test/test_task/test_spr.wdl

@@ -4,28 +4,46 @@ import 'compare_md5sum.wdl' as compare_md5sum
 
 workflow test_spr {
     input {
-        String pe_ta
-        String se_ta
+        File pe_ta
+        File se_ta
 
-        String ref_pe_ta_pr1
-        String ref_pe_ta_pr2
-        String ref_se_ta_pr1
-        String ref_se_ta_pr2
+        File ref_pe_ta_pr1
+        File ref_pe_ta_pr2
+        File ref_se_ta_pr1
+        File ref_se_ta_pr2
+        File ref_pe_seed_10_ta_pr1
+        File ref_pe_seed_10_ta_pr2
+        File ref_se_seed_10_ta_pr1
+        File ref_se_seed_10_ta_pr2
     }
     Float spr_mem_factor = 0.0
     Float spr_disk_factor = 6.0
 
     call chip.spr as pe_spr { input :
         ta = pe_ta,
         paired_end = true,
-
+        pseudoreplication_random_seed = 0,
         mem_factor = spr_mem_factor,
         disk_factor = spr_disk_factor,
     }    
     call chip.spr as se_spr { input :
         ta = se_ta,
         paired_end = false,
-
+        pseudoreplication_random_seed = 0,
+        mem_factor = spr_mem_factor,
+        disk_factor = spr_disk_factor,
+    }
+    call chip.spr as pe_spr_seed_10 { input :
+        ta = pe_ta,
+        paired_end = true,
+        pseudoreplication_random_seed = 10,
+        mem_factor = spr_mem_factor,
+        disk_factor = spr_disk_factor,
+    }
+    call chip.spr as se_spr_seed_10 { input :
+        ta = se_ta,
+        paired_end = false,
+        pseudoreplication_random_seed = 10,
         mem_factor = spr_mem_factor,
         disk_factor = spr_disk_factor,
     }
@@ -36,18 +54,30 @@ workflow test_spr {
             'pe_spr_pr2',
             'se_spr_pr1',
             'se_spr_pr2',
+            'pe_spr_seed_10_pr1',
+            'pe_spr_seed_10_pr2',
+            'se_spr_seed_10_pr1',
+            'se_spr_seed_10_pr2',
         ],
         files = [
             pe_spr.ta_pr1,
             pe_spr.ta_pr2,
             se_spr.ta_pr1,
             se_spr.ta_pr2,
+            pe_spr_seed_10.ta_pr1,
+            pe_spr_seed_10.ta_pr2,
+            se_spr_seed_10.ta_pr1,
+            se_spr_seed_10.ta_pr2,
         ],
         ref_files = [
             ref_pe_ta_pr1,
             ref_pe_ta_pr2,
             ref_se_ta_pr1,
             ref_se_ta_pr2,
+            ref_pe_seed_10_ta_pr1,
+            ref_pe_seed_10_ta_pr2,
+            ref_se_seed_10_ta_pr1,
+            ref_se_seed_10_ta_pr2,
         ],
     }
 }

docs/input.md

@@ -225,6 +225,7 @@ Parameter|Type | Description
 Parameter|Default|Description
 ---------|-------|-----------
 `chip.filter_chrs` | `[]` (empty array of string) | Array of chromosome names to be filtered out from a final (filtered/nodup) BAM. No chromosomes are filtered out by default.
+`chip.pseudoreplication_random_seed` | `0` | Random seed (positive integer) used for pseudo-replication (shuffling reads in TAG-ALIGN and then split it into two). If `0` then TAG-ALIGN file's size (in bytes) is used for random seed.
 
 ## Resource parameters
 

example_input_json/template.full.json

@@ -35,6 +35,7 @@
     "chip.subsample_reads" : 0,
     "chip.ctl_subsample_reads" : 0,
     "chip.xcor_subsample_reads" : 15000000,
+    "chip.pseudoreplication_random_seed" : 0,
 
     "chip.xcor_trim_bp" : 50,
     "chip.use_filt_pe_ta_for_xcor" : false,

src/encode_task_spr.py

@@ -18,6 +18,11 @@ def parse_arguments():
                         help='Path for TAGALIGN file.')
     parser.add_argument('--paired-end', action="store_true",
                         help='Paired-end TAGALIGN.')
+    parser.add_argument('--pseudoreplication-random-seed',
+                        type=int, default=0,
+                        help='Set it to 0 to use file\'s size (in bytes) as random seed.'
+                             'Otherwise this seed will be used for GNU shuf --random-source=sha256(seed).'
+                             'It is useful when random seed based on input file size does not work.')
     parser.add_argument('--out-dir', default='', type=str,
                         help='Output directory.')
     parser.add_argument('--log-level', default='INFO',
@@ -32,7 +37,7 @@ def parse_arguments():
     return args
 
 
-def spr_se(ta, out_dir):
+def spr_se(ta, pseudoreplication_random_seed, out_dir):
     prefix = os.path.join(out_dir,
                           os.path.basename(strip_ext_ta(ta)))
     tmp_pr1 = '{}.00'.format(prefix)
@@ -41,35 +46,42 @@ def spr_se(ta, out_dir):
     ta_pr2 = '{}.pr2.tagAlign.gz'.format(prefix)
     nlines = int((get_num_lines(ta)+1)/2)
 
+    if pseudoreplication_random_seed == 0:
+        random_seed = run_shell_cmd('zcat -f {ta} | wc -c'.format(ta=ta))
+        log.info(
+            'Using input file\'s size {random_seed} as random seed for pseudoreplication.'.format(
+                random_seed=random_seed,
+            )
+        )
+    else:
+        random_seed = pseudoreplication_random_seed
+        log.info(
+            'Using a fixed integer {random_seed} as random seed for pseudoreplication.'.format(
+                random_seed=random_seed,
+            )
+        )
+
     # bash-only
-    cmd1 = 'zcat {} | shuf --random-source=<(openssl enc '
-    cmd1 += '-aes-256-ctr -pass pass:$(zcat -f {} | wc -c) '
-    cmd1 += '-nosalt </dev/zero 2>/dev/null) | '
-    cmd1 += 'split -d -l {} - {}.'
-    cmd1 = cmd1.format(
-        ta,
-        ta,
-        nlines,
-        prefix)
-    run_shell_cmd(cmd1)
-
-    cmd2 = 'gzip -nc {} > {}'
-    cmd2 = cmd2.format(
-        tmp_pr1,
-        ta_pr1)
-    run_shell_cmd(cmd2)
-
-    cmd3 = 'gzip -nc {} > {}'
-    cmd3 = cmd3.format(
-        tmp_pr2,
-        ta_pr2)
-    run_shell_cmd(cmd3)
+    run_shell_cmd(
+        'zcat {ta} | shuf --random-source=<(openssl enc '
+        '-aes-256-ctr -pass pass:{random_seed} '
+        '-nosalt </dev/zero 2>/dev/null) | '
+        'split -d -l {nlines} - {prefix}.'.format(
+            ta=ta,
+            random_seed=random_seed,
+            nlines=nlines,
+            prefix=prefix,
+        )
+    )
+
+    run_shell_cmd('gzip -nc {tmp_pr1} > {ta_pr1}'.format(tmp_pr1=tmp_pr1, ta_pr1=ta_pr1))
+    run_shell_cmd('gzip -nc {tmp_pr2} > {ta_pr2}'.format(tmp_pr2=tmp_pr2, ta_pr2=ta_pr2))
 
     rm_f([tmp_pr1, tmp_pr2])
     return ta_pr1, ta_pr2
 
 
-def spr_pe(ta, out_dir):
+def spr_pe(ta, pseudoreplication_random_seed, out_dir):
     prefix = os.path.join(out_dir,
                           os.path.basename(strip_ext_ta(ta)))
     tmp_pr1 = '{}.00'.format(prefix)
@@ -78,40 +90,57 @@ def spr_pe(ta, out_dir):
     ta_pr2 = '{}.pr2.tagAlign.gz'.format(prefix)
     nlines = int((get_num_lines(ta)/2+1)/2)
 
+    if pseudoreplication_random_seed == 0:
+        random_seed = run_shell_cmd('zcat -f {ta} | wc -c'.format(ta=ta))
+        log.info(
+            'Using input file\'s size {random_seed} as random seed for pseudoreplication.'.format(
+                random_seed=random_seed,
+            )
+        )
+    else:
+        random_seed = pseudoreplication_random_seed
+        log.info(
+            'Using a fixed integer {random_seed} as random seed for pseudoreplication.'.format(
+                random_seed=random_seed,
+            )
+        )
+
     # bash-only
-    cmd1 = 'zcat -f {} | sed \'N;s/\\n/\\t/\' | '
-    cmd1 += 'shuf --random-source=<(openssl enc -aes-256-ctr '
-    cmd1 += '-pass pass:$(zcat -f {} | wc -c) '
-    cmd1 += '-nosalt </dev/zero 2>/dev/null) | '
-    cmd1 += 'split -d -l {} - {}.'
-    cmd1 = cmd1.format(
-        ta,
-        ta,
-        nlines,
-        prefix)
-    run_shell_cmd(cmd1)
-
-    cmd2 = 'zcat -f {} | '
-    cmd2 += 'awk \'BEGIN{{OFS="\\t"}} '
-    cmd2 += '{{printf "%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n'
-    cmd2 += '%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n",'
-    cmd2 += '$1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}}\' | '
-    cmd2 += 'gzip -nc > {}'
-    cmd2 = cmd2.format(
-        tmp_pr1,
-        ta_pr1)
-    run_shell_cmd(cmd2)
-
-    cmd3 = 'zcat -f {} | '
-    cmd3 += 'awk \'BEGIN{{OFS="\\t"}} '
-    cmd3 += '{{printf "%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n'
-    cmd3 += '%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n",'
-    cmd3 += '$1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}}\' | '
-    cmd3 += 'gzip -nc > {}'
-    cmd3 = cmd3.format(
-        tmp_pr2,
-        ta_pr2)
-    run_shell_cmd(cmd3)
+    run_shell_cmd(
+        'zcat -f {ta} | sed \'N;s/\\n/\\t/\' | '
+        'shuf --random-source=<(openssl enc -aes-256-ctr '
+        '-pass pass:{random_seed} -nosalt </dev/zero 2>/dev/null) | '
+        'split -d -l {nlines} - {prefix}.'.format(
+            ta=ta,
+            random_seed=random_seed,
+            nlines=nlines,
+            prefix=prefix,
+        )
+    )
+
+    run_shell_cmd(
+        'zcat -f {tmp_pr1} | '
+        'awk \'BEGIN{{OFS="\\t"}} '
+        '{{printf "%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n'
+        '%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n",'
+        '$1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}}\' | '
+        'gzip -nc > {ta_pr1}'.format(
+            tmp_pr1=tmp_pr1,
+            ta_pr1=ta_pr1,
+        )
+    )
+
+    run_shell_cmd(
+        'zcat -f {tmp_pr2} | '
+        'awk \'BEGIN{{OFS="\\t"}} '
+        '{{printf "%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n'
+        '%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n",'
+        '$1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}}\' | '
+        'gzip -nc > {ta_pr2}'.format(
+            tmp_pr2=tmp_pr2,
+            ta_pr2=ta_pr2,
+        )
+    )
 
     rm_f([tmp_pr1, tmp_pr2])
     return ta_pr1, ta_pr2
@@ -125,9 +154,13 @@ def main():
 
     log.info('Making self-pseudo replicates...')
     if args.paired_end:
-        ta_pr1, ta_pr2 = spr_pe(args.ta, args.out_dir)
+        ta_pr1, ta_pr2 = spr_pe(
+            args.ta, args.pseudoreplication_random_seed, args.out_dir,
+        )
     else:
-        ta_pr1, ta_pr2 = spr_se(args.ta, args.out_dir)
+        ta_pr1, ta_pr2 = spr_se(
+            args.ta, args.pseudoreplication_random_seed, args.out_dir,
+        )
 
     log.info('List all files in output directory...')
     ls_l(args.out_dir)