<p>I think you have a design problem here. A better way would be perform this operation in your <em>workflow</em> block. This would allow you to feed in the tumor and normal BAMs (and avoid accessing files outside of the working directory as previously discussed). You can use the <a href="https://www.nextflow.io/docs/latest/operator.html#map" rel="nofollow noreferrer"><code>map</code></a> operator to extract the patient or case id from the sample names and the <a href="https://www.nextflow.io/docs/latest/operator.html#filter" rel="nofollow noreferrer"><code>filter</code></a> operator to separate the samples into two channels. You could also use the <a href="https://www.nextflow.io/docs/latest/operator.html#branch" rel="nofollow noreferrer"><code>branch</code></a> operator for this like in the example below. The <a href="https://www.nextflow.io/docs/latest/operator.html#join" rel="nofollow noreferrer"><code>join</code></a> operator can then be used to join together the items emitted in a consistent way. For example:</p>
<pre><code>workflow {

    ref_fasta = file( params.ref_fasta )
    targets = file( params.targets )

    ...

    applybqsr.out.recal_bam
        .branch { sample_id, bam -&gt;

            def case_id = sample_id.substring(0, sample_id.lastIndexOf('_'))

            tumor: sample_id.toUpperCase().endsWith('_T')
                return tuple( case_id, sample_id, bam )

            normal: sample_id.toUpperCase().endsWith('_N')
                return tuple( case_id, sample_id, bam )

        }
        .set { cases_ch }

    joined_ch = cases_ch.tumor.join( cases_ch.normal )

    mutect2( joined_ch, ref_fasta, targets )
}
</code></pre>
<pre><code>process mutect2 {

    tag { sample_id }

    module 'gatk/4.1.3.0'

    input:
    tuple val(case_id), val(tumor), path(tumor_bam), val(normal), path(normal_bam)
    path reference_fasta
    path intervals_bed

    output:
    tuple val(case_id), path(&quot;${case_id}.mutect2.vcf&quot;)

    script:
    def avail_mem = task.memory ? task.memory.toGiga() : 0
    def java_options = [
        avail_mem ? &quot;-Xmx${avail_mem}G&quot; : &quot;&quot;,
        &quot;-XX:+UseSerialGC&quot;,
    ]

    &quot;&quot;&quot;
    gatk \\
        --java-options &quot;<span class="math-container">${java_options.join(' ')}" \\
    Mutect2 \\
        -R "$</span>{reference_fasta}&quot; \\
        -L &quot;<span class="math-container">${intervals_bed}" \\
        -I "$</span>{tumor_bam}&quot; \\
        -I &quot;<span class="math-container">${normal_bam}" \\
        -normal "$</span>{normal}&quot; \\
        -tumor &quot;<span class="math-container">${tumor}" \\
        -O "$</span>{case_id}&quot;.mutect2.vcf&quot;
    &quot;&quot;&quot;
}
</code></pre>
<p>There's also likely an issue with how <em>applybqsr</em> is invoked. If you require <a href="https://www.nextflow.io/docs/latest/process.html#process-multiple-input-channels" rel="nofollow noreferrer">multiple input channels</a>, most of the time what you'll want is one <a href="https://www.nextflow.io/docs/latest/channel.html#queue-channel" rel="nofollow noreferrer">queue channel</a> and one or more <a href="https://www.nextflow.io/docs/latest/channel.html#value-channel" rel="nofollow noreferrer">value channels</a>. I suspect the inputs for <em>applybqsr</em> also need to be joined into a single queue channel.</p>

CC BY-SA 4.0
2023-07-13T13:39:35.497Z
I really don't understand what you did with `.branch` in applybqsr. Nextflow is so mysterious. How are you defining `normal:` and `tumor:` without define keyword? Also, I have samples as: MM_tiny_n_R1.fastq.gz MM_tiny_t_R1.fastq.gz Thank you for your reply.
Death Metal
CC BY-SA 4.0
2023-07-13T13:41:16.590Z
Also, where is `cases.tumor` stored/created?
Death Metal
CC BY-SA 4.0
2023-07-13T14:16:19.967Z
@DeathMetal Sorry that should have been `cases_ch.tumor`. Fixed now. The boolean expressions we define in our branch closure can be given unique labels, which we can then use as named channels.
Steve
CC BY-SA 4.0
2023-07-13T14:19:44.750Z
@DeathMetal With your current glob pattern, the `fromFilePairs` factory method should give you the sample names: `MM_tiny_n` and `MM_tiny_t`. I assumed your _FASTP_ and dowstream processes also return tuples containing the sample names/ids, but I might be wrong..
Steve
CC BY-SA 4.0
2023-07-13T18:32:55.077Z
This is complicated for me, sorry. Would you mind explaining  in `tumor: sample_id.endsWith('_T')  return tuple( case_id, sample_id, bam )` what is **case_id**? Second, how does tuple has three elements in it?
Death Metal
CC BY-SA 4.0
2023-07-13T18:33:34.980Z
If you get a chance, I've a new question at https://bioinformatics.stackexchange.com/questions/21310/nextflow-missing-value-declared-as-output-parameter-output-variable-defined
Death Metal
CC BY-SA 4.0
2023-07-14T00:55:49.507Z
@DeathMetal The _case_id_ is really just the patient or case or donor from which the samples are derived. For example, `MM_tiny` would be the case or donor id here. Tuples can contain any number of elements. The elements don't even need to be of the same type necessarily.
Steve
CC BY-SA 4.0
2023-07-13T01:28:23.290Z
2023-07-14T00:51:33.487Z
2023-07-14T00:51:33.487Z
Steve
<p>I'm developing pipeline in nextflow it uses tools: fastp, bwa mem index, bwamem, gatk mark duplicates, gatk setupnmd, gatk applybqsr, gatk recalibrate. I do this using FASTQ.gz of normal and tumor DNA samples.<br />
I'm interested to use mutect2, lancet and strelka in nextsteps. However, I am unable to write code in the pipeline that makes sense. I use loads of shell scripting in mutect2 to use tumor-normal.<br />
To expand, I've input in FASTQ as tiny_n, tiny_t, where _n is normal and _t is tumor. I run pipeline until applybqsr for each individually. However, as I reach mutect2, I check if sample contains <strong>_n</strong> or otherwise, I then create an empty file for the output as required in the mutect2 process.</p>
<p>How do I overcome creating dummy/empty file?</p>
<p>Please see below code of <strong>main.nf</strong>:</p>
<pre><code>params.raw_files=&quot;/sc/arion/projects/tiny/all_tumors_normals/tiny_*{R1,R2}.fastq&quot;

include { FASTP} from './fastp_process.nf'
//include {bwa_index} from './index_process.nf'
include { align_bwa_mem} from './bwamem_process_already_index.nf'
include { gatk_markduplicates} from './gatk_markduplicates_process.nf'
include { recalibrate_bam } from './recalibratebam_process.nf'
include {setupnmdtags} from './setupnmdtags_process.nf'
include { applybqsr } from './applybqsr_process.nf'
include { mutect2 } from './mutect2_process.nf'
workflow {

        read_pairs_ch = Channel.fromFilePairs(params.raw_files) 

        FASTP(read_pairs_ch) 

        align_bwa_mem(FASTP.out.reads) 
        gatk_markduplicates(align_bwa_mem.out.sorted_bams)
        setupnmdtags(gatk_markduplicates.out.sorted_bam_mark_duplicates)
        recalibrate_bam(setupnmdtags.out.setupnmdtags_bam)

applybqsr(gatk_markduplicates.out.sorted_bam_mark_duplicates,recalibrate_bam.out.recalib_table)


mutect2(applybqsr.out.recal_bam)
}
</code></pre>
<p>below is code of mutect2 to handle only one sample, run for normal, create dummy/empty file for tumor sample, <strong>mutect2_process.nf</strong>:</p>
<pre><code>process mutect2 {
        maxForks 3
        debug true
        tag &quot;${sample_id}&quot;
 memory '15 GB'

        input:
        tuple val(sample_id),path(bqsrbam)
        
        output:
        tuple val(sample_id),path(&quot;${sample_id}.mutect2out.vcf&quot;)
        

        script:
                         def avail_mem = task.memory ? task.memory.toGiga() : 0
                def java_options = [
        avail_mem ? &quot;-Xmx${avail_mem}G&quot; : &quot;&quot;,
        &quot;-XX:+UseSerialGC&quot;,
    ]

        &quot;&quot;&quot;
        echo &quot;${sample_id}\n&quot;

        grep_value=\\\$(echo &quot;${sample_id}|  | grep &quot;_n&quot; &quot;)

        if [[ ! -z &quot;\\\${grep_value}&quot; ]]
        then 

                identifier=\\\$(echo ${sample_id} | awk -F&quot;_&quot; '{print \\\$1}')
                echo &quot;\\\${identifier}&quot;_t.sorted.markdup.recal.bam&quot;&quot;

                ml gatk/4.1.3.0
gatk --java-options &quot;<span class="math-container">${java_options.join(' ')}" Mutect2 -L $</span>params.target -I &quot;<span class="math-container">${params.outdir}/applybqsr/""\\\${identifier}"_t.sorted.markdup.recal.bam"" -I  "$</span>{params.outdir}/applybqsr/&quot;&quot;\\\${identifier}&quot;_n.sorted.markdup.recal.bam&quot;&quot; -normal \\\${identifier}&quot;_n&quot; -tumor \\\${identifier}&quot;_t&quot; --output <span class="math-container">${sample_id}".mutect2out.vcf" -R $</span>params.hg38genome

        else 
                echo &quot;Nothing to do\n&quot;
                touch ${sample_id}&quot;.mutect2out.vcf&quot;
 
        fi

        &quot;&quot;&quot;
        }


workflow.onComplete { 
        log.info ( workflow.success ? &quot;Done mutect2!&quot; : &quot;Oops .. something went wrong in mutect2&quot; )
}
</code></pre>
<p>There is a load of unix/shell code that does trivial checks. It runs gatk mutect2 oly when _n is found.
How do I overcome: output defined <code>tuple val(sample_id),path(&quot;${sample_id}.mutect2out.vcf&quot;)</code> ?<br />
Is there a better way to do this?</p>

<nextflow>
nextflow: tumor normal sample - how to make code organic
# nextflow: tumor normal sample - how to make code organic

(Asked by: username_0 on 2023-07-11)

I'm developing pipeline in nextflow it uses tools: fastp, bwa mem index, bwamem, gatk mark duplicates, gatk setupnmd, gatk applybqsr, gatk recalibrate. I do this using FASTQ.gz of normal and tumor DNA samples.
I'm interested to use mutect2, lancet and strelka in nextsteps. However, I am unable to write code in the pipeline that makes sense. I use loads of shell scripting in mutect2 to use tumor-normal.
To expand, I've input in FASTQ as tiny_n, tiny_t, where _n is normal and _t is tumor. I run pipeline until applybqsr for each individually. However, as I reach mutect2, I check if sample contains _n or otherwise, I then create an empty file for the output as required in the mutect2 process.
How do I overcome creating dummy/empty file?
Please see below code of main.nf:
params.raw_files="/sc/arion/projects/tiny/all_tumors_normals/tiny_*{R1,R2}.fastq"

include { FASTP} from './fastp_process.nf'
//include {bwa_index} from './index_process.nf'
include { align_bwa_mem} from './bwamem_process_already_index.nf'
include { gatk_markduplicates} from './gatk_markduplicates_process.nf'
include { recalibrate_bam } from './recalibratebam_process.nf'
include {setupnmdtags} from './setupnmdtags_process.nf'
include { applybqsr } from './applybqsr_process.nf'
include { mutect2 } from './mutect2_process.nf'
workflow {

        read_pairs_ch = Channel.fromFilePairs(params.raw_files) 

        FASTP(read_pairs_ch) 

        align_bwa_mem(FASTP.out.reads) 
        gatk_markduplicates(align_bwa_mem.out.sorted_bams)
        setupnmdtags(gatk_markduplicates.out.sorted_bam_mark_duplicates)
        recalibrate_bam(setupnmdtags.out.setupnmdtags_bam)

applybqsr(gatk_markduplicates.out.sorted_bam_mark_duplicates,recalibrate_bam.out.recalib_table)


mutect2(applybqsr.out.recal_bam)
}

below is code of mutect2 to handle only one sample, run for normal, create dummy/empty file for tumor sample, mutect2_process.nf:
process mutect2 {
        maxForks 3
        debug true
        tag "${sample_id}"
 memory '15 GB'

        input:
        tuple val(sample_id),path(bqsrbam)
        
        output:
        tuple val(sample_id),path("${sample_id}.mutect2out.vcf")
        

        script:
                         def avail_mem = task.memory ? task.memory.toGiga() : 0
                def java_options = [
        avail_mem ? "-Xmx${avail_mem}G" : "",
        "-XX:+UseSerialGC",
    ]

        """
        echo "${sample_id}\n"

        grep_value=\\\$(echo "${sample_id}|  | grep "_n" ")

        if [[ ! -z "\\\${grep_value}" ]]
        then 

                identifier=\\\$(echo ${sample_id} | awk -F"_" '{print \\\$1}')
                echo "\\\${identifier}"_t.sorted.markdup.recal.bam""

                ml gatk/4.1.3.0
gatk --java-options "${java_options.join(' ')}" Mutect2 -L $params.target -I "${params.outdir}/applybqsr/""\\\${identifier}"_t.sorted.markdup.recal.bam"" -I  "${params.outdir}/applybqsr/""\\\${identifier}"_n.sorted.markdup.recal.bam"" -normal \\\${identifier}"_n" -tumor \\\${identifier}"_t" --output ${sample_id}".mutect2out.vcf" -R $params.hg38genome

        else 
                echo "Nothing to do\n"
                touch ${sample_id}".mutect2out.vcf"
 
        fi

        """
        }


workflow.onComplete { 
        log.info ( workflow.success ? "Done mutect2!" : "Oops .. something went wrong in mutect2" )
}

There is a load of unix/shell code that does trivial checks. It runs gatk mutect2 oly when _n is found.
How do I overcome: output defined tuple val(sample_id),path("${sample_id}.mutect2out.vcf") ?
Is there a better way to do this?

**Tags:** nextflow

---

## 1 answers

### Answer by username_1 on 2023-07-13. Score: 3 | Accepted answer

I think you have a design problem here. A better way would be perform this operation in your workflow block. This would allow you to feed in the tumor and normal BAMs (and avoid accessing files outside of the working directory as previously discussed). You can use the map operator to extract the patient or case id from the sample names and the filter operator to separate the samples into two channels. You could also use the branch operator for this like in the example below. The join operator can then be used to join together the items emitted in a consistent way. For example:
workflow {

    ref_fasta = file( params.ref_fasta )
    targets = file( params.targets )

    ...

    applybqsr.out.recal_bam
        .branch { sample_id, bam ->

            def case_id = sample_id.substring(0, sample_id.lastIndexOf('_'))

            tumor: sample_id.toUpperCase().endsWith('_T')
                return tuple( case_id, sample_id, bam )

            normal: sample_id.toUpperCase().endsWith('_N')
                return tuple( case_id, sample_id, bam )

        }
        .set { cases_ch }

    joined_ch = cases_ch.tumor.join( cases_ch.normal )

    mutect2( joined_ch, ref_fasta, targets )
}

process mutect2 {

    tag { sample_id }

    module 'gatk/4.1.3.0'

    input:
    tuple val(case_id), val(tumor), path(tumor_bam), val(normal), path(normal_bam)
    path reference_fasta
    path intervals_bed

    output:
    tuple val(case_id), path("${case_id}.mutect2.vcf")

    script:
    def avail_mem = task.memory ? task.memory.toGiga() : 0
    def java_options = [
        avail_mem ? "-Xmx${avail_mem}G" : "",
        "-XX:+UseSerialGC",
    ]

    """
    gatk \\
        --java-options "${java_options.join(' ')}" \\
    Mutect2 \\
        -R "${reference_fasta}" \\
        -L "${intervals_bed}" \\
        -I "${tumor_bam}" \\
        -I "${normal_bam}" \\
        -normal "${normal}" \\
        -tumor "${tumor}" \\
        -O "${case_id}".mutect2.vcf"
    """
}

There's also likely an issue with how applybqsr is invoked. If you require multiple input channels, most of the time what you'll want is one queue channel and one or more value channels. I suspect the inputs for applybqsr also need to be joined into a single queue channel.

> **Comments:**
> **username_0**: I really don't understand what you did with `.branch` in applybqsr. Nextflow is so mysterious. How are you defining `normal:` and `tumor:` without define keyword? Also, I have samples as: MM_tiny_n_R1.fastq.gz MM_tiny_t_R1.fastq.gz Thank you for your reply.  
> **username_0**: Also, where is `cases.tumor` stored/created?  
> **username_1**: @DeathMetal Sorry that should have been `cases_ch.tumor`. Fixed now. The boolean expressions we define in our branch closure can be given unique labels, which we can then use as named channels.  
> **username_1**: @DeathMetal With your current glob pattern, the `fromFilePairs` factory method should give you the sample names: `MM_tiny_n` and `MM_tiny_t`. I assumed your _FASTP_ and dowstream processes also return tuples containing the sample names/ids, but I might be wrong..  
> **username_0**: This is complicated for me, sorry. Would you mind explaining  in `tumor: sample_id.endsWith('_T')  return tuple( case_id, sample_id, bam )` what is **case_id**? Second, how does tuple has three elements in it?  
> **username_0**: If you get a chance, I've a new question at https://bioinformatics.stackexchange.com/questions/21310/nextflow-missing-value-declared-as-output-parameter-output-variable-defined  
> **username_1**: @DeathMetal The _case_id_ is really just the patient or case or donor from which the samples are derived. For example, `MM_tiny` would be the case or donor id here. Tuples can contain any number of elements. The elements don't even need to be of the same type necessarily.  

---

