diff --git a/docs/CHANGELOG.rst b/docs/CHANGELOG.rst index d23d221ec..61cc28eb9 100644 --- a/docs/CHANGELOG.rst +++ b/docs/CHANGELOG.rst @@ -16,6 +16,9 @@ Added Changed ------- +- Deprecate the old yaml based ATAC-seq workflow in favour of the Python one + (previously in beta) implemented in the process ``workflow-atac-seq`` + and bump the version of the new workflow to 4.0.0 Fixed ----- diff --git a/resolwe_bio/processes/workflows/atac_seq_beta.py b/resolwe_bio/processes/workflows/atac_seq.py similarity index 99% rename from resolwe_bio/processes/workflows/atac_seq_beta.py rename to resolwe_bio/processes/workflows/atac_seq.py index 611ebef01..953b8724a 100644 --- a/resolwe_bio/processes/workflows/atac_seq_beta.py +++ b/resolwe_bio/processes/workflows/atac_seq.py @@ -14,7 +14,7 @@ class WorkflowATACSeq(Process): - """Beta ATAC-Seq workflow. + """ATAC-Seq workflow. This ATAC-seq pipeline follows the official ENCODE DCC pipeline with additional modifications. @@ -31,13 +31,13 @@ class WorkflowATACSeq(Process): Finally, a bigWig file is produced along with the ChipQC and MultiQC reports. """ - slug = "workflow-atac-seq-beta" - name = "ATAC-Seq (Beta)" + slug = "workflow-atac-seq" + name = "ATAC-Seq" requirements = { "expression-engine": "jinja", } data_name = "{{ reads|name|default('?') }}" - version = "0.1.0" + version = "4.0.0" entity = {"type": "sample"} process_type = "data:workflow:atacseq" category = "ATAC-seq" @@ -78,7 +78,7 @@ class DownsamplingOptions: ) two_pass = BooleanField( label="2-pass mode", - default=False, + default=True, disabled="!downsampling_options.downsample_reads", description="Enable two-pass mode when down-sampling. " "Two-pass mode is twice as slow but with much reduced memory.", diff --git a/resolwe_bio/processes/workflows/atac_seq.yml b/resolwe_bio/processes/workflows/atac_seq.yml deleted file mode 100644 index 47efd4096..000000000 --- a/resolwe_bio/processes/workflows/atac_seq.yml +++ /dev/null @@ -1,492 +0,0 @@ -- slug: workflow-atac-seq - name: ATAC-Seq - data_name: "{{ reads|name|default('?') }}" - entity: - type: sample - requirements: - expression-engine: jinja - version: 3.1.1 - type: data:workflow:atacseq - category: ATAC-seq - description: | - This ATAC-seq pipeline closely follows the official ENCODE DCC pipeline. It is comprised of - three steps; alignment, pre-peakcall QC, and calling peaks (with post-peakcall QC). - - First, reads are aligned to a genome using - [Bowtie2](http://bowtie-bio.sourceforge.net/index.shtml) aligner. Next, pre-peakcall QC - metrics are calculated. QC report contains ENCODE 3 proposed QC metrics -- - [NRF](https://www.encodeproject.org/data-standards/terms/), - [PBC bottlenecking coefficients, NSC, and RSC](https://genome.ucsc.edu/ENCODE/qualityMetrics.html#chipSeq). - Finally, the peaks are called using [MACS2](https://github.com/taoliu/MACS/). - The post-peakcall QC report includes additional QC metrics -- number of peaks, - fraction of reads in peaks (FRiP), number of reads in peaks, and if promoter - regions BED file is provided, number of reads in promoter regions, fraction of - reads in promoter regions, number of peaks in promoter regions, and fraction of reads in promoter regions. - input: - - name: reads - label: Select sample(s) - type: data:reads:fastq - - name: genome - label: Genome - type: data:index:bowtie2 - - name: promoter - label: Promoter regions BED file - type: data:bed - required: false - description: | - BED file containing promoter regions (TSS+-1000 bp for example). Needed to get the number - of peaks and reads mapped to promoter regions. - - name: alignment - label: Alignment (Bowtie2) - group: - - name: mode - label: Alignment mode - type: basic:string - description: | - End to end: Bowtie 2 requires that the entire read align from one end to the other, - without any trimming (or "soft clipping") of characters from either end. - Local: Bowtie 2 does not require that the entire read align from one end to the other. - Rather, some characters may be omitted ("soft clipped") from the ends in order to - achieve the greatest possible alignment score. - default: --local - choices: - - label: end to end mode - value: --end-to-end - - label: local - value: --local - - name: speed - label: Speed vs. Sensitivity - type: basic:string - default: --sensitive - choices: - - label: Very fast - value: --very-fast - - label: Fast - value: --fast - - label: Sensitive - value: --sensitive - - label: Very sensitive - value: --very-sensitive - - name: PE_options - label: Paired end alignment options - group: - - name: use_se - label: Map as single-ended (for paired-end reads only) - type: basic:boolean - default: false - description: | - If this option is selected paired-end reads will be mapped as single-ended and - other paired-end options are ignored. - - name: discordantly - label: Report discordantly matched read - type: basic:boolean - default: true - description: | - If both mates have unique alignments, but the alignments do not match paired-end - expectations (orientation and relative distance) then alignment will be reported. - Useful for detecting structural variations. - - name: rep_se - label: Report single ended - type: basic:boolean - default: true - description: | - If paired alignment can not be found Bowtie2 tries to find alignments for the - individual mates. - - name: minins - label: Minimal distance - type: basic:integer - description: | - The minimum fragment length for valid paired-end alignments. 0 imposes no minimum. - default: 0 - - name: maxins - label: Maximal distance - type: basic:integer - description: | - The maximum fragment length for valid paired-end alignments. - default: 2000 - - name: start_trimming - label: Initial trimming - group: - - name: trim_5 - label: Bases to trim from 5' - type: basic:integer - description: | - Number of bases to trim from from 5' (left) end of each read before alignment. - default: 0 - - name: trim_3 - label: Bases to trim from 3' - type: basic:integer - description: | - Number of bases to trim from from 3' (right) end of each read before alignment - default: 0 - - name: trimming - label: Iterative trimming - description: | - Trimming options only apply to single-end reads and paired-end reads mapped as - single-ended. - group: - - name: trim_iter - label: Iterations - type: basic:integer - description: | - Number of iterations. - default: 0 - - name: trim_nucl - label: Bases to trim - type: basic:integer - description: | - Number of bases to trim from 3' end in each iteration. - default: 2 - - name: reporting - label: Reporting - group: - - name: rep_mode - label: Report mode - type: basic:string - description: | - Default mode: search for multiple alignments, report the best one; - -k mode: search for one or more alignments, report each; - -a mode: search for and report all alignments - default: def - choices: - - label: Default mode - value: def - - label: -k mode - value: k - - label: -a mode (very slow) - value: a - - name: k_reports - label: Number of reports (for -k mode only) - type: basic:integer - description: | - Searches for at most X distinct, valid alignments for each read. The search - terminates when it can't find more distinct valid alignments, or when it finds X, - whichever happens first. - default: 5 - - name: prepeakqc_settings - label: Pre-peak QC settings - group: - - name: q_threshold - label: Quality filtering threshold - type: basic:integer - default: 30 - - name: n_sub - label: Number of reads to subsample - type: basic:integer - default: 25000000 - - name: tn5 - label: Tn5 shifting - type: basic:boolean - default: true - description: | - Tn5 transposon shifting. Shift reads on "+" strand by 4 bp - and reads on "-" strand by 5 bp. - - name: shift - label: User-defined cross-correlation peak strandshift - type: basic:integer - default: 0 - description: | - If defined, SPP tool will not try to estimate fragment length but will use the given value - as fragment length. - - name: settings - label: MACS2 settings - group: - - name: tagalign - label: Use tagAlign files - type: basic:boolean - default: true - description: | - Use filtered tagAlign files as case (treatment) and control - (background) samples. If extsize parameter is not set, run MACS - using input's estimated fragment length. - - name: duplicates - label: Number of duplicates - type: basic:string - required: false - hidden: "settings.tagalign" - choices: - - label: "1" - value: "1" - - label: "auto" - value: "auto" - - label: "all" - value: "all" - description: | - It controls the MACS behavior towards duplicate tags at the exact same location -- the - same coordination and the same strand. The 'auto' option makes MACS calculate the - maximum tags at the exact same location based on binomal distribution using 1e-5 as - pvalue cutoff and the 'all' option keeps all the tags. If an integer is given, at most - this number of tags will be kept at the same location. The default is to keep one tag - at the same location. - - name: duplicates_prepeak - label: Number of duplicates - type: basic:string - hidden: "!settings.tagalign" - required: false - default: all - choices: - - label: "1" - value: "1" - - label: "auto" - value: "auto" - - label: "all" - value: "all" - description: | - It controls the MACS behavior towards duplicate tags at the exact same location -- the - same coordination and the same strand. The 'auto' option makes MACS calculate the - maximum tags at the exact same location based on binomal distribution using 1e-5 as - pvalue cutoff and the 'all' option keeps all the tags. If an integer is given, at most - this number of tags will be kept at the same location. The default is to keep one tag - at the same location. - - name: qvalue - label: Q-value cutoff - type: basic:decimal - required: false - disabled: "settings.pvalue && settings.pvalue_prepeak" - description: | - The q-value (minimum FDR) cutoff to call significant regions. Q-values - are calculated from p-values using Benjamini-Hochberg procedure. - - name: pvalue - label: P-value cutoff - type: basic:decimal - required: false - disabled: "settings.qvalue" - hidden: "settings.tagalign" - description: | - The p-value cutoff. If specified, MACS2 will use p-value instead of q-value cutoff. - - name: pvalue_prepeak - label: P-value cutoff - type: basic:decimal - disabled: "settings.qvalue" - hidden: "!settings.tagalign || settings.qvalue" - default: 0.01 - description: | - The p-value cutoff. If specified, MACS2 will use p-value instead of q-value cutoff. - - name: cap_num - label: Cap number of peaks by taking top N peaks - type: basic:integer - default: 300000 - disabled: "settings.broad" - description: | - To keep all peaks set value to 0. - - name: mfold_lower - label: MFOLD range (lower limit) - type: basic:integer - required: false - description: | - This parameter is used to select the regions within MFOLD range of high-confidence - enrichment ratio against background to build model. The regions must be lower than - upper limit, and higher than the lower limit of fold enrichment. DEFAULT:10,30 means - using all regions not too low (>10) and not too high (<30) to build paired-peaks - model. If MACS can not find more than 100 regions to build model, it will use the - --extsize parameter to continue the peak detection ONLY if --fix-bimodal is set. - - name: mfold_upper - label: MFOLD range (upper limit) - type: basic:integer - required: false - description: | - This parameter is used to select the regions within MFOLD range of high-confidence - enrichment ratio against background to build model. The regions must be lower than - upper limit, and higher than the lower limit of fold enrichment. DEFAULT:10,30 means - using all regions not too low (>10) and not too high (<30) to build paired-peaks - model. If MACS can not find more than 100 regions to build model, it will use the - --extsize parameter to continue the peak detection ONLY if --fix-bimodal is set. - - name: slocal - label: Small local region - type: basic:integer - required: false - description: | - Slocal and llocal parameters control which two levels of regions will be checked - around the peak regions to calculate the maximum lambda as local lambda. By default, - MACS considers 1000 bp for small local region (--slocal), and 10000 bp for large local - region (--llocal) which captures the bias from a long range effect like an open - chromatin domain. You can tweak these according to your project. Remember that if the - region is set too small, a sharp spike in the input data may kill the significant - peak. - - name: llocal - label: Large local region - type: basic:integer - required: false - description: | - Slocal and llocal parameters control which two levels of regions will be checked - around the peak regions to calculate the maximum lambda as local lambda. By default, - MACS considers 1000 bp for small local region (--slocal), and 10000 bp for large local - region (--llocal) which captures the bias from a long range effect like an open - chromatin domain. You can tweak these according to your project. Remember that if the - region is set too small, a sharp spike in the input data may kill the significant - peak. - - name: extsize - label: extsize - type: basic:integer - default: 150 - description: | - While '--nomodel' is set, MACS uses this parameter to extend reads in 5'->3' direction - to fix-sized fragments. For example, if the size of binding region for your - transcription factor is 200 bp, and you want to bypass the model building by MACS, - this parameter can be set as 200. This option is only valid when --nomodel is set or - when MACS fails to build model and --fix-bimodal is on. - - name: shift - label: Shift - type: basic:integer - default: -75 - description: | - Note, this is NOT the legacy --shiftsize option which is replaced by --extsize! You - can set an arbitrary shift in bp here. Please Use discretion while setting it other - than default value (0). When --nomodel is set, MACS will use this value to move - cutting ends (5') then apply --extsize from 5' to 3' direction to extend them to - fragments. When this value is negative, ends will be moved toward 3'->5' direction, - otherwise 5'->3' direction. Recommended to keep it as default 0 for ChIP-Seq datasets, - or -1 * half of EXTSIZE together with --extsize option for detecting enriched cutting - loci such as certain DNAseI-Seq datasets. Note, you can't set values other than 0 if - format is BAMPE for paired-end data. Default is 0. - - name: band_width - label: Band width - type: basic:integer - required: false - description: | - The band width which is used to scan the genome ONLY for model building. You can set - this parameter as the sonication fragment size expected from wet experiment. The - previous side effect on the peak detection process has been removed. So this parameter - only affects the model building. - - name: nolambda - label: Use backgroud lambda as local lambda - type: basic:boolean - default: false - description: | - With this flag on, MACS will use the background lambda as local lambda. This means - MACS will not consider the local bias at peak candidate regions. - - name: fix_bimodal - label: Turn on the auto paired-peak model process - type: basic:boolean - default: false - description: | - Turn on the auto paired-peak model process. If it's set, when MACS failed - to build paired model, it will use the nomodel settings, the '--extsize' parameter - to extend each tag. If set, MACS will be terminated if paired-peak model has failed. - - name: nomodel - label: Bypass building the shifting model - type: basic:boolean - hidden: "settings.tagalign" - default: false - description: | - While on, MACS will bypass building the shifting model. - - name: nomodel_prepeak - label: Bypass building the shifting model - type: basic:boolean - default: true - hidden: "!settings.tagalign" - description: | - While on, MACS will bypass building the shifting model. - - name: down_sample - label: Down-sample - type: basic:boolean - default: false - description: | - When set to true, random sampling method will scale down the bigger sample. By default, MACS - uses linear scaling. This option will make the results unstable and irreproducible - since each time, random reads would be selected, especially the numbers (pileup, - pvalue, qvalue) would change. - - name: bedgraph - label: Save fragment pileup and control lambda - type: basic:boolean - default: true - description: | - If this flag is on, MACS will store the fragment pileup, control lambda, -log10pvalue - and -log10qvalue scores in bedGraph files. The bedGraph files will be stored in - current directory named NAME+'_treat_pileup.bdg' for treatment data, - NAME+'_control_lambda.bdg' for local lambda values from control, - NAME+'_treat_pvalue.bdg' for Poisson pvalue scores (in -log10(pvalue) form), and - NAME+'_treat_qvalue.bdg' for q-value scores from Benjamini-Hochberg-Yekutieli - procedure. - - name: spmr - label: Save signal per million reads for fragment pileup profiles - type: basic:boolean - default: true - disabled: "settings.bedgraph === false" - - name: call_summits - label: Call summits - type: basic:boolean - default: true - description: | - MACS will now reanalyze the shape of signal profile (p or q-score depending on cutoff - setting) to deconvolve subpeaks within each peak called from general procedure. It's - highly recommended to detect adjacent binding events. While used, the output subpeaks - of a big peak region will have the same peak boundaries, and different scores and peak - summit positions. - - name: broad - label: Composite broad regions - type: basic:boolean - default: false - disabled: "settings.call_summits === true" - description: | - When this flag is on, MACS will try to composite broad regions in BED12 (a - gene-model-like format) by putting nearby highly enriched regions into a broad region - with loose cutoff. The broad region is controlled by another cutoff through - --broad-cutoff. The maximum length of broad region length is 4 times of d from MACS. - - name: broad_cutoff - label: Broad cutoff - type: basic:decimal - required: false - disabled: "settings.call_summits === true || settings.broad !== true" - description: | - Cutoff for broad region. This option is not available unless --broad is set. If -p is - set, this is a p-value cutoff, otherwise, it's a q-value cutoff. DEFAULT = 0.1 - run: - language: workflow - program: - - id: alignment - run: alignment-bowtie2 - input: - genome: '{{ input.genome }}' - reads: '{{ input.reads }}' - mode: '{{ input.alignment.mode }}' - speed: '{{ input.alignment.speed }}' - PE_options: - use_se: '{{ input.alignment.PE_options.use_se }}' - discordantly: '{{ input.alignment.PE_options.discordantly }}' - rep_se: '{{ input.alignment.PE_options.rep_se }}' - minins: '{{ input.alignment.PE_options.minins }}' - maxins: '{{ input.alignment.PE_options.maxins }}' - start_trimming: - trim_5: '{{ input.alignment.start_trimming.trim_5 }}' - trim_3: '{{ input.alignment.start_trimming.trim_3 }}' - trimming: - trim_iter: '{{ input.alignment.trimming.trim_iter }}' - trim_nucl: '{{ input.alignment.trimming.trim_nucl }}' - reporting: - rep_mode: '{{ input.alignment.reporting.rep_mode }}' - k_reports: '{{ input.alignment.reporting.k_reports }}' - - id: macs2 - run: macs2-callpeak - input: - case: '{{ steps.alignment }}' - promoter: '{{ input.promoter }}' - tagalign: '{{ input.settings.tagalign }}' - prepeakqc_settings: - q_threshold: '{{ input.prepeakqc_settings.q_threshold }}' - n_sub: '{{ input.prepeakqc_settings.n_sub }}' - tn5: '{{ input.prepeakqc_settings.tn5 }}' - shift: '{{ input.prepeakqc_settings.shift }}' - settings: - duplicates: '{{ input.settings.duplicates }}' - duplicates_prepeak: '{{ input.settings.duplicates_prepeak }}' - qvalue: '{{ input.settings.qvalue }}' - pvalue: '{{ input.settings.pvalue }}' - pvalue_prepeak: '{{ input.settings.pvalue_prepeak }}' - cap_num: '{{ input.settings.cap_num }}' - mfold_lower: '{{ input.settings.mfold_lower }}' - mfold_upper: '{{ input.settings.mfold_upper }}' - slocal: '{{ input.settings.slocal }}' - llocal: '{{ input.settings.llocal }}' - extsize: '{{ input.settings.extsize }}' - shift: '{{ input.settings.shift }}' - nolambda: '{{ input.settings.nolambda }}' - fix_bimodal: '{{ input.settings.fix_bimodal }}' - nomodel: '{{ input.settings.nomodel }}' - nomodel_prepeak: '{{ input.settings.nomodel_prepeak }}' - down_sample: '{{ input.settings.down_sample }}' - bedgraph: '{{ input.settings.bedgraph }}' - spmr: '{{ input.settings.spmr }}' - call_summits: '{{ input.settings.call_summits }}' - broad: '{{ input.settings.broad }}' - broad_cutoff: '{{ input.settings.broad_cutoff }}' diff --git a/resolwe_bio/tests/workflows/test_atac_seq.py b/resolwe_bio/tests/workflows/test_atac_seq.py index 2c734e591..4ea6165fb 100644 --- a/resolwe_bio/tests/workflows/test_atac_seq.py +++ b/resolwe_bio/tests/workflows/test_atac_seq.py @@ -27,37 +27,6 @@ def test_atac_seq_workflow(self): self.run_process( "workflow-atac-seq", - { - "reads": reads.id, - "genome": bowtie2_index.id, - "promoter": promoters.id, - }, - ) - for data in Data.objects.all(): - self.assertStatus(data, Data.STATUS_DONE) - - @tag_process("workflow-atac-seq-beta") - def test_beta_atac_seq_workflow(self): - with self.preparation_stage(): - inputs = { - "src": "mm10_chr17_44M-45M.fa.gz", - "species": "Mus musculus", - "build": "mm10", - } - ref_seq = self.run_process("upload-fasta-nucl", inputs) - bowtie2_index = self.run_process("bowtie2-index", {"ref_seq": ref_seq.id}) - reads = self.prepare_paired_reads( - mate1=["atac_R1.fastq.gz"], mate2=["atac_R2.fastq.gz"] - ) - inputs = { - "src": "promoters_mm10_chr17_subregion.bed", - "species": "Mus musculus", - "build": "mm10", - } - promoters = self.run_process("upload-bed", inputs) - - self.run_process( - "workflow-atac-seq-beta", { "reads": reads.id, "genome": bowtie2_index.id,