Merge pull request #34 from luslab/dev

{

    "Exclude": ["README.md", "\\.bed$", "conf/igenomes.config"]

    "Exclude": ["README.md", "docs/output.md", "docs/usage.md", "\\.bed$", "conf/igenomes.config"]

}

## Introduction

**nf-core/cutandrun** is a bioinformatics best-practise analysis pipeline for CUT&Run and CUT&Tag sequencing data analysis to study protein-DNA interactions and epigenomic profiling.

**nf-core/cutandrun** is a best-practice bioinformatic analysis pipeline for CUT&Run and CUT&Tag experimental protocols that where developed to study protein-DNA interactions and epigenomic profiling.

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 uses Docker/Singularity containers making installation trivial and results highly reproducible. The [Nextflow DSL2](https://www.nextflow.io/docs/latest/dsl2.html) implementation of this pipeline uses one container per process which makes it much easier to maintain and update software dependencies. Where possible, these processes have been submitted to and installed from [nf-core/modules](https://github.com/nf-core/modules) in order to make them available to all nf-core pipelines, and to everyone within the Nextflow community!

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 is capable of using docker/singularity containers making installation trivial and results highly reproducible. The [Nextflow DSL2](https://www.nextflow.io/docs/latest/dsl2.html) implementation of this pipeline uses one container per process which makes it easier to maintain and update software dependencies. Where possible, these processes have been submitted to and installed from [nf-core/modules](https://github.com/nf-core/modules).

<!-- TODO nf-core: Add full-sized test dataset and amend the paragraph below if applicable -->

On release, automated continuous integration tests run the pipeline on a full-sized dataset on the AWS cloud infrastructure. This ensures that the pipeline runs on AWS, has sensible resource allocation defaults set to run on real-world datasets, and permits the persistent storage of results to benchmark between pipeline releases and other analysis sources. The results obtained from the full-sized test can be viewed on the [nf-core website](https://nf-co.re/cutandrun/results).

The pipeline has been developed with continuous integration (CI) in mind. nf-core code and module linting as well as a battery of over 100 unit and integration tests run on pull request to the main repository and on release of the pipeline. On official release, automated CI tests run the pipeline on a full-sized dataset on the AWS cloud infrastructure. This ensures that the pipeline runs on AWS, has sensible resource allocation defaults set to run on real-world datasets, and permits the persistent storage of results to benchmark between pipeline releases and other analysis sources. The results obtained from the full-sized test can be viewed on the [nf-core website](https://nf-co.re/cutandrun/results).

## Pipeline summary

1. Merge re-sequenced FastQ files ([`cat`](http://www.linfo.org/cat.html))

2. Read QC ([`FastQC`](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))

3. Adapter and quality trimming ([`Trim Galore!`](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/))

4. Alignment to both target and spike-in genomes ([`Bowtie 2`](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml))

5. Filter on quality, sort and index alignments ([`SAMtools`](https://sourceforge.net/projects/samtools/files/samtools/))

6. Duplicate read marking ([`picard MarkDuplicates`](https://broadinstitute.github.io/picard/))

7. Create bedGraph files ([`BEDTools`](https://github.com/arq5x/bedtools2/)

8. Create bigWig coverage files ([`bedGraphToBigWig`](http://hgdownload.soe.ucsc.edu/admin/exe/))

9. Peak calling specifically tailored for low background noise ([`SEACR`](https://github.com/FredHutch/SEACR))

10. Quality control and analysis:

    1. Alignment, fragment length and peak analysis and replicate reproducibility ([`Python`](https://www.python.org/))

    2. Differential peak analysis ([`DESeq2`](https://bioconductor.org/packages/release/bioc/html/DESeq2.html))

    3. Heatmap peak analysis ([`deepTools`](https://github.com/deeptools/deepTools/))

11. Genome browser session ([`IGV`](https://software.broadinstitute.org/software/igv/))

12. Present QC for raw read, alignment and duplicate reads ([`MultiQC`](http://multiqc.info/))

1. Check input files

2. Merge re-sequenced FastQ files ([`cat`](http://www.linfo.org/cat.html))

3. Read QC ([`FastQC`](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))

4. Adapter and quality trimming ([`Trim Galore!`](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/))

5. Alignment to both target and spike-in genomes ([`Bowtie 2`](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml))

6. Filter on quality, sort and index alignments ([`samtools`](https://sourceforge.net/projects/samtools/files/samtools/))

7. Duplicate read marking ([`picard`](https://broadinstitute.github.io/picard/))

8. Create bedGraph files ([`bedtools`](https://github.com/arq5x/bedtools2/)

9. Create bigWig coverage files ([`bedGraphToBigWig`](http://hgdownload.soe.ucsc.edu/admin/exe/))

10. Peak calling specifically tailored for low background noise experiments ([`SEACR`](https://github.com/FredHutch/SEACR))

11. Consensus peak merging and reporting ([`bedtools`](https://github.com/arq5x/bedtools2/))

12. Quality control and analysis:

    1. Alignment, fragment length and peak analysis and replicate reproducibility ([`python`](https://www.python.org/))

    2. Heatmap peak analysis ([`deepTools`](https://github.com/deeptools/deepTools/))

13. Genome browser session ([`IGV`](https://software.broadinstitute.org/software/igv/))

14. Present QC for raw read, alignment and duplicate reads ([`MultiQC`](http://multiqc.info/))

## Quick Start

        nextflow run nf-core/cutandrun \

            -profile <docker/singularity/podman/conda/institute> \

            --input samplesheet.csv \

            --genome GRCh37

            --genome GRCh38

        ```

See [usage docs](https://nf-co.re/cutandrun/usage) for all of the available options when running the pipeline.

## Samplesheet input

You will need to create a samplesheet file with information about the samples in your experiment before running the pipeline. Use this parameter to specify its location. It has to be a comma-separated file with 4 columns, and a header row as shown in the examples below.

You will need to create a samplesheet file with information about the samples in your experiment before running the pipeline. Use this parameter to specify its location. It has to be a comma-separated file with the correct data structure as shown in the examples below.

```bash

--input '[path to samplesheet file]'

--input <path to samplesheet file>

```

### Multiple replicates

The `group` identifier is the same when you have multiple biological replicates from the same experimental group, just increment the `replicate` identifier appropriately. A special case for `group` is if you have non-specific IgG antibody control data that can be used for normalising your experimental CUT&Run (OR CUT&Tag) data. In this case, the `group` name for the IgG control data _must_ be set to `igg`. Below is an example for a single target group in triplicate, complemented by an IgG control duplicate:

```bash

group,replicate,fastq_1,fastq_2

target,1,H3K27me3_S1_L001_R1.fastq.gz,H3K27me3_S1_L001_R2.fastq.gz

target,2,H3K27me3_S2_L001_R1.fastq.gz,H3K27me3_S2_L001_R2.fastq.gz

target,3,H3K27me3_S3_L001_R1.fastq.gz,H3K27me3_S3_L001_R2.fastq.gz

igg,1,IGG_S1_L001_R1.fastq.gz,IGG_S1_L001_R2.fastq.gz

igg,2,IGG_S2_L001_R1.fastq.gz,IGG_S2_L001_R2.fastq.gz

```

There are 4 use-cases for various combinations of experimental and IgG control replicate numbers that are note-worthy:

* One-to-one:

Across all experimental groups and IgG control, the number of replicates are the same and numbered in a uniform fashion. In this case, IgG normalisation at the peak-calling level is done for matching replicate numbers.

* Many experimental replicates, one IgG replicate:

Each experimental replicate will be normalised by the given single IgG replicate.

* Equal numbers of replicates across experimental groups and more than one IgG replicate:

If the user puts forward groups with equal amounts of uniformly numbered replicates, but different to the multiple IgG replicates, then a warning will show informing the user that just the first IgG replicate will be used for normalisation.

* Varying experimental and IgG replicates:

If the user puts forward groups with varying, ununiform experimental replicate numbers and more than one IgG replicate, an error will be thrown up and the pipeline will hault. If you wish to merge the IgG control data in this case, see next steps.

It is _recommended_ to have an IgG control for normalising your experimental data and this is the default action for the pipeline. However, if you run the pipeline without IgG control data you must supply `--igg_control false`.

### Multiple runs of the same library

The `group` and `replicate` identifiers are the same when you have re-sequenced the same sample more than once (e.g. to increase sequencing depth), or if you would like to merge technical replicates. The pipeline will concatenate the raw reads before alignment. Below is an example for two samples, one experimental and one control, sequenced across multiple lanes:

```bash

group,replicate,fastq_1,fastq_2

target,1,H3K27me3_S1_L001_R1.fastq.gz,H3K27me3_S1_L001_R2.fastq.gz

target,1,H3K27me3_S1_L002_R1.fastq.gz,H3K27me3_S1_L002_R2.fastq.gz

igg,1,IGG_S1_L001_R1.fastq.gz,IGG_S1_L001_R2.fastq.gz

igg,1,IGG_S1_L002_R1.fastq.gz,IGG_S1_L002_R2.fastq.gz

```

### Full design

A final design file may look something like the one below. This is for one experimental group in triplicate where the last replicate of the `treatment` group has been sequenced twice, another experimental group in duplicate, and one IgG control group.

An example sample sheet structure is shown below. This defines two target experimental groups for the histone marks h3k27me3 and h3k4me3 with two biological replicates per group. Each antibody target also had an IgG control performed alongside. The two IgG experiments are configured as biological replicates with the `igg` group keyword with a unique control group assignment. The target experiments are then assigned to the igg control group using the `control_group` column.

```bash

group,replicate,fastq_1,fastq_2

h3k27me3,1,H3K27me3_S1_L001_R1.fastq.gz,H3K27me3_S1_L001_R2.fastq.gz

h3k27me3,2,H3K27me3_S2_L001_R1.fastq.gz,H3K27me3_S2_L001_R2.fastq.gz

h3k27me3,3,H3K27me3_S3_L001_R1.fastq.gz,H3K27me3_S3_L001_R2.fastq.gz

h3k27me3,3,H3K27me3_S3_L002_R1.fastq.gz,H3K27me3_S3_L002_R2.fastq.gz

h3k4me3,1,H3K4me3_S1_L001_R1.fastq.gz,H3K4me3_S1_L001_R2.fastq.gz

h3k4me3,2,H3K4me3_S2_L001_R1.fastq.gz,H3K4me3_S2_L001_R2.fastq.gz

igg,1,IGG_S1_L001_R1.fastq.gz,IGG_S1_L001_R2.fastq.gz

igg,2,IGG_S2_L001_R1.fastq.gz,IGG_S2_L001_R2.fastq.gz

group,replicate,control_group,fastq_1,fastq_2

h3k27me3,1,1,H3K27ME3_S1_L001_R1.fastq.gz,H3K27ME3_S1_L001_R2.fastq.gz

h3k27me3,2,1,H3K27ME3_S2_L001_R1.fastq.gz,H3K27ME3_S2_L001_R2.fastq.gz

h3k4me3,1,2,H3K4ME3_S1_L001_R1.fastq.gz,H3K4ME3_S1_L001_R2.fastq.gz

h3k4me3,2,2,H3K4ME3_S2_L001_R1.fastq.gz,H3K4MES_S2_L001_R2.fastq.gz

igg,1,1,IGG_S1_L001_R1.fastq.gz,IGG_S1_L001_R2.fastq.gz

igg,2,2,IGG_S2_L001_R1.fastq.gz,IGG_S2_L001_R2.fastq.gz

```

| Column         | Description                                                                                                 |

|----------------|-------------------------------------------------------------------------------------------------------------|

| `group`        | Group identifier for sample. This will be identical for replicate samples from the same experimental group. |

| `replicate`    | Integer representing replicate number.                                                                      |

| `control_group`    | Integer representing the IgG control group the target is assigned to.                                                                      |

| `fastq_1`      | Full path to FastQ file for read 1. File has to be zipped and have the extension ".fastq.gz" or ".fq.gz".   |

| `fastq_2`      | Full path to FastQ file for read 2. File has to be zipped and have the extension ".fastq.gz" or ".fq.gz".   |

An [example samplesheet](../assets/samplesheet.csv) has been provided with the pipeline.

### Multiple replicates and IgG controls

To assign biological replicates to the same group use the same group identifier but increment the `replicate` column appropriately. To merge multiple fastq files from the same library or to merge technical replicate data, use identical `group` AND `replicate` column identifiers. The pipeline will automatically merge these samples and treat them as one biological replicate. An example of this type of sample sheet configuration is shown below where each biological replicate has two fastq files:

```bash

group,replicate,control_group,fastq_1,fastq_2

h3k27me3,1,1,H3K27ME3_S1_L001_R1.fastq.gz,H3K27ME3_S1_L001_R2.fastq.gz

h3k27me3,1,1,H3K27ME3_S2_L001_R1.fastq.gz,H3K27ME3_S2_L001_R2.fastq.gz

h3k27me3,2,1,H3K27ME3_S3_L001_R1.fastq.gz,H3K27ME3_S3_L001_R2.fastq.gz

h3k27me3,2,1,H3K27ME3_S4_L001_R1.fastq.gz,H3K27ME3_S4_L001_R2.fastq.gz

igg,1,1,IGG_S1_L001_R1.fastq.gz,IGG_S1_L001_R2.fastq.gz

```

To add any IgG controls that were processed, use the `igg` keyword in the `group` column and increment as with the target samples. The IgG control will be used to normalising your experimental CUT&Run (OR CUT&Tag) data. It is _recommended_ to have an IgG control for normalising your experimental data and this is the default action for the pipeline. However, if you run the pipeline without IgG control data you must supply `--igg_control false`.

## Running the pipeline

The typical command for running the pipeline is as follows:

```bash

nextflow run nf-core/cutandrun --input samplesheet.csv --genome GRCh37 -profile docker

nextflow run nf-core/cutandrun --input samplesheet.csv --genome GRCh38 -profile docker

```

This will launch the pipeline with the `docker` configuration profile. See below for more information about profiles.

Several generic profiles are bundled with the pipeline which instruct the pipeline to use software packaged using different methods (Docker, Singularity, Podman, Shifter, Charliecloud, Conda) - see below. When using Biocontainers, most of these software packaging methods pull Docker containers from quay.io e.g [FastQC](https://quay.io/repository/biocontainers/fastqc) except for Singularity which directly downloads Singularity images via https hosted by the [Galaxy project](https://depot.galaxyproject.org/singularity/) and Conda which downloads and installs software locally from [Bioconda](https://bioconda.github.io/).

> We highly recommend the use of Docker or Singularity containers for full pipeline reproducibility, however when this is not possible, Conda is also supported.

> We recommend the use of Docker or Singularity containers for full pipeline reproducibility, however when this is not possible, Conda is also supported.

The pipeline also dynamically loads configurations from [https://github.com/nf-core/configs](https://github.com/nf-core/configs) when it runs, making multiple config profiles for various institutional clusters available at run time. For more information and to see if your system is available in these configs please see the [nf-core/configs documentation](https://github.com/nf-core/configs#documentation).

    igg_control                = true

    peak_threshold             = 0.05

    replicate_threshold        = 1

    //skip_peak_calling          = false

    only_peak_calling          = false

    // Reporting and Visualisation