This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Is it possible to run haplotypecaller in gnu parallel ?

I have a bash script about HaplotypeCaller.sh. Is it possible to run this in gnu parallel ?

cat HaplotypeCaller.sh | parallel

#!/bin/bash

for i in *.bam; do 
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I $i -O $i.vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list; done
gnu parallel haplotypecaller parallel gatk
• 5,391 views • link
written by DareDevil DareDevil

Something like this should work

find . -name '*.bam' | rev | cut -c4- | parallel -I{} gatk --java-options "-Xmx40g" HaplotypeCaller \
    -R /media/gatk/Homo_sapiens_assembly38.fasta \
    -I {}.bam \
    -O {}.vcf.gz \
    --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz \
    -L /media/gatk/target/ag_V6_list.interval_list
• 0 views • link
written by 4galaxy77 4galaxy77

I get following error while running above command

Error was: .detros_1C100/..bam with exception: Cannot read non-existent file:
• 0 views • link
written by DareDevil DareDevil

@ 4galaxy77 rev | cut -c4- removes the .bam extension. But {.} in parallel removes the extension.

• 0 views • link
written by cpad0112 cpad0112

see if following prints exact commands you want to run and then remove --dry-run option

parallel --dry-run 'gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I {} -O {.}.vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list' ::: *.bam

make sure that system has enough resources to execute parallel commands. Other wise, limit memory / threads/ cpus.

• 0 views • link
written by cpad0112 cpad0112

3 answers 

#!/bin/bash

do_one() {
  gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I "$1" -O "$1".vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list
}
export -f do_one
# test that do_one works on a single bam file. Then:
parallel do_one ::: *.bam

If you do not want the output to be named foo.bam.vcf.gz but instead foo.vcf.gz:

#!/bin/bash

do_one() {
  gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I "$1" -O "$2" --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list
}
export -f do_one
# test that do_one works on a single bam file:
#   do_one foo.bam foo.vcf.gz
# Then:
parallel do_one {} {.} ::: *.bam
• 0 views • link
written by ole.tange ole.tange

This gives me an error export: Illegal option -f

• 0 views • link
written by DareDevil DareDevil

What is your shell? Something other than bash? In any case, see my other comment. You don't need a function at all.

• 0 views • link
written by Malcolm.Cook Malcolm.Cook

Are you writing exactly what I wrote? Because if you change #!/bin/bash to #!/bin/sh then it will not work.

• 0 views • link
written by ole.tange ole.tange

Ya I have taken the same code provided by you. It gives an error export: Illegal option -f

• 0 views • link
written by DareDevil DareDevil

As written above, the output files will have .bam.vcf.gz extension, where retaining the .bam was probably not desirable, and is why the OP's pipeline includes | cut -c4-.

I would not use a bash function in this case, and would rather write:

parallel gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I {} -O {.}.vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list ::: *.bam

Taking advantage of parallel's {.} variable which is the file without its last extension (in this case, trimming the trailing .bam).

• 0 views • link
written by Malcolm.Cook Malcolm.Cook

Cheers !! this works!!

• 0 views • link
written by DareDevil DareDevil

Update: When I executed your code one by one it started running. But When I saved your code in a file as parallel_haplo.sh it does not work.

sh parallel_haplo.sh It gives the error: parallel_haplo.sh : 6: export: Illegal option -f

It works when I run as bash parallel_haplo.sh. I don't know what made the difference

• 0 views • link
written by DareDevil DareDevil

nextflow (not tested)

params.ref="/media/gatk/Homo_sapiens_assembly38.fasta"
params.bams=""
params.dbsnp="/media/gatk/dbsnp_138.hg38.vcf.gz"
params.intervals="/media/gatk/target/ag_V6_list.interval_list"

Channel.fromPath(params.bams).splitCsv(header: false,sep:',',strip:true).map{T->file(T[0])}.set{bamfiles}

process hapCaller {
tag "${bam.name}"
memory "40g"
input:
   file bam from bamfiles
output:
   file("${bam.getBaseName()}.vcf.gz") into vcf
script:
"""
gatk --java-options " -Xmx${task.memory.giga}g -Djava.io.tmpdir=." HaplotypeCaller \
    -R ${params.ref} \
    -I {bam.toRealPath()} \
    -O "${bam.getBaseName()}.vcf.gz" \
    --dbsnp "${params.dbsnp}" \
    -L "${params.intervals}"
"""
}
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

The answer given by Malcolm.Cook, ole.tange and cpad0112 works well

Method 1:

parallel gatk --java-options "-Xmx40g" HaplotypeCaller \
-R /media/gatk/Homo_sapiens_assembly38.fasta -I {} -O {.}.vcf.gz \
 --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz \
-L /media/gatk/target/ag_V6_list.interval_list ::: *.bam

Method 2:

#!/bin/bash

do_one() {
  gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I "$1" -O "$1".vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list
}
export -f do_one
# test that do_one works on a single bam file. Then:
parallel do_one ::: *.bam
• 0 views • link
written by DareDevil DareDevil

Log in to answer this question.

More posts like this