How to handle output from a python script that gives different outputs depending upon the condition - Nextflow

Question

I have created three Nextflow processes: One checks the strandedness of the fastq file and outputs it to a .txt file. Second one executes a python script with different commands depending upon the condition. And the third one aligns them.

/*
 * the `CHECK_STRANDEDNESS` process takes SAM file 
 * and a BED file
 * and outputs a a txt file containing information about strandedness. 
 */
 
process CHECK_STRANDEDNESS {

    input:
    path "${pair_id}.sam"
    path "${gtf.baseName}.bed"

    output:
    path "${pair_id}_strand_check.txt"

    script:
    """
    infer_experiment.py -i "${pair_id}.sam" -r "${gtf.baseName}.bed" > "${pair_id}_strand_check.txt"
    """
}


/*
 * the `FINDTAIL` process takes "${pair_id}_strand_check.txt" file  
 * and executes commands according to the reads strandedness. 
 * The output generated is FASTA files.  
 */
 
process FINDTAIL {

    input:
    path "${pair_id}_strand_check.txt"
    path tuple val(pair_id), path(reads)

    output:
    

    script:
    """
    #!/usr/bin/env python3

    import subprocess
    import pandas as pd
    import os


    result = pd.read_csv("${pair_id}_strand_check.txt", sep="\r\n", header=None, engine='python')

    failed = float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', ''))
    fwd = float(result.iloc[2,0].replace('Fraction of reads explained by "1++,1--,2+-,2-+": ', ''))
    rev = float(result.iloc[3,0].replace('Fraction of reads explained by "1+-,1-+,2++,2--": ', ''))
    fwd_percent = fwd/(fwd+rev)
    rev_percent = rev/(fwd+rev)

    if float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', '')) > 0.50:
        cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
        cmd2 = "findtail_v1.01 --input_file ${reads[1]}  --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
        print(cmd1)
        subprocess.call(cmd1, shell=True)
        print(cmd2)
        subprocess.call(cmd2, shell=True)

        if fwd_percent > 0.9:
            #Over 90% of reads explained by "1++,1--,2+-,2-+
            #Data is likely FR/fr-secondstrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif rev_percent > 0.9:
            # Over 90% of reads explained by "1+-,1-+,2++,2--
            # Data is likely RF/fr-firststrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif max(fwd_percent, rev_percent) < 0.6:
            #Under 60% of reads explained by one direction
            #Data is likely unstranded
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.sl.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.pr.fasta
            cmd3 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.sl.fasta
            cmd4 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.pr.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
            print(cmd3)
            subprocess.call(cmd3, shell=True)
            print(cmd4)
            subprocess.call(cmd4, shell=True)

        else:
            #if strand couldnt be detected
            # we assume it is first-stranded
            # since most illumina reads are
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
        """
}

/*
 * define the `ALIGN` process that outputs unmapped files in FASTA format
 * It will take the FASTA files creates by the `FINDTAIL` process 
 * and try to align them INDIVIDUALLY.
 */
 
process ALIGN {

    input:
    path index
    path "tail_${pair_id}*.fasta"

    output:
    path "unmapped_tail_${pair_id}*.fasta"

    script:
    """
    bowtie2 -x $index -f "tail_${pair_id}*.fasta" --un "unmapped_tail_${pair_id}*.fasta"
    """
}

After creating a "${pair_id}_strand_check.txt" file, I input it to the FINDTAIL process that finds reads with polyadenylated ends. It takes fastq files individually, so I have to run it individually for each file.

If fastq file pair is first stranded or second stranded, I executed findtail twice,once for R1 and once for R2, which creates two fasta files. But if it is unstranded, I have to run findtail four times, twice for each file. Thus, depending on the paired-end reads' strandedness, I have different number of inputs. How do I define the output of this process?

In the next process I align these files individually to get unmapped reads.

This issue has been eating my brain away for 5 days. So please let me know if something is unclear.

Best wishes!

Answer 1

You should be able to just use a glob pattern that accounts for all possible output files, regardless of the specified conditions. As long as a single file matches the glob pattern when the process completes (successfully), Nextflow won't complain. I think this should suffice, but please let me know in the comments if it doesn't:

output:
tuple val(pair_id), path("tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")