HaplotypeCaller Process Only Processes One File

[zbaida-achraf]

Hi, I’ve created this pipeline for detecting SNPs.

Important

The pipeline works correctly for the majority of the processes (such as Bowtie2, SAMtoBAM, AddReadGroup, etc.), but when it reaches the HaplotypeCaller process, it only processes one file and does not proceed with the others. The HaplotypeCaller seems to stop after handling one sample (SRR4420293), and the subsequent samples (SRR4420294 and others) are not processed.

#!/usr/bin/env nextflow

nextflow.enable.dsl=2

// Define parameters with default values
params.samplesheet = "./samplesheet.csv"
params.index_bowtie = "/data/SNP_scenario/requirements/index_bowtie/genome"
params.genome_fasta = "/data/SNP_scenario/requirements/GRCh38.primary_assembly.genome.fa"
params.genome_dict = "/data/SNP_scenario/requirements/GRCh38.primary_assembly.genome.dict"
params.dbSNP = "/data/SNP_scenario/requirements/00-All_chr.vcf.gz"
params.threads = 32
params.outdir = "results"

// Input validation with detailed error messages
def checkPathExists(path, description) {
    if (!file(path).exists()) {
        error "ERROR: ${description} not found at: ${path}"
    }
}

// Input validation
checkPathExists(params.samplesheet, "Sample sheet file")
checkPathExists(params.genome_fasta, "Reference genome")
checkPathExists(params.dbSNP, "dbSNP file")

// Ensure dbSNP index file exists
checkPathExists("${params.dbSNP}.tbi", "dbSNP index file")

// Create channels for reference files
genome_fasta = Channel.fromPath(params.genome_fasta, checkIfExists: true)
dbsnp = Channel.fromPath(params.dbSNP, checkIfExists: true)
dbsnp_idx = Channel.fromPath("${params.dbSNP}.tbi", checkIfExists: true)
fasta_index = Channel.fromPath("${params.genome_fasta}.fai", checkIfExists: true)
genome_dict = Channel.fromPath(params.genome_dict, checkIfExists: true)


// Modified channel creation to match your samplesheet format
samplesheet_ch = Channel
    .fromPath(params.samplesheet)
    .splitCsv(header: true)
    .map { row ->
        def fastq1 = file(row.fastq_1)  // Changed from fastq1 to fastq_1
        def fastq2 = file(row.fastq_2)  // Changed from fastq2 to fastq_2
        if (!fastq1.exists()) error "FASTQ R1 file not found: ${row.fastq_1}"
        if (!fastq2.exists()) error "FASTQ R2 file not found: ${row.fastq_2}"
        tuple(row.sample_name, fastq1, fastq2)
    }

// Define processes
process Cutadapt {
    tag "${sample_name}"
    publishDir "${params.outdir}/trimmed", mode: 'copy'  // Output to trimmed directory
    
    input:
    tuple val(sample_name), path(fastq1), path(fastq2)

    output:
    tuple val(sample_name), path("${sample_name}_1.trimmed.fastq.gz"), path("${sample_name}_2.trimmed.fastq.gz"), emit: trimmed_reads

    script:
    """
    echo "FASTQ1: ${fastq1}"
    echo "FASTQ2: ${fastq2}"
    cutadapt -a AACTGTAGGCACCATCAAT -o ${sample_name}_1.trimmed.fastq.gz -p ${sample_name}_2.trimmed.fastq.gz ${fastq1} ${fastq2} -q 20 -m 16 -j ${params.threads} > ${sample_name}_cutadapt.log
    """
}

process Bowtie2 {
    tag "${sample_name}"
    publishDir "${params.outdir}/aligned", mode: 'copy'
    
    input:
    tuple val(sample_name), path(fastq1), path(fastq2)

    output:
    tuple val(sample_name), path("${sample_name}.sam"), emit: aligned_sam
    path "${sample_name}.bowtie2.log", emit: logs

    script:
    """
    echo "Using Bowtie2 index: ${params.index_bowtie}"
    echo "FASTQ1: ${fastq1}"
    echo "FASTQ2: ${fastq2}"
    bowtie2 -x ${params.index_bowtie} \
        -1 ${fastq1} \
        -2 ${fastq2} \
        -S ${sample_name}.sam \
        -p ${params.threads} 2> ${sample_name}.bowtie2.log
    """
}

process SAMtoBAM {
    tag "${sample_name}"
    publishDir "${params.outdir}/bam", mode: 'copy'  // Output to bam directory
    
    input:
    tuple val(sample_name), path(sam)

    output:
    tuple val(sample_name), path("${sample_name}.sorted.bam"), path("${sample_name}.sorted.bam.bai"), emit: sorted_bam

    script:
    """
    samtools view -@ ${params.threads} -bS ${sam} > ${sample_name}.bam
    samtools sort -@ ${params.threads} ${sample_name}.bam -o ${sample_name}.sorted.bam
    samtools index -@ ${params.threads} ${sample_name}.sorted.bam
    """
}

process AddReadGroup {
    tag "${sample_name}"
    publishDir "${params.outdir}/bam_rg", mode: 'copy'  // Output to bam_rg directory
    
    input:
    tuple val(sample_name), path(bam), path(bai)

    output:
    tuple val(sample_name), path("${sample_name}.rg.bam"), path("${sample_name}.rg.bam.bai"), emit: rg_bam

    script:
    """
    samtools addreplacerg -@ ${params.threads} \
        -r '@RG\\tID:${sample_name}\\tSM:${sample_name}\\tPL:illumina' \
        -o ${sample_name}.rg.bam \
        ${bam}
    samtools index -@ ${params.threads} ${sample_name}.rg.bam
    """
}

process HaplotypeCaller {
    tag "${sample_name}"
    publishDir "${params.outdir}/vcf", mode: 'copy'  // Output to vcf directory
    
    input:
    tuple val(sample_name), path(bam), path(bai)
    path genome_fasta

    output:
    tuple val(sample_name), path("${sample_name}.vcf"), path("${sample_name}.vcf.idx"), emit: raw_vcf

    script:
    """
    ln -s ${params.genome_fasta}.fai .
    ln -s ${params.genome_dict} .

    gatk HaplotypeCaller \
        -R ${genome_fasta} \
        -I ${bam} \
        -O ${sample_name}.vcf \
        --native-pair-hmm-threads ${params.threads}
    """
}

process SelectVariants {
    tag "${sample_name}"
    publishDir "${params.outdir}/vcf_snps", mode: 'copy'  // Output to vcf_snps directory
    
    input:
    tuple val(sample_name), path(vcf), path(idx)
    path genome_fasta

    output:
    tuple val(sample_name), path("${sample_name}.snps.vcf"), path("${sample_name}.snps.vcf.idx"), emit: snp_vcf

    script:
    """
    ln -s ${params.genome_fasta}.fai .
    ln -s ${params.genome_dict} .
   
    gatk SelectVariants \
        -R ${genome_fasta} \
        -V ${vcf} \
        --select-type-to-include SNP \
        -O ${sample_name}.snps.vcf
    """
}

process VariantFiltration {
    tag "${sample_name}"
    publishDir "${params.outdir}/vcf_filtered", mode: 'copy'  // Output to vcf_filtered directory
    
    input:
    tuple val(sample_name), path(vcf), path(idx)
    path genome_fasta

    output:
    tuple val(sample_name), path("${sample_name}.filtered.vcf"), path("${sample_name}.filtered.vcf.idx"), emit: filtered_vcf

    script:
    """
    ln -s ${params.genome_fasta}.fai .
    ln -s ${params.genome_dict} .
    
    gatk VariantFiltration \
        -R ${genome_fasta} \
        -V ${vcf} \
        --filter-expression "QUAL < 30.0" \
        --filter-name "low_qual" \
        -O ${sample_name}.filtered.vcf
    """
}

process ExcludeFilteredVariants {
    tag "${sample_name}"
    publishDir "${params.outdir}/vcf_final", mode: 'copy'  // Output to vcf_final directory
    
    input:
    tuple val(sample_name), path(vcf), path(idx)
    path genome_fasta

    output:
    tuple val(sample_name), path("${sample_name}.final.vcf"), path("${sample_name}.final.vcf.idx"), emit: final_vcf

    script:
    """
    ln -s ${params.genome_fasta}.fai .
    ln -s ${params.genome_dict} .

    gatk SelectVariants \
        -R ${genome_fasta} \
        -V ${vcf} \
        --exclude-filtered \
        -O ${sample_name}.final.vcf
    """
}

process VariantAnnotator {
    tag "${sample_name}"
    publishDir "${params.outdir}/vcf_annotated", mode: 'copy'  // Output to vcf_annotated directory
    
    input:
    tuple val(sample_name), path(vcf), path(idx)
    path genome_fasta
    path dbsnp
    path dbsnp_idx

    output:
    tuple val(sample_name), path("${sample_name}.annotated.vcf"), path("${sample_name}.annotated.vcf.idx")

    script:
    """
    ln -s ${params.genome_fasta}.fai .
    ln -s ${params.genome_dict} .
   
    gatk VariantAnnotator \
        -R ${genome_fasta} \
        -V ${vcf} \
        -O ${sample_name}.annotated.vcf \
        --dbsnp ${dbsnp}
    """
}


process FilterSNPs {
    tag "${sample_name}"
    publishDir "${params.outdir}/vcf_filtered", mode: 'copy'  // Output to vcf_filtered directory

    input:
    tuple val(sample_name), path(vcf), path(idx)

    output:
    tuple val(sample_name), path("${sample_name}.filtered_snps.vcf"), emit: filtered_snps_vcf

    script:
    """
    bcftools view -m2 -M2 -v snps ${vcf} -o ${sample_name}.filtered_snps.vcf --threads ${params.threads}
    """
}

process ConvertToPlink {
    tag "${sample_name}"
    publishDir "${params.outdir}/plink", mode: 'copy'  // Output to plink directory

    input:
    tuple val(sample_name), path(vcf)

    output:
    tuple val(sample_name), path("${sample_name}.bed"), path("${sample_name}.bim"), path("${sample_name}.fam"), emit: plink_files

    script:
    """
    plink2 --vcf ${vcf} \
        --allow-extra-chr \
        --not-chr X \
        --make-bed \
        --out ${sample_name} \
        --threads ${params.threads}
    """
}

workflow {
    // Validate input files existence
    checkPathExists(params.samplesheet, "Sample sheet file")
    checkPathExists(params.genome_fasta, "Reference genome")
    checkPathExists(params.dbSNP, "dbSNP file")
    
    // Main workflow execution
    Cutadapt(samplesheet_ch)
    Bowtie2(Cutadapt.out.trimmed_reads)
    SAMtoBAM(Bowtie2.out.aligned_sam)
    AddReadGroup(SAMtoBAM.out.sorted_bam)
    HaplotypeCaller(AddReadGroup.out.rg_bam, genome_fasta)
    SelectVariants(HaplotypeCaller.out.raw_vcf, genome_fasta)
    VariantFiltration(SelectVariants.out.snp_vcf, genome_fasta)
    ExcludeFilteredVariants(VariantFiltration.out.filtered_vcf, genome_fasta)
    VariantAnnotator(
        ExcludeFilteredVariants.out.final_vcf,
        genome_fasta,
        dbsnp,
        dbsnp_idx
    )
    FilterSNPs(VariantAnnotator.out)
    ConvertToPlink(FilterSNPs.out.filtered_snps_vcf)
}

Results:

Launching pipeline.nf [loving_bhabha] DSL2 - revision: b0f6a0cd2c Connection to 3.77.126.214 closed. executor > local (15) DELL Precision@DESKTOP-PH6B8O3 MINGW64 /c/Perso/VM [100%] 2 of 2 ✔

$ e/324cd2] process > Bowtie2 (SRR4420294) [100%] 2 of 2 ✔

[c8/b8c876] process > SAMtoBAM (SRR4420294) [100%] 2 of 2 ✔

[f2/9ca3d0] process > AddReadGroup (SRR4420294) [100%] 2 of 2 ✔

[26/5e0941] process > HaplotypeCaller (SRR4420293) [100%] 1 of 1 ✔

[36/71db6d] process > SelectVariants (SRR4420293) [100%] 1 of 1 ✔

[cb/b190eb] process > VariantFiltration (SRR4420293) [100%] 1 of 1 ✔

[c6/e05d44] process > ExcludeFilteredVariants (SRR4420293) [100%] 1 of 1 ✔

[e6/b38f1d] process > VariantAnnotator (SRR4420293) [100%] 1 of 1 ✔

[0f/c3c7c8] process > FilterSNPs (SRR4420293) [100%] 1 of 1 ✔

[41/b84233] process > ConvertToPlink (SRR4420293) [100%] 1 of 1 ✔

Completed at: 29-Jan-2025 15:27:28
Duration : 1h 10m 17s
CPU hours : 2.3
Succeeded : 15

[mcolpus]

Just a guess, but you might need the genome reference channels to be value channels. Try adding a .first() to the channels after // Create channels for reference files

[zbaida-achraf]

Even added that he still has the same problem.

// Create channels for reference files
genome_fasta = Channel.fromPath(params.genome_fasta, checkIfExists: true).first()
dbsnp = Channel.fromPath(params.dbSNP, checkIfExists: true).first()
dbsnp_idx = Channel.fromPath("${params.dbSNP}.tbi", checkIfExists: true).first()
fasta_index = Channel.fromPath("${params.genome_fasta}.fai", checkIfExists: true).first()
genome_dict = Channel.fromPath(params.genome_dict, checkIfExists: true).first()