3
$\begingroup$

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!

$\endgroup$
2
  • 1
    $\begingroup$ From the (now removed) edits, the updated path glob pattern could just be: path("tail_${pair_id}_{1,2}.{pr,sl}.fasta"). If it's still not working, check to see what files are actually being created. They might be being created in another directory or not at all. $\endgroup$
    – Steve
    Dec 25, 2022 at 12:23
  • 1
    $\begingroup$ I ran it the same thing again and it worked somehow. So gathered it was fault on my end and remove the edit. $\endgroup$
    – pubsurfted
    Dec 26, 2022 at 6:20

1 Answer 1

4
$\begingroup$

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")
$\endgroup$
1
  • 1
    $\begingroup$ Thank you so much! Merry Christmas and Happy Holidays! $\endgroup$
    – pubsurfted
    Dec 24, 2022 at 10:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.