diff --git a/README.md b/README.md
index ce231bc2d9517a8db52b47206db25ea3a71370f9..cf702bf337b100b48ab98405e67cb3b1e7c855d9 100644
--- a/README.md
+++ b/README.md
@@ -12,8 +12,9 @@ https://img.shields.io/badge/singularity-available-7E4C74.svg)
### Introduction
This pipeline is based on the [HiC-Pro workflow](https://github.com/nservant/HiC-Pro).
-It was designed to process Hi-C data from raw fastq files (paired-end Illumina data) to normalized contact maps. The current version supports digestion protocols.
-Support for other protocols is ongoing.
+It was designed to process Hi-C data from raw fastq files (paired-end Illumina data) to normalized contact maps.
+The current version supports most protocols, including digestion protocols as well as protocols that do not require restriction enzymes such as DNase Hi-C.
+In practice, this workflow was successfully applied to many data-sets including dilution Hi-C, in situ Hi-C, DNase Hi-C, Micro-C, capture-C, capture Hi-C or HiChip data.
The pipeline is built using [Nextflow](https://www.nextflow.io), a workflow tool to run tasks across multiple compute infrastructures in a very portable manner. It comes with docker / singularity containers making installation trivial and results highly reproducible.
diff --git a/bin/mapped_2hic_dnase.py b/bin/mapped_2hic_dnase.py
new file mode 100755
index 0000000000000000000000000000000000000000..36c5a605d0001de3775bb70e7934d06be7145797
--- /dev/null
+++ b/bin/mapped_2hic_dnase.py
@@ -0,0 +1,464 @@
+#!/usr/bin/env python
+
+# HiC-Pro
+# Copyleft 2015 Institut Curie
+# Author(s): Nicolas Servant, Eric Viara
+# Contact: nicolas.servant@curie.fr
+# This software is distributed without any guarantee under the terms of the
+# GNU General
+# Public License, either Version 2, June 1991 or Version 3, June 2007.
+
+"""
+Script to keep only valid pairs when no restriction enzyme are used (i.e. DNAse or Micro-HiC)
+"""
+
+import getopt
+import sys
+import os
+import re
+import pysam
+
+
+def usage():
+ """Usage function"""
+ print "Usage : python mapped_2hic_dnase.py"
+ print "-r/--mappedReadsFile <BAM/SAM file of mapped reads>"
+ print "[-o/--outputDir] <Output directory. Default is current directory>"
+ print "[-d/--minCisDist] <Minimum distance between intrachromosomal contact to consider>"
+ print "[-g/--gtag] <Genotype tag. If specified, this tag will be reported in the valid pairs output for allele specific classification>"
+ print "[-a/--all] <Write all additional output files, with information about the discarded reads (self-circle, dangling end, etc.)>"
+ print "[-v/--verbose] <Verbose>"
+ print "[-h/--help] <Help>"
+ return
+
+
+def get_args():
+ """Get argument"""
+ try:
+ opts, args = getopt.getopt(
+ sys.argv[1:],
+ "r:o:d:g:avh",
+ ["mappedReadsFile=",
+ "outputDir=", "minDist=", "gatg", "all", "verbose", "help"])
+ except getopt.GetoptError:
+ usage()
+ sys.exit(-1)
+ return opts
+
+
+def get_read_strand(read):
+ """
+ Conversion of read position to naive strand representation
+
+ Parameters
+ ----------
+ read : list
+ list of aligned reads
+ """
+ strand = "+"
+ if read.is_reverse:
+ strand = "-"
+ return strand
+
+
+def get_read_pos(read, st="start"):
+ """
+ Return the read position (zero-based) used for the intersection with
+ the restriction fragment
+
+ The 5' end is not a good choice for the reverse reads (which contain part
+ of the restriction site, and thus overlap the next restriction fragment)
+ Using the left-most position (5' for forward, 3' for reverse) or the
+ middle of the read should work but the middle of the reads might be more
+ safe
+
+ Parameters
+ -----------
+ read : list
+ list of aligned reads
+ """
+ if st == "middle":
+ pos = read.pos + int(read.alen/2)
+ elif st =="start":
+ pos = get_read_start(read)
+ elif st == "left":
+ pos = read.pos
+
+ return pos
+
+
+def get_read_start(read):
+ """
+ Return the 5' end of the read
+ """
+ if read.is_reverse:
+ pos = read.pos + read.alen -1
+ else:
+ pos = read.pos
+ return pos
+
+
+def get_ordered_reads(read1, read2):
+ """
+ Reorient reads
+
+ The sequencing is usually not oriented. Reorient the reads so that r1 is
+ always before r2
+
+ read1 = [AlignedRead]
+ read2 = [AlignedRead]
+ """
+ if read1.tid == read2.tid:
+ if get_read_pos(read1) < get_read_pos(read2):
+ r1 = read1
+ r2 = read2
+ else:
+ r1 = read2
+ r2 = read1
+ else:
+ if read1.tid < read2.tid:
+ r1 = read1
+ r2 = read2
+ else:
+ r1 = read2
+ r2 = read1
+
+ return r1, r2
+
+
+def isIntraChrom(read1, read2):
+ """
+ Return true is the reads pair is intrachromosomal
+
+ read1 : [AlignedRead]
+ read2 : [AlignedRead]
+
+ """
+ if read1.tid == read2.tid:
+ return True
+ else:
+ return False
+
+
+def get_valid_orientation(read1, read2):
+ """
+ Both reads are expected to be on the different restriction fragments
+
+ Check the orientation of reads ->-> / <-<- / -><- / <-->
+
+ read1 : [AlignedRead]
+ read2 : [AlignedRead]
+
+ """
+ # Get oriented reads
+ r1, r2 = get_ordered_reads(read1, read2)
+
+ direction = None
+ if get_read_strand(r1) == "+" and get_read_strand(r2) == "+":
+ direction = "FF"
+ elif get_read_strand(r1) == "-" and get_read_strand(r2) == "-":
+ direction = "RR"
+ elif get_read_strand(r1) == "+" and get_read_strand(r2) == "-":
+ direction = "FR"
+ elif get_read_strand(r1) == "-" and get_read_strand(r2) == "+":
+ direction = "RF"
+
+ return direction
+
+
+def get_cis_dist(read1, read2):
+ """
+ Calculte the size of the DNA fragment library
+
+ read1 : [AlignedRead]
+ read2 : [AlignedRead]
+
+ """
+ # Get oriented reads
+ ##r1, r2 = get_ordered_reads(read1, read2)
+ dist = None
+ if not r1.is_unmapped and not r2.is_unmapped:
+ ## Contact distances can be calculated for intrachromosomal reads only
+ if isIntraChrom(read1, read2):
+ r1pos = get_read_pos(read1)
+ r2pos = get_read_pos(read2)
+ dist = abs(r1pos - r2pos)
+ return dist
+
+
+def get_read_tag(read, tag):
+ for t in read.tags:
+ if t[0] == tag:
+ return t[1]
+ return None
+
+
+if __name__ == "__main__":
+ # Read command line arguments
+ opts = get_args()
+ verbose = False
+ allOutput = False
+ minInsertSize = None
+ maxInsertSize = None
+ minDist = None
+ outputDir = "."
+ gtag = None
+
+ if len(opts) == 0:
+ usage()
+ sys.exit()
+
+ for opt, arg in opts:
+ if opt in ("-h", "--help"):
+ usage()
+ sys.exit()
+ elif opt in ("-r", "--mappedReadsFile"):
+ mappedReadsFile = arg
+ elif opt in ("-o", "--outputDir"):
+ outputDir = arg
+ elif opt in ("-d", "--minCisDist"):
+ minDist = arg
+ elif opt in ("-g", "--gtag"):
+ gtag = arg
+ elif opt in ("-a", "--all"):
+ allOutput = True
+ elif opt in ("-v", "--verbose"):
+ verbose = True
+ else:
+ assert False, "unhandled option"
+
+ # Verbose mode
+ if verbose:
+ print "## overlapMapped2HiCFragments.py"
+ print "## mappedReadsFile=", mappedReadsFile
+ print "## minCisDist=", minDist
+ print "## allOuput=", allOutput
+ print "## verbose=", verbose, "\n"
+
+ # Initialize variables
+ reads_counter = 0
+ valid_counter = 0
+ valid_counter_FF = 0
+ valid_counter_RR = 0
+ valid_counter_FR = 0
+ valid_counter_RF = 0
+ single_counter = 0
+ dump_counter = 0
+ filt_counter = 0
+
+ # AS counter
+ G1G1_ascounter = 0
+ G2G2_ascounter = 0
+ G1U_ascounter = 0
+ UG1_ascounter = 0
+ G2U_ascounter = 0
+ UG2_ascounter = 0
+ G1G2_ascounter = 0
+ G2G1_ascounter = 0
+ UU_ascounter = 0
+ CF_ascounter = 0
+
+ baseReadsFile = os.path.basename(mappedReadsFile)
+ baseReadsFile = re.sub(r'\.bam$|\.sam$', '', baseReadsFile)
+
+ # Open handlers for output files
+ handle_valid = open(outputDir + '/' + baseReadsFile + '.validPairs', 'w')
+
+ if allOutput:
+ handle_dump = open(outputDir + '/' + baseReadsFile + '.DumpPairs', 'w')
+ handle_single = open(outputDir + '/' + baseReadsFile + '.SinglePairs','w')
+ handle_filt = open(outputDir + '/' + baseReadsFile + '.FiltPairs','w')
+
+ # Read the SAM/BAM file
+ if verbose:
+ print "## Opening SAM/BAM file '", mappedReadsFile, "'..."
+ samfile = pysam.Samfile(mappedReadsFile, "rb")
+
+ # Reads are 0-based too (for both SAM and BAM format)
+ # Loop on all reads
+ for read in samfile.fetch(until_eof=True):
+ reads_counter += 1
+ cur_handler = None
+ interactionType = None
+ htag = ""
+
+ # First mate
+ if read.is_read1:
+ r1 = read
+ if not r1.is_unmapped:
+ r1_chrom = samfile.getrname(r1.tid)
+ else:
+ r1_chrom = None
+
+ # Second mate
+ elif read.is_read2:
+ r2 = read
+ if not r2.is_unmapped:
+ r2_chrom = samfile.getrname(r2.tid)
+ else:
+ r2_chrom = None
+
+ if isIntraChrom(r1,r2):
+ dist = get_cis_dist(r1, r2)
+ else:
+ dist = None
+
+ # Check singleton
+ if r1.is_unmapped or r2.is_unmapped:
+ interactionType = "SI"
+ single_counter += 1
+ cur_handler = handle_single if allOutput else None
+
+ # Check Distance criteria - Filter
+ if (minDist is not None and dist is not None and dist < int(minDist)):
+ interactionType = "FILT"
+ filt_counter += 1
+ cur_handler = handle_filt if allOutput else None
+
+ # By default pair is valid
+ if interactionType == None:
+ interactionType = "VI"
+ valid_counter += 1
+ cur_handler = handle_valid
+ validType = get_valid_orientation(r1, r2)
+ if validType == "RR":
+ valid_counter_RR += 1
+ elif validType == "FF":
+ valid_counter_FF += 1
+ elif validType == "FR":
+ valid_counter_FR += 1
+ elif validType == "RF":
+ valid_counter_RF += 1
+ else:
+ interactionType = "DUMP"
+ dump_counter += 1
+ cur_handler = handle_dump if allOutput else None
+
+
+
+ # Split valid pairs based on XA tag
+ if gtag is not None:
+ r1as = get_read_tag(r1, gtag)
+ r2as = get_read_tag(r2, gtag)
+
+ if r1as == 1 and r2as == 1:
+ G1G1_ascounter += 1
+ elif r1as == 2 and r2as == 2:
+ G2G2_ascounter += 1
+ elif r1as == 1 and r2as == 0:
+ G1U_ascounter += 1
+ elif r1as == 0 and r2as == 1:
+ UG1_ascounter += 1
+ elif r1as == 2 and r2as == 0:
+ G2U_ascounter += 1
+ elif r1as == 0 and r2as == 2:
+ UG2_ascounter += 1
+ elif r1as == 1 and r2as == 2:
+ G1G2_ascounter += 1
+ elif r1as == 2 and r2as == 1:
+ G2G1_ascounter += 1
+ elif r1as == 3 or r2as == 3:
+ CF_ascounter += 1
+ else:
+ UU_ascounter += 1
+
+
+ if cur_handler is not None:
+ if not r1.is_unmapped and not r2.is_unmapped:
+
+ ##reorient reads to ease duplicates removal
+ or1, or2 = get_ordered_reads(r1, r2)
+ or1_chrom = samfile.getrname(or1.tid)
+ or2_chrom = samfile.getrname(or2.tid)
+
+ ##reset as tag now that the reads are oriented
+ r1as = get_read_tag(or1, gtag)
+ r2as = get_read_tag(or2, gtag)
+ if gtag is not None:
+ htag = str(r1as)+"-"+str(r2as)
+
+ cur_handler.write(
+ or1.qname + "\t" +
+ or1_chrom + "\t" +
+ str(get_read_pos(or1)+1) + "\t" +
+ str(get_read_strand(or1)) + "\t" +
+ or2_chrom + "\t" +
+ str(get_read_pos(or2)+1) + "\t" +
+ str(get_read_strand(or2)) + "\t" +
+ "NA" + "\t" + ##dist
+ "NA" + "\t" + ##resfrag1
+ "NA" + "\t" + ##resfrag2
+ str(or1.mapping_quality) + "\t" +
+ str(or2.mapping_quality) + "\t" +
+ str(htag) + "\n")
+
+ elif r2.is_unmapped and not r1.is_unmapped:
+ cur_handler.write(
+ r1.qname + "\t" +
+ r1_chrom + "\t" +
+ str(get_read_pos(r1)+1) + "\t" +
+ str(get_read_strand(r1)) + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ str(r1.mapping_quality) + "\t" +
+ "*" + "\n")
+ elif r1.is_unmapped and not r2.is_unmapped:
+ cur_handler.write(
+ r2.qname + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ r2_chrom + "\t" +
+ str(get_read_pos(r2)+1) + "\t" +
+ str(get_read_strand(r2)) + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ "*" + "\t" +
+ str(r2.mapping_quality) + "\n")
+
+ if (reads_counter % 100000 == 0 and verbose):
+ print "##", reads_counter
+
+ # Close handler
+ handle_valid.close()
+ if allOutput:
+ handle_dump.close()
+ handle_single.close()
+ handle_filt.close()
+
+ # Write stats file
+ handle_stat = open(outputDir + '/' + baseReadsFile + '.RSstat', 'w')
+ handle_stat.write("## Hi-C processing - no restriction fragments\n")
+ handle_stat.write("Valid_interaction_pairs\t" + str(valid_counter) + "\n")
+ handle_stat.write(
+ "Valid_interaction_pairs_FF\t" + str(valid_counter_FF) + "\n")
+ handle_stat.write(
+ "Valid_interaction_pairs_RR\t" + str(valid_counter_RR) + "\n")
+ handle_stat.write(
+ "Valid_interaction_pairs_RF\t" + str(valid_counter_RF) + "\n")
+ handle_stat.write(
+ "Valid_interaction_pairs_FR\t" + str(valid_counter_FR) + "\n")
+ handle_stat.write("Single-end_pairs\t" + str(single_counter) + "\n")
+ handle_stat.write("Filtered_pairs\t" + str(filt_counter) + "\n")
+ handle_stat.write("Dumped_pairs\t" + str(dump_counter) + "\n")
+
+ ## Write AS report
+ if gtag is not None:
+ handle_stat.write("## ======================================\n")
+ handle_stat.write("## Allele specific information\n")
+ handle_stat.write("Valid_pairs_from_ref_genome_(1-1)\t" + str(G1G1_ascounter) + "\n")
+ handle_stat.write("Valid_pairs_from_ref_genome_with_one_unassigned_mate_(0-1/1-0)\t" + str(UG1_ascounter+G1U_ascounter) + "\n")
+ handle_stat.write("Valid_pairs_from_alt_genome_(2-2)\t" + str(G2G2_ascounter) + "\n")
+ handle_stat.write("Valid_pairs_from_alt_genome_with_one_unassigned_mate_(0-2/2-0)\t" + str(UG2_ascounter+G2U_ascounter) + "\n")
+ handle_stat.write("Valid_pairs_from_alt_and_ref_genome_(1-2/2-1)\t" + str(G1G2_ascounter+G2G1_ascounter) + "\n")
+ handle_stat.write("Valid_pairs_with_both_unassigned_mated_(0-0)\t" + str(UU_ascounter) + "\n")
+ handle_stat.write("Valid_pairs_with_at_least_one_conflicting_mate_(3-)\t" + str(CF_ascounter) + "\n")
+
+ handle_stat.close()
+
+
diff --git a/conf/hicpro.config b/conf/hicpro.config
index b4eac51dacad0bb55f4695d9bfcc0a29551d4211..0a2c9b9e0db09f4f9861ba353b84a534820aba38 100644
--- a/conf/hicpro.config
+++ b/conf/hicpro.config
@@ -21,7 +21,7 @@ params {
min_restriction_fragment_size =
max_restriction_fragment_size =
min_insert_size =
- max_insert_size =
+ max_insert_size =
// Hi-C Processing
min_cis_dist =
@@ -29,7 +29,7 @@ params {
rm_multi = true
rm_dup = true
- bins_size = '1000000,500000'
+ bin_size = '1000000,500000'
ice_max_iter = 100
ice_filer_low_count_perc = 0.02
diff --git a/docs/usage.md b/docs/usage.md
index 15405acf307a67d179bad24da88ea1dc1a7aa600..853c38414b6e53090c3d9b0f19e849a0972b1243 100644
--- a/docs/usage.md
+++ b/docs/usage.md
@@ -34,7 +34,9 @@
* [`--max_restriction_fragment_size`](#--max_restriction_fragment_size)
* [`--min_insert_size`](#--min_insert_size)
* [`--max_insert_size`](#--max_insert_size)
- * [Hi-C Processing](#hi-c-processing)
+ * [DNase Hi-C](#dnase-hi-c)
+ * [`--dnase`](#--dnase)
+ * [Hi-C Processing](#hi-c-processing)
* [`--min_cis_dist`](#--min_cis_dist)
* [`--rm_singleton`](#--rm_singleton)
* [`--rm_dup`](#--rm_dup)
@@ -149,7 +151,7 @@ Please note the following requirements:
If left unspecified, a default pattern is used: `data/*{1,2}.fastq.gz`
-## Reference genomes
+## Reference genomes and annotation files
The pipeline config files come bundled with paths to the illumina iGenomes reference index files. If running with docker or AWS, the configuration is set up to use the [AWS-iGenomes](https://ewels.github.io/AWS-iGenomes/) resource.
@@ -223,7 +225,7 @@ If not specified, this file will be automatically created by the pipeline. In th
```
```bash
---bwt2_index '[path to chromosome size file]'
+--chromosome_size '[path to chromosome size file]'
```
### `--restriction_fragments`
@@ -251,7 +253,7 @@ Note that the `--restriction_site` parameter is mandatory to create this file.
The following options are defined in the `hicpro.config` file, and can be updated either using a custom configuration file (see `-c` option) or using command line parameter.
-## Reads mapping
+### Reads mapping
The reads mapping is currently based on the two-steps strategy implemented in the HiC-pro pipeline. The idea is to first align reads from end-to-end.
Reads that do not aligned are then trimmed at the ligation site, and their 5' end is re-aligned to the reference genome.
@@ -281,7 +283,7 @@ Minimum mapping quality. Reads with lower quality are discarded. Default: 10
--min_mapq '[Minimum quality value]'
```
-## Digestion Hi-C
+### Digestion Hi-C
#### `--restriction_site`
@@ -340,7 +342,18 @@ Maximum reads insert size. Longer 3C products are discarded. Default: ''
--max_insert_size '[numeric]'
```
-## Hi-C processing
+### DNAse Hi-C
+
+#### `--dnase`
+
+In DNAse Hi-C mode, all options related to digestion Hi-C (see previous section) are ignored.
+In this case, it is highly recommanded to use the `--min_cis_dist` parameter to remove spurious ligation products.
+
+```bash
+--dnase'
+```
+
+### Hi-C processing
#### `--min_cis_dist`
@@ -376,7 +389,7 @@ If specified, reads that aligned multiple times on the genome are discarded. Not
## Genome-wide contact maps
-#### `--bins_size`
+#### `--bin_size`
Resolution of contact maps to generate (space separated). Default:'1000000,500000'
diff --git a/main.nf b/main.nf
index 8ac8bd57b2e2bb1b2674d2e0fb457a7ae935f877..85ac949ebe87fa5431e470f7e3d7b58715aacb64 100644
--- a/main.nf
+++ b/main.nf
@@ -36,7 +36,7 @@ def helpMessage() {
-profile Configuration profile to use. Can use multiple (comma separated)
Available: conda, docker, singularity, awsbatch, test and more.
- References If not specified in the configuration file or you wish to overwrite any of the references.
+ References: If not specified in the configuration file or you wish to overwrite any of the references.
--genome Name of iGenomes reference
--bwt2_index Path to Bowtie2 index
--fasta Path to Fasta reference
@@ -50,11 +50,13 @@ def helpMessage() {
--restriction_site Cutting motif(s) of restriction enzyme(s) (comma separated)
--ligation_site Ligation motifs to trim (comma separated)
-
--min_restriction_fragment_size Minimum size of restriction fragments to consider
--max_restriction_framgnet_size Maximum size of restriction fragmants to consider
--min_insert_size Minimum insert size of mapped reads to consider
--max_insert_size Maximum insert size of mapped reads to consider
+
+ --dnase Run DNase Hi-C mode. All options related to restriction fragments are not considered
+
--min_cis_dist Minimum intra-chromosomal distance to consider
--rm_singleton Remove singleton reads
--rm_multi Remove multi-mapped reads
@@ -221,7 +223,7 @@ else {
}
// Resolutions for contact maps
-map_res = Channel.from( params.bins_size.tokenize(',') )
+map_res = Channel.from( params.bin_size.tokenize(',') )
// Stage config files
ch_multiqc_config = Channel.fromPath(params.multiqc_config)
@@ -358,7 +360,7 @@ if(!params.chromosome_size && params.fasta){
}
}
-if(!params.restriction_fragments && params.fasta){
+if(!params.restriction_fragments && params.fasta && !params.dnase){
process getRestrictionFragments {
tag "$fasta [${params.restriction_site}]"
publishDir path: { params.saveReference ? "${params.outdir}/reference_genome" : params.outdir },
@@ -527,39 +529,66 @@ process combine_mapped_files{
* STEP2 - DETECT VALID PAIRS
*/
-
-process get_valid_interaction{
- tag "$sample"
- publishDir "${params.outdir}/hic_results/data", mode: 'copy',
+if (!params.dnase){
+ process get_valid_interaction{
+ tag "$sample"
+ publishDir "${params.outdir}/hic_results/data", mode: 'copy',
saveAs: {filename -> filename.indexOf("*stat") > 0 ? "stats/$filename" : "$filename"}
- input:
- set val(sample), file(pe_bam) from paired_bam
- file frag_file from res_frag_file.collect()
-
- output:
- set val(sample), file("*.validPairs") into valid_pairs
- set val(sample), file("*.validPairs") into valid_pairs_4cool
- set val(sample), file("*RSstat") into all_rsstat
-
- script:
-
- if (params.splitFastq){
- sample = sample.toString() - ~/(\.[0-9]+)$/
- }
-
- def opts = ""
- if ("$params.min_cis_dist".isInteger()) opts="${opts} -d ${params.min_cis_dist}"
- if ("$params.min_insert_size".isInteger()) opts="${opts} -s ${params.min_insert_size}"
- if ("$params.max_insert_size".isInteger()) opts="${opts} -l ${params.max_insert_size}"
- if ("$params.min_restriction_fragment_size".isInteger()) opts="${opts} -t ${params.min_restriction_fragment_size}"
- if ("$params.max_restriction_fragment_size".isInteger()) opts="${opts} -m ${params.max_restriction_fragment_size}"
+ input:
+ set val(sample), file(pe_bam) from paired_bam
+ file frag_file from res_frag_file.collect()
+
+ output:
+ set val(sample), file("*.validPairs") into valid_pairs
+ set val(sample), file("*.validPairs") into valid_pairs_4cool
+ set val(sample), file("*RSstat") into all_rsstat
+
+ script:
+ if (params.splitFastq){
+ sample = sample.toString() - ~/(\.[0-9]+)$/
+ }
+
+ def opts = ""
+ if ("$params.min_cis_dist".isInteger()) opts="${opts} -d ${params.min_cis_dist}"
+ if ("$params.min_insert_size".isInteger()) opts="${opts} -s ${params.min_insert_size}"
+ if ("$params.max_insert_size".isInteger()) opts="${opts} -l ${params.max_insert_size}"
+ if ("$params.min_restriction_fragment_size".isInteger()) opts="${opts} -t ${params.min_restriction_fragment_size}"
+ if ("$params.max_restriction_fragment_size".isInteger()) opts="${opts} -m ${params.max_restriction_fragment_size}"
+
+ """
+ mapped_2hic_fragments.py -f ${frag_file} -r ${pe_bam} ${opts}
+ """
+ }
+}
+else{
+ process get_valid_interaction_dnase{
+ tag "$sample"
+ publishDir "${params.outdir}/hic_results/data", mode: 'copy',
+ saveAs: {filename -> filename.indexOf("*stat") > 0 ? "stats/$filename" : "$filename"}
- """
- mapped_2hic_fragments.py -f ${frag_file} -r ${pe_bam} ${opts}
- """
+ input:
+ set val(sample), file(pe_bam) from paired_bam
+
+ output:
+ set val(sample), file("*.validPairs") into valid_pairs
+ set val(sample), file("*.validPairs") into valid_pairs_4cool
+ set val(sample), file("*RSstat") into all_rsstat
+
+ script:
+ if (params.splitFastq){
+ sample = sample.toString() - ~/(\.[0-9]+)$/
+ }
+
+ def opts = ""
+ if ("$params.min_cis_dist".isInteger()) opts="${opts} -d ${params.min_cis_dist}"
+ """
+ mapped_2hic_dnase.py -r ${pe_bam} ${opts}
+ """
+ }
}
+
/*
* STEP3 - BUILD MATRIX
*/
diff --git a/nextflow.config b/nextflow.config
index 5cd6dfec6888a0727fdfc7c46862313328d60c3b..32486aab168ccff7647220cd75a61f72f695fbc5 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -21,6 +21,7 @@ params {
restriction_fragments = false
skip_cool = false
skip_multiqc = false
+ dnase = false
// Boilerplate options
name = false