Merge branch 'develop'

Author: Maarten-vd-Sande

.gitignore

@@ -4,6 +4,7 @@ slurm*.out
 report.html
 samples.tsv
 config.yaml
+seq2science.*.log
 build/
 *.egg*
 .idea/

CHANGELOG.md

@@ -8,6 +8,32 @@ All changed fall under either one of these types: `Added`, `Changed`, `Deprecate
 
 ## [Unreleased]
 
+### Added
+
+- added CLI functionality to the deseq2.R script (try it with `Rscript /path/to/deseq2.R --help`!)
+- --force flag to seq2science init to automatically overwrite existing samples.tsv and config.yaml
+- local fastqs with Illumina's '_100' are now recognized
+- added the workflow explanation to the multiqc report
+
+### Changed
+
+- config checks: all keys converted to lower case & duplicate keys throw an exception
+- MultiQC updated to v1.10
+- Link to seq2science log instead of snakemake log in final message
+
+### Fixed
+
+- Issue when filtering a combination of single-end and paired-end reads on template length
+- explain functionality testing
+- scATAC can properly use SE fastqs
+- scRNA can use fqexts other than R1/R2
+- fastq renaming works again
+- added missing schemas to extended docs
+
+### Fixed
+
+- Bug with edgeR.upperquartile normalization. Now makes everything NaN, so pipeline finishes succesfully.
+
 ## [0.5.0] - 2021-03-03
 
 Version 0.5.0 brings many quality of life improvements, such as seq2science automatically inferring what needs to be re-run when changing the samples.tsv and/or the config.yaml, differential peak analysis for chip/atac workflows and tab-completion!

docs/content/DESeq2.md

@@ -2,20 +2,21 @@
 Seq2science outputs counts matrices for each assembly in any bulk-sequencing workflow, which can optionally can be used for differential gene/peak analysis with [DESeq2](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8).
 To do so, add one or more contrast design(s) in the `config.yaml` file, with matching identifiers in the `samples.tsv` file.
 
+After completing the workflow, rerunning Seq2science with new contrasts will only perform the new analyses.
+Alternatively, you can call the DESeq2 script from the command line.
+To get started with the command line script, try out `Rscript /path/to/seq2science/scripts/deseq2.R --help`.
+
 This section details how the contrast designs are created.
 Examples are given [below](./DESeq2.html#deseq2-contrast-designs), and each `config.yaml` contains a commented-out example contrast design at the bottom.
 
-Note: (additional) design contrasts can be added at any time.
-After completing the workflow, rerunning Seq2science with new contrasts will only perform the new analyses.
-
 ##### Overview of the DESeq2 method
 DESeq2 automatically performs library bias correction when loading your data, and batch correction is performed if it is included in the contrast design.
 After calculating differentially expressed genes/peaks, a multiple testing procedure is applied, which is either the Benjamini-Hochberg procedure (the default) or Independent Hypothesis Weighing.
 The False Discovery Rate cutoff is set by alpha, which is 0.1 by default.
 Finally, count values are log transformed and shrunk (by default using the apeglm method).
 These defaults can be changed in the [config.yaml](./schemas.html#deseq2), under the `deseq2` variables using the `multiple_testing_procedure`, `alpha_value` and `shrinkage_estimator` options respectively.
 
-For more information, check out the steps in this [vignette](www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html), which Seq2science follows.
+For more information, check out the steps in this [vignette](https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html), which Seq2science follows.
 
 ##### Output
 For each contrast design, the list of *all* genes/peaks is saved to file, with analysis results for expressed genes. Briefly, these include:

docs/content/faq.md

@@ -52,4 +52,4 @@ There are differences between these methods however, and caution is advised with
 Additionally, Seq2science employs strict quality checks by default.
 You could try reducing the `min_mapping_quality` (for HTSeq and FeatureCounts) or change the parameters to the quantifier.
 
-The biggest difference in expressed genes found is in switching to or from Salmon.
+The biggest differences in expressed genes can be found by switching gene annotation, followed by switching quantifier (to or from Salmon).

docs/content/workflows.rst

@@ -11,7 +11,7 @@ These workflows generally overlap a lot, so once you understand the logic of one
    Downloading fastqs <workflows/download_fastq.md>
    Alignment <workflows/alignment.md>
    ATAC-seq <workflows/atac_seq.md>
-   RNA-seq <workflows/rna_seq.md>
    ChIP-seq <workflows/chip_seq.md>
+   RNA-seq <workflows/rna_seq.md>
    scATAC-seq <workflows/scatac_seq.md>
    scRNA-seq <workflows/scrna_seq.md>

docs/content/workflows/alignment.md

@@ -8,26 +8,28 @@ Aligning samples has never been easier!
 
 #### Downloading of sample(s)
 Depending on whether the samples you start seq2science with is your own data, public data, or a mix, the pipeline might start with downloading samples.
-Take a look at the [downloading_fastq](./download_fastq.html) workflow for extensive documentation about downloading of public samples.
+You control which samples are used in the [samples.tsv](#filling-out-the-samples-tsv).
+Background on public data can be found [here](./download_fastq.html#download-sra-file).
 
 #### Read trimming
-The pipeline starts by trimming the reads with Trim galore or Fastp.
-The trimmer will automatically trims the low quality 3' ends of reads, and removes short reads.
-After the quality trimming it automatically detects which adapter was used, and trims it.
+The pipeline starts by trimming the reads with Trim Galore! or Fastp (the default).
+The trimmer will automatically trim the low quality 3' ends of reads, and removes short reads.
+During the quality trimming it automatically detects which adapter was used, and trims this as well.
 Trimming parameters for the pipeline can be set in the configuration.
 
 #### Alignment & Sorting
 After trimming the reads are aligned against an assembly. Currently we support *bowtie2*, *bwa*, *bwa-mem2*, *hisat2* and *STAR* as aligners. Choosing which aligner is as easy as setting the *aligner* variable in the `config.yaml`, for example: `aligner: bwa`. Sensible defaults have been set for every aligner, but can be overwritten for either (or both) the indexing and alignment by specifying them in the `config.yaml`:
 ```
 aligner: 
-  star:
-    index: '--limitGenomeGenerateRAM 55000000000`  # give STAR 55 GB RAM during indexing (default is 31 GB)
-    align: '--readFilesCommand gunzip -c'          # unzips fastqs before aligning
+  bwa-mem:
+    index: '-a bwtsw'
+    align: '-M'
 ```
 
-The pipeline will check if the assembly you specified is present in the *genome_dir*, and otherwise will download it for you through [genomepy](https://github.com/vanheeringen-lab/genomepy). All these aligners require an index to be formed first for each assembly, but don't worry, the pipeline does this for you. 
+The pipeline will check if the assembly you specified is present in the *genome_dir*, and otherwise will download it for you through [genomepy](https://github.com/vanheeringen-lab/genomepy).
+All these aligners require an index to be formed first for each assembly, but don't worry, the pipeline does this for you.
 
-The outputted alignment.bam is immediately sorted by either samtools or sambamba (*bam_sorter*) either in *queryname* or *coordinate* (default) order. Take a look at our [Choosing an aligner](https://github.com/vanheeringen-lab/snakemake-workflows/wiki/3.-Alignment#choosing-an-appropriate-aligner) section for tips which aligner to use.
+The outputted alignment.bam is immediately sorted by either samtools or sambamba (*bam_sorter*) either in *queryname* or *coordinate* (default) order.
 
 #### Custom assembly extensions
 The genome and/or gene annotation can be extended with custom files, such as ERCC spike-ins for scRNA-seq. To do so, add `custom_genome_extension: path/to/spike_in.fa` and `custom_annotation_extension: path/to/spike_in.gtf` to the config.
@@ -47,7 +49,9 @@ A UCSC compatible trackhub can be generated for this workflow.
 See the [trackhub page](../results.html#trackhub)<!-- @IGNORE PREVIOUS: link --> for more information!
 
 ### Filling out the samples.tsv
-Before running a workflow you will have to specify which samples you want to run the workflow on. Each workflow starts with a `samples.tsv` as an example, and you should adapt it to your specific needs. As an example, the `samples.tsv` could look something like this:
+Before running a workflow you will have to specify which samples you want to run the workflow on.
+Each workflow starts with a `samples.tsv` as an example, and you should adapt it to your specific needs.
+As an example, the `samples.tsv` could look something like this:
 ```
 sample    assembly    technical_replicate    descriptive_name
 GSM123    GRCh38      heart_1      heart_merged
@@ -58,13 +62,30 @@ GSM890    danRer11    stage_9      stage_9
 ```
 
 #### Sample column
-This column is necessary for all workflows, not just the atac-seq workflow. If you use the pipeline on public data this should be the name of the accession (e.g. GSM2837484). If you use the pipeline on local data this should be the *basename* of the file without the *extension(s)*. For instance, `/home/user/myfastqs/sample1.fastq.gz` would be `sample1` (for single-ended data). For paired-ended data `/home/user/myfastqs/sample2_R1.fastq.gz` and `/home/user/myfastqs/sample2_R2.fastq.gz` would be `sample2`.
+If you use the pipeline on **public data** this should be the name of the accession (e.g. GSM2837484).
+Accepted formats start with "GSM", "SRR", "SRX", "DRR", "DRX", "ERR" or "ERX".
+
+If you use the pipeline on **local data** this should be the *basename* of the file without the *extension(s)*.
+For example:
+- `/home/user/myfastqs/sample1.fastq.gz` -------> `sample1` for single-ended data
+- `/home/user/myfastqs/sample2_R1.fastq.gz` ┬> `sample2` for paired-ended data <br> `/home/user/myfastqs/sample2_R2.fastq.gz` ┘
+
+For **local data**, some fastq files may have slightly different naming formats.
+For instance, Illumina may produce a sample named `sample3_S1_L001_R1_001.fastq.gz` (and the `R2` fastq).
+Seq2science will attempt to recognize these files based on the sample name `sample3`.
+
+For **both local and public data**, identifiers used to recognize fastq files are the fastq read extensions (`R1` and `R2` by default) and the fastq suffix (`fastq` by default).
+The directory where seq2science will store (or look for) fastqs is determined by the `fastq_dir` config option.
+In the example above, the `fastq_dir` should be set to `/home/user/myfastqs`.
+These setting can be changed in the `config.yaml`.
 
 #### Assembly column
-This column is necessary for all workflows, except the *downloading samples* workflow. Here you simply add the name of the assembly you want your samples aligned against and the workflow will download it for you.
+Here you simply add the name of the assembly you want your samples aligned against and the workflow will download it for you.
 
 #### Descriptive_name column
-The descriptive_name column is used for the trackhub and multiqc report. In the trackhub your tracks will be called after the descriptive name, and in the multiqc report there will be a button to rename your samples after this column. The descriptive name can not contain '-' characters, but underscores '_' are allowed.
+The descriptive\_name column is used for the trackhub and multiqc report.
+In the trackhub your tracks will be called after the descriptive name, and in the multiqc report there will be a button to rename your samples after this column.
+The descriptive name can not contain '-' characters, but underscores '\_' are allowed.
 
 #### Technical_replicate column
 Technical replicates, or any fastq file you may wish to merge, are set using the `technical_replicate` column in the samples.tsv file.

docs/content/workflows/atac_seq.md

@@ -8,15 +8,16 @@ Running an ATAC-seq analysis has never been easier!
 
 #### Downloading of sample(s)
 Depending on whether the samples you start seq2science with is your own data, public data, or a mix, the pipeline might start with downloading samples.
-Take a look at the [downloading_fastq](./download_fastq.html) workflow for extensive documentation about downloading of public samples.
+You control which samples are used in the [samples.tsv](#filling-out-the-samples-tsv).
+Background on public data can be found [here](./download_fastq.html#download-sra-file).
 
 #### Downloading and indexing of assembly(s)
 Depending on whether the assembly and its index you align your samples against already exist seq2science will start with downloading of the assembly through [genomepy](https://github.com/vanheeringen-lab/genomepy).
 
 #### Read trimming
-The pipeline starts by trimming the reads with Trim galore or Fastp.
-The trimmer will automatically trims the low quality 3' ends of reads, and removes short reads.
-After the quality trimming it automatically detects which adapter was used, and trims it.
+The pipeline starts by trimming the reads with Trim Galore! or Fastp (the default).
+The trimmer will automatically trim the low quality 3' ends of reads, and removes short reads.
+During the quality trimming it automatically detects which adapter was used, and trims this as well.
 Trimming parameters for the pipeline can be set in the configuration.
 
 #### Alignment
@@ -87,16 +88,33 @@ GSM890    danRer11    stage_9      stage_9             GSM234
 ```
 
 #### Sample column
-This column is necessary for all workflows, not just the ATAC-seq workflow. If you use the pipeline on public data this should be the name of the accession (e.g. GSM2837484), if you use the pipeline on local data this should be the *basename* of the file without the *extension*. For instance `/home/user/myfastqs/sample1.fastq.gz` would be `sample1`.
+If you use the pipeline on **public data** this should be the name of the accession (e.g. GSM2837484).
+Accepted formats start with "GSM", "SRR", "SRX", "DRR", "DRX", "ERR" or "ERX".
+
+If you use the pipeline on **local data** this should be the *basename* of the file without the *extension(s)*.
+For example:
+- `/home/user/myfastqs/sample1.fastq.gz` -------> `sample1` for single-ended data
+- `/home/user/myfastqs/sample2_R1.fastq.gz` ┬> `sample2` for paired-ended data <br> `/home/user/myfastqs/sample2_R2.fastq.gz` ┘
+
+For **local data**, some fastq files may have slightly different naming formats.
+For instance, Illumina may produce a sample named `sample3_S1_L001_R1_001.fastq.gz` (and the `R2` fastq).
+Seq2science will attempt to recognize these files based on the sample name `sample3`.
+
+For **both local and public data**, identifiers used to recognize fastq files are the fastq read extensions (`R1` and `R2` by default) and the fastq suffix (`fastq` by default).
+The directory where seq2science will store (or look for) fastqs is determined by the `fastq_dir` config option.
+In the example above, the `fastq_dir` should be set to `/home/user/myfastqs`.
+These setting can be changed in the `config.yaml`.
 
 #### Assembly column
-This column is necessary for all workflows, except the *downloading samples* workflow. Here you simply add the name of the assembly you want your samples aligned against and the workflow will download it for you.
+Here you simply add the name of the assembly you want your samples aligned against and the workflow will download it for you.
 
 #### Control column
 In the `control` column you can (optionally) add the "sample name" of the input control, although it is not common to generate this for ATAC-seq. 
 
 #### Descriptive_name column
-The descriptive_name column is used for the trackhub and multiqc report. In the trackhub your tracks will be called after the descriptive name, and in the multiqc report there will be a button to rename your samples after this column. The descriptive name can not contain '-' characters, but underscores '_' are allowed.
+The descriptive\_name column is used for the trackhub and multiqc report.
+In the trackhub your tracks will be called after the descriptive name, and in the multiqc report there will be a button to rename your samples after this column.
+The descriptive name can not contain '-' characters, but underscores '\_' are allowed.
 
 #### Technical_replicate column
 Technical replicates, or any fastq file you may wish to merge, are set using the `technical_replicate` column in the samples.tsv file.