diff --git a/.travis.yml b/.travis.yml
index 966de4b291e58f20cafbedcf791057f2789dbd8a..2eaea150d6ecf198773030802425b4a34edd7219 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -11,9 +11,9 @@ before_install:
# PRs to master are only ok if coming from dev branch
- '[ $TRAVIS_PULL_REQUEST = "false" ] || [ $TRAVIS_BRANCH != "master" ] || ([ $TRAVIS_PULL_REQUEST_SLUG = $TRAVIS_REPO_SLUG ] && [ $TRAVIS_PULL_REQUEST_BRANCH = "dev" ])'
# Pull the docker image first so the test doesn't wait for this
- - docker pull nservant/nf-core-hic:dev
+ - docker pull nfcore/hic:dev
# Fake the tag locally so that the pipeline runs properly
- - docker tag nservant/nf-core-hic:dev nservant/nf-core-hic:latest
+ - docker tag nfcore/hic:dev nfcore/hic:dev
install:
# Install Nextflow
diff --git a/Dockerfile b/Dockerfile
index 4d28385b74a7fd65cebe0848cd39097a9518864a..490ee3637819dd1ba8028a9d57e978fb794b1f2a 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -2,6 +2,9 @@ FROM nfcore/base
LABEL authors="Nicolas Servant" \
description="Docker image containing all requirements for nf-core/hic pipeline"
+## Install gcc for pip iced install
+RUN apt-get update && apt-get install -y gcc g++ && apt-get clean -y
+
COPY environment.yml /
RUN conda env create -f /environment.yml && conda clean -a
ENV PATH /opt/conda/envs/nf-core-hic-1.0dev/bin:$PATH
diff --git a/README.md b/README.md
index 192fa0de6b9cfd727673b531668abdfdb38c3e2d..cf702bf337b100b48ab98405e67cb3b1e7c855d9 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,5 @@
-# nf-core/hic
+# 
+
**Analysis of Chromosome Conformation Capture data (Hi-C)**
[](https://travis-ci.org/nf-core/hic)
@@ -10,8 +11,21 @@
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 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.
+### Pipeline summary
+1. Mapping using a two steps strategy to rescue reads spanning the ligation sites (bowtie2)
+2. Detection of valid interaction products
+3. Duplicates removal
+4. Create genome-wide contact maps at various resolution
+5. Contact maps normalization using the ICE algorithm (iced)
+6. Quality controls and report (MultiQC)
+7. Addition export for visualisation and downstream analysis (cooler)
### Documentation
The nf-core/hic pipeline comes with documentation about the pipeline, found in the `docs/` directory:
@@ -25,7 +39,5 @@ The nf-core/hic pipeline comes with documentation about the pipeline, found in t
4. [Output and how to interpret the results](docs/output.md)
5. [Troubleshooting](docs/troubleshooting.md)
-<!-- TODO nf-core: Add a brief overview of what the pipeline does and how it works -->
-
### Credits
nf-core/hic was originally written by Nicolas Servant.
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/base.config b/conf/base.config
index 98c386132177729a50a6937e716ae11533476491..7f99f28907f6c88829f0cc094de6672a4bacd873 100644
--- a/conf/base.config
+++ b/conf/base.config
@@ -11,7 +11,7 @@
process {
- container = params.container
+ container = process.container
cpus = { check_max( 2, 'cpus' ) }
memory = { check_max( 8.GB * task.attempt, 'memory' ) }
@@ -22,20 +22,18 @@ process {
maxErrors = '-1'
// Process-specific resource requirements
-
withName:makeBowtie2Index {
cpus = { check_max( 1, 'cpus' ) }
memory = { check_max( 10.GB * task.attempt, 'memory' ) }
time = { check_max( 12.h * task.attempt, 'time' ) }
}
-
withName:bowtie2_end_to_end {
- cpus = { check_max( 2, 'cpus' ) }
+ cpus = { check_max( 4, 'cpus' ) }
memory = { check_max( 4.GB * task.attempt, 'memory' ) }
time = { check_max( 5.h * task.attempt, 'time' ) }
}
withName:bowtie2_on_trimmed_reads {
- cpus = { check_max( 2, 'cpus' ) }
+ cpus = { check_max( 4, 'cpus' ) }
memory = { check_max( 4.GB * task.attempt, 'memory' ) }
time = { check_max( 5.h * task.attempt, 'time' ) }
}
@@ -59,12 +57,12 @@ process {
memory = { check_max( 4.GB * task.attempt, 'memory' ) }
time = { check_max( 5.h * task.attempt, 'time' ) }
}
- withName:build_contact_maps {
+ withName:build_contact_maps {
cpus = { check_max( 1, 'cpus' ) }
memory = { check_max( 6.GB * task.attempt, 'memory' ) }
time = { check_max( 5.h * task.attempt, 'time' ) }
}
-withName:run_ice {
+ withName:run_ice {
cpus = { check_max( 1, 'cpus' ) }
memory = { check_max( 10.GB * task.attempt, 'memory' ) }
time = { check_max( 5.h * task.attempt, 'time' ) }
diff --git a/conf/hicpro.config b/conf/hicpro.config
index d969960bcc10e78ace7187f521863dbbe8ea62b8..0a2c9b9e0db09f4f9861ba353b84a534820aba38 100644
--- a/conf/hicpro.config
+++ b/conf/hicpro.config
@@ -13,15 +13,15 @@ params {
splitFastq = false
bwt2_opts_end2end = '--very-sensitive -L 30 --score-min L,-0.6,-0.2 --end-to-end --reorder'
bwt2_opts_trimmed = '--very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder'
- min_mapq = 0
+ min_mapq = 10
// Digestion Hi-C
restriction_site = 'A^AGCTT'
ligation_site = 'AAGCTAGCTT'
- min_restriction_fragment_size = 100
- max_restriction_fragment_size = 100000
- min_insert_size = 100
- max_insert_size = 600
+ min_restriction_fragment_size =
+ max_restriction_fragment_size =
+ min_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/conf/test.config b/conf/test.config
index 1dd09920231bdd5da123402760d5c858d5d78b92..b4fd1845c65349aa6a58a82dc033b38d6bf76815 100644
--- a/conf/test.config
+++ b/conf/test.config
@@ -9,20 +9,26 @@
params {
- config_profile_name = 'Hi-C test data from Dixon et al. (2012)'
+ config_profile_name = 'Hi-C test data from Schalbetter et al. (2017)'
config_profile_description = 'Minimal test dataset to check pipeline function'
// Limit resources so that this can run on Travis
- max_cpus = 4
- max_memory = 6.GB
+ max_cpus = 2
+ max_memory = 4.GB
max_time = 1.h
// Input data
readPaths = [
- ['SRR400264_00', ['https://github.com/nf-core/test-datasets/raw/hic/SRR400264_00_R1.fastq.gz', 'https://github.com/nf-core/test-datasets/raw/hic/SRR400264_00_R2.fastq.gz']],
- ['SRR400264_01', ['https://github.com/nf-core/test-datasets/raw/hic/SRR400264_01_R1.fastq.gz', 'https://github.com/nf-core/test-datasets/raw/hic/SRR400264_01_R2.fastq.gz']]
- ]
+ ['SRR4292758_00', ['https://github.com/nf-core/test-datasets/raw/hic/data/SRR4292758_00_R1.fastq.gz', 'https://github.com/nf-core/test-datasets/raw/hic/data/SRR4292758_00_R2.fastq.gz']]
+ ]
// Annotations
- genome = 'GRCh37'
+ fasta = 'https://github.com/nf-core/test-datasets/raw/hic/reference/W303_SGD_2015_JRIU00000000.fsa'
+ restriction_site = 'A^AGCTT'
+ ligation_site = 'AAGCTAGCTT'
+ min_mapq = 0
+
+ // Options
+ skip_cool = true
}
+
diff --git a/docs/configuration/local.md b/docs/configuration/local.md
index 657422ffdaed246d92580b427e7e56cdabbf5cab..9cd485e2cd60b7670e242f4a399fde8e15bdaf92 100644
--- a/docs/configuration/local.md
+++ b/docs/configuration/local.md
@@ -11,7 +11,7 @@ First, install docker on your system: [Docker Installation Instructions](https:/
Then, simply run the analysis pipeline:
```bash
-nextflow run nf-core/hic -profile docker --genome '<genome ID>' --design '<path to your design file>'
+nextflow run nf-core/hic -profile docker --genome '<genome ID>'
```
Nextflow will recognise `nf-core/hic` and download the pipeline from GitHub. The `-profile docker` configuration lists the [nf-core/hic](https://hub.docker.com/r/nfcore/hic/) image that we have created and is hosted at dockerhub, and this is downloaded.
diff --git a/docs/images/nfcore-hic_logo.png b/docs/images/nfcore-hic_logo.png
new file mode 100644
index 0000000000000000000000000000000000000000..713ca9c500139f90bc9f3ad8e09d8c59b1f634b4
Binary files /dev/null and b/docs/images/nfcore-hic_logo.png differ
diff --git a/docs/output.md b/docs/output.md
index 9de7067b987fe10ef7435c135bbd58545f13ecd2..518ac60f545d4f87051dfec5ace6f972cd93d65b 100644
--- a/docs/output.md
+++ b/docs/output.md
@@ -2,31 +2,117 @@
This document describes the output produced by the pipeline. Most of the plots are taken from the MultiQC report, which summarises results at the end of the pipeline.
-<!-- TODO nf-core: Write this documentation describing your workflow's output -->
-
## Pipeline overview
The pipeline is built using [Nextflow](https://www.nextflow.io/)
and processes data using the following steps:
-* [FastQC](#fastqc) - read quality control
-* [MultiQC](#multiqc) - aggregate report, describing results of the whole pipeline
+* [Reads alignment](#reads-alignment)
+* [Valid pairs detection](#valid-pairs-detection)
+* [Duplicates removal](#duplicates-removal)
+* [Contact maps](#contact-maps)
+* [MultiQC](#multiqc) - aggregate report and quality controls, describing results of the whole pipeline
+* [Export](#exprot) - additionnal export for compatibility with downstream analysis tool and visualization
+
+The current version is mainly based on the [HiC-Pro](https://github.com/nservant/HiC-Pro) pipeline.
+For details about the workflow, see [Servant et al. 2015](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0831-x)
+
+## Reads alignment
+
+Using Hi-C data, each reads mate has to be independantly aligned on the reference genome.
+The current workflow implements a two steps mapping strategy. First, the reads are aligned using an end-to-end aligner.
+Second, reads spanning the ligation junction are trimmmed from their 3' end, and aligned back on the genome.
+Aligned reads for both fragment mates are then paired in a single paired-end BAM file.
+Singletons are discarded, and multi-hits are filtered according to the configuration parameters (`--rm-multi`).
+Note that if the `--dnase` mode is activated, HiC-Pro will skip the second mapping step.
+
+**Output directory: `results/mapping`**
+
+* `*bwt2pairs.bam` - final BAM file with aligned paired data
+* `*.pairstat` - mapping statistics
+
+if `--saveAlignedIntermediates` is specified, additional mapping file results are available ;
+
+* `*.bam` - Aligned reads (R1 and R2) from end-to-end alignment
+* `*_unmap.fastq` - Unmapped reads after end-to-end alignment
+* `*_trimmed.fastq` - Trimmed reads after end-to-end alignment
+* `*_trimmed.bam` - Alignment of trimmed reads
+* `*bwt2merged.bam` - merged BAM file after the two-steps alignment
+* `*.mapstat` - mapping statistics per read mate
+
+Usually, a high fraction of reads is expected to be aligned on the genome (80-90%). Among them, we usually observed a few percent (around 10%) of step 2 aligned reads. Those reads are chimeric fragments for which we detect a ligation junction. An abnormal level of chimeric reads can reflect a ligation issue during the library preparation.
+The fraction of singleton or multi-hits depends on the genome complexity and the fraction of unmapped reads. The fraction of singleton is usually close to the sum of unmapped R1 and R2 reads, as it is unlikely that both mates from the same pair were unmapped.
+
+## Valid pairs detection
+
+Each aligned reads can be assigned to one restriction fragment according to the reference genome and the digestion protocol.
+
+Invalid pairs are classified as follow:
+* Dangling end, i.e. unligated fragments (both reads mapped on the same restriction fragment)
+* Self circles, i.e. fragments ligated on themselves (both reads mapped on the same restriction fragment in inverted orientation)
+* Religation, i.e. ligation of juxtaposed fragments
+* Filtered pairs, i.e. any pairs that do not match the filtering criteria on inserts size, restriction fragments size
+* Dumped pairs, i.e. any pairs for which we were not able to reconstruct the ligation product.
+
+Only valid pairs involving two different restriction fragments are used to build the contact maps.
+Duplicated valid pairs associated to PCR artefacts are discarded (see `--rm_dup`.
-## FastQC
-[FastQC](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) gives general quality metrics about your reads. It provides information about the quality score distribution across your reads, the per base sequence content (%T/A/G/C). You get information about adapter contamination and other overrepresented sequences.
+In case of Hi-C protocols that do not require a restriction enzyme such as DNase Hi-C or micro Hi-C, the assignment to a restriction is not possible (see `--dnase`).
+Short range interactions that are likely to be spurious ligation products can thus be discarded using the `--min_cis_dist` parameter.
-For further reading and documentation see the [FastQC help](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/).
+* `*.validPairs` - List of valid ligation products
+* `*RSstat` - Statitics of number of read pairs falling in each category
-> **NB:** The FastQC plots displayed in the MultiQC report shows _untrimmed_ reads. They may contain adapter sequence and potentially regions with low quality. To see how your reads look after trimming, look at the FastQC reports in the `trim_galore` directory.
+The validPairs are stored using a simple tab-delimited text format ;
-**Output directory: `results/fastqc`**
+```
+read name / chr_reads1 / pos_reads1 / strand_reads1 / chr_reads2 / pos_reads2 / strand_reads2 / fragment_size / res frag name R1 / res frag R2 / mapping qual R1 / mapping qual R2 [/ allele_specific_tag]
+```
-* `sample_fastqc.html`
- * FastQC report, containing quality metrics for your untrimmed raw fastq files
-* `zips/sample_fastqc.zip`
- * zip file containing the FastQC report, tab-delimited data file and plot images
+The ligation efficiency can be assessed using the filtering of valid and invalid pairs. As the ligation is a random process, 25% of each valid ligation class is expected. In the same way, a high level of dangling-end or self-circle read pairs is associated with a low quality experiment, and reveals a problem during the digestion, fill-in or ligation steps.
+In the context of Hi-C protocol without restriction enzyme, this analysis step is skipped. The aligned pairs are therefore directly used to generate the contact maps. A filter of the short range contact (typically <1kb) is recommanded as this pairs are likely to be self ligation products.
+
+## Duplicates removal
+
+Note that validPairs file are generated per reads chunck.
+These files are then merged in the allValidPairs file, and duplicates are removed if the `--rm_dup` parameter is used.
+
+* `*allValidPairs` - combined valid pairs from all read chunks
+* `*mergestat` - statistics about duplicates removal and valid pairs information
+
+Additional quality controls such as fragment size distribution can be extracted from the list of valid interaction products.
+We usually expect to see a distribution centered around 300 pb which correspond to the paired-end insert size commonly used.
+The fraction of dplicates is also presented. A high level of duplication indicates a poor molecular complexity and a potential PCR bias.
+Finaly, an important metric is to look at the fraction of intra and inter-chromosomal interactions, as well as long range (>20kb) versus short range (<20kb) intra-chromosomal interactions.
+
+## Contact maps
+
+Intra et inter-chromosomal contact maps are build for all specified resolutions.
+The genome is splitted into bins of equal size. Each valid interaction is associated with the genomic bins to generate the raw maps.
+In addition, Hi-C data can contain several sources of biases which has to be corrected.
+The current workflow uses the [ìced](https://github.com/hiclib/iced) and [Varoquaux and Servant, 2018](http://joss.theoj.org/papers/10.21105/joss.01286) python package which proposes a fast implementation of the original ICE normalization algorithm (Imakaev et al. 2012), making the assumption of equal visibility of each fragment.
+
+* `*.matrix` - genome-wide contact maps
+* `*_iced.matrix` - genome-wide iced contact maps
+
+The contact maps are generated for all specified resolution (see `--bin_size` argument)
+A contact map is defined by :
+* A list of genomic intervals related to the specified resolution (BED format).
+* A matrix, stored as standard triplet sparse format (i.e. list format).
+
+Based on the observation that a contact map is symmetric and usually sparse, only non-zero values are stored for half of the matrix. The user can specified if the 'upper', 'lower' or 'complete' matrix has to be stored. The 'asis' option allows to store the contacts as they are observed from the valid pairs files.
+
+```
+ A B 10
+ A C 23
+ B C 24
+ (...)
+```
+
+This format is memory efficient, and is compatible with several software for downstream analysis.
## MultiQC
+
[MultiQC](http://multiqc.info) is a visualisation tool that generates a single HTML report summarising all samples in your project. Most of the pipeline QC results are visualised in the report and further statistics are available in within the report data directory.
The pipeline has special steps which allow the software versions used to be reported in the MultiQC output for future traceability.
diff --git a/docs/troubleshooting.md b/docs/troubleshooting.md
index 927d1b712a15d979f05fc3d74f1c16ddfd9d457c..e6772eb34bd66f12b8477547a1c3cc250d34f33d 100644
--- a/docs/troubleshooting.md
+++ b/docs/troubleshooting.md
@@ -1,7 +1,5 @@
# nf-core/hic: Troubleshooting
-<!-- TODO nf-core: Change this documentation if these parameters/errors are not relevant for your workflow -->
-
## Input files not found
If only no file, only one input file , or only read one and not read two is picked up then something is wrong with your input file declaration
diff --git a/docs/usage.md b/docs/usage.md
index 95d863f688afbf78f738249aade28bd711deb60b..853c38414b6e53090c3d9b0f19e849a0972b1243 100644
--- a/docs/usage.md
+++ b/docs/usage.md
@@ -19,6 +19,38 @@
* [`--genome`](#--genome)
* [`--fasta`](#--fasta)
* [`--igenomesIgnore`](#--igenomesignore)
+ * [`--bwt2_index`](#--bwt2_index)
+ * [`--chromosome_size`](#--chromosome_size)
+ * [`--restriction_fragments`](#--restriction_fragments)
+* [Hi-C specific options](#hi-c-specific-options)
+ * [Reads mapping](#reads-mapping)
+ * [`--bwt2_opts_end2end`](#--bwt2_opts_end2end)
+ * [`--bwt2_opts_trimmed`](#--bwt2_opts_trimmed)
+ * [`--min_mapq`](#--min_mapq)
+ * [Digestion Hi-C](#digestion-hi-c)
+ * [`--restriction_site`](#--restriction_site)
+ * [`--ligation_site`](#--ligation_site)
+ * [`--min_restriction_fragment_size`](#--min_restriction_fragment_size)
+ * [`--max_restriction_fragment_size`](#--max_restriction_fragment_size)
+ * [`--min_insert_size`](#--min_insert_size)
+ * [`--max_insert_size`](#--max_insert_size)
+ * [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)
+ * [`--rm_multi`](#--rm_multi)
+ * [Genome-wide contact maps](#genome-wide-contact-maps)
+ * [`--bins_size`](#--bins_size)
+ * [`--ice_max_iter`](#--ice_max_iter)
+ * [`--ice_filer_low_count_perc`](#--ice_filer_low_count_perc)
+ * [`--ice_filer_high_count_perc`](#--ice_filer_high_count_perc)
+ * [`--ice_eps`](#--ice_eps)
+ * [Inputs/Outputs](#inputs-outputs)
+ * [`--splitFastq`](#--splitFastq)
+ * [`--saveReference`](#--saveReference)
+ * [`--saveAlignedIntermediates`](#--saveAlignedIntermediates)
* [Job resources](#job-resources)
* [Automatic resubmission](#automatic-resubmission)
* [Custom resource requests](#custom-resource-requests)
@@ -48,11 +80,11 @@ It is recommended to limit the Nextflow Java virtual machines memory. We recomme
```bash
NXF_OPTS='-Xms1g -Xmx4g'
```
-<!-- TODO nf-core: Document required command line parameters to run the pipeline-->
+
## Running the pipeline
The typical command for running the pipeline is as follows:
```bash
-nextflow run nf-core/hic --reads '*_R{1,2}.fastq.gz' -profile docker
+nextflow run nf-core/hic --reads '*_R{1,2}.fastq.gz' -genome GRCh37 -profile docker
```
This will launch the pipeline with the `docker` configuration profile. See below for more information about profiles.
@@ -119,17 +151,7 @@ Please note the following requirements:
If left unspecified, a default pattern is used: `data/*{1,2}.fastq.gz`
-### `--singleEnd`
-By default, the pipeline expects paired-end data. If you have single-end data, you need to specify `--singleEnd` on the command line when you launch the pipeline. A normal glob pattern, enclosed in quotation marks, can then be used for `--reads`. For example:
-
-```bash
---singleEnd --reads '*.fastq'
-```
-
-It is not possible to run a mixture of single-end and paired-end files in one run.
-
-
-## 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.
@@ -153,19 +175,19 @@ Note that you can use the same configuration setup to save sets of reference fil
The syntax for this reference configuration is as follows:
-<!-- TODO nf-core: Update reference genome example according to what is needed -->
+
```nextflow
params {
genomes {
'GRCh37' {
- fasta = '<path to the genome fasta file>' // Used if no star index given
+ fasta = '<path to the genome fasta file>' // Used if no annotations are given
+ bowtie2 = '<path to bowtie2 index files>'
}
// Any number of additional genomes, key is used with --genome
}
}
```
-<!-- TODO nf-core: Describe reference path flags -->
### `--fasta`
If you prefer, you can specify the full path to your reference genome when you run the pipeline:
@@ -176,6 +198,264 @@ If you prefer, you can specify the full path to your reference genome when you r
### `--igenomesIgnore`
Do not load `igenomes.config` when running the pipeline. You may choose this option if you observe clashes between custom parameters and those supplied in `igenomes.config`.
+### `--bwt2_index`
+
+The bowtie2 indexes are required to run the Hi-C pipeline. If the `--bwt2_index` is not specified, the pipeline will either use the igenome bowtie2 indexes (see `--genome` option) or build the indexes on-the-fly (see `--fasta` option)
+
+```bash
+--bwt2_index '[path to bowtie2 index (with basename)]'
+```
+
+### `--chromosome_size`
+
+The Hi-C pipeline will also requires a two-columns text file with the chromosome name and its size (tab separated).
+If not specified, this file will be automatically created by the pipeline. In the latter case, the `--fasta` reference genome has to be specified.
+```
+ chr1 249250621
+ chr2 243199373
+ chr3 198022430
+ chr4 191154276
+ chr5 180915260
+ chr6 171115067
+ chr7 159138663
+ chr8 146364022
+ chr9 141213431
+ chr10 135534747
+ (...)
+```
+
+```bash
+--chromosome_size '[path to chromosome size file]'
+```
+
+### `--restriction_fragments`
+
+Finally, Hi-C experiments based on restriction enzyme digestion requires a BED file with coordinates of restriction fragments.
+
+```
+ chr1 0 16007 HIC_chr1_1 0 +
+ chr1 16007 24571 HIC_chr1_2 0 +
+ chr1 24571 27981 HIC_chr1_3 0 +
+ chr1 27981 30429 HIC_chr1_4 0 +
+ chr1 30429 32153 HIC_chr1_5 0 +
+ chr1 32153 32774 HIC_chr1_6 0 +
+ chr1 32774 37752 HIC_chr1_7 0 +
+ chr1 37752 38369 HIC_chr1_8 0 +
+ chr1 38369 38791 HIC_chr1_9 0 +
+ chr1 38791 39255 HIC_chr1_10 0 +
+ (...)
+```
+
+If not specified, this file will be automatically created by the pipline. In this case, the `--fasta` reference genome will be used.
+Note that the `--restriction_site` parameter is mandatory to create this file.
+
+## Hi-C specific options
+
+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
+
+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.
+Note that the default option are quite stringent, and can be updated according to the reads quality or the reference genome.
+
+#### `--bwt2_opts_end2end`
+
+Bowtie2 alignment option for end-to-end mapping. Default: '--very-sensitive -L 30 --score-min L,-0.6,-0.2 --end-to-end --reorder'
+
+```bash
+--bwt2_opts_end2end '[Options for bowtie2 step1 mapping on full reads]'
+```
+
+#### `--bwt2_opts_trimmed`
+
+Bowtie2 alignment option for trimmed reads mapping (step 2). Default: '--very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder'
+
+```bash
+--bwt2_opts_trimmed '[Options for bowtie2 step2 mapping on trimmed reads]'
+```
+
+#### `--min_mapq`
+
+Minimum mapping quality. Reads with lower quality are discarded. Default: 10
+
+```bash
+--min_mapq '[Minimum quality value]'
+```
+
+### Digestion Hi-C
+
+#### `--restriction_site`
+
+Restriction motif(s) for Hi-C digestion protocol. The restriction motif(s) is(are) used to generate the list of restriction fragments.
+The precise cutting site of the restriction enzyme has to be specified using the '^' character. Default: 'A^AGCTT'
+Here are a few examples:
+* MboI: '^GATC'
+* DpnII: '^GATC'
+* BglII: 'A^GATCT'
+* HindIII: 'A^AGCTT'
+
+Note that multiples restriction motifs can be provided (comma-separated).
+
+```bash
+--restriction_size '[Cutting motif]'
+```
+
+#### `--ligation_site`
+
+Ligation motif after reads ligation. This motif is used for reads trimming and depends on the fill in strategy.
+Note that multiple ligation sites can be specified. Default: 'AAGCTAGCTT'
+
+```bash
+--ligation_site '[Ligation motif]'
+```
+
+#### `--min_restriction_fragment_size`
+
+Minimum size of restriction fragments to consider for the Hi-C processing. Default: ''
+
+```bash
+--min_restriction_fragment_size '[numeric]'
+```
+
+#### `--max_restriction_fragment_size`
+
+Maximum size of restriction fragments to consider for the Hi-C processing. Default: ''
+
+```bash
+--max_restriction_fragment_size '[numeric]'
+```
+
+#### `--min_insert_size`
+
+Minimum reads insert size. Shorter 3C products are discarded. Default: ''
+
+```bash
+--min_insert_size '[numeric]'
+```
+
+#### `--max_insert_size`
+
+Maximum reads insert size. Longer 3C products are discarded. Default: ''
+
+```bash
+--max_insert_size '[numeric]'
+```
+
+### 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`
+
+Filter short range contact below the specified distance. Mainly useful for DNase Hi-C. Default: ''
+
+```bash
+--min_cis_dist '[numeric]'
+```
+
+#### `--rm_singleton`
+
+If specified, singleton reads are discarded at the mapping step.
+
+```bash
+--rm_singleton
+```
+
+#### `--rm_dup`
+
+If specified, duplicates reads are discarded before building contact maps.
+
+```bash
+--rm_dup
+```
+
+#### `--rm_multi`
+
+If specified, reads that aligned multiple times on the genome are discarded. Note the default mapping options are based on random hit assignment, meaning that only one position is kept per read.
+
+```bash
+--rm_multi
+```
+
+## Genome-wide contact maps
+
+#### `--bin_size`
+
+Resolution of contact maps to generate (space separated). Default:'1000000,500000'
+
+```bash
+--bins_size '[numeric]'
+```
+
+#### `--ice_max_iter`
+
+Maximum number of iteration for ICE normalization. Default: 100
+
+```bash
+--ice_max_iter '[numeric]'
+```
+
+#### `--ice_filer_low_count_perc`
+
+Define which pourcentage of bins with low counts should be force to zero. Default: 0.02
+
+```bash
+--ice_filter_low_count_perc '[numeric]'
+```
+
+#### `--ice_filer_high_count_perc`
+
+Define which pourcentage of bins with low counts should be discarded before normalization. Default: 0
+
+```bash
+--ice_filter_high_count_perc '[numeric]'
+```
+
+#### `--ice_eps`
+
+The relative increment in the results before declaring convergence for ICE normalization. Default: 0.1
+
+```bash
+--ice_eps '[numeric]'
+```
+
+## Inputs/Outputs
+
+#### `--splitFastq`
+
+By default, the nf-core Hi-C pipeline expects one read pairs per sample. However, for large Hi-C data processing single fastq files can be very time consuming.
+The `--splitFastq` option allows to automatically split input read pairs into chunks of reads. In this case, all chunks will be processed in parallel and merged before generating the contact maps, thus leading to a significant increase of processing performance.
+
+```bash
+--splitFastq '[Number of reads per chunk]'
+```
+
+#### `--saveReference`
+
+If specified, annotation files automatically generated from the `--fasta` file are exported in the results folder. Default: false
+
+```
+--saveReference
+```
+
+#### `--saveAlignedIntermediates`
+
+If specified, all intermediate mapping files are saved and exported in the results folder. Default: false
+
+```
+--saveReference
+```
+
## Job resources
### Automatic resubmission
Each step in the pipeline has a default set of requirements for number of CPUs, memory and time. For most of the steps in the pipeline, if the job exits with an error code of `143` (exceeded requested resources) it will automatically resubmit with higher requests (2 x original, then 3 x original). If it still fails after three times then the pipeline is stopped.
@@ -198,8 +478,6 @@ Please make sure to also set the `-w/--work-dir` and `--outdir` parameters to a
## Other command line parameters
-<!-- TODO nf-core: Describe any other command line flags here -->
-
### `--outdir`
The output directory where the results will be saved.
diff --git a/environment.yml b/environment.yml
index b5cee2f9016661f493d9e1fc04f9a62d3ee31516..afddba2ef81dd47d399b156c582c4641a6e28dc3 100644
--- a/environment.yml
+++ b/environment.yml
@@ -4,12 +4,16 @@ channels:
- bioconda
- defaults
dependencies:
- - python=2.7.11
+ - python=2.7.13
+ - pip
- conda-forge::scipy=1.0.1
- conda-forge::numpy=1.9.3
+ - conda-forge::r-markdown=0.8
- bcbio::bx-python=0.7.3
- bioconda::pysam=0.14.1
- cooler=0.8.3
- bowtie2=2.3.5
- samtools=1.9
- multiqc=1.6
+ - pip:
+ - iced==0.4.2
\ No newline at end of file
diff --git a/main.nf b/main.nf
index 53798f4aeca7cc2266f93189294c5b4719566dae..5550fb1d156943fb8cdf46f6ae09a64b347a463a 100644
--- a/main.nf
+++ b/main.nf
@@ -33,11 +33,11 @@ def helpMessage() {
Mandatory arguments:
--reads Path to input data (must be surrounded with quotes)
- --genome Name of iGenomes reference
-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
--chromosome_size Path to chromosome size file
@@ -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
@@ -72,6 +74,10 @@ def helpMessage() {
--email Set this parameter to your e-mail address to get a summary e-mail with details of the run sent to you when the workflow exits
-name Name for the pipeline run. If not specified, Nextflow will automatically generate a random mnemonic.
+ Step options:
+ --skip_cool Skip generation of cool files
+ --skip_multiQC Skip MultiQC
+
AWSBatch options:
--awsqueue The AWSBatch JobQueue that needs to be set when running on AWSBatch
--awsregion The AWS Region for your AWS Batch job to run on
@@ -93,6 +99,11 @@ if (params.genomes && params.genome && !params.genomes.containsKey(params.genome
exit 1, "The provided genome '${params.genome}' is not available in the iGenomes file. Currently the available genomes are ${params.genomes.keySet().join(", ")}"
}
+// Check Digestion or DNase Hi-C mode
+if (!params.dnase && !params.ligation_site) {
+ exit 1, "Ligation motif not found. For DNase Hi-C, please use '--dnase' option"
+}
+
// Reference index path configuration
params.bwt2_index = params.genome ? params.genomes[ params.genome ].bowtie2 ?: false : false
params.fasta = params.genome ? params.genomes[ params.genome ].fasta ?: false : false
@@ -138,8 +149,7 @@ if (params.readPaths){
.from( params.readPaths )
.map { row -> [ row[0], [file(row[1][0]), file(row[1][1])]] }
.separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0], a[1][0]), tuple(a[0], a[1][1])] }
- .println()
-}else{
+ }else{
raw_reads = Channel.create()
raw_reads_2 = Channel.create()
@@ -218,7 +228,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)
@@ -239,31 +249,36 @@ log.info """=======================================================
nf-core/hic v${workflow.manifest.version}"
======================================================="""
def summary = [:]
-summary['Pipeline Name'] = 'nf-core/hic'
+summary['Pipeline Name'] = 'nf-core/hic'
summary['Pipeline Version'] = workflow.manifest.version
-summary['Run Name'] = custom_runName ?: workflow.runName
-
-summary['Reads'] = params.reads
-summary['Fasta Ref'] = params.fasta
-
-
-summary['Max Memory'] = params.max_memory
-summary['Max CPUs'] = params.max_cpus
-summary['Max Time'] = params.max_time
-summary['Output dir'] = params.outdir
-summary['Working dir'] = workflow.workDir
+summary['Run Name'] = custom_runName ?: workflow.runName
+
+summary['Reads'] = params.reads
+summary['splitFastq'] = params.splitFastq
+summary['Fasta Ref'] = params.fasta
+summary['Ligation Motif'] = params.ligation_site
+summary['DNase Mode'] = params.dnase
+summary['Remove Dup'] = params.rm_dup
+summary['Maps resolution'] = params.bin_size
+
+summary['Max Memory'] = params.max_memory
+summary['Max CPUs'] = params.max_cpus
+summary['Max Time'] = params.max_time
+summary['Output dir'] = params.outdir
+summary['Working dir'] = workflow.workDir
summary['Container Engine'] = workflow.containerEngine
-if(workflow.containerEngine) summary['Container'] = workflow.container
-summary['Current home'] = "$HOME"
-summary['Current user'] = "$USER"
-summary['Current path'] = "$PWD"
-summary['Working dir'] = workflow.workDir
-summary['Output dir'] = params.outdir
-summary['Script dir'] = workflow.projectDir
-summary['Config Profile'] = workflow.profile
+if(workflow.containerEngine)
+ summary['Container'] = workflow.container
+summary['Current home'] = "$HOME"
+summary['Current user'] = "$USER"
+summary['Current path'] = "$PWD"
+summary['Working dir'] = workflow.workDir
+summary['Output dir'] = params.outdir
+summary['Script dir'] = workflow.projectDir
+summary['Config Profile'] = workflow.profile
if(workflow.profile == 'awsbatch'){
- summary['AWS Region'] = params.awsregion
- summary['AWS Queue'] = params.awsqueue
+ summary['AWS Region'] = params.awsregion
+ summary['AWS Queue'] = params.awsqueue
}
if(params.email) summary['E-mail Address'] = params.email
log.info summary.collect { k,v -> "${k.padRight(15)}: $v" }.join("\n")
@@ -355,7 +370,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 },
@@ -398,14 +413,26 @@ process bowtie2_end_to_end {
script:
prefix = reads.toString() - ~/(\.fq)?(\.fastq)?(\.gz)?$/
def bwt2_opts = params.bwt2_opts_end2end
- """
- bowtie2 --rg-id BMG --rg SM:${prefix} \\
+
+ if (!params.dnase){
+ """
+ bowtie2 --rg-id BMG --rg SM:${prefix} \\
${bwt2_opts} \\
-p ${task.cpus} \\
-x ${index}/${bwt2_base} \\
--un ${prefix}_unmap.fastq \\
-U ${reads} | samtools view -F 4 -bS - > ${prefix}.bam
- """
+ """
+ }else{
+ """
+ bowtie2 --rg-id BMG --rg SM:${prefix} \\
+ ${bwt2_opts} \\
+ -p ${task.cpus} \\
+ -x ${index}/${bwt2_base} \\
+ --un ${prefix}_unmap.fastq \\
+ -U ${reads} > ${prefix}.bam
+ """
+ }
}
process trim_reads {
@@ -413,6 +440,9 @@ process trim_reads {
publishDir path: { params.saveAlignedIntermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.saveAlignedIntermediates ? it : null }, mode: 'copy'
+ when:
+ !params.dnase
+
input:
set val(prefix), file(reads) from unmapped_end_to_end
@@ -432,6 +462,9 @@ process bowtie2_on_trimmed_reads {
publishDir path: { params.saveAlignedIntermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.saveAlignedIntermediates ? it : null }, mode: 'copy'
+ when:
+ !params.dnase
+
input:
set val(prefix), file(reads) from trimmed_reads
file index from bwt2_index_trim.collect()
@@ -450,47 +483,80 @@ process bowtie2_on_trimmed_reads {
"""
}
-process merge_mapping_steps{
- tag "$sample = $bam1 + $bam2"
- publishDir path: { params.saveAlignedIntermediates ? "${params.outdir}/mapping" : params.outdir },
+if (!params.dnase){
+ process merge_mapping_steps{
+ tag "$sample = $bam1 + $bam2"
+ publishDir path: { params.saveAlignedIntermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.saveAlignedIntermediates ? it : null }, mode: 'copy'
- input:
- set val(prefix), file(bam1), file(bam2) from end_to_end_bam.join( trimmed_bam )
+ input:
+ set val(prefix), file(bam1), file(bam2) from end_to_end_bam.join( trimmed_bam )
- output:
- set val(sample), file("${prefix}_bwt2merged.bam") into bwt2_merged_bam
- set val(oname), file("${prefix}.mapstat") into all_mapstat
+ output:
+ set val(sample), file("${prefix}_bwt2merged.bam") into bwt2_merged_bam
+ set val(oname), file("${prefix}.mapstat") into all_mapstat
- script:
- sample = prefix.toString() - ~/(_R1|_R2|_val_1|_val_2)/
- tag = prefix.toString() =~/_R1|_val_1/ ? "R1" : "R2"
- oname = prefix.toString() - ~/(\.[0-9]+)$/
+ script:
+ sample = prefix.toString() - ~/(_R1|_R2|_val_1|_val_2)/
+ tag = prefix.toString() =~/_R1|_val_1/ ? "R1" : "R2"
+ oname = prefix.toString() - ~/(\.[0-9]+)$/
- """
- samtools merge -@ ${task.cpus} \\
+ """
+ samtools merge -@ ${task.cpus} \\
-f ${prefix}_bwt2merged.bam \\
${bam1} ${bam2}
- samtools sort -@ ${task.cpus} -m 800M \\
+ samtools sort -@ ${task.cpus} -m 800M \\
-n -T /tmp/ \\
-o ${prefix}_bwt2merged.sorted.bam \\
${prefix}_bwt2merged.bam
- mv ${prefix}_bwt2merged.sorted.bam ${prefix}_bwt2merged.bam
-
- echo "## ${prefix}" > ${prefix}.mapstat
- echo -n "total_${tag}\t" >> ${prefix}.mapstat
- samtools view -c ${prefix}_bwt2merged.bam >> ${prefix}.mapstat
- echo -n "mapped_${tag}\t" >> ${prefix}.mapstat
- samtools view -c -F 4 ${prefix}_bwt2merged.bam >> ${prefix}.mapstat
- echo -n "global_${tag}\t" >> ${prefix}.mapstat
- samtools view -c -F 4 ${bam1} >> ${prefix}.mapstat
- echo -n "local_${tag}\t" >> ${prefix}.mapstat
- samtools view -c -F 4 ${bam2} >> ${prefix}.mapstat
- """
+ mv ${prefix}_bwt2merged.sorted.bam ${prefix}_bwt2merged.bam
+
+ echo "## ${prefix}" > ${prefix}.mapstat
+ echo -n "total_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c ${prefix}_bwt2merged.bam >> ${prefix}.mapstat
+ echo -n "mapped_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c -F 4 ${prefix}_bwt2merged.bam >> ${prefix}.mapstat
+ echo -n "global_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c -F 4 ${bam1} >> ${prefix}.mapstat
+ echo -n "local_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c -F 4 ${bam2} >> ${prefix}.mapstat
+ """
+ }
+}else{
+ process dnase_mapping_stats{
+ tag "$sample = $bam1"
+ publishDir path: { params.saveAlignedIntermediates ? "${params.outdir}/mapping" : params.outdir },
+ saveAs: { params.saveAlignedIntermediates ? it : null }, mode: 'copy'
+
+ input:
+ set val(prefix), file(bam1) from end_to_end_bam
+
+ output:
+ set val(sample), file(bam1) into bwt2_merged_bam
+ set val(oname), file("${prefix}.mapstat") into all_mapstat
+
+ script:
+ sample = prefix.toString() - ~/(_R1|_R2|_val_1|_val_2)/
+ tag = prefix.toString() =~/_R1|_val_1/ ? "R1" : "R2"
+ oname = prefix.toString() - ~/(\.[0-9]+)$/
+
+ """
+ echo "## ${prefix}" > ${prefix}.mapstat
+ echo -n "total_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c ${bam1} >> ${prefix}.mapstat
+ echo -n "mapped_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c -F 4 ${bam1} >> ${prefix}.mapstat
+ echo -n "global_${tag}\t" >> ${prefix}.mapstat
+ samtools view -c -F 4 ${bam1} >> ${prefix}.mapstat
+ echo -n "local_${tag}\t0" >> ${prefix}.mapstat
+ """
+ }
}
+
+
process combine_mapped_files{
tag "$sample = $r1_prefix + $r2_prefix"
publishDir "${params.outdir}/mapping", mode: 'copy',
@@ -524,39 +590,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
*/
@@ -674,6 +767,9 @@ process generate_cool{
tag "$sample"
publishDir "${params.outdir}/export/cool", mode: 'copy'
+ when:
+ !params.skip_cool
+
input:
set val(sample), file(vpairs) from all_valid_pairs_4cool
file chrsize from chromosome_size_cool.collect()
@@ -694,15 +790,18 @@ process generate_cool{
process multiqc {
publishDir "${params.outdir}/MultiQC", mode: 'copy'
+ when:
+ !params.skip_multiqc
+
input:
- file multiqc_config from ch_multiqc_config
- file ('input_*/*') from all_mstats.concat(all_mergestat).collect()
- file ('software_versions/*') from software_versions_yaml
- file workflow_summary from create_workflow_summary(summary)
+ file multiqc_config from ch_multiqc_config
+ file ('input_*/*') from all_mstats.concat(all_mergestat).collect()
+ file ('software_versions/*') from software_versions_yaml
+ file workflow_summary from create_workflow_summary(summary)
output:
- file "*multiqc_report.html" into multiqc_report
- file "*_data"
+ file "*multiqc_report.html" into multiqc_report
+ file "*_data"
script:
rtitle = custom_runName ? "--title \"$custom_runName\"" : ''
diff --git a/nextflow.config b/nextflow.config
index 0075f7f7e82e408308915d5c82bb219bef2eff6f..32486aab168ccff7647220cd75a61f72f695fbc5 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -11,11 +11,6 @@
// Global default params, used in configs
params {
- // Container slug. Stable releases should specify release tag!
- // Developmental code should specify :latest
- //container = 'nfcore/hic:latest'
- container = 'nservant/nf-core-hic:latest'
-
// Workflow flags
// TODO nf-core: Specify your pipeline's command line flags
reads = "*{1,2}.fastq.gz"
@@ -24,7 +19,9 @@ params {
readPaths = false
chromosome_size = false
restriction_fragments = false
-
+ skip_cool = false
+ skip_multiqc = false
+ dnase = false
// Boilerplate options
name = false
@@ -41,6 +38,7 @@ params {
custom_config_version = 'master'
}
+process.container = 'nfcore/hic:dev'
// Load base.config by default for all pipelines
includeConfig 'conf/base.config'
@@ -56,14 +54,8 @@ profiles {
awsbatch { includeConfig 'conf/awsbatch.config' }
conda { process.conda = "$baseDir/environment.yml" }
debug { process.beforeScript = 'echo $HOSTNAME' }
- docker {
- docker.enabled = true
- process.container = params.container
- }
- singularity {
- singularity.enabled = true
- process.container = {"shub://${params.container.replace('nfcore', 'nf-core')}"}
- }
+ docker { docker.enabled = true }
+ singularity { singularity.enabled = true }
test { includeConfig 'conf/test.config' }
}