Sample with no reads

### Operating System

Mac

### Other Linux

_No response_

### Workflow Version

v1.7.0

### Workflow Execution

EPI2ME Desktop (Local)

### Other workflow execution

_No response_

### EPI2ME Version

v5.2.5

### CLI command run

_No response_

### Workflow Execution - CLI Execution Profile

None

### What happened?

Hi I am using the wf-transcriptomes but with slight modifications.
Currently I have sequencing runs where I have only used some barcodes so for the unused barcodes I have some erroneously assigned reads.

When I run the workflow for all the barcodes detected during base calling (located in the same directory) I get an error saying that a specific sample does not have any reads after QC using pychopper and thus the pipeline stops.

Can you help me modify your original pipeline so it does not need BAM files for all samples to function?
I tried to indicate that the output in each process is optional but still the pipeline stops.

Thanks!,
Antonio

P.S.: Under these lines you can find my modified version of the "main.nf" file.

#!/usr/bin/env nextflow

/* This workflow is a adapted from two previous pipeline written in Snakemake:
- https://github.com/nanoporetech/pipeline-nanopore-ref-isoforms
*/

import groovy.json.JsonBuilder;
import nextflow.util.BlankSeparatedList;
import java.util.ArrayList;
nextflow.enable.dsl = 2

include { fastq_ingress; xam_ingress } from './lib/ingress'
include { configure_igv } from './lib/common'
include { reference_assembly } from './subworkflows/reference_assembly'
include { differential_expression } from './subworkflows/differential_expression'

OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")

process getVersions {
    label "isoforms"
    cpus 1
    memory "2 GB"
    output:
        path "versions.txt"
    script:
    """
    python -c "import pysam; print(f'pysam,{pysam.__version__}')" >> versions.txt
    python -c "import pandas; print(f'pandas,{pandas.__version__}')" >> versions.txt
    python -c "import sklearn; print(f'scikit-learn,{sklearn.__version__}')" >> versions.txt
    fastcat --version | sed 's/^/fastcat,/' >> versions.txt
    minimap2 --version | sed 's/^/minimap2,/' >> versions.txt
    samtools --version | head -n 1 | sed 's/ /,/' >> versions.txt
    bedtools --version | head -n 1 | sed 's/ /,/' >> versions.txt
    python -c "import pychopper; print(f'pychopper,{pychopper.__version__}')" >> versions.txt
    gffread --version | sed 's/^/gffread,/' >> versions.txt
    seqkit version | head -n 1 | sed 's/ /,/' >> versions.txt
    stringtie --version | sed 's/^/stringtie,/' >> versions.txt
    gffcompare --version | head -n 1 | sed 's/ /,/' >> versions.txt
    """
}


process getParams {
    label "isoforms"
    cpus 1
    memory "2 GB"
    output:
        path "params.json"
    script:
        def paramsJSON = new JsonBuilder(params).toPrettyString()
    """
    # Output nextflow params object to JSON
    echo '$paramsJSON' > params.json
    """
}



process decompress_ref {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path compressed_ref
    output:
        path "${compressed_ref.baseName}", emit: decompressed_ref
    """
    gzip -df ${compressed_ref}
    """
}

process validate_ref_annotation {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path "annotation.gtf"
        path "reference.fasta"
    output: 
        stdout
    // Checks for overlap between seq_id column in annotation gtf and fasta reference ID's
    // If no overlap is found exit
    // Partial overlap (eg. user supplies genes/contigs of interest in annotation but the genome sequence) - warning
    script:
    """
    grep -v '^#' annotation.gtf | cut -f1  | sort -u > seq_ids.txt
    awk '/^>/ {print substr(\$1,2)}' reference.fasta | sort -u > ref_ids.txt
    matches=\$(comm -12 seq_ids.txt ref_ids.txt)
    only_in_annotation=\$(comm -23 seq_ids.txt ref_ids.txt)
    only_in_reference=\$(comm -13 seq_ids.txt ref_ids.txt)
    if [[ -z "\$matches" ]]; then
        echo "
        ERROR: Seqid mismatch found between the provided ref_annotation (GTF/GFF)
        file and ref_genome (FASTA).
        For the reference guided differential expression subworkflow they must overlap.
        " >&2
        echo "Annotation ID examples:"
        head -n 5 seq_ids.txt
        echo "Reference ID examples:"
        head -n 5 ref_ids.txt
        echo "We recommend getting both files from the same source.
        eg. both from Ensembl or both from NCBI.
        Alternatively provide a pre-computed transcriptome using the ref_transcriptome parameter
        See the README for more details on which inputs are supported."
        exit 78 
    fi
    if [[ -n "\$only_in_annotation" ]]; then
        echo "Warning: Some sequence IDs are only present in the reference annotation and not the
             reference genome so will not be used in downstream analysis eg."
        echo "\$only_in_annotation" | head -n 5
    fi
    if [[ -n "\$only_in_reference" ]]; then
        echo "Warning: Some FASTA reference IDs are only present in the reference genome
            and not the reference annotation so will not be used in downstream analysis eg."
        echo "\$only_in_reference" | head -n 5
    fi
    """
}


process decompress_annotation {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path compressed_annotation
    output:
        path "${compressed_annotation.baseName}"
    """
    gzip -df ${compressed_annotation}
    """
}



process decompress_transcriptome {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path "compressed_ref.gz"
    output:
        path "compressed_ref", emit: decompressed_ref
    """
    gzip -df "compressed_ref.gz"
    """
}


// Remove empty transcript ID fields
process preprocess_ref_annotation {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path ref_annotation
    output:
        path "amended.${ref_annotation}"
    """
    sed -i -e 's/transcript_id "";//g' ${ref_annotation}
    mv ${ref_annotation} "amended.${ref_annotation}"
    """
}

// Just keep transcript ID for each transcriptome fasta
process preprocess_ref_transcriptome {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path "ref_transcriptome"
    output:
        path "amended.${ref_transcriptome}"
    """
    sed -i -e 's/|.*//' ${ref_transcriptome}
    mv ${ref_transcriptome} "amended.${ref_transcriptome}"
    """
}



process preprocess_reads {
    /*
    Concatenate reads from a sample directory.
    Optionally classify, trim, and orient cDNA reads using pychopper
    */

    label "isoforms"
    cpus params.threads
    memory "2 GB"
    input:
        tuple val(meta), path('seqs.fastq.gz')
    output:
        tuple val("${meta.alias}"),
              path("${meta.alias}_pychopper_output/${meta.alias}_full_length_reads.fastq"),
              emit: full_len_reads
        tuple val("${meta.alias}"),
              path("${meta.alias}_pychopper_output/"),
              emit: pychopper_output
        path("${meta.alias}_pychopper_output/pychopper.tsv"),
              emit: report
    script:
        def cdna_kit = params.cdna_kit.split("-")[-1]
	    def extra_params = params.pychopper_opts ?: ''
        """
        pychopper -t ${params.threads} -u unclassified_reads.fastq -l len_fail_reads.fastq -w rescue_reads.fastq -K qc_fail_reads.fastq -k ${cdna_kit} -m ${params.pychopper_backend} ${extra_params} 'seqs.fastq.gz' ${meta.alias}_full_length_reads.fastq
        workflow-glue generate_pychopper_stats --data pychopper.tsv --output .

        # Add sample id column
        sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${meta.alias}/" pychopper.tsv > tmp
        mv tmp pychopper.tsv


        mkdir "${meta.alias}_pychopper_output/"
        shopt -s extglob  # Allow extended pattern matching so we can exclude files from the mv
        mv !("${meta.alias}_pychopper_output"|seqs.fastq.gz) "${meta.alias}_pychopper_output/"
        """
}

process build_minimap_index{
    /*
    Build minimap index from reference genome
    */
    label "isoforms"
    cpus params.threads
    memory "31 GB"

    input:
        path reference
    output:
        path "genome_index.mmi", emit: index
    script:
    """
    minimap2 -t ${params.threads} ${params.minimap2_index_opts} -I 1000G -d "genome_index.mmi" ${reference}
    """
}

process split_bam{
    /*
    Partition BAM file into loci or bundles with `params.bundle_min_reads` minimum size
    If no splitting required, just create single symbolic link to a single bundle.

    */

    label 'isoforms'
    cpus params.threads
    memory "15 GB"

    input:
        tuple val(sample_id), path(bam)
    output:
        tuple val(sample_id), path('*.bam'), emit: bundles, optional: true
    script:
    """
    n=`samtools view -c $bam`
    echo "Number of reads before split_bam: \$n"
    if [[ n -lt 1 ]]
    then
        echo 'There are no reads mapping for $sample_id. Skipping $sample_id!'
        exit 0
    fi

    re='^[0-9]+\$'

    if [[ $params.bundle_min_reads =~ \$re ]]
    then
        echo "Bundling up the bams"
        seqkit bam -j ${params.threads} -N ${params.bundle_min_reads} ${bam} -o  bam_bundles/
        let i=1
        for b in bam_bundles/*.bam; do
            echo \$b
            newname="${sample_id}_batch_\${i}.bam"
            mv \$b \$newname
           ((i++))
        done
    else
        echo 'no bundling'
        ln -s ${bam} ${sample_id}_batch_1.bam
    fi
    """
}


process assemble_transcripts{
    /*
    Assemble transcripts using stringtie.
    Take aligned reads in bam format that may be a chunk of a larger alignment file.
    Optionally use reference annotation to guide assembly.

    Output gff annotation files in a tuple with `sample_id` for combining into samples later in the pipeline.
    */
    label 'isoforms'
    cpus params.threads
    memory "2 GB"

    input:
        tuple val(sample_id), path(bam), path(ref_annotation)
        val use_ref_ann
    output:
        tuple val(sample_id), path('*.gff'), emit: gff_bundles, optional: true
    script:
        def G_FLAG = use_ref_ann == false ? '' : "-G ${ref_annotation}"
        def prefix =  bam.name.split(/\./)[0]

    """
    stringtie --rf ${G_FLAG} -L -v -p ${task.cpus} ${params.stringtie_opts} \
    -o  ${prefix}.gff -l ${prefix} ${bam}
     """
}


process merge_gff_bundles{
    /*
    Merge gff bundles into a single gff file per sample, and get summary statistics
    */
    label 'isoforms'
    cpus params.threads
    memory "2 GB"

    input:
        tuple val(sample_id), path ('gff_bundles/annotation*.gff')
    output:
        tuple val(sample_id), path("${sample_id}.gff"), emit: gff, optional: true
        tuple val(sample_id), path("transcriptome_summary.pickle"), emit: summary, optional: true
    script:
    def merged_gff = "${sample_id}.gff"
    """
    echo '##gff-version 2' >> $merged_gff;
    echo '#pipeline-nanopore-isoforms: stringtie' >> $merged_gff;

    find -L gff_bundles -type f -name "*.gff" \
        -exec awk '!/^#/ {print}' {} \\; >> "${sample_id}.gff"
    if ! [ -s "${sample_id}.gff" ]; then
        echo "No transcripts found for ${sample_id}"
        # This is unlikely to ever happen, but if it does, we should fail the workflow.
        exit 0
    fi

    workflow-glue summarise_gff \
        $merged_gff \
        $sample_id \
        transcriptome_summary.pickle
    """
}

process run_gffcompare{
    /*
    Compare query and reference annotations.
    If ref_annotation is an optional file, just make an empty directory to satisfy
    the requirements of the downstream processes.
    */

    label 'isoforms'
    cpus 1
    memory "2 GB"
    input:
       tuple val(sample_id), path(query_annotation)
       path ref_annotation
    output:
        tuple val(sample_id), path("${sample_id}"), emit: gffcmp_dir, optional: true
        path ("${sample_id}_annotated.gtf"), emit: gtf, optional: true
        tuple val(sample_id), path("${sample_id}_transcripts_table.tsv"), optional: true,
            emit: isoforms_table
    script:
    def out_dir = "${sample_id}"
    """
        mkdir $out_dir
        echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation"

        gffcompare -o ${out_dir}/str_merged -r ${ref_annotation} \
            ${params.gffcompare_opts} ${query_annotation}

        mv *.tmap "${out_dir}"
        mv *.refmap "${out_dir}"
        cp "${out_dir}/str_merged.annotated.gtf" "${sample_id}_annotated.gtf"

        workflow-glue parse_gffcompare \
            --sample_id "${sample_id}" \
            --gffcompare_dir "${out_dir}" \
            --isoform_table_out "${sample_id}_transcripts_table.tsv" \
            --tracking $out_dir/str_merged.tracking \
            --annotation ${ref_annotation}
        """
}

process get_transcriptome{
        /*
        Write out a transcriptome file based on the query gff annotations.
        */
        label 'isoforms'
        cpus 1
        memory "2 GB"
        input:
            tuple val(sample_id), path("transcripts.gff"), path(gffcompare_dir), path("reference.fa")
        output:
            tuple val(sample_id), path("*transcriptome.fas"), emit: transcriptome, optional: true
        script:
        def transcriptome = "${sample_id}_transcriptome.fas"
        def merged_transcriptome = "${sample_id}_merged_transcriptome.fas"
        // if no ref_annotation gffcmp_dir will be optional file
        // so skip getting transcriptome FASTA from the annotated files.
        if (params.ref_annotation){
        """
        gffread -F -g reference.fa -w ${merged_transcriptome} $gffcompare_dir/str_merged.annotated.gtf
        """
        } else {
        """
        gffread -g reference.fa -w ${transcriptome} "transcripts.gff"
        """
        }
}

process merge_transcriptomes {
    // Merge the transcriptomes from all samples
    label 'isoforms'
    cpus 2
    memory "2 GB"
    input:
        path "query_annotations/*"
        path ref_annotation
        path ref_genome
    output:
        path "final_non_redundant_transcriptome.fasta", emit: fasta
        path "stringtie.gtf", emit: gtf
    """
    stringtie --merge -G "${ref_annotation}" -p ${task.cpus} -o stringtie.gtf query_annotations/*
    gffread -g "${ref_genome}" -w "final_non_redundant_transcriptome.fasta" "stringtie.gtf"
    """
}

process makeReport {

    label "wf_common"
    cpus 2
    memory "4 GB"
    publishDir "${params.out_dir}", mode: 'copy', pattern: "wf-transcriptomes-report.html"
    input:
        val metadata
        path stats, stageAs: "stats_*"
        path versions
        val wf_version
        path "params.json"
        path "transcriptome_aln_stats/*"
        path pychopper, stageAs: "pychopper_report/*"
        path aln_stats, stageAs: "aln_stats/*"
        path gffcmp_dir, stageAs: "gffcmp_dir/*"
        path gff_annotation, stageAs: "gff_annotation/*"
        path de_report, stageAs: "de_report/*"
        path isoforms_table, stageAs: "isoforms_table/*"
        path transcriptome_summary, stageAs: "transcriptome_summary/summary_*.pkl"

    output:
        path ("wf-transcriptomes-*.html"), emit: report
        path ("results_dge.tsv"), emit: results_dge, optional: true
        path ("unfiltered_tpm_transcript_counts.tsv"), emit: tpm, optional: true
        path ("unfiltered_transcript_counts_with_genes.tsv"), emit: unfiltered, optional: true
        path ("filtered_transcript_counts_with_genes.tsv"), emit: filtered, optional: true
        path ("all_gene_counts.tsv"), emit: gene_counts, optional: true
    script:
        String report_name = "wf-transcriptomes-report.html"
        String metadata = new JsonBuilder(metadata).toPrettyString()
        String gff_opts = gff_annotation.fileName.name == OPTIONAL_FILE.name ? "" :  "--gff_annotation gff_annotation/"
        String de_report_opts = de_report.fileName.name == OPTIONAL_FILE.name ? "" : "--de_report de_report/ --de_stats transcriptome_aln_stats/"
        String gffcmp_opts = gffcmp_dir.fileName.name == OPTIONAL_FILE.name  ? "" :  "--gffcompare_dir gffcmp_dir/"
        String aln_stats_opts = aln_stats.fileName.name == OPTIONAL_FILE.name ? "" :  "--alignment_stats aln_stats/"
        String pychop_opts = pychopper.fileName.name == OPTIONAL_FILE.name ? "" :  "--pychop_report pychopper_report/"
        String iso_table_opts = isoforms_table.fileName.name == OPTIONAL_FILE.name ? "" :  "--isoform_table isoforms_table/"
        String tr_summary_opts = transcriptome_summary.fileName.name == OPTIONAL_FILE.name ? "" :  "--transcriptome_summary transcriptome_summary/"
    """
    echo '${metadata}' > metadata.json
    workflow-glue report \
        --report $report_name \
        --versions $versions \
        --wf_version $wf_version \
        --params params.json \
        $aln_stats_opts \
        $pychop_opts \
        --stats $stats \
        --metadata metadata.json \
        $gff_opts \
        $iso_table_opts \
        $gffcmp_opts \
        --isoform_table_nrows ${params.isoform_table_nrows} \
        $de_report_opts \
        $tr_summary_opts
    """
}


// Creates a new directory named after the sample alias and moves the fastcat results
// into it.
process collectFastqIngressResultsInDir {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        // 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 publish_results {
    // publish inputs to output directory
    label "isoforms"
    cpus 1
    memory "2 GB"
    publishDir (
        params.out_dir,
        mode: "copy",
        saveAs: { dirname ? "$dirname/$fname" : fname }
    )
    input:
        tuple path(fname), val(dirname)
    output:
        path fname
    """
    """
}


// Check ref_annotation transcript strand column for "." if in de_analysis mode
process check_annotation_strand {
    label "isoforms"
    cpus 1
    memory "2 GB"
    input:
        path "ref_annotation.gtf"
    output:
        tuple stdout, path("ref_annotation.gtf")   
    """
    awk '{if (\$3=="transcript" && \$7 != "+" && \$7 != "-") print \$3, \$7}' "ref_annotation.gtf"
    """
}


// Process to create the faidx index
process faidx {
    // If the input file is gzipped, we need to emit the indexes for the input gzip file
    // only. Therefore, this become redundant to be emitted as it won't be used by the
    // IGV configuration, but only by internal processes together with the decompressed
    // FASTA file. To avoid unnecessary emissions, we enable only if the input file is
    // decompressed.
    publishDir "${params.out_dir}/igv_reference", mode: 'copy', pattern: "*", enabled: !params.ref_genome.toLowerCase().endsWith("gz")
    label "wf_common"
    cpus 1
    memory 4.GB
    input:
        path(ref)
    output:
        path("${ref}.fai")
    script:
    """
    samtools faidx ${ref}
    """
}


// Process to create the faidx indexes for a gzipped reference
process gz_faidx {
    publishDir "${params.out_dir}/igv_reference", mode: 'copy', pattern: "*"
    label "wf_common"
    cpus 1
    memory 4.GB
    // If a user provides a non-bgzipped file, the process won't
    // generate the indexes. We should tolerate that, still avoid emitting
    // the reference and simply have a broken IGV file.
    // The gzi is not required to operate the workflow, so we actually tolerate any failure.
    errorStrategy 'ignore'
    input:
        path(ref)
    output:
        tuple path("${ref}.fai"), path("${ref}.gzi")
    script:
    """
    samtools faidx ${ref}
    """
}



// workflow module
workflow pipeline {
    take:
        reads
        ref_genome
        ref_annotation
        ref_transcriptome
        use_ref_ann
    main:
        if (params.ref_genome && file(params.ref_genome).extension == "gz") {
            // gzipped ref not supported by some downstream tools
            // easier to just decompress and pass it around.
            ref_genome = decompress_ref(ref_genome)
        }else {
            ref_genome = Channel.fromPath(ref_genome)
        }
        if (params.ref_annotation && file(params.ref_annotation).extension == "gz") {
            // gzipped ref not supported by some downstream tools
            // easier to just decompress and pass it around.
            decompress_annot= decompress_annotation(ref_annotation)
            ref_annotation = preprocess_ref_annotation(decompress_annot)
        }else {
            ref_annotation = preprocess_ref_annotation(ref_annotation)
        }
        
        fastq_ingress_results = reads
        | collectFastqIngressResultsInDir
        // fastq_ingress doesn't have the index; add one extra null for compatibility.
        // We do not use variable name as assigning variable name with a tuple
        // not matching (e.g. meta, bam, bai, stats <- [meta, bam, stats]) causes
        // the workflow to crash.
        reads = reads
        .map{
            it.size() == 4 ? it : [it[0], it[1], null, it[2]]
        }
        map_sample_ids_cls = {it ->
        /* Harmonize tuples
        output:
            tuple val(sample_id), path('*.gff')
        When there are multiple paths, will emit:
            [sample_id, [path, path ..]]
        when there's a single path, this:
            [sample_id, path]
        This closure makes both cases:
            [[sample_id, path][sample_id, path]].
        */
            if (it[1].getClass() != java.util.ArrayList){
                // If only one path, `it` will be [sample_id, path]
                return [it]
            }
            l = [];
            for (x in it[1]){
                l.add(tuple(it[0], x))
            }
            return l
        }

     
        results = Channel.empty()
        // Define BAM output Directory
        String publish_prefix_bams = "BAMS"
        software_versions = getVersions()
        workflow_params = getParams()
        input_reads = reads.map{ meta, samples, index, stats -> [meta, samples]}
        sample_ids = input_reads.flatMap({meta,samples -> meta.alias})

        if (!params.direct_rna){
            preprocess_reads(input_reads)
            full_len_reads = preprocess_reads.out.full_len_reads
            pychopper_report = preprocess_reads.out.report.collectFile(keepHeader: true)
            pychopper_results_dir = preprocess_reads.out.pychopper_output.map{ it -> it[1]}
            results = results.concat(pychopper_results_dir)
        }
        else{
            full_len_reads = input_reads.map{ meta, reads -> [meta.alias, reads]}
            pychopper_report = OPTIONAL_FILE
        }
        
        if (params.transcriptome_source != "precomputed"){
            build_minimap_index(ref_genome)
            log.info("Doing reference based transcript analysis")
            assembly = reference_assembly(build_minimap_index.out.index, ref_genome, full_len_reads, publish_prefix_bams)
        
            assembly_stats = assembly.stats_1_supp.map { it[1] }.collect()
            // New lines to store also the info from other subsets not included in the original script
            //assembly_stats_1_supp = assembly.stats_1_supp.map { it[1] }.collect()
            //assembly_stats_unmapped = assembly.stats_unmapped.map { it[1] }.collect()
            //assembly_stats_lowMAPQ = assembly.stats_lowMAPQ.map { it[1] }.collect()

            split_bam(assembly.bam_1_supp.map {sample_id, bam, bai -> [sample_id, bam]})
            // New lines to store also the info from other subsets not included in the original script
            //split_bam(assembly.bam_1_supp.map {sample_id, bam, bai -> [sample_id, bam]})
            //split_bam(assembly.bam_unmapped.map {sample_id, bam, bai -> [sample_id, bam]})
            //split_bam(assembly.bam_lowMAPQ.map {sample_id, bam, bai -> [sample_id, bam]})

            assemble_transcripts(split_bam.out.bundles.flatMap(map_sample_ids_cls).combine(ref_annotation),use_ref_ann)

            merge_gff_bundles(assemble_transcripts.out.gff_bundles.groupTuple())
            transcriptome_summary = merge_gff_bundles.out.summary.map {it[1]}.collect()

            // only run gffcompare if ref annotation provided. Otherwise create optional files and channels
            if (params.ref_annotation){
                run_gffcompare(merge_gff_bundles.out.gff, ref_annotation)
                gff_compare_dir = run_gffcompare.out.gffcmp_dir
                gff_compare = run_gffcompare.out.gffcmp_dir.map{ it -> it[1]}.collect()
                isoforms_table = run_gffcompare.out.isoforms_table.map{ it -> it[1]}.collect()
                // create per sample gff tuples with gff compare directories
                gff_tuple = merge_gff_bundles.out.gff
                .join(gff_compare_dir)
            } else {
                // create per sample gff tuples with optional files as no ref_annotation
                optional_channel = Channel.fromPath("$projectDir/data/OPTIONAL_FILE")
                gff_tuple = merge_gff_bundles.out.gff.combine(optional_channel)
                gff_compare = OPTIONAL_FILE
                isoforms_table = OPTIONAL_FILE
            }
            // For reference based assembly, there is only one reference
            // So map this reference to all sample_ids
            seq_for_transcriptome_build = sample_ids.flatten().combine(ref_genome)
            get_transcriptome(
                gff_tuple
                .join(seq_for_transcriptome_build))


            merge_gff = merge_gff_bundles.out.gff.map{ it -> it[1]}.collect()
        }
        else{
            gff_compare = OPTIONAL_FILE
            isoforms_table = OPTIONAL_FILE
            merge_gff = OPTIONAL_FILE
            assembly_stats = OPTIONAL_FILE
            transcriptome_summary = OPTIONAL_FILE
            use_ref_ann = false
        }
        if (params.de_analysis){
            sample_sheet = file(params.sample_sheet, type:"file")
            // check ref annotation contains only + or - strand as DE analysis will error on .
            check_annotation_strand(ref_annotation).map { stdoutput, annotation ->
            // check if there was an error message
            if (stdoutput) error "In ref_annotation, transcript features must have a strand of either '+' or '-'."
                    stdoutput
                }
            if (!params.ref_transcriptome){
                validate_ref_annotation(ref_annotation, ref_genome).map { stdoutput ->
                    if (stdoutput) {
                        log.warn(stdoutput)
                    }
                }
                merge_transcriptomes(run_gffcompare.output.gtf.collect(), ref_annotation, ref_genome)
                transcriptome = merge_transcriptomes.out.fasta
                gtf = merge_transcriptomes.out.gtf
            }
            else {
                transcriptome =  Channel.fromPath(ref_transcriptome)
                if (file(params.ref_transcriptome).extension == "gz") {
                    transcriptome = decompress_transcriptome(ref_transcriptome)
                }
                transcriptome = preprocess_ref_transcriptome(transcriptome)
                gtf = ref_annotation
            }
            de = differential_expression(transcriptome, full_len_reads.map{ sample_id, reads -> [[alias:sample_id], reads]}, sample_sheet, gtf)
            de_report = de.all_de
            de_outputs = de.de_outputs
            de_alignment_stats = de.de_alignment_stats
        } else{
            de_report = OPTIONAL_FILE
            de_alignment_stats = OPTIONAL_FILE
        }

        // get metadata and stats files, keeping them ordered (could do with transpose I suppose)
        reads.multiMap{ meta, path, index, stats ->
            meta: meta
            stats: stats
        }.set { for_report }
        metadata = for_report.meta.collect()
        stats = for_report.stats.collect()

        makeReport(
            metadata,
            stats,
            software_versions,
            workflow.manifest.version,
            workflow_params,
            de_alignment_stats,
            pychopper_report,
            assembly_stats,
            gff_compare,
            merge_gff,
            de_report,
            isoforms_table,
            transcriptome_summary)

       report = makeReport.out.report

       results = results.concat(report)
       
       if (use_ref_ann){
            results = run_gffcompare.output.gffcmp_dir.concat(
                      assembly.stats,
                      run_gffcompare.out.isoforms_table,
                      get_transcriptome.out.transcriptome.flatMap(map_sample_ids_cls))
                      .map {it -> it[1]}
                      .concat(results)

       }
    
        if (!use_ref_ann && params.transcriptome_source == "reference-guided"){
            results =  assembly.stats.concat(
                       get_transcriptome.out.transcriptome.flatMap(map_sample_ids_cls))
                      .map {it -> it[1]}
                      .concat(results)

        }

       results = results.map{ [it, null] }.concat(fastq_ingress_results.map { [it, "fastq_ingress_results"] })

        if (params.de_analysis){
           de_results = report.concat(
            transcriptome, de_outputs.flatten(),
            makeReport.out.results_dge,  makeReport.out.tpm,
            makeReport.out.filtered,  makeReport.out.unfiltered,
            makeReport.out.gene_counts)
            // Output de_analysis results in the dedicated directory.
            results = results.concat(de_results.map{ [it, "de_analysis"] })
        }


        results.concat(workflow_params.map{ [it, null]})
        // IGV config
        if (params.transcriptome_source == "precomputed" && params.igv){
            log.warn("IGV configuration does not work if transcriptome sources is set to `precomputed`.")
        }
        if (params.transcriptome_source != "precomputed" && params.igv){
            is_compressed = file("${params.ref_genome}").extension == "gz"
            String publish_ref = "igv_reference"
            reference_genome = Channel.fromPath("${params.ref_genome}")
            igv_ref = reference_genome | flatten | map { it -> "${it.toUriString()}" }
            if (is_compressed){
                // Define indexes names.
                String input_fai_index = "${params.ref_genome}.fai"
                String input_gzi_index = "${params.ref_genome}.gzi"
                // Check whether the input gzref is indexed. If so, pass these as indexes.
                // Otherwise, generate the gzip + fai indexes for the compressed reference.
                if (file(input_fai_index).exists() && file(input_gzi_index).exists()){
                    gzindexes = Channel.fromPath(input_fai_index)
                    | mix(
                        Channel.fromPath(input_gzi_index)
                    )
                    gz_igv = gzindexes | flatten | map { it -> "${it.toUriString()}" }
                } else {
                    gz_igv = gz_faidx(Channel.fromPath("${params.ref_genome}"))
                    | flatten
                    | map  { it -> "$publish_ref/${it.Name}" }
                    gz_igv | ifEmpty{
                        if (params.containsKey("igv") && params.igv){
                            log.warn """\
                                The input reference is compressed but not with bgzip, which is required to create an index.
                                The workflow will proceed but it will not be possible to load the reference in the IGV Viewer.
                                To use the IGV Viewer, provide an uncompressed, or bgzip compressed version of the input reference next time you run the workflow.
                                """.stripIndent()
                        }
                    }
                }
            } else {
                gzindexes = Channel.empty()
                gz_igv = Channel.empty()
            }

            // Generate fai index if the file is either compressed, or if fai doesn't exists
            if (!is_compressed && file("${params.ref_genome}.fai").exists()){
                ref_idx = Channel.fromPath("${params.ref_genome}.fai")
                igv_index = ref_idx | flatten | map { it -> "${it.toUriString()}" }
            } else {
                ref_idx = faidx(reference_genome)
                igv_index = ref_idx | map { it -> "$publish_ref/${it.Name}" }
            }

            // get list of file names
            // Absolute paths required for directories
            igv_files = reads
            | map { meta, sample, index, stats -> meta.alias }
            | toSortedList
            | map { list -> list.collect{
                [
                 "$publish_prefix_bams/${it}_reads_aln_sorted.bam",
                 "$publish_prefix_bams/${it}_reads_aln_sorted.bam.bai"
                ]
            } }
            | concat (igv_ref)
            | flatten
            | concat ( igv_index)
            | concat (gz_igv)
            | flatten
            | collectFile(name: "file-names.txt", newLine: true, sort: false)

            // configure IGV
            igv_conf = configure_igv(
                igv_files,
                Channel.of(null), // igv locus
                [displayMode: "SQUISHED", colorBy: "strand"], // bam extra opts
                Channel.of(null), // vcf extra opts
                )

            results = results.concat(igv_conf.map{ [it, null]})
      }

       
    emit:
        results
}

// entrypoint workflow
WorkflowMain.initialise(workflow, params, log)
workflow {

    Pinguscript.ping_start(nextflow, workflow, params)

    error = null

    if (params.containsValue("jaffal_refBase")) {
        error = "JAFFAL fusion detection has been removed from this workflow."
    }
    if (params.containsKey("minimap_index_opts")) {
        error = "`--minimap_index_opts` parameter is deprecated. Use parameter `--minimap2_index_opts` instead."
    }

    if (params.transcriptome_source == "precomputed" && !params.ref_transcriptome){
        error = "As transcriptome source parameter is precomputed you must include a ref_transcriptome parameter"
    }
    if (params.transcriptome_source == "reference-guided" && !params.ref_genome){
        error = "As transcriptome source is reference guided you must include a ref_genome parameter"
    }
    if (params.ref_genome){
        ref_genome = file(params.ref_genome, type: "file")
        if (!ref_genome.exists()) {
            error = "--ref_genome: File doesn't exist, check path."
        }
    }else {
        ref_genome = OPTIONAL_FILE
    }
    if (params.containsValue("denovo")) {
        error = "Denovo transcriptome source is no longer supported. Please use the reference-guided or precomputed options."
    }

    if (params.ref_annotation){
        ref_annotation = file(params.ref_annotation, type: "file")

        if (!ref_annotation.exists()) {
            error = "--ref_annotation: File doesn't exist, check path."
        }
        use_ref_ann = true
    }else{
        ref_annotation= OPTIONAL_FILE
        use_ref_ann = false
    }
    ref_transcriptome = OPTIONAL_FILE
    if (params.ref_transcriptome){
        log.info("Reference Transcriptome provided will be used for differential expression.")
        ref_transcriptome = file(params.ref_transcriptome, type:"file")
    }
    if (params.de_analysis){
        if (!params.ref_annotation){
            error = "When running in --de_analysis mode you must provide a reference annotation."
        }
        if (!params.sample_sheet){
            error = "You must provide a sample_sheet with at least alias and condition columns."
        }
        if (params.containsKey("condition_sheet")) {
        error = "Condition sheets have been deprecated. Please add a 'condition' column to your sample sheet instead. Check the quickstart for more information."
        }
    } else{
        if (!params.ref_annotation){
            log.info("Warning: As no --ref_annotation was provided, the output transcripts will not be annotated.")
        }
    }
    if (error){
        throw new Exception(error)
    }

    if (params.fastq) {
        samples = fastq_ingress([
            "input":params.fastq,
            "sample":params.sample,
            "sample_sheet":params.sample_sheet,
            "analyse_unclassified":params.analyse_unclassified,
            "stats": true,
            "fastcat_extra_args": "",
            "per_read_stats": true])
    } else {
        samples = xam_ingress([
            "input":params.bam,
            "sample":params.sample,
            "sample_sheet":params.sample_sheet,
            "analyse_unclassified":params.analyse_unclassified,
            "keep_unaligned": true,
            "return_fastq": true,
            "stats": true,
            "per_read_stats": true])
    }

    pipeline(samples, ref_genome, ref_annotation, ref_transcriptome, use_ref_ann)
    publish_results(pipeline.out.results)
}

workflow.onComplete {
    Pinguscript.ping_complete(nextflow, workflow, params)
}
workflow.onError {
    Pinguscript.ping_error(nextflow, workflow, params)
}


### Relevant log output

```shell
.
```

### Application activity log entry

```shell

```

### Were you able to successfully run the latest version of the workflow with the demo data?

yes

### Other demo data information

```shell

```

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Sample with no reads #177

Operating System

Other Linux

Workflow Version

Workflow Execution

Other workflow execution

EPI2ME Version

CLI command run

Workflow Execution - CLI Execution Profile

What happened?

Relevant log output

Application activity log entry

Were you able to successfully run the latest version of the workflow with the demo data?

Other demo data information

Metadata

Assignees

Labels

Type

Fields

Projects

Milestone

Relationships

Development

Uh oh!

Sample with no reads #177

Description

Operating System

Other Linux

Workflow Version

Workflow Execution

Other workflow execution

EPI2ME Version

CLI command run

Workflow Execution - CLI Execution Profile

What happened?

Relevant log output

Application activity log entry

Were you able to successfully run the latest version of the workflow with the demo data?

Other demo data information

Metadata

Metadata

Assignees

Labels

Type

Fields

Projects

Milestone

Relationships

Development

Issue actions