diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6be2eab8616b01e97266a50938accaf430f00609..677f8cc18bbac6047efdf249bb3a7b800dbbe5cd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,9 +5,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
## v1.3.0dev - 2020-11-01
+* New `--digestion` parameter to automatically set the restriction_site and ligation_site motifs
+* New `--keep_multi` and `keep_dup` options. Default: false
* Template update for nf-core/tools v1.11
* Minor fix to summary log messages in pipeline header
+### `Fixed`
+
+* Fix recurrent bug in input file extension (#86)
+* Fix bug in `--bin_size` parameter (#85)
+* `min_mapq` is ignored if `--keep_multi` is used
+
+### Deprecated
+
+* `--rm_dup` and `rm_multi` are replaced by `--keep_dup` and `--keep_multi`
+
## v1.2.2 - 2020-09-02
### `Added`
diff --git a/bin/scrape_software_versions.py b/bin/scrape_software_versions.py
index 9f5650db54daa74c59c3967500a399ff55534f29..f882f23325768249c6f6f0122f0d238bd9a7df37 100755
--- a/bin/scrape_software_versions.py
+++ b/bin/scrape_software_versions.py
@@ -7,9 +7,9 @@ import re
regexes = {
'nf-core/hic': ['v_pipeline.txt', r"(\S+)"],
'Nextflow': ['v_nextflow.txt', r"(\S+)"],
- 'Bowtie2': ['v_bowtie2.txt', r"Bowtie2 v(\S+)"],
- 'Python': ['v_python.txt', r"Python v(\S+)"],
- 'Samtools': ['v_samtools.txt', r"Samtools v(\S+)"],
+ 'Bowtie2': ['v_bowtie2.txt', r"bowtie2-align-s version (\S+)"],
+ 'Python': ['v_python.txt', r"Python (\S+)"],
+ 'Samtools': ['v_samtools.txt', r"samtools (\S+)"],
'MultiQC': ['v_multiqc.txt', r"multiqc, version (\S+)"],
}
results = OrderedDict()
diff --git a/conf/test.config b/conf/test.config
index 2ab8e57eda3d8fddb90c6d7ed3ddf2c0fd0672ca..5b3ac0e1395c89032ebfe04e8688e62d1abb660f 100644
--- a/conf/test.config
+++ b/conf/test.config
@@ -24,14 +24,8 @@ config_profile_name = 'Hi-C test data from Schalbetter et al. (2017)'
// Annotations
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 = 2
- rm_dup = true
- rm_singleton = true
- rm_multi = true
-
+ digestion = 'hindiii'
+ min_mapq = 20
min_restriction_fragment_size = 100
max_restriction_fragment_size = 100000
min_insert_size = 100
diff --git a/docs/usage.md b/docs/usage.md
index 31eabe9fb8044a14de7319c5824c18a53b06936e..be630a4298976e67bec49659f6dc6bc38ac86cc6 100644
--- a/docs/usage.md
+++ b/docs/usage.md
@@ -188,6 +188,30 @@ We recommend adding the following line to your environment to limit this
NXF_OPTS='-Xms1g -Xmx4g'
```
+## Use case
+
+### Hi-C digestion protocol
+
+Here is an command line example for standard DpnII digestion protocols.
+Alignment will be performed on the `mm10` genome with default paramters.
+Multi-hits will not be considered and duplicates will be removed.
+Note that by default, no filters are applied on DNA and restriction fragment sizes.
+
+```bash
+nextflow run main.nf --input './*_R{1,2}.fastq.gz' --genome 'mm10' --digestion 'dnpii'
+```
+
+### DNase Hi-C protocol
+
+Here is an command line example for DNase protocol.
+Alignment will be performed on the `mm10` genome with default paramters.
+Multi-hits will not be considered and duplicates will be removed.
+Contacts involving fragments separated by less than 1000bp will be discarded.
+
+```bash
+nextflow run main.nf --input './*_R{1,2}.fastq.gz' --genome 'mm10' --dnase --min_cis 1000
+```
+
## Inputs
### `--input`
@@ -207,16 +231,7 @@ notation to specify read pairs.
If left unspecified, a default pattern is used: `data/*{1,2}.fastq.gz`
-By default, the pipeline expects paired-end data. If you have single-end data,
-you need to specify `--single_end` on the command line when you launch the pipeline.
-A normal glob pattern, enclosed in quotation marks, can then be used for `--input`.
-For example:
-
-```bash
---single_end --input '*.fastq'
-```
-
-It is not possible to run a mixture of single-end and paired-end files in one run.
+Note that the Hi-C data analysis requires paired-end data.
## Reference genomes
@@ -302,7 +317,7 @@ 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
+The following options are defined in the `nextflow.config` file, and can be
updated either using a custom configuration file (see `-c` option) or using
command line parameter.
@@ -345,17 +360,28 @@ Minimum mapping quality. Reads with lower quality are discarded. Default: 10
### Digestion Hi-C
+#### `--digestion`
+
+This parameter allows to automatically set the `--restriction_site` and
+`--ligation_site` parameter according to the restriction enzyme you used.
+Available keywords are 'hindiii', 'dpnii', 'mboi', 'arima'.
+
+```bash
+--digestion 'hindiii'
+```
+
#### `--restriction_site`
-Restriction motif(s) for Hi-C digestion protocol. The restriction motif(s)
-is(are) used to generate the list of restriction fragments.
+If the restriction enzyme is not available through the `--digestion`
+parameter, you can also defined manually the 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
* ARIMA kit: ^GATC,G^ANTC
@@ -383,7 +409,7 @@ Exemple of the ARIMA kit: GATCGATC,GANTGATC,GANTANTC,GATCANTC
#### `--min_restriction_fragment_size`
Minimum size of restriction fragments to consider for the Hi-C processing.
-Default: ''
+Default: '0' - no filter
```bash
--min_restriction_fragment_size '[numeric]'
@@ -392,7 +418,7 @@ Default: ''
#### `--max_restriction_fragment_size`
Maximum size of restriction fragments to consider for the Hi-C processing.
-Default: ''
+Default: '0' - no filter
```bash
--max_restriction_fragment_size '[numeric]'
@@ -401,7 +427,7 @@ Default: ''
#### `--min_insert_size`
Minimum reads insert size. Shorter 3C products are discarded.
-Default: ''
+Default: '0' - no filter
```bash
--min_insert_size '[numeric]'
@@ -410,7 +436,7 @@ Default: ''
#### `--max_insert_size`
Maximum reads insert size. Longer 3C products are discarded.
-Default: ''
+Default: '0' - no filter
```bash
--max_insert_size '[numeric]'
@@ -434,36 +460,29 @@ to remove spurious ligation products.
#### `--min_cis_dist`
Filter short range contact below the specified distance.
-Mainly useful for DNase Hi-C. Default: ''
+Mainly useful for DNase Hi-C. Default: '0'
```bash
--min_cis_dist '[numeric]'
```
-#### `--rm_singleton`
-
-If specified, singleton reads are discarded at the mapping step.
-
-```bash
---rm_singleton
-```
-
-#### `--rm_dup`
+#### `--keep_dups`
-If specified, duplicates reads are discarded before building contact maps.
+If specified, duplicates reads are not discarded before building contact maps.
```bash
---rm_dup
+--keep_dups
```
-#### `--rm_multi`
+#### `--keep_multi`
-If specified, reads that aligned multiple times on the genome are discarded.
+If specified, reads that aligned multiple times on the genome are not discarded.
Note the default mapping options are based on random hit assignment, meaning
that only one position is kept per read.
+Note that in this case the `--min_mapq` parameter is ignored.
```bash
---rm_multi
+--keep_multi
```
## Genome-wide contact maps
diff --git a/main.nf b/main.nf
index 53976cae18ae30e055c3c75ae1a041269bfe7117..46b3d9b9920442e0cd23a7fd8717adddddda06c5 100644
--- a/main.nf
+++ b/main.nf
@@ -22,48 +22,50 @@ def helpMessage() {
Mandatory arguments:
--input [file] Path to input data (must be surrounded with quotes)
+ --genome [str] Name of iGenomes reference
-profile [str] 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.
- --genome [str] Name of iGenomes reference
--bwt2_index [file] Path to Bowtie2 index
--fasta [file] Path to Fasta reference
+
+ Digestion Hi-C If not specified in the configuration file or you wish to set up specific digestion protocol
+ --digestion [str] Digestion Hi-C. Name of restriction enzyme used for digestion pre-configuration. Default: 'hindiii'
+ --ligation_site [str] Ligation motifs to trim (comma separated) if not available in --digestion. Default: false
+ --restriction_site [str] Cutting motif(s) of restriction enzyme(s) (comma separated) if not available in --digestion. Default: false
--chromosome_size [file] Path to chromosome size file
--restriction_fragments [file] Path to restriction fragment file (bed)
--save_reference [bool] Save reference genome to output folder. Default: False
+ DNase Hi-C
+ --dnase [bool] Run DNase Hi-C mode. All options related to restriction fragments are not considered. Default: False
+ --min_cis_dist [int] Minimum intra-chromosomal distance to consider. Default: None
+
Alignments
- --split_fastq [bool] Split fastq files in reads chunks to speed up computation. Default: false
- --fastq_chunks_size [int] Size of read chunks if split_fastq is true. Default: 20000000
- --save_aligned_intermediates [bool] Save intermediates alignment files. Default: False
--bwt2_opts_end2end [str] Options for bowtie2 end-to-end mappinf (first mapping step). See hic.config for default.
--bwt2_opts_trimmed [str] Options for bowtie2 mapping after ligation site trimming. See hic.config for default.
--min_mapq [int] Minimum mapping quality values to consider. Default: 10
- --restriction_site [str] Cutting motif(s) of restriction enzyme(s) (comma separated). Default: 'A^AGCTT'
- --ligation_site [str] Ligation motifs to trim (comma separated). Default: 'AAGCTAGCTT'
- --rm_singleton [bool] Remove singleton reads. Default: true
- --rm_multi [bool] Remove multi-mapped reads. Default: true
- --rm_dup [bool] Remove duplicates. Default: true
-
- Contacts calling
- --min_restriction_fragment_size [int] Minimum size of restriction fragments to consider. Default: 0
- --max_restriction_fragment_size [int] Maximum size of restriction fragments to consider. Default: 0
- --min_insert_size [int] Minimum insert size of mapped reads to consider. Default: 0
- --max_insert_size [int] Maximum insert size of mapped reads to consider. Default: 0
- --save_interaction_bam [bool] Save BAM file with interaction tags (dangling-end, self-circle, etc.). Default: False
+ --keep_multi [bool] Keep multi-mapped reads (--min_mapq is ignored). Default: false
+ --keep_dups [bool] Keep duplicates. Default: false
+ --save_aligned_intermediates [bool] Save intermediates alignment files. Default: False
+ --split_fastq [bool] Split fastq files in reads chunks to speed up computation. Default: false
+ --fastq_chunks_size [int] Size of read chunks if split_fastq is true. Default: 20000000
- --dnase [bool] Run DNase Hi-C mode. All options related to restriction fragments are not considered. Default: False
- --min_cis_dist [int] Minimum intra-chromosomal distance to consider. Default: 0
+ Valid Pairs Detection
+ --min_restriction_fragment_size [int] Minimum size of restriction fragments to consider. Default: None
+ --max_restriction_fragment_size [int] Maximum size of restriction fragments to consider. Default: None
+ --min_insert_size [int] Minimum insert size of mapped reads to consider. Default: None
+ --max_insert_size [int] Maximum insert size of mapped reads to consider. Default: None
+ --save_interaction_bam [bool] Save BAM file with interaction tags (dangling-end, self-circle, etc.). Default: False
Contact maps
- --bin_size [int] Bin size for contact maps (comma separated). Default: '1000000,500000'
+ --bin_size [str] Bin size for contact maps (comma separated). Default: '1000000,500000'
--ice_max_iter [int] Maximum number of iteration for ICE normalization. Default: 100
--ice_filter_low_count_perc [float] Percentage of low counts columns/rows to filter before ICE normalization. Default: 0.02
--ice_filter_high_count_perc [float] Percentage of high counts columns/rows to filter before ICE normalization. Default: 0
--ice_eps [float] Convergence criteria for ICE normalization. Default: 0.1
-
Workflow
--skip_maps [bool] Skip generation of contact maps. Useful for capture-C. Default: False
--skip_ice [bool] Skip ICE normalization. Default: False
@@ -100,9 +102,15 @@ 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(", ")}"
}
+if (params.digest && params.digestion && !params.digest.containsKey(params.digestion)) {
+ exit 1, "Unknown digestion protocol. Currently, the available digestion options are ${params.digest.keySet().join(", ")}. Please set manually the '--restriction_site' and '--ligation_site' parameters."
+}
+params.restriction_site = params.digestion ? params.digest[ params.digestion ].restriction_site ?: false : false
+params.ligation_site = params.digestion ? params.digest[ params.digestion ].ligation_site ?: false : false
+
// 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"
+ exit 1, "Ligation motif not found. Please either use the `--digestion` parameters or specify the `--restriction_site` and `--ligation_site`. For DNase Hi-C, please use '--dnase' option"
}
// Reference index path configuration
@@ -148,27 +156,37 @@ if (params.input_paths){
Channel
.from( params.input_paths )
.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])] }
- }else{
+ .separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0] + "_R1", a[1][0]), tuple(a[0] + "_R2", a[1][1])] }
+}else{
raw_reads = Channel.create()
raw_reads_2 = Channel.create()
- Channel
- .fromFilePairs( params.input )
- .separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0], a[1][0]), tuple(a[0], a[1][1])] }
+ if ( params.split_fastq ){
+ Channel
+ .fromFilePairs( params.input, flat:true )
+ .splitFastq( by: params.fastq_chunks_size, pe:true, file: true, compress:true)
+ .separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0] + "_R1", a[1]), tuple(a[0] + "_R2", a[2])] }
+ }else{
+ Channel
+ .fromFilePairs( params.input )
+ .separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0] + "_R1", a[1][0]), tuple(a[0] + "_R2", a[1][1])] }
+ }
}
-// SPlit fastq files
-// https://www.nextflow.io/docs/latest/operator.html#splitfastq
-
-if ( params.split_fastq ){
- raw_reads_full = raw_reads.concat( raw_reads_2 )
- raw_reads = raw_reads_full.splitFastq( by: params.fastq_chunks_size, file: true)
- }else{
- raw_reads = raw_reads.concat( raw_reads_2 ).dump(tag: "data")
+// Update sample name if splitFastq is used
+def updateSampleName(x) {
+ if ((matcher = x[1] =~ /\s*(\.[\d]+).fastq.gz/)) {
+ res = matcher[0][1]
+ }
+ return [x[0] + res, x[1]]
}
+if (params.split_fastq ){
+ raw_reads = raw_reads.concat( raw_reads_2 ).map{it -> updateSampleName(it)}.dump(tag:'input')
+}else{
+ raw_reads = raw_reads.concat( raw_reads_2 ).dump(tag:'input')
+}
/*
* Other input channels
@@ -222,12 +240,12 @@ else if ( params.fasta && params.restriction_site ){
.ifEmpty { exit 1, "Restriction fragments: Fasta file not found: ${params.fasta}" }
.set { fasta_for_resfrag }
}
-else {
+else if (! params.dnase) {
exit 1, "No restriction fragments file specified!"
}
// Resolutions for contact maps
-map_res = Channel.from( params.bin_size.tokenize(',') )
+map_res = Channel.from( params.bin_size.toString() ).splitCsv().flatten()
/**********************************************************
* SET UP LOGS
@@ -243,17 +261,21 @@ summary['splitFastq'] = params.split_fastq
if (params.split_fastq)
summary['Read chunks Size'] = params.fastq_chunks_size
summary['Fasta Ref'] = params.fasta
-summary['Restriction Motif']= params.restriction_site
-summary['Ligation Motif'] = params.ligation_site
-summary['DNase Mode'] = params.dnase
-summary['Remove Dup'] = params.rm_dup
-summary['Remove MultiHits'] = params.rm_multi
+if (params.restriction_site){
+ summary['Digestion'] = params.digestion
+ summary['Restriction Motif']= params.restriction_site
+ summary['Ligation Motif'] = params.ligation_site
+ summary['Min Fragment Size']= params.min_restriction_fragment_size
+ summary['Max Fragment Size']= params.max_restriction_fragment_size
+ summary['Min Insert Size'] = params.min_insert_size
+ summary['Max Insert Size'] = params.max_insert_size
+}else{
+ summary['DNase Mode'] = params.dnase
+ summary['Min CIS dist'] = params.min_cis_dist
+}
summary['Min MAPQ'] = params.min_mapq
-summary['Min Fragment Size']= params.min_restriction_fragment_size
-summary['Max Fragment Size']= params.max_restriction_fragment_size
-summary['Min Insert Size'] = params.min_insert_size
-summary['Max Insert Size'] = params.max_insert_size
-summary['Min CIS dist'] = params.min_cis_dist
+summary['Keep Duplicates'] = params.keep_dups ? 'Yes' : 'No'
+summary['Keep Multihits'] = params.keep_multi ? 'Yes' : 'No'
summary['Maps resolution'] = params.bin_size
summary['Max Resources'] = "$params.max_memory memory, $params.max_cpus cpus, $params.max_time time per job"
if (workflow.containerEngine) summary['Container'] = "$workflow.containerEngine - $workflow.container"
@@ -330,10 +352,10 @@ def create_workflow_summary(summary) {
def yaml_file = workDir.resolve('workflow_summary_mqc.yaml')
yaml_file.text = """
- id: 'nf-core-chipseq-summary'
+ id: 'nf-core-hic-summary'
description: " - this information is collected when the pipeline is started."
- section_name: 'nf-core/chipseq Workflow Summary'
- section_href: 'https://github.com/nf-core/chipseq'
+ section_name: 'nf-core/hic Workflow Summary'
+ section_href: 'https://github.com/nf-core/hic'
plot_type: 'html'
data: |
<dl class=\"dl-horizontal\">
@@ -423,7 +445,7 @@ if(!params.restriction_fragments && params.fasta && !params.dnase){
*/
process bowtie2_end_to_end {
- tag "$prefix"
+ tag "$sample"
label 'process_medium'
publishDir path: { params.save_aligned_intermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.save_aligned_intermediates ? it : null }, mode: params.publish_dir_mode
@@ -433,13 +455,12 @@ process bowtie2_end_to_end {
file index from bwt2_index_end2end.collect()
output:
- set val(prefix), file("${prefix}_unmap.fastq") into unmapped_end_to_end
- set val(prefix), file("${prefix}.bam") into end_to_end_bam
+ set val(sample), file("${prefix}_unmap.fastq") into unmapped_end_to_end
+ set val(sample), file("${prefix}.bam") into end_to_end_bam
script:
prefix = reads.toString() - ~/(\.fq)?(\.fastq)?(\.gz)?$/
def bwt2_opts = params.bwt2_opts_end2end
-
if (!params.dnase){
"""
bowtie2 --rg-id BMG --rg SM:${prefix} \\
@@ -462,7 +483,7 @@ process bowtie2_end_to_end {
}
process trim_reads {
- tag "$prefix"
+ tag "$sample"
label 'process_low'
publishDir path: { params.save_aligned_intermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.save_aligned_intermediates ? it : null }, mode: params.publish_dir_mode
@@ -471,12 +492,13 @@ process trim_reads {
!params.dnase
input:
- set val(prefix), file(reads) from unmapped_end_to_end
+ set val(sample), file(reads) from unmapped_end_to_end
output:
- set val(prefix), file("${prefix}_trimmed.fastq") into trimmed_reads
+ set val(sample), file("${prefix}_trimmed.fastq") into trimmed_reads
script:
+ prefix = reads.toString() - ~/(\.fq)?(\.fastq)?(\.gz)?$/
"""
cutsite_trimming --fastq $reads \\
--cutsite ${params.ligation_site} \\
@@ -485,7 +507,7 @@ process trim_reads {
}
process bowtie2_on_trimmed_reads {
- tag "$prefix"
+ tag "$sample"
label 'process_medium'
publishDir path: { params.save_aligned_intermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.save_aligned_intermediates ? it : null }, mode: params.publish_dir_mode
@@ -494,11 +516,11 @@ process bowtie2_on_trimmed_reads {
!params.dnase
input:
- set val(prefix), file(reads) from trimmed_reads
+ set val(sample), file(reads) from trimmed_reads
file index from bwt2_index_trim.collect()
output:
- set val(prefix), file("${prefix}_trimmed.bam") into trimmed_bam
+ set val(sample), file("${prefix}_trimmed.bam") into trimmed_bam
script:
prefix = reads.toString() - ~/(_trimmed)?(\.fq)?(\.fastq)?(\.gz)?$/
@@ -512,22 +534,24 @@ process bowtie2_on_trimmed_reads {
}
if (!params.dnase){
- process merge_mapping_steps{
- tag "$sample = $bam1 + $bam2"
+ process bowtie2_merge_mapping_steps{
+ tag "$prefix = $bam1 + $bam2"
label 'process_medium'
publishDir path: { params.save_aligned_intermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.save_aligned_intermediates ? it : null }, mode: params.publish_dir_mode
input:
- set val(prefix), file(bam1), file(bam2) from end_to_end_bam.join( trimmed_bam )
+ set val(prefix), file(bam1), file(bam2) from end_to_end_bam.join( trimmed_bam ).dump(tag:'merge')
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|_1$|_2)/
- tag = prefix.toString() =~/_R1|_val_1|_1/ ? "R1" : "R2"
+ //sample = prefix.toString() - ~/(_R1|_R2|_val_1|_val_2|_1|_2)/
+ sample = prefix.toString() - ~/(_R1|_R2)/
+ //tag = prefix.toString() =~/_R1|_val_1|_1/ ? "R1" : "R2"
+ tag = prefix.toString() =~/_R1/ ? "R1" : "R2"
oname = prefix.toString() - ~/(\.[0-9]+)$/
"""
samtools merge -@ ${task.cpus} \\
@@ -554,46 +578,48 @@ if (!params.dnase){
}
}else{
process dnase_mapping_stats{
- tag "$sample = $bam1"
+ tag "$sample = $bam"
label 'process_medium'
publishDir path: { params.save_aligned_intermediates ? "${params.outdir}/mapping" : params.outdir },
saveAs: { params.save_aligned_intermediates ? it : null }, mode: params.publish_dir_mode
input:
- set val(prefix), file(bam1) from end_to_end_bam
+ set val(prefix), file(bam) from end_to_end_bam
output:
- set val(sample), file(bam1) into bwt2_merged_bam
+ set val(sample), file(bam) into bwt2_merged_bam
set val(oname), file("${prefix}.mapstat") into all_mapstat
script:
- sample = prefix.toString() - ~/(_R1|_R2|_val_1|_val_2|_1|_2)/
- tag = prefix.toString() =~/_R1|_val_1|_1/ ? "R1" : "R2"
+ //sample = prefix.toString() - ~/(_R1|_R2|_val_1|_val_2|_1|_2)/
+ sample = prefix.toString() - ~/(_R1|_R2)/
+ //tag = prefix.toString() =~/_R1|_val_1|_1/ ? "R1" : "R2"
+ tag = prefix.toString() =~/_R1/ ? "R1" : "R2"
oname = prefix.toString() - ~/(\.[0-9]+)$/
"""
echo "## ${prefix}" > ${prefix}.mapstat
echo -n "total_${tag}\t" >> ${prefix}.mapstat
- samtools view -c ${bam1} >> ${prefix}.mapstat
+ samtools view -c ${bam} >> ${prefix}.mapstat
echo -n "mapped_${tag}\t" >> ${prefix}.mapstat
- samtools view -c -F 4 ${bam1} >> ${prefix}.mapstat
+ samtools view -c -F 4 ${bam} >> ${prefix}.mapstat
echo -n "global_${tag}\t" >> ${prefix}.mapstat
- samtools view -c -F 4 ${bam1} >> ${prefix}.mapstat
+ samtools view -c -F 4 ${bam} >> ${prefix}.mapstat
echo -n "local_${tag}\t0" >> ${prefix}.mapstat
"""
}
}
-process combine_mapped_files{
+process combine_mates{
tag "$sample = $r1_prefix + $r2_prefix"
label 'process_low'
publishDir "${params.outdir}/mapping", mode: params.publish_dir_mode,
saveAs: {filename -> filename.indexOf(".pairstat") > 0 ? "stats/$filename" : "$filename"}
input:
- set val(sample), file(aligned_bam) from bwt2_merged_bam.groupTuple()
+ set val(sample), file(aligned_bam) from bwt2_merged_bam.groupTuple().dump(tag:'mates')
output:
- set val(sample), file("${sample}_bwt2pairs.bam") into paired_bam
+ set val(oname), file("${sample}_bwt2pairs.bam") into paired_bam
set val(oname), file("*.pairstat") into all_pairstat
script:
@@ -604,9 +630,11 @@ process combine_mapped_files{
oname = sample.toString() - ~/(\.[0-9]+)$/
def opts = "-t"
- opts = params.rm_singleton ? "${opts}" : "--single ${opts}"
- opts = params.rm_multi ? "${opts}" : "--multi ${opts}"
- if ("$params.min_mapq".isInteger()) opts="${opts} -q ${params.min_mapq}"
+ if (params.keep_multi) {
+ opts="${opts} --multi"
+ }else if (params.min_mapq){
+ opts="${opts} -q ${params.min_mapq}"
+ }
"""
mergeSAM.py -f ${r1_bam} -r ${r2_bam} -o ${sample}_bwt2pairs.bam ${opts}
"""
@@ -615,7 +643,7 @@ process combine_mapped_files{
/*
* STEP2 - DETECT VALID PAIRS
-*/
+ */
if (!params.dnase){
process get_valid_interaction{
@@ -697,7 +725,7 @@ process remove_duplicates {
saveAs: {filename -> filename.indexOf("*stat") > 0 ? "stats/$sample/$filename" : "$filename"}
input:
- set val(sample), file(vpairs) from valid_pairs.groupTuple()
+ set val(sample), file(vpairs) from valid_pairs.groupTuple().dump(tag:'final')
output:
set val(sample), file("*.allValidPairs") into all_valid_pairs
@@ -705,7 +733,7 @@ process remove_duplicates {
file("stats/") into all_mergestat
script:
- if ( params.rm_dup ){
+ if ( ! params.keep_dups ){
"""
mkdir -p stats/${sample}
diff --git a/nextflow.config b/nextflow.config
index bd5235b8ae579131bbc8dde2a62f6c3d044aa476..6e69e46d4f8b51058abd1e92f373a2307294f19c 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -8,51 +8,71 @@
// Global default params, used in configs
params {
- // Workflow flags
+ // Inputs / outputs
genome = false
input = "data/*{1,2}.fastq.gz"
- single_end = false
-
outdir = './results'
genome = false
input_paths = false
- split_fastq = false
- fastq_chunks_size = 20000000
chromosome_size = false
restriction_fragments = false
- skip_maps = false
- skip_ice = false
- skip_cool = false
- skip_multiqc = false
save_reference = false
- save_interaction_bam = false
- save_aligned_intermediates = 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 = 10
-
+
// Digestion Hi-C
- restriction_site = 'A^AGCTT'
- ligation_site = 'AAGCTAGCTT'
+ digestion = false
+ digest {
+ 'hindiii'{
+ restriction_site='A^AGCTT'
+ ligation_site='AAGCTAGCTT'
+ }
+ 'mboi' {
+ restriction_site='^GATC'
+ ligation_site='GATCGATC'
+ }
+ 'dpnii' {
+ restriction_site='^GATC'
+ ligation_site='GATCGATC'
+ }
+ 'arima' {
+ restriction_site='^GATC,G^ANT'
+ ligation_site='GATCGATC,GATCGANT,GANTGATC,GANTGANT'
+ }
+ }
min_restriction_fragment_size = 0
max_restriction_fragment_size = 0
min_insert_size = 0
max_insert_size = 0
+
+ // Dnase Hi-C
dnase = false
min_cis_dist = 0
- rm_dup = true
- rm_singleton = true
- rm_multi = true
+
+ // Mapping
+ split_fastq = false
+ fastq_chunks_size = 20000000
+ save_interaction_bam = false
+ save_aligned_intermediates = 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'
+ keep_dups = false
+ keep_multi = false
+ min_mapq = 10
+
+ // Contact maps
bin_size = '1000000,500000'
ice_max_iter = 100
ice_filer_low_count_perc = 0.02
ice_filer_high_count_perc = 0
ice_eps = 0.1
-
- publish_dir_mode = 'copy'
+ // Workflow
+ skip_maps = false
+ skip_ice = false
+ skip_cool = false
+ skip_multiqc = false
+
// Boilerplate options
+ publish_dir_mode = 'copy'
multiqc_config = false
name = false
email = false
diff --git a/nextflow_schema.json b/nextflow_schema.json
index 84bb558fd02002db37df4755daa2dff36d09d89f..c0733c3263055245da9996ae775173c488c5f03e 100644
--- a/nextflow_schema.json
+++ b/nextflow_schema.json
@@ -5,11 +5,11 @@
"description": "Analysis of Chromosome Conformation Capture data (Hi-C)",
"type": "object",
"definitions": {
- "input_output_options": {
- "title": "Input/output options",
+ "mandatory_arguments": {
+ "title": "Mandatory arguments",
"type": "object",
"fa_icon": "fas fa-terminal",
- "description": "Define where the pipeline should find input data and save output data.",
+ "description": "Mandatory arguments to run the pipeline",
"required": [
"input"
],
@@ -26,22 +26,11 @@
"description": "Input FastQ files for test only",
"default": "undefined"
},
- "split_fastq": {
- "type": "boolean",
- "description": "Split the reads into chunks before running the pipelne",
- "fa_icon": "fas fa-dna",
- "default": "false"
- },
- "fastq_chunks_size":{
- "type": "integer",
- "description": "Read number per chunks if split_fastq is used",
- "default": "20000000"
- },
- "single_end": {
- "type": "boolean",
- "description": "Specifies that the input is single-end reads.",
- "fa_icon": "fas fa-align-center",
- "help_text": "By default, the pipeline expects paired-end data. If you have single-end data, you need to specify `--single_end` on the command line when you launch the pipeline. A normal glob pattern, enclosed in quotation marks, can then be used for `--input`. For example:\n\n```bash\n--single_end --input '*.fastq'\n```\n\nIt is not possible to run a mixture of single-end and paired-end files in one run."
+ "genome": {
+ "type": "string",
+ "description": "Name of iGenomes reference.",
+ "fa_icon": "fas fa-book",
+ "help_text": "If using a reference genome configured in the pipeline using iGenomes, use this parameter to give the ID for the reference. This is then used to build the full paths for all required reference genome files e.g. `--genome GRCh38`.\n\nSee the [nf-core website docs](https://nf-co.re/usage/reference_genomes) for more details."
},
"outdir": {
"type": "string",
@@ -64,12 +53,6 @@
"fa_icon": "fas fa-dna",
"description": "Options for the reference genome indices used to align reads.",
"properties": {
- "genome": {
- "type": "string",
- "description": "Name of iGenomes reference.",
- "fa_icon": "fas fa-book",
- "help_text": "If using a reference genome configured in the pipeline using iGenomes, use this parameter to give the ID for the reference. This is then used to build the full paths for all required reference genome files e.g. `--genome GRCh38`.\n\nSee the [nf-core website docs](https://nf-co.re/usage/reference_genomes) for more details."
- },
"fasta": {
"type": "string",
"fa_icon": "fas fa-font",
@@ -94,6 +77,29 @@
"type": "string",
"description": "Full path to directory containing Bowtie index including base name. i.e. `/path/to/index/base`.",
"fa_icon": "far fa-file-alt"
+ }
+ }
+ },
+ "digestion_hi_c": {
+ "title": "Digestion Hi-C",
+ "type": "object",
+ "description": "Parameters for protocols based on restriction enzyme",
+ "default": "",
+ "properties": {
+ "digestion": {
+ "type": "string",
+ "default": "hindiii",
+ "description": "Name of restriction enzyme to automatically set the restriction_site and ligation_site options"
+ },
+ "restriction_site": {
+ "type": "string",
+ "default": "'A^AGCTT'",
+ "description": "Restriction motifs used during digestion. Several motifs (comma separated) can be provided."
+ },
+ "ligation_site": {
+ "type": "string",
+ "default": "'AAGCTAGCTT",
+ "description": "Expected motif after DNA ligation. Several motifs (comma separated) can be provided."
},
"chromosome_size": {
"type": "string",
@@ -115,41 +121,40 @@
}
}
},
- "data_processing_options": {
- "title": "Data processing",
+ "dnase_hi_c": {
+ "title": "DNAse Hi-C",
"type": "object",
- "description": "Parameters for Hi-C data processing",
+ "description": "Parameters for protocols based on DNAse digestion",
"default": "",
- "fa_icon": "fas fa-bahai",
"properties": {
"dnase": {
"type": "boolean",
"description": "For Hi-C protocols which are not based on enzyme digestion such as DNase Hi-C"
},
- "restriction_site": {
- "type": "string",
- "default": "'A^AGCTT'",
- "description": "Restriction motifs used during digestion. Several motifs (comma separated) can be provided."
- },
- "ligation_site": {
- "type": "string",
- "default": "'AAGCTAGCTT",
- "description": "Expected motif after DNA ligation. Several motifs (comma separated) can be provided."
- },
- "rm_dup": {
- "type": "boolean",
- "description": "Remove duplicates",
- "default": true
- },
- "rm_multi": {
+ "min_cis_dist": {
+ "type": "integer",
+ "default": "O",
+ "description": "Minimum distance between loci to consider. Useful for --dnase mode to remove spurious ligation products. Only values > 0 are considered"
+ }
+ }
+ },
+ "alignments": {
+ "title": "Alignments",
+ "type": "object",
+ "description": "Parameters for reads aligments",
+ "default": "",
+ "fa_icon": "fas fa-bahai",
+ "properties": {
+ "split_fastq": {
"type": "boolean",
- "description": "Remove multi-mapped reads",
- "default": true
+ "description": "Split the reads into chunks before running the pipelne",
+ "fa_icon": "fas fa-dna",
+ "default": "false"
},
- "rm_singleton": {
- "type": "boolean",
- "description": "Remove singleton",
- "default": true
+ "fastq_chunks_size": {
+ "type": "integer",
+ "description": "Read number per chunks if split_fastq is used",
+ "default": "20000000"
},
"min_mapq": {
"type": "integer",
@@ -166,27 +171,28 @@
"default": "'--very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder'",
"description": "Option for trimmed reads mapping"
},
- "save_interaction_bam": {
- "type": "boolean",
- "description": "Save a BAM file where all reads are flagged by their interaction classes"
- },
"save_aligned_intermediates": {
"type": "boolean",
"description": "Save all BAM files during two-steps mapping"
}
}
},
- "contacts_calling_options": {
- "title": "Contacts calling",
+ "valid_pairs_detection": {
+ "title": "Valid Pairs Detection",
"type": "object",
"description": "Options to call significant interactions",
"default": "",
"fa_icon": "fas fa-signature",
"properties": {
- "min_cis_dist": {
- "type": "integer",
- "default": "O",
- "description": "Minimum distance between loci to consider. Useful for --dnase mode to remove spurious ligation products. Only values > 0 are considered"
+ "keep_dups": {
+ "type": "string",
+ "description": "Keep duplicated reads",
+ "default": "False"
+ },
+ "keep_multi": {
+ "type": "string",
+ "description": "Keep multi-aligned reads",
+ "default": "False"
},
"max_insert_size": {
"type": "integer",
@@ -207,6 +213,10 @@
"type": "integer",
"default": "0",
"description": "Minimum restriction fragment size to consider. Only values > 0 are considered"
+ },
+ "save_interaction_bam": {
+ "type": "boolean",
+ "description": "Save a BAM file where all reads are flagged by their interaction classes"
}
}
},
@@ -435,16 +445,22 @@
},
"allOf": [
{
- "$ref": "#/definitions/input_output_options"
+ "$ref": "#/definitions/mandatory_arguments"
},
{
"$ref": "#/definitions/reference_genome_options"
},
{
- "$ref": "#/definitions/data_processing_options"
+ "$ref": "#/definitions/digestion_hi_c"
+ },
+ {
+ "$ref": "#/definitions/dnase_hi_c"
+ },
+ {
+ "$ref": "#/definitions/alignments"
},
{
- "$ref": "#/definitions/contacts_calling_options"
+ "$ref": "#/definitions/valid_pairs_detection"
},
{
"$ref": "#/definitions/contact_maps_options"
@@ -462,4 +478,4 @@
"$ref": "#/definitions/institutional_config_options"
}
]
-}
+}
\ No newline at end of file