nextflow: unable to find index

Question

I've a pipeline with many processes.
One of the steps is to run lancet tool. Lancet tool runs per chromosome (1-22, x,y). The input in the lancet step is output from the applybqsr (gatk).

There are tumor and normal bam generated along with the index. human genome FASTA has a index. When I run lancet through the main.nf I run into an error

ERROR: index not found for BAM file

I've a config file as:

params {

        hg38genome="/sc/arion/reference_genome/human/hg38/Homo_sapiens_assembly38.fasta.gz"
        outdir = '/sc/arion/output_pipeline_nextflow/test_tiny_datasets/'
        dbsnp="/sc/arion/projects/Daphni2_Pipeline_resources/resourceDBsnp146/dbsnp_146.hg38.vcf.gz"
        target="/sc/arion/projects/Twist_Exome_Target_hg38.padded.interval_list.bed/Twist_Exome_Target_hg38.padded.interval_list.bed"

}
        withName: applybqsr {
        publishDir = [
                [  
                path: "${params.outdir}/applybqsr/normal/",
                mode: 'copy',
                pattern:'*normal*'
                ],
                [
                path: "${params.outdir}/applybqsr/tumor/",
                mode: 'copy',
                pattern:'*tumor*'
                ]       
        ]
        }
        //applybqsr ends        

        
        withName: lancet {
        publishDir = [  
                path: "${params.outdir}/lancet/",
                mode: 'copy'
        ]
        }
        //lancet ends

I have applybqsr as:

process applybqsr{

        maxForks 3
        debug true

        memory '6 GB'

        input:
        tuple val(sample_id_tumor), path(markbam_tumor)
        tuple val(sample_id_tumor), path(recaltable_tumor)
        tuple val(sample_id_normal), path(markbam_normal)
        tuple val(sample_id_normal), path(recaltable_normal)

        output:
        tuple val(sample_id_tumor), path("${sample_id_tumor}.sorted.markdup.recal.bam"), emit: recal_bam_tumor
        tuple val(sample_id_tumor), path("${sample_id_tumor}.sorted.markdup.recal.bai"),  emit: recal_bai_tumor
        tuple val(sample_id_normal), path("${sample_id_normal}.sorted.markdup.recal.bam"), emit: recal_bam_normal
        tuple val(sample_id_normal), path("${sample_id_normal}.sorted.markdup.recal.bai"),  emit: recal_bai_normal
        
        script:
        """
        ml gatk/4.1.3.0 samtools/1.9
        samtools index ${markbam_tumor}
        samtools index ${markbam_normal}

        gatk ApplyBQSR \\
        -R $params.hg38genome \\
        -I ${markbam_tumor} \\
        --bqsr-recal-file ${recaltable_tumor} \\
        -O ${sample_id_tumor}.sorted.markdup.recal.bam \\
        -L $params.target \\
        --create-output-bam-index true \\
        --java-options -Xmx${task.memory.toGiga()}g

        gatk ApplyBQSR \\
        -R $params.hg38genome \\
        -I ${markbam_normal} \\
        --bqsr-recal-file ${recaltable_normal} \\
        -O ${sample_id_normal}.sorted.markdup.recal.bam \\
        -L $params.target \\
        --create-output-bam-index true \\
        --java-options -Xmx${task.memory.toGiga()}g


        """

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

I've lancet.nf as:

process lancet {
        maxForks 3
        debug true

        input:
        tuple val(sample_id_tumor),path(bqsrbam_tumor)
        tuple val(sample_id_normal),path(bqsrbam_normal)
        val chromosome

        output:
        path "CHR${chromosome}_${sample_name}.lancet.vcf"

        script:

        sample_name=sample_id_normal.split('_normal')[0]

"""
ml lancet

lancet --tumor ${bqsrbam_tumor} --normal ${bqsrbam_normal} --ref $params.hg38genome --reg "chr"$chromosome --num-threads 8 > "CHR"${chromosome}"_${sample_name}.lancet.vcf" 
"""

}

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

I have main.nf as:

params.outdir_index="/sc/arion/output_pipeline_nextflow/test_tiny_datasets/bwa_index"
params.raw_files_normal="/sc/arion/tiny/all_tumors_normals/*_normal_*{R1,R2}.fastq"
params.raw_files_tumor="/sc/arion/tiny/all_tumors_normals/*_tumor_*{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 {setupnmdtags} from './setupnmdtags_process.nf'
include { recalibrate_bam } from './recalibratebam_process.nf'
include { applybqsr } from './applybqsr_process.nf'
include { mutect2 } from './mutect2_process.nf'
include { lancet } from './lancet_process.nf'

workflow {

        normal_read_pairs=Channel.fromFilePairs(params.raw_files_normal)
        tumor_read_pairs=Channel.fromFilePairs(params.raw_files_tumor)
        //FASTP   align_bwa_mem        gatk_markduplicates         setupnmdtags        recalibrate_bam 
        
        applybqsr(gatk_markduplicates.out.sorted_bam_mark_duplicates_tumor,recalibrate_bam.out.recalib_table_tumor,gatk_markduplicates.out.sorted_bam_mark_duplicates_normal,recalibrate_bam.out.recalib_table_normal)

        mutect2(applybqsr.out.recal_bam_tumor,applybqsr.out.recal_bam_normal)
        
        chr=Channel.of( 1..22, 'X', 'Y' )
        lancet(applybqsr.out.recal_bam_tumor,applybqsr.out.recal_bam_normal,chr)
}

Lancet runs fine if I run it hard coded into a process code block with same input (tumor-normal bam) and human genome FASTA. How do I enable index files to be staged into the lancet process generated one step before?

Answer 1

I am reluctant to go through this wall of code, but generally it is advisable to always keep bam files and its index in a tuple so the index is staged with its bam file.

For example a process that creates the bam and now emits a tuple with bam and bai file:

(...)
output:
tuple path("${id}.bam"), path("${id}.bam.bai"), emit: tuple_bam

and then a process that receives such a tuple_bam:

(...)
input:
tuple path(bam), path(bai)

The first entry is the bam, the second the index.