Nextflow 生物信息学 snp 调用初学者问题：一个进程的一对文件输出作为另一个进程的输入

Question

我刚刚开始学习nextflow。我一直在关注此处提供的 rna seq 教程：https://training.nextflow.io/basic_training/rnaseq_pipeline/#quality-control 但有一些事情我无法在自己的工作流程中重新实现。我正在慢慢尝试从我经常使用的管道构建 SNP 调用工作流程，只是作为我学习如何使用 nextflow 的测试。

这是我的代码，它使用在nexflow.config文件中指定的两个容器。

params.reads = "/home/*_R{1,2}*.gz"
params.outdir = "./"
params.ref_genome = "/home/GCF_000001735.4_TAIR10.1_genomic.fasta"

log.info """\
    S N P   C A L L I N G   P I P E L I N E
    ===================================
    reference    : ${params.ref_genome}
    reads        : ${params.reads}
    outdir       : ${params.outdir}
    """
    .stripIndent(true)


/*
 * Processes
 */
process trimSequences {
    tag "Fastp on $sample_id"
    publishDir params.outdir, mode: 'copy'
    cpus 4
    debug true

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


    output:
    tuple val(sample_id), path("${sample_id}_{1,2}_trimmed.fastq.gz")

    script:
    """
    fastp -w $task.cpus -i ${reads[0]} -I ${reads[1]} -o ${sample_id}_1_trimmed.fastq.gz -O ${sample_id}_2_trimmed.fastq.gz
    """

}



process bwaIndex {
    publishDir params.outdir, mode: 'copy'

    input:
    path reference

    output:
    file "${reference.baseName}.fasta.amb"
    file "${reference.baseName}.fasta.ann"
    file "${reference.baseName}.fasta.bwt"
    file "${reference.baseName}.fasta.pac"
    file "${reference.baseName}.fasta.sa"

    script:
    """
    bwa index -a bwtsw $reference
    """
}

process bwaMap {
    publishDir params.outdir, mode: 'copy'
    cpus 4
    input:
    path reference
    tuple val(trimmed_reads), path(trimmed_ch)

    output:
    file "final.sam"

    script:
    """
    bwa mem -t $task.cpus $reference ${trimmed_reads[0]} ${trimmed_reads[1]} > final.sam
   """
}


/*
 * Define the workflow
 */
workflow {

    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs_ch }

   // Process read pairs through trimSequences
    trimmed_ch = trimSequences(read_pairs_ch)
    bwaIndex(params.ref_genome)
    bwaMap(params.ref_genome, trimmed_ch)
}

因此第一个过程“trimSequences”有效。这基本上需要目录 /home/ 中的样本的一对 fastq 文件，其中包含多个样本的 fastq 文件对，并修剪每一对，为输出提供我想要的正确输出名称。因此，对于每对输入文件，我都会生成一对输出文件。然后我将它们复制到当前目录中。

第二个过程，bwa索引起作用（bwaIndex）。但我认为我的代码一点也不优雅。在此过程结束时，我还复制输出（在当前工作目录中为我的参考而生成的 bwa 索引文件）。索引文件自动生成 通过函数 bwa index -a bwtsw，通过将后缀附加到参考 fasta 名称，我需要它们以这种格式用于下一步，与参考 fasta 位于同一位置。

主要问题是第三个进程，实际的映射进程，bwaMap。首先，它应该只在 bwaIndex 完成后启动。但是我不知道如何强制执行，因为 bwa mem 只需要引用的 fasta 文件作为输入（与 bwaIndex 进程相同），然后它根据与引用相同的名称加上后缀在同一目录中搜索索引文件。这也是为什么在最后的 bwaIndex 中我将当前目录中的索引文件与引用一起复制的原因。所以 bwaMap 找不到参考索引文件。

另一个问题是 bwaMap 没有正确获取 trimSequences 的输出。我在工作流程部分或 bwaMap 过程的输入和输出部分中指定错误。我尝试了多种方法，但找不到解决方案。当我在当前目录中有索引（手动将它们复制到那里）时，bwaMap进程失败，因为它没有正确地将修剪后的fastq对作为输入（表明没有正确读取修剪后的读取对，trimSequences的输出）。

然后我想做的就是简单地使用 bwa 来映射每对修剪的读数（每个样本的修剪 r1 和修剪 r2），trimSequences 的输出到参考基因组。

错误：

Caused by:
 Process `bwaMap (1)` terminated with an error exit status (1)

Command executed:

 bwa mem -t 1 GCF_000001735.4_TAIR10.1_genomic.fasta R A > final.sam

Command exit status:
 1

Command output:
 (empty)

Command error:
 [E::bwa_idx_load_from_disk] fail to locate the index files

Answer 1

所以 bwaIndex 问题——我希望 bwaMap 在开始之前等待 bwaIndex 完成，这很容易解决，这不是我发帖的原因，但确实造成了问题。我只是将 bwaIndex 的输出添加为流程规范中 bwaMap 的输入。然后在工作流程部分中，我简单地保存了 bwaIndex 的输出，并将其作为 bwaMap 的输入。

发布的主要问题和原因是我无法理解为什么 bwaMap 无法正确读取 trimSequences 的输出。我知道它没有正确读取，因为每次它抛出错误时都会报告 bwa mem 命令运行，并且我可以看到在错误中报告的命令中没有正确读取修剪后的读取名称。例如，在我在问题中报告的错误中，我有“R”和“A”，而不是修剪后的读取全名/路径。问题在于我在 bwaMap 进程中指定输入的方式。我修复了它，如下面的代码所示。现在一切正常了，但可能是以一种不太优雅的方式。仍然欢迎任何有关如何改进此语法的建议。我仍在寻找解决这个问题的方法，不确定我是否完全理解何时使用路径、文件和各种输入类型，但会到达那里。 我的解决方案：

#!/usr/bin/env nextflow

params.reads = "/home/*_R{1,2}*.gz"
params.outdir = "./"
params.ref_genome = "/home/GCF_000001735.4_TAIR10.1_genomic.fasta"

log.info """\
    S N P   C A L L I N G   P I P E L I N E
    ===================================
    reference    : ${params.ref_genome}
    reads        : ${params.reads}
    outdir       : ${params.outdir}
    """
    .stripIndent(true)


/*
 * Processes
 */
process trimSequences {
    tag "Fastp on $sample_id"
    cpus 4


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


    output:
    tuple val(sample_id), path("${sample_id}_{1,2}_trimmed.fastq.gz")

    script:
    """
    fastp -w $task.cpus -i ${reads[0]} -I ${reads[1]} -o ${sample_id}_1_trimmed.fastq.gz -O ${sample_id}_2_trimmed.fastq.gz
    """

}



process bwaIndex {

    input:
    path reference

    output:
    file "${reference.baseName}.fasta.amb"
    file "${reference.baseName}.fasta.ann"
    file "${reference.baseName}.fasta.bwt"
    file "${reference.baseName}.fasta.pac"
    file "${reference.baseName}.fasta.sa"

    script:
    """
    bwa index -a bwtsw $reference
    """
}

process bwaMap {

    publishDir params.outdir, mode: 'copy'
    cpus 4

    input:
    path reference
    file "${reference.baseName}.fasta.amb"
    file "${reference.baseName}.fasta.ann"
    file "${reference.baseName}.fasta.bwt"
    file "${reference.baseName}.fasta.pac"
    file "${reference.baseName}.fasta.sa"
    tuple val(sample_id), path(trimmed_reads)

    output:
    file "final.sam"

    script:
    """
    bwa mem -t $task.cpus $reference ${trimmed_reads[0]} ${trimmed_reads[1]} > final.sam
    """
}



/*
 * Define the workflow
 */
workflow {

    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs_ch }


    trimmed_ch = trimSequences(read_pairs_ch)
    bwa_index = bwaIndex(params.ref_genome)
    bwaMap(params.ref_genome, bwa_index, trimmed_ch)
}