From 8053448d613bad1da9ccb297434b87f4f26ff715 Mon Sep 17 00:00:00 2001 From: Sarah Griffiths Date: Tue, 4 Apr 2023 11:54:35 +0000 Subject: [PATCH] CW-1882 CSV --- .pre-commit-config.yaml | 2 +- CHANGELOG.md | 5 + bin/workflow_glue/report.py | 17 +- lib/fastqingress.nf | 826 ++++++++++++------------ main.nf | 131 ++-- nextflow.config | 2 +- nextflow_schema.json | 4 +- subworkflows/differential_expression.nf | 8 +- 8 files changed, 526 insertions(+), 469 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f5d1036..1d1d82d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -18,7 +18,7 @@ repos: additional_dependencies: - epi2melabs - repo: https://github.com/pycqa/flake8 - rev: 3.7.9 + rev: 5.0.4 hooks: - id: flake8 pass_filenames: false diff --git a/CHANGELOG.md b/CHANGELOG.md index 3aec571..544f11d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,11 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [unreleased] +### Updated +- Condition sheet parameter description fixed to CSV +- Update fastqingress + ## [v0.1.9] ### Updated - Simplify JAFFAL docs diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index ee75d2d..b04ae6c 100755 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -31,7 +31,7 @@ def argparser(): """Argument parser for entrypoint.""" parser = wf_parser("report") parser.add_argument("--report", help="Report output file") - parser.add_argument("--summaries", nargs='+', help="Read summary file.") + parser.add_argument("--stats", help="Read stats file.") parser.add_argument( "--versions", required=True, help="directory containing CSVs containing name,version.") @@ -821,17 +821,16 @@ def load_sample_data(files, sample_ids, read_func=None): return df_ -def seq_stats_tabs(report, sample_ids, stats): +def seq_stats_tabs(report, stats): """Make tabs of sequence summaries by sample.""" tabs = [] - for id_, summ in sorted(zip(sample_ids, stats)): - df_sum = pd.read_csv(summ, index_col=False, sep='\t') - rlp = read_length_plot(df_sum) - rqp = read_quality_plot(df_sum) + df_all = pd.read_csv(stats, sep="\t") + for sample_id, df_sample in df_all.groupby('sample_name'): + rlp = read_length_plot(df_sample) + rqp = read_quality_plot(df_sample) grid = gridplot( [rlp, rqp], ncols=2, sizing_mode="stretch_width") - - tabs.append(Panel(child=grid, title=id_)) + tabs.append(Panel(child=grid, title=sample_id)) section = report.add_section() section.markdown(""" ### Sequence summaries""") @@ -908,7 +907,7 @@ def main(args): revision=args.revision, commit=args.commit) # QC - seq_stats_tabs(report, args.sample_ids, args.summaries) + seq_stats_tabs(report, args.stats) if args.alignment_stats is not None: df_aln_stats = load_sample_data(args.alignment_stats, sample_ids) diff --git a/lib/fastqingress.nf b/lib/fastqingress.nf index 13988c2..e24048e 100644 --- a/lib/fastqingress.nf +++ b/lib/fastqingress.nf @@ -1,411 +1,437 @@ +import java.nio.file.NoSuchFileException + import ArgumentParser -// Downstream tooling assumes FASTQ files are nicely organised into directories. -// In the case where a single FASTQ file has been input and the parent directory -// contains other valid FASTQ, we will create a directory in the work area to -// hold it instead. We stage the file in with `copy` (rather than `link`) -// to ensure that when the new dir is mounted to containers downstream it does -// not contain a symlink that cannot be read. -// See CW-1154 -process isolateSingleFile { - label params.process_label - stageInMode 'copy' - cpus 1 - input: - file reads - output: - path "$reads.simpleName" - script: - def name = reads.simpleName - """ - mkdir $name - mv $reads $name - """ -} - - -process checkSampleSheet { - label params.process_label - cpus 1 - input: - file "sample_sheet.txt" - output: - file "samples.txt" - """ - workflow-glue check_sample_sheet sample_sheet.txt samples.txt - """ -} +EXTENSIONS = ["fastq", "fastq.gz", "fq", "fq.gz"] /** - * Compare number of samples in samplesheet - * with the number of barcoded dirs found and - * print warnings + * Take a map of input arguments, find valid inputs, and return a channel + * with elements of `[metamap, seqs.fastq.gz, path-to-fastcat-stats]`. + * The last item is `null` if `fastcat` was not run. It is only run on directories + * containing more than one FASTQ file or when `fastcat_stats: true`. * - * - * @param number of samples in sample sheet - * @param number of barcoded directories - * @return null - */ - -def compareSampleSheetFastq(int sample_sheet_count, int valid_dir_count) -{ - - if (sample_sheet_count != valid_dir_count) { - log.warn "The number of samplesheet entries ({}) does not match the number of barcoded directories ({})", sample_sheet_count, valid_dir_count - } - -} - -/** - * Take an input file and sample name to return a channel with - * a single named sample. If the input file is in a directory with other valid - * input files (or other directories containing valid files), a copy of it will - * be made to the working directory using the isolateSingleFile process. - * - * - * @param input_file Single fastq file - * @param sample_name Name to give the sample - * @return Channel of tuples (path, map(sample_id, type, barcode)) - */ -def handle_single_file(input_file, sample_name) -{ - singleFile = Channel.fromPath(input_file) - ArrayList valid_files_in_dir = find_fastq(input_file.parent, true) - if (valid_files_in_dir.size() == 1) { - // Avoid a stageInMode copy if the parent directory contains only one valid FASTQ anyway - return singleFile.map { it -> tuple(it.parent, create_metamap([sample_id:sample_name ?: it.simpleName])) } - } - else { - // Isolate the file via copy with isolateSingleFile - sample = isolateSingleFile(singleFile) - return sample.map { it -> tuple(it, create_metamap([sample_id:sample_name ?: it.simpleName])) } - } -} - - -/** - * Find fastq data using various globs. Wrapper around Nextflow `file` - * method. - * - * @param pattern file object corresponding to top level input folder. - * @param search_subdirs boolean flag to search subdirectories of pattern - * @return list of files. - */ - -def find_fastq(pattern, search_subdirs) -{ - ArrayList files = [] - ArrayList extensions = ["fastq", "fastq.gz", "fq", "fq.gz"] - for (ext in extensions) { - if (search_subdirs) { - files += file(pattern.resolve("**.${ext}"), type: 'file') - } - else { - files += file(pattern.resolve("*.${ext}"), type: 'file') - } - } - return files -} - - -/** - * Take an input directory return the barcode and non barcode - * sub directories contained within. - * - * - * @param input_directory Top level input folder to locate sub directories - * @param unclassified Keep unclassified directory - * - * @return A list containing sublists of barcode and non_barcode sub directories - */ -def get_subdirectories(input_directory, unclassified) -{ - barcode_dirs = file(input_directory.resolve("barcode*"), type: 'dir', maxdepth: 1) - all_dirs = file(input_directory.resolve("*"), type: 'dir', maxdepth: 1) - if (!unclassified) { - all_dirs.removeIf(it -> it.SimpleName.toLowerCase() == "unclassified") - } - non_barcoded = (all_dirs + barcode_dirs) - all_dirs.intersect(barcode_dirs) - - return [barcode_dirs, non_barcoded] -} - - -/** - * Load a sample sheet into a Nextflow channel to map barcodes - * to sample names. - * - * @param samples CSV file according to MinKNOW sample sheet specification - * @return A Nextflow Channel of tuples (barcode, sample name, sample type) - */ -def get_sample_sheet(sample_sheet) -{ - log.info "Checking sample sheet." - sample_sheet = file(sample_sheet); - is_file = sample_sheet.isFile() - - if (!is_file) { - log.error "`--samples` is not a file." - exit 1 - } - - return checkSampleSheet(sample_sheet) - .splitCsv(header: true) - .map { row -> tuple( - row.barcode, - row.sample_id, - row.type ? row.type : 'test_sample') - } -} - - -/** - * Take a list of input directories and return directories which are - * valid, i.e. contains only .fastq(.gz) files. - * - * - * @param input_dirs List of barcoded directories (barcodeXX,mydir...) - * @return List of valid directories - */ -def get_valid_directories(input_dirs) -{ - valid_dirs = [] - no_fastq_dirs = [] - invalid_files_dirs = [] - for (d in input_dirs) { - valid = true - fastq = find_fastq(d, false) - all_files = file(d.resolve("*"), type: 'file', maxdepth: 1) - non_fastq = ( all_files + fastq ) - all_files.intersect(fastq) - - if (non_fastq) { - valid = false - invalid_files_dirs << d - } - if (!fastq) { - valid = false - no_fastq_dirs << d - } - if (valid) { - valid_dirs << d - } - } - if (valid_dirs.size() == 0) { - log.error "None of the directories given contain .fastq(.gz) files." - exit 1 - } - if (no_fastq_dirs.size() > 0) { - log.warn "Excluding directories not containing .fastq(.gz) files:" - for (d in no_fastq_dirs) { - log.warn " - ${d}" - } - } - if (invalid_files_dirs.size() > 0) { - log.warn "Excluding directories containing non .fastq(.gz) files:" - for (d in invalid_files_dirs) { - log.warn " - ${d}" - } - } - return valid_dirs -} - - -/** - * Take an input directory and sample name to return a channel - * with a single named sample. - * - * - * @param input_directory Directory of fastq files - * @param sample_name Name to give the sample - * @return Channel of tuples (path, map(sample_id, type, barcode)) - */ -def handle_flat_dir(input_directory, sample_name) -{ - valid_dirs= get_valid_directories([ file(input_directory) ]) - return Channel.fromPath(valid_dirs) - .map { it -> tuple(it, create_metamap([sample_id:sample_name ?: it.baseName])) } - -} - - -/** - * Take a list of barcode directories and a sample sheet to return - * a channel of named samples. - * - * - * @param barcoded_dirs List of barcoded directories (barcodeXX,...) - * @param sample_sheet List of tuples mapping barcode to sample name - * or a simple string for non-multiplexed data. - * @param min_barcode Minimum barcode to accept. - * @param max_barcode Maximum (inclusive) barcode to accept. - * @return Channel of tuples (path, map(sample_id, type, barcode)) - */ -def handle_barcoded_dirs(barcoded_dirs, sample_sheet, min_barcode, max_barcode) -{ - valid_dirs = get_valid_directories(barcoded_dirs) - // link sample names to barcode through sample sheet - if (!sample_sheet) { - sample_sheet = Channel - .fromPath(valid_dirs) - .filter(~/.*barcode[0-9]{1,3}$/) // up to 192 - .filter { barcode_in_range(it, min_barcode, max_barcode) } - .map { path -> tuple(path.baseName, path.baseName, 'test_sample')} - } else { - - // return warning if there is a discrepancy between the samplesheet and barcode dirs - - // unclassfied will never be in the sample_sheet so remove - non_unclassified = valid_dirs - non_unclassified -= 'unclassified' - - barcode_dirs_found = non_unclassified.size() - - int count = 0 - - // We do this instead of .count() because valid_dirs is a list and - // sample_sheet is a channel - the channel is only populated after - // checkSampleSheet is complete and so if you compare without - // waiting for that then the comparisson fails - sample_sheet_entries = sample_sheet.subscribe onNext: { count++ }, onComplete: { compareSampleSheetFastq(count,barcode_dirs_found) } - - } - - return Channel - .fromPath(valid_dirs) - .filter(~/.*barcode[0-9]{1,3}$/) // up to 192 - .filter { barcode_in_range(it, min_barcode, max_barcode) } - .map { path -> tuple(path.baseName, path) } - .join(sample_sheet) - .map { barcode, path, sample, type -> tuple(path, create_metamap([sample_id:sample, type:type, barcode:barcode])) } -} - - -/** - * Determine if a barcode path is within a required numeric range - * - * @param path barcoded directory (barcodeXX). - * @param min_barcode Minimum barcode to accept. - * @param max_barcode Maximum (inclusive) barcode to accept. - */ -def barcode_in_range(path, min_barcode, max_barcode) -{ - pattern = ~/barcode(\d+)/ - matcher = "${path}" =~ pattern - def value = null - try{ - value = matcher[0][1].toInteger() - }catch(ArrayIndexOutOfBoundsException ex){ - print("${path} is not a barcoded directory") - } - valid = ((value >= min_barcode) && (value <= max_barcode)) - return valid -} - - -/** - * Take a list of non-barcode directories to return a channel - * of named samples. Samples are named by directory baseName. - * - * - * @param non_barcoded_dirs List of directories (mydir,...) - * @return Channel of tuples (path, map(sample_id, type, barcode)) - */ -def handle_non_barcoded_dirs(non_barcoded_dirs) -{ - valid_dirs = get_valid_directories(non_barcoded_dirs) - return Channel.fromPath(valid_dirs) - .map { path -> tuple(path, create_metamap([sample_id:path.baseName])) } -} - -def create_metamap(Map arguments) { - def parser = new ArgumentParser( - args:["sample_id"], - kwargs:[ - "type": "test_sample", - "barcode": null, - ], - name:"create_metamap", - ) - return parser.parse_args(arguments) -} - -/** - * Take an input (file or directory) and return a channel of - * named samples. - * - * @param input Top level input file or folder to locate fastq data. - * @param sample string to name single sample data. - * @param sample_sheet Path to sample sheet CSV file. - * @param min_barcode Minimum barcode to accept. - * @param max_barcode Maximum (inclusive) barcode to accept. - * @param unclassified Keep unclassified reads. - * - * @return Channel of tuples (path, map(sample_id, type, barcode)) + * @param arguments: map with arguments containing + * - "input": path to either: (i) input FASTQ file, (ii) top-level directory containing + * FASTQ files, (iii) directory containing sub-directories which contain FASTQ + * files + * - "sample": string to name single sample + * - "sample_sheet": path to CSV sample sheet + * - "analyse_unclassified": boolean whether to keep unclassified reads + * - "fastcat_stats": boolean whether to write the `fastcat` stats + * @return Channel of `[Map(alias, barcode, type, ...), Path, Path|null]`. + * The first element is a map with metadata, the second is the path to the + * `.fastq.gz` file with the (potentially concatenated) sequences and the third is + * the path to the directory with the fastcat statistics (or `null` if `fastcat` + * wasn't run). */ def fastq_ingress(Map arguments) { - def parser = new ArgumentParser( - args:["input"], - kwargs:[ - "sample":null, "sample_sheet":null, - "min_barcode":0, "max_barcode":Integer.MAX_VALUE, - "unclassified":false], - name:"fastq_ingress") - Map margs = parser.parse_args(arguments) - - - log.info "Checking fastq input." - input = file(margs.input) - - // Handle file input - if (input.isFile()) { - // Assume sample is a string at this point - log.info "Single file input detected." - if (margs.sample_sheet) { - log.warn "Warning: `--sample_sheet` given but single file input found. Ignoring." - } - return handle_single_file(input, margs.sample) + // check arguments + Map margs = parse_arguments(arguments) + // define the channel for holding the inputs [metamap, input_path]. It will be + // either filled by `watchPath` (only emitting files) or by the data of the three + // input types (single file or dir with fastq or subdirs with fastq). + def ch_input + // handle `watchPath` case + if (margs["watch_path"]) { + ch_input = watch_path(margs) + } else { + // create a channel with the inputs (single file / dir with fastq / subdirs + // with fastq) + ch_input = get_valid_inputs(margs) } - - // Handle directory input - if (input.isDirectory()) { - // Get barcoded and non barcoded subdirectories - (barcoded, non_barcoded) = get_subdirectories(input, margs.unclassified) - - // Case 03: If no subdirectories, handle the single dir - if (!barcoded && !non_barcoded) { - log.info "Single directory input detected." - if (margs.sample_sheet) { - log.warn "`--sample_sheet` given but single non-barcode directory found. Ignoring." - } - return handle_flat_dir(input, margs.sample) - } - - if (margs.sample) { - log.warn "`--sample` given but multiple directories found, ignoring." - } - - // Case 01, 02, 04: Handle barcoded and non_barcoded dirs - // Handle barcoded folders - barcoded_samples = Channel.empty() - if (barcoded) { - log.info "Barcoded directories detected." - sample_sheet = null - if (margs.sample_sheet) { - sample_sheet = get_sample_sheet(margs.sample_sheet) - } - barcoded_samples = handle_barcoded_dirs(barcoded, sample_sheet, margs.min_barcode, margs.max_barcode) - } - - non_barcoded_samples = Channel.empty() - if (non_barcoded) { - log.info "Non barcoded directories detected." - if (!barcoded && margs.sample_sheet) { - log.warn "Warning: `--sample_sheet` given but no barcode directories found." - } - non_barcoded_samples = handle_non_barcoded_dirs(non_barcoded) - } - - return barcoded_samples.mix(non_barcoded_samples) + // `ch_input` might contain elements of `[metamap, null]` if there were entries in + // the sample sheet for which no FASTQ files were found. We put these into an extra + // channel and combine with the result channel before returning. + ch_input = ch_input.branch { meta, path -> + reads_found: path as boolean + no_reads_found: true } + def ch_result + if (margs.fastcat_stats) { + // run fastcat regardless of input type + ch_result = fastcat(ch_input.reads_found, margs["fastcat_extra_args"]) + } else { + // the fastcat stats were not requested --> run fastcat only on directories with + // more than one FASTQ file (and not on single files or directories with a + // single file) + def ch_branched = ch_input.reads_found.map {meta, path -> + // find directories with only a single FASTQ file and "unwrap" the file + if (path.isDirectory()) { + List fq_files = get_fq_files_in_dir(path) + if (fq_files.size() == 1) { + path = fq_files[0] + } + } + [meta, path] + } .branch { meta, path -> + // now there can only be two cases: + // (i) single FASTQ file (pass to `move_or_compress` later) + // (ii) dir with multiple fastq files (pass to `fastcat` later) + single_file: path.isFile() + dir_with_fastq_files: true + } + // call the respective processes on both branches and return + ch_result = fastcat( + ch_branched.dir_with_fastq_files, margs["fastcat_extra_args"] + ).concat( + ch_branched.single_file | move_or_compress | map { + meta, path -> [meta, path, null] + } + ) + } + return ch_result.concat(ch_input.no_reads_found.map { [*it, null] }) +} + + +/** + * Run `watchPath` on the input directory and return a channel [metamap, path-to-fastq]. + * The meta data is taken from the sample sheet in case one was provided. Otherwise it + * only contains the `alias` (either `margs["sample"]` or the name of the parent + * directory of the file). + * + * @param margs: map with parsed input arguments + * @return: Channel of [metamap, path-to-fastq] + */ +def watch_path(Map margs) { + // we have two cases to consider: (i) files being generated in the top-level + // directory and (ii) files being generated in sub-directories. If we find files of + // both kinds, throw an error. + Path input + try { + input = file(margs.input, checkIfExists: true) + } catch (NoSuchFileException e) { + error "Input path $margs.input does not exist." + } + if (input.isFile()) { + error "Input ($input) must be a directory when using `watch_path`." + } + // get existing FASTQ files first (look for relevant files in the top-level dir and + // all sub-dirs) + def ch_existing_input = Channel.fromPath(input) + | concat(Channel.fromPath("$input/*", type: 'dir')) + | map { get_fq_files_in_dir(it) } + | flatten + // now get channel with files found by `watchPath` + def ch_watched = Channel.watchPath("$input/**").until { it.name.startsWith('STOP') } + // only keep FASTQ files + | filter { + for (ext in EXTENSIONS) { + if (it.name.endsWith(ext)) return true + } + return false + } + // merge the channels + ch_watched = ch_existing_input | concat(ch_watched) + // check if input is as expected; start by throwing an error when finding files in + // top-level dir and sub-directories + String prev_input_type + ch_watched + | map { + String input_type = (it.parent == input) ? "top-level" : "sub-dir" + if (prev_input_type && (input_type != prev_input_type)) { + error "`watchPath` found FASTQ files in the top-level directory " + + "as well as in sub-directories." + } + // if file is in a sub-dir, make sure it's not a sub-sub-dir + if ((input_type == "sub-dir") && (it.parent.parent != input)) { + error "`watchPath` found a FASTQ file more than one level of " + + "sub-directories deep ('$it')." + } + // we also don't want files in the top-level dir when we got a sample sheet + if ((input_type == "top-level") && margs["sample_sheet"]) { + error "`watchPath` found files in top-level directory even though a " + + "sample sheet was provided ('${margs["sample_sheet"]}')." + } + prev_input_type = input_type + } + if (margs.sample_sheet) { + // add metadata from sample sheet (we can't use join here since it does not work + // with repeated keys; we therefore need to transform the sample sheet data into + // a map with the barcodes as keys) + def ch_sample_sheet = get_sample_sheet(file(margs.sample_sheet)) + | collect + | map { it.collectEntries { [(it["barcode"]): it] } } + // now we can use this channel to annotate all files with the corresponding info + // from the sample sheet + ch_watched = ch_watched + | combine(ch_sample_sheet) + | map { file_path, sample_sheet_map -> + String barcode = file_path.parent.name + Map meta = sample_sheet_map[barcode] + // throw error if the barcode was not in the sample sheet + if (!meta) { + error "Sub-directory $barcode was not found in the sample sheet." + } + [meta, file_path] + } + } else { + ch_watched = ch_watched + | map { + // This file could be in the top-level dir or a sub-dir. In the first case + // check if a sample name was provided. In the second case, the alias is + // always the name of the sub-dir. + String alias + if (it.parent == input) { + // top-level dir + alias = margs["sample"] ?: it.parent.name + } else { + // sub-dir + alias = it.parent.name + } + [create_metamap([alias: alias]), it] + } + } + return ch_watched +} + + +process move_or_compress { + label params.process_label + cpus params.threads + input: + tuple val(meta), path(input) + output: + tuple val(meta), path("seqs.fastq.gz") + script: + String out = "seqs.fastq.gz" + if (input.name.endsWith('.gz')) { + // we need to take into account that the file could already be named + // "seqs.fastq.gz" in which case `mv` would fail + """ + [ "$input" == "$out" ] || mv $input $out + """ + } else { + """ + cat $input | bgzip -@ $task.cpus > $out + """ + } +} + + +process fastcat { + label params.process_label + cpus params.threads + input: + tuple val(meta), path(input) + val extra_args + output: + tuple val(meta), path("seqs.fastq.gz"), path("fastcat_stats") + script: + String out = "seqs.fastq.gz" + String fastcat_stats_outdir = "fastcat_stats" + """ + mkdir $fastcat_stats_outdir + fastcat \ + -s ${meta["alias"]} \ + -r $fastcat_stats_outdir/per-read-stats.tsv \ + -f $fastcat_stats_outdir/per-file-stats.tsv \ + $extra_args \ + $input \ + | bgzip -@ $task.cpus > $out + """ +} + + +/** + * Parse input arguments for `fastq_ingress`. + * + * @param arguments: map with input arguments (see `fastq_ingress` for details) + * @return: map of parsed arguments + */ +Map parse_arguments(Map arguments) { + ArgumentParser parser = new ArgumentParser( + args:["input"], + kwargs:["sample": null, + "sample_sheet": null, + "analyse_unclassified": false, + "fastcat_stats": false, + "fastcat_extra_args": "", + "watch_path": false], + name: "fastq_ingress") + return parser.parse_args(arguments) +} + + +/** + * Find valid inputs based on the input type. + * + * @param margs: parsed arguments (see `fastq_ingress` for details) + * @return: channel of `[metamap, input-path]`; `input-path` can be the path to + * a single FASTQ file or to a directory containing FASTQ files + */ +def get_valid_inputs(Map margs){ + log.info "Checking fastq input." + Path input + try { + input = file(margs.input, checkIfExists: true) + } catch (NoSuchFileException e) { + error "Input path $margs.input does not exist." + } + // declare resulting input channel and other variables needed in the outer scope + def ch_input + ArrayList sub_dirs_with_fastq_files + // handle case of `input` being a single file + if (input.isFile()) { + // the `fastcat` process can deal with directories or single file inputs + ch_input = Channel.of( + [create_metamap([alias: margs["sample"] ?: input.simpleName]), input]) + } else if (input.isDirectory()) { + // input is a directory --> we accept two cases: (i) a top-level directory with + // fastq files and no sub-directories or (ii) a directory with one layer of + // sub-directories containing fastq files + boolean dir_has_fastq_files = get_fq_files_in_dir(input) + // find potential sub-directories (and sub-dirs with FASTQ files; note that + // these lists can be empty) + ArrayList sub_dirs = file(input.resolve('*'), type: "dir") + sub_dirs_with_fastq_files = sub_dirs.findAll { get_fq_files_in_dir(it) } + // deal with first case (top-lvl dir with FASTQ files and no sub-directories + // containing FASTQ files) + if (dir_has_fastq_files) { + if (sub_dirs_with_fastq_files) { + error "Input directory '$input' cannot contain FASTQ " + + "files and sub-directories with FASTQ files." + } + ch_input = Channel.of( + [create_metamap([alias: margs["sample"] ?: input.baseName]), input]) + } else { + // deal with the second case (sub-directories with fastq data) --> first + // check whether we actually found sub-directories + if (!sub_dirs_with_fastq_files) { + error "Input directory '$input' must contain either FASTQ files " + + "or sub-directories containing FASTQ files." + } + // make sure that there are no sub-sub-directories with FASTQ files and that + // the sub-directories actually contain fastq files) + if (sub_dirs.any { + ArrayList subsubdirs = file(it.resolve('*'), type: "dir") + subsubdirs.any { get_fq_files_in_dir(it) } + }) { + error "Input directory '$input' cannot contain more " + + "than one level of sub-directories with FASTQ files." + } + // remove directories called 'unclassified' unless otherwise specified + if (!margs.analyse_unclassified) { + sub_dirs_with_fastq_files = sub_dirs_with_fastq_files.findAll { + it.baseName != "unclassified" + } + } + // filter based on sample sheet in case one was provided + if (margs.sample_sheet) { + // get channel of entries in the sample sheet + def ch_sample_sheet = get_sample_sheet(file(margs.sample_sheet)) + // get the union of both channels (missing values will be replaced with + // `null`) + def ch_union = Channel.fromPath(sub_dirs_with_fastq_files).map { + [it.baseName, it] + }.join(ch_sample_sheet.map{[it.barcode, it]}, remainder: true) + // after joining the channels, there are three possible cases: + // (i) valid input path and sample sheet entry are both present + // (ii) there is a sample sheet entry but no corresponding input dir + // --> we'll emit `[metamap-from-sample-sheet-entry, null]` + // (iii) there is a valid path, but the sample sheet entry is missing + // --> drop this entry and print a warning to the log + ch_input = ch_union.map {barcode, path, sample_sheet_entry -> + if (sample_sheet_entry) { + [create_metamap(sample_sheet_entry), path] + } else { + log.warn "Input directory '$barcode' was found, but sample " + + "sheet '$margs.sample_sheet' has no such entry." + } + } + } else { + ch_input = Channel.fromPath(sub_dirs_with_fastq_files).map { + [create_metamap([alias: it.baseName, barcode: it.baseName]), it] + } + } + } + } else { + error "Input $input appears to be neither a file nor a directory." + } + // a sample sheet only makes sense in the case of a directory with + // sub-directories + if (margs.sample_sheet && !sub_dirs_with_fastq_files) { + error "Sample sheet was provided, but input does not contain " + + "sub-directories with FASTQ files." + } + return ch_input +} + + +/** + * Create a map that contains at least these keys: `[alias, barcode, type]`. + * `alias` is required, `barcode` and `type` are filled with default values if + * missing. Additional entries are allowed. + * + * @param kwargs: map with input parameters; must contain `alias` + * @return: map(alias, barcode, type, ...) + */ +Map create_metamap(Map arguments) { + ArgumentParser parser = new ArgumentParser( + args: ["alias"], + kwargs: [ + "barcode": null, + "type": "test_sample", + ], + name: "create_metamap", + ) + return parser.parse_known_args(arguments) +} + + +/** + * Get the fastq files in the directory (non-recursive). + * + * @param dir: path to the target directory + * @return: list of found fastq files + */ +ArrayList get_fq_files_in_dir(Path dir) { + return EXTENSIONS.collect { file(dir.resolve("*.$it"), type: "file") } .flatten() +} + + +/** + * Check the sample sheet and return a channel with its rows if it is valid. + * + * @param sample_sheet: path to the sample sheet CSV + * @return: channel of maps (with values in sample sheet header as keys) + */ +def get_sample_sheet(Path sample_sheet) { + // If `validate_sample_sheet` does not return an error message, we can assume that + // the sample sheet is valid and parse it. However, because of Nextflow's + // asynchronous magic, we might emit values from `.splitCSV()` before the + // error-checking closure finishes. This is no big deal, but undesired nonetheless + // as the error message might be overwritten by the traces of new nextflow processes + // in STDOUT. Thus, we use the somewhat clunky construct with `concat` and `last` + // below. This lets the CSV channel only start to emit once the error checking is + // done. + ch_err = validate_sample_sheet(sample_sheet).map { + // check if there was an error message + if (it) error "Invalid sample sheet: ${it}." + it + } + // concat the channel holding the path to the sample sheet to `ch_err` and call + // `.last()` to make sure that the error-checking closure above executes before + // emitting values from the CSV + return ch_err.concat(Channel.fromPath(sample_sheet)).last().splitCsv( + header: true, quote: '"' + ) +} + + +/** + * Python script for validating a sample sheet. The script will write messages + * to STDOUT if the sample sheet is invalid. In case there are no issues, no + * message is emitted. + * + * @param: path to sample sheet CSV + * @return: string (optional) + */ +process validate_sample_sheet { + label params.process_label + input: path csv + output: stdout + """ + workflow-glue check_sample_sheet $csv + """ } diff --git a/main.nf b/main.nf index ab4b7a1..dc8d283 100644 --- a/main.nf +++ b/main.nf @@ -16,24 +16,7 @@ include { denovo_assembly } from './subworkflows/denovo_assembly' include { gene_fusions } from './subworkflows/JAFFAL/gene_fusions' include { differential_expression } from './subworkflows/differential_expression' - - -process summariseConcatReads { - // concatenate fastq and fastq.gz in a dir write stats - - label "isoforms" - cpus 1 - input: - tuple path(directory), val(meta) - output: - tuple val(meta.sample_id), path("${meta.sample_id}.fastq"), emit: input_reads - tuple val(meta.sample_id), path('*.stats'), emit: summary - script: - """ - - fastcat -s ${meta.sample_id} -r ${meta.sample_id}.stats -x ${directory} > ${meta.sample_id}.fastq - """ -} +OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE") process getVersions { label "isoforms" @@ -68,8 +51,6 @@ process getParams { path "params.json" script: def paramsJSON = new JsonBuilder(params).toPrettyString() - println('test') - println(params.workDir) """ # Output nextflow params object to JSON echo '$paramsJSON' > params.json @@ -86,19 +67,19 @@ process preprocess_reads { cpus 4 input: - tuple val(sample_id), path(input_reads) + tuple val(meta), path(input_reads) output: - tuple val(sample_id), path("${sample_id}_full_length_reads.fastq"), emit: full_len_reads + tuple val("${meta.alias}"), path("${meta.alias}_full_length_reads.fastq"), emit: full_len_reads path '*.tsv', emit: report script: """ - pychopper -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${sample_id}_full_length_reads.fastq - mv pychopper.tsv ${sample_id}_pychopper.tsv - workflow-glue generate_pychopper_stats --data ${sample_id}_pychopper.tsv --output . + pychopper -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${meta.alias}_full_length_reads.fastq + mv pychopper.tsv ${meta.alias}_pychopper.tsv + workflow-glue generate_pychopper_stats --data ${meta.alias}_pychopper.tsv --output . # Add sample id column - sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${sample_id}/" ${sample_id}_pychopper.tsv > tmp - mv tmp ${sample_id}_pychopper.tsv + sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${meta.alias}/" ${meta.alias}_pychopper.tsv > tmp + mv tmp ${meta.alias}_pychopper.tsv """ } @@ -312,7 +293,7 @@ process makeReport { path "pychopper_report/*" path"jaffal_csv/*" val sample_ids - path seq_summaries + path per_read_stats path "aln_stats/*" path gffcmp_dir path "gff_annotation/*" @@ -359,7 +340,7 @@ process makeReport { \$OPT_ALN \ \$OPT_PC_REPORT \ --sample_ids $sids \ - --summaries $seq_summaries \ + --stats $per_read_stats \ \$OPT_GFF \ --isoform_table_nrows $params.isoform_table_nrows \ \$OPT_JAFFAL_CSV \ @@ -369,22 +350,56 @@ process makeReport { """ } -// See https://github.com/nextflow-io/nextflow/issues/1636 -// This is the only way to publish files from a workflow whilst -// decoupling the publish from the process steps. -process output { - // publish inputs to output directory - publishDir "${params.out_dir}", mode: 'copy', pattern: "*" + +// Creates a new directory named after the sample alias and moves the fastcat results +// into it. +process collectFastqIngressResultsInDir { label "isoforms" input: - path fname + // both the fastcat seqs as well as stats might be `OPTIONAL_FILE` --> stage in + // different sub-directories to avoid name collisions + tuple val(meta), path(concat_seqs, stageAs: "seqs/*"), path(fastcat_stats, + stageAs: "stats/*") + output: + // use sub-dir to avoid name clashes (in the unlikely event of a sample alias + // being `seq` or `stats`) + path "out/*" + script: + String outdir = "out/${meta["alias"]}" + String metaJson = new JsonBuilder(meta).toPrettyString() + String concat_seqs = \ + (concat_seqs.fileName.name == OPTIONAL_FILE.name) ? "" : concat_seqs + String fastcat_stats = \ + (fastcat_stats.fileName.name == OPTIONAL_FILE.name) ? "" : fastcat_stats + """ + mkdir -p $outdir + echo '$metaJson' > metamap.json + mv metamap.json $concat_seqs $fastcat_stats $outdir + """ +} + + +// See https://github.com/nextflow-io/nextflow/issues/1636. This is the only way to +// publish files from a workflow whilst decoupling the publish from the process steps. +// The process takes a tuple containing the filename and the name of a sub-directory to +// put the file into. If the latter is `null`, puts it into the top-level directory. +process output { + // publish inputs to output directory + label "isoforms" + publishDir ( + params.out_dir, + mode: "copy", + saveAs: { dirname ? "$dirname/$fname" : fname } + ) + input: + tuple path(fname), val(dirname) output: path fname """ - echo "Writing output files" """ } + // workflow module workflow pipeline { take: @@ -397,6 +412,10 @@ workflow pipeline { condition_sheet ref_transcriptome main: + fastq_ingress_results = reads + // replace `null` with path to optional file + | map { [ it[0], it[1] ?: OPTIONAL_FILE, it[2] ?: OPTIONAL_FILE ] } + | collectFastqIngressResultsInDir map_sample_ids_cls = {it -> /* Harmonize tuples output: @@ -419,19 +438,23 @@ workflow pipeline { return l } - summariseConcatReads(reads) - sample_ids = summariseConcatReads.out.summary.flatMap({it -> it[0]}) + software_versions = getVersions() workflow_params = getParams() + input_reads = reads.map{ meta, samples, stats -> [meta, samples]} + sample_ids = input_reads.flatMap({meta,samples -> meta.alias}) + stats = reads.map { + it[2] ? it[2].resolve('per-read-stats.tsv') : null + } if (!params.direct_rna){ - preprocess_reads(summariseConcatReads.out.input_reads) + preprocess_reads(input_reads) full_len_reads = preprocess_reads.out.full_len_reads pychopper_report = preprocess_reads.out.report.collectFile(keepHeader: true) } else{ - full_len_reads = summariseConcatReads.out.input_reads + full_len_reads = input_reads.map{ meta, reads -> [meta.alias, reads]} pychopper_report = file("$projectDir/data/OPTIONAL_FILE") } if (params.transcriptome_source != "precomputed"){ @@ -446,7 +469,7 @@ workflow pipeline { assembly = reference_assembly(build_minimap_index.out.index, ref_genome, full_len_reads) } assembly_stats = assembly.stats.map{ it -> it[1]}.collect() - + split_bam(assembly.bam) assemble_transcripts(split_bam.out.bundles.flatMap(map_sample_ids_cls), ref_annotation) @@ -465,7 +488,6 @@ workflow pipeline { // So map this reference to all sample_ids seq_for_transcriptome_build = sample_ids.flatten().combine(Channel.fromPath(params.ref_genome)) } - get_transcriptome( merge_gff_bundles.out.gff .join(run_gffcompare.out.gffcmp_dir) @@ -505,8 +527,9 @@ workflow pipeline { check_condition_sheet = check_match.splitCsv(header: true).map{ row -> tuple( row.sample_id) } - check_condition_sheet.join(summariseConcatReads.out.input_reads, failOnMismatch: true) - de = differential_expression(transcriptome, summariseConcatReads.out.input_reads, condition_sheet, gtf) + join_reads = input_reads.map{ meta, reads -> [meta.alias, reads]} + check_condition_sheet.join(join_reads, failOnMismatch: true) + de = differential_expression(transcriptome, input_reads, condition_sheet, gtf) de_report = de.all_de count_transcripts_file = de.count_transcripts dtu_plots = de.dtu_plots @@ -521,8 +544,8 @@ workflow pipeline { workflow_params, pychopper_report, jaffal_out, - summariseConcatReads.out.summary.map{it->it[0]}.collect(), - summariseConcatReads.out.summary.map{it->it[1]}.collect(), + input_reads.map{ meta, fastq -> meta.alias}.collect(), + stats, assembly_stats, gff_compare, merge_gff, @@ -530,7 +553,7 @@ workflow pipeline { count_transcripts_file) report = makeReport.out.report - + results = results.concat(makeReport.out.report) if (use_ref_ann){ @@ -575,6 +598,7 @@ workflow pipeline { results = results.concat(de.dtu_plots, de_outputs) } + results = fastq_ingress_results.map { [it, "fastq_ingress_results"] }.concat(results.map{ [it, null]}) emit: results telemetry = workflow_params @@ -652,10 +676,13 @@ workflow { if (error){ throw new Exception(error) }else{ - reads = fastq_ingress([ - "input":params.fastq, - "sample":params.sample, - "sample_sheet":params.sample_sheet]) + reads = samples = fastq_ingress([ + "input":params.fastq, + "sample":params.sample, + "sample_sheet":params.sample_sheet, + "analyse_unclassified":params.analyse_unclassified, + "fastcat_stats": true, + "fastcat_extra_args": ""]) pipeline(reads, ref_genome, ref_annotation, jaffal_refBase, params.jaffal_genome, params.jaffal_annotation, diff --git a/nextflow.config b/nextflow.config index 4177444..5f6b688 100644 --- a/nextflow.config +++ b/nextflow.config @@ -105,7 +105,7 @@ params { "--jaffal_annotation genCode22" ] agent = null - container_sha = "sha3d3c83523695550f398cbe095551b1192de5085a" + container_sha = "mr131_sha203915eb4b4dd444cb2e845d0b9f7814e26b7b5c" } } diff --git a/nextflow_schema.json b/nextflow_schema.json index e30b4a8..d6b4c02 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -219,9 +219,9 @@ "condition_sheet": { "type": "string", "format": "file-path", - "description": "TSV file with sample_id, condition", + "description": "CSV file with sample_id, condition", "default": "null", - "help_text": "The condition sheet should be a headed TSV file with two columns sample_id,condition. Should be at least 3 repeats for each condition." + "help_text": "The condition sheet should be a headed CSV file with two columns sample_id,condition. Should be at least 3 repeats for each condition." }, "min_gene_expr": { "type": "integer", diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index 915f5df..8a5998d 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -118,16 +118,16 @@ process map_transcriptome{ cpus params.threads input: - tuple val(sample_id), path (fastq_reads) + tuple val(meta), path (fastq_reads) file index file transcript_reference output: - tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam + tuple val("${meta.alias}"), path("${meta.alias}_reads_aln_sorted.bam"), emit: bam """ minimap2 -t ${task.cpus} -ax splice -uf -p 1.0 "${index}" "${fastq_reads}" \ | samtools view -Sb > "output.bam" - samtools sort -@ ${task.cpus} "output.bam" -o "${sample_id}_reads_aln_sorted.bam" - samtools index "${sample_id}_reads_aln_sorted.bam" + samtools sort -@ ${task.cpus} "output.bam" -o "${meta.alias}_reads_aln_sorted.bam" + samtools index "${meta.alias}_reads_aln_sorted.bam" """ }