Using shell command in Nextflow DSL=2 for BWA mem aligner

Question

[rohitsatyam102]

Hi Everyone!!

I am trying to write a nextflow script for bwa alignment, however, I am facing problem in adding read group in bwa mem command. The sample code and error is given below:

1. Using shell block

nextflow.enable.dsl=2
params.raw = "data/*{1,2}.fastq.gz"
params.ref = "results/00_indexes/bwaidx"
params.outdir="results/04_alignments"

process BWAMEM{
    publishDir "$params.outdir/bwa.alignments"
    input:
    tuple val(sample_id), path(reads)
    output:
    path '${sample_id}.sorted.bam'
    shell:
    '''
    id=$(zcat !{reads[0]} | head -n 1 | cut -f 3-4 -d":" | sed 's/@//')
    echo "$id"
    bwa mem -M -R "$(echo "@RG\\tID:${id}\\tSM:!{sample_id}\\tPL:ILLUMINA")" -t 8 !{params.ref} !{reads[0]} !{reads[1]} | samtools sort -@8 -o !{sample_id}.sorted.bam -
    '''
}

reads_ch = Channel.fromFilePairs(params.raw, checkIfExists: true )

workflow {
  BWAMEM(reads_ch)
}

################# Error ###########################
Error executing process > 'BWAMEM (1)'

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

Command executed:

  id=$(zcat S9_1.fastq.gz | head -n 1 | cut -f 3-4 -d":" | sed 's/@//')
  echo "$id"
  bwa mem -M -R "$(echo "@RG\tID:${id}\tSM:S9\tPL:ILLUMINA")" -t 8 results/00_indexes/bwaidx S9_1.fastq.gz S9_2.fastq.gz | samtools sort -@8 -o S9.sorted.bam -

Command exit status:
  1

Command output:
  000000000-JLT9H:1

Command error:
  [E::bwa_idx_load_from_disk] fail to locate the index files
  samtools sort: failed to read header from "-"

Work dir:
  /home/subudhak/Documents/nxtflow/work/73/7f89163a91a197851e8b84284c192b

Tip: you can replicate the issue by changing to the process work dir and entering the command `bash .command.run`

Even usage of \ didn't help.

Answer 1

[Pierre Lindenbaum]

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

use a full path for params.ref : /full/path/to/results/00_indexes/bwaidx

[rohitsatyam102]

Also why do you think that $id isn't getting substituted by anything in read group?

[rohitsatyam102]

Thanks. Adding $baseDir worked. Here is the final code:

nextflow.enable.dsl=2
params.raw = "data/*{1,2}.fastq.gz"
params.ref = "$baseDir/results/00_indexes/bwaidx/sequence.fasta"
params.outdir="results/04_alignments"

process BWAMEM{
    publishDir "$params.outdir/bwa.alignments"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "${sample_id}.sorted.bam"
    shell:
    '''
    id=$(zcat !{reads[0]} | head -n 1 | cut -f 3-4 -d":" | sed 's/@//')
    echo "$id"
    bwa mem -M -R "$(echo "@RG\\tID:${id}\\tSM:!{sample_id}\\tPL:ILLUMINA")" -t 8 !{params.ref} !{reads[0]} !{reads[1]} | samtools sort -@8 -o !{sample_id}.sorted.bam -
    '''
}

reads_ch = Channel.fromFilePairs(params.raw, checkIfExists: true )

workflow {
  BWAMEM(reads_ch)
}

Answer 2

[colindaven]

If I were you, I'd write the complicated RG stuff in a separate script section of your process.

That way, you're only dealing with the complexities of groovy variables, rather than groovy and bash vars (which I find confusing).

You can then just println $yourVar to print without the complexities of the workflow.

example script section:

input:
file bam
file bai


output:
file "plots"
path "*.calmd_cov_window.txt", emit: window_txt


script:
prefix = bam.name.toString().tokenize('.').get(0)
name = bam

"""
cp -R ${params.WOCHENENDE_DIR}/plots/ .
cp -R ${params.WOCHENENDE_DIR}/scripts/ .
cp scripts/*.sh .
bash runbatch_sambamba_depth.sh


"""

[rohitsatyam102]

Thanks for the suggestion colindaven but I am newbie to nextflow. I have yet to learn what does .tokenize and emit does. I will revisit your solution, when I learn advance stuff.