From d83807f38d391b5bda48d0bfdb6d8c31127d4892 Mon Sep 17 00:00:00 2001 From: Sarah Griffiths Date: Mon, 16 Sep 2024 09:26:12 +0000 Subject: [PATCH] Tidy de analysis --- CHANGELOG.md | 1 + bin/de_analysis.R | 3 +- main.nf | 4 +-- nextflow.config | 13 ++++---- subworkflows/differential_expression.nf | 40 ++++++++++--------------- 5 files changed, 25 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9ea7221..2e6153c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - If required for IGV, reference indexes are output in to a `igv_reference` directory ### Changed - BAMS are output in to a BAMS directory. +- Reconcile with template 5.2.6. ## [v1.3.0] ### Removed diff --git a/bin/de_analysis.R b/bin/de_analysis.R index b76c106..cb31755 100755 --- a/bin/de_analysis.R +++ b/bin/de_analysis.R @@ -9,13 +9,14 @@ min_samps_gene_expr <- as.numeric(args[2]) min_samps_feature_expr <- as.numeric(args[3]) min_gene_expr <- as.numeric(args[4]) min_feature_expr <- as.numeric(args[5]) +sample_sheet <- args[6] cat("Loading counts, conditions and parameters.\n") cts <- as.matrix(read.csv("all_counts.tsv", sep="\t", row.names="Reference", stringsAsFactors=FALSE)) # Set up sample data frame: #changed this to sample_id -coldata <- read.csv("sample_sheet.csv", row.names="alias", sep=",", stringsAsFactors=TRUE) +coldata <- read.csv(sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE) coldata$sample_id <- rownames(coldata) # check if control condition exists, sets as reference diff --git a/main.nf b/main.nf index 78a5795..8bd590d 100644 --- a/main.nf +++ b/main.nf @@ -717,9 +717,7 @@ workflow pipeline { de = differential_expression(transcriptome, input_reads, sample_sheet, gtf) de_report = de.all_de count_transcripts_file = de.count_transcripts - dtu_plots = de.dtu_plots de_outputs = de.de_outputs - counts = de.counts } else{ de_report = OPTIONAL_FILE count_transcripts_file = OPTIONAL_FILE @@ -763,7 +761,7 @@ workflow pipeline { if (params.de_analysis){ de_results = report.concat( - transcriptome, de_outputs.flatten(), counts.flatten(), + transcriptome, de_outputs.flatten(), makeReport.out.results_dge, makeReport.out.tpm, makeReport.out.filtered, makeReport.out.unfiltered, makeReport.out.gene_counts) diff --git a/nextflow.config b/nextflow.config index 136cacd..c13c030 100644 --- a/nextflow.config +++ b/nextflow.config @@ -20,7 +20,6 @@ params { threads = 4 // Thresholds for viewing isoforms in report table isoform_table_nrows = 5000 - out_dir = "output" sample = null @@ -42,7 +41,7 @@ params { pychopper_opts = null pychopper_backend = "edlib" cdna_kit = "SQK-PCS109" - + // Extra option passed to minimap2 when generating index minimap2_index_opts = "-k14" @@ -79,9 +78,9 @@ params { de_analysis = false ref_transcriptome = null min_samps_gene_expr = 3 - min_samps_feature_expr = 1 + min_samps_feature_expr = 1 min_gene_expr = 10 - min_feature_expr = 3 + min_feature_expr = 3 wf { @@ -95,8 +94,8 @@ params { "--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'", ] agent = null - container_sha = "shafb1e2372e1535f0b42891ed2c68ffdac2ca1d658" - common_sha = "shae58638742cf84dbeeec683ba24bcdee67f64b986" + container_sha = "shad8671ea3a8ed52f2c0f40355e8eb5c6f00d2cbda" + common_sha = "shad28e55140f75a68f59bbecc74e880aeab16ab158" } } @@ -165,7 +164,7 @@ profiles { } withLabel:wf_common { container = "${params.aws_image_prefix}-wf-common:${params.wf.common_sha}" - } + } shell = ['/bin/bash', '-euo', 'pipefail'] } } diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index ea40be5..f21d1b5 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -69,16 +69,21 @@ process deAnalysis { output: path "de_analysis/results_dtu_stageR.tsv", emit: stageR path "merged/filtered_transcript_counts_with_genes.tsv", emit: flt_counts - path "de_analysis/unfiltered_transcript_counts_with_genes.tsv", emit: unflt_counts path "merged/all_gene_counts.tsv", emit: gene_counts + path "de_analysis/unfiltered_transcript_counts_with_genes.tsv", emit: unflt_counts path "de_analysis/results_dge.tsv", emit: dge path "de_analysis/results_dexseq.tsv", emit: dexseq - path "de_analysis", emit: de_analysis + path "de_analysis/results_dge.pdf", emit: dge_pdf + path "de_analysis/results_dge.tsv", emit: dge_tsv + path "de_analysis/results_dtu_gene.tsv", emit: dtu_gene + path "de_analysis/results_dtu_transcript.tsv", emit: dtu_transcript + path "de_analysis/results_dtu_stageR.tsv", emit: dtu_stageR + path "de_analysis/results_dtu.pdf", emit: dtu_pdf path "de_analysis/cpm_gene_counts.tsv", emit: cpm """ mkdir merged mkdir de_analysis - de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr + de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr "sample_sheet.csv" """ } @@ -91,23 +96,10 @@ process plotResults { path "filtered_transcript_counts_with_genes.tsv" path "results_dtu_stageR.tsv" path "sample_sheet.tsv" - path de_analysis output: - path "de_analysis/dtu_plots.pdf", emit: dtu_plots - path "sample_sheet.tsv", emit: sample_sheet_csv - // Output all DE files for use in report process - path "de_analysis/*", emit: stageR - // Output selected DE files to be output in out_dir - path "de_analysis/results_dge.pdf", emit: dge_pdf - path "de_analysis/results_dge.tsv", emit: dge_tsv - path "de_analysis/results_dtu_gene.tsv", emit: dtu_gene - path "de_analysis/results_dtu_transcript.tsv", emit: dtu_transcript - path "de_analysis/results_dtu_stageR.tsv", emit: dtu_stageR - path "de_analysis/results_dtu.pdf", emit: dtu_pdf + path "dtu_plots.pdf", emit: dtu_plots """ plot_dtu_results.R - # output plots to common analysis output directory - cp dtu_plots.pdf de_analysis/dtu_plots.pdf """ } @@ -159,6 +151,7 @@ workflow differential_expression { sample_sheet ref_annotation main: + sample_sheet = Channel.fromPath(sample_sheet) checkSampleSheetCondition(sample_sheet) t_index = build_minimap_index_transcriptome(ref_transcriptome) mapped = map_transcriptome(full_len_reads.combine(t_index) @@ -167,19 +160,16 @@ workflow differential_expression { merged = mergeCounts(count_transcripts.out.counts.collect()) merged_TPM = mergeTPM(count_transcripts.out.counts.collect()) analysis = deAnalysis(sample_sheet, merged, ref_annotation) - plotResults(analysis.flt_counts, analysis.stageR, sample_sheet, analysis.de_analysis) + plotResults(analysis.flt_counts, analysis.stageR, sample_sheet) // Concat files required for making the report de_report = analysis.flt_counts.concat(analysis.gene_counts, analysis.dge, analysis.dexseq, - analysis.stageR, plotResults.out.sample_sheet_csv, merged, ref_annotation, merged_TPM, analysis.unflt_counts).collect() - // Concat files required to be output to user - de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, plotResults.out.dge_pdf, plotResults.out.dge_tsv, - plotResults.out.dtu_gene, plotResults.out.dtu_transcript, plotResults.out.dtu_stageR, plotResults.out.dtu_pdf).collect() + analysis.stageR, sample_sheet, merged, ref_annotation, merged_TPM, analysis.unflt_counts).collect() + // Concat files required to be output to user without any changes + de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, analysis.dge_pdf, analysis.dge_tsv, + analysis.dtu_gene, analysis.dtu_transcript, analysis.dtu_stageR, analysis.dtu_pdf, analysis.flt_counts, analysis.gene_counts, merged_TPM).collect() count_transcripts_file = count_transcripts.out.seqkit_stats.collect() - all_counts = merged_TPM.concat(analysis.flt_counts, analysis.gene_counts) emit: all_de = de_report count_transcripts = count_transcripts_file - dtu_plots = plotResults.out.dtu_plots de_outputs = de_outputs_concat - counts = all_counts }