Nextflow中bwa_map仅处理单个双端样本的问题排查与解决

刚接触Nextflow编程，运行流程时发现fastqc_trimmed步骤可处理所有双端样本，但bwa_map步骤仅处理一个样本。以下是相关代码、输入信息及解决方案：

流程代码 #

params.reads="/path/to/data/*_{1,2}_subsample.fastq"
params.outdir = "/path/to/results"
params.datadir = "/path/to/data"
params.genome = "/path/to/data/genome.fasta"

reads_ch = channel.fromFilePairs( params.reads, checkIfExists: true )
genome_ch = Channel.fromPath(params.genome, checkIfExists: true )

process fastqc_raw {

    tag "FASTQC on $sample_id"
    conda 'bioconda::fastqc=0.11.9'
    publishDir "$params.outdir/Raw_data_total/fastQC/$sample_id"

    input:
        tuple val(sample_id), path(reads)

    output:
      file("*.{html,zip}")

    script:
    """
    fastqc -t 4 ${reads}
    """
}

process trim_raw {

    tag "trimmomatic on $sample_id"
    conda 'bioconda::trimmomatic=0.39'
    publishDir "$params.outdir/Cleaned_data/$sample_id"

    input:
        tuple val(sample_id), path(reads)

    output:
        tuple val(sample_id), file("*.trimmed_Q20.fastq.gz")

    script:
    """
    trimmomatic PE -phred33 -threads 4 ${reads[0]} ${reads[1]} ${sample_id}_1.trimmed_Q20.fastq.gz output_forward_unpaired.fq.gz ${sample_id}_2.trimmed_Q20.fastq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:/trinity/shared/apps/Trimmomatic-0.39/adapters/Total_adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:75
    """
}

process fastqc_trimmed {

    tag "FASTQC on trimmed $sample_id"
    conda 'bioconda::fastqc=0.11.9'
    publishDir "$params.outdir/Cleaned_data/$sample_id"

    input:
        tuple val(sample_id), path(reads)

    output:
      file("*.{html,zip}")

    script:
    """
    fastqc -t 4 ${reads}
    """
}

process bwa_build_bwt {

    tag "bwa-mem2 indexing on ${genome}"
    conda 'bioconda::bwa-mem2=2.2.1'
    publishDir "$params.datadir"

    input:
        val genome

    output:
      file("*.0123")
      file("*.bwt.2bit.64")
      file("*.amb")
      file("*.ann")
      file("*.pac")

    script:
    def prefix = task.ext.prefix ?: "${genome.baseName}"
    """
    bwa-mem2 index -p ${prefix}.fasta ${genome}
    """

}

process samtools_faidx {

    tag "samtools faidx on ${genome}"
    conda 'bioconda::samtools=1.15.1'
    publishDir "$params.datadir"

    input:
        val genome

    output:
      file("*.fai")

    script:
    def prefix = task.ext.prefix ?: "${genome.baseName}"
    """
    samtools faidx --fai-idx ${prefix}.fasta.fai ${genome}
    """

}

// ne traite qu'un échantillon voir pourquoi

process bwa_map {

    tag "BWA map $sample_id"
    conda 'bioconda::bwa-mem2=2.2.1 bioconda::samtools=1.15.1'
    publishDir "$params.outdir/Mapped_bam/"

    input:
        tuple val(sample_id), path(reads)
        val bwa_build_bwt
        val samtools_faidx
        val genome

    output:
      tuple val(sample_id), file("*_BWA.bam")

    script:
    rg="@RG\tID:${sample_id}\tPL:ILLUMINA\tSM:${sample_id}"
    """
    bwa-mem2 mem -t 4 -R '${rg}' ${genome} ${reads[0]} ${reads[1]} | samtools view -Sb - > ${sample_id}_BWA.bam
    """
}

workflow {
    fastqc_raw(reads_ch)
    trim_raw(reads_ch)
    fastqc_trimmed(trim_raw.out)
    bwa_build_bwt(genome_ch)
    samtools_faidx(genome_ch)
    bwa_map(trim_raw.out, bwa_build_bwt.out[0], samtools_faidx.out, genome_ch)

}

bwa_map输入信息 #

genome_ch : /path/to/data/genome.fasta
bwa_build_bwt.out[0] : /path/to/genome.fasta.0123
samtools_faidx.out : /path/to/genome.fasta.fai

trim_raw.out :

[ERR8014964, [/path/to/ERR8014964_1.trimmed_Q20.fastq.gz, /path/to/ERR8014964_2.trimmed_Q20.fastq.gz]]
[ERR8014965, [/path/to/ERR8014965_1.trimmed_Q20.fastq.gz, /path/to/ERR8014965_2.trimmed_Q20.fastq.gz]]

问题原因与解决方案 #

原因 #

普通Nextflow channel只能被消费一次，而genome_ch被bwa_build_bwt、samtools_faidx和bwa_map三个流程重复调用，导致bwa_map仅能获取到genome_ch的最后一次消费结果，因此只处理了一个样本。

解决方法 #

将genome_ch转换为value channel，通过添加.first()方法实现。value channel存储单个值，可被多个流程或操作符重复消费。修改后的代码如下：

Channel
    .fromPath(params.genome,
              checkIfExists: true)
    .first()
    .set { genome_ch }

内容的提问来源于stack exchange，提问作者user3383071