how to get one vcf file for freebayes from many bam files

Question

[evelyn]

Hello,

I want to generate one vcf file out of many bam files. I used this command:

RUN=${SLURM_ARRAY_TASK_ID}
LIST_OF_BAMS=/path/list_of_bam.txt
freebayes  -f  /path/genome.fasta -L ${LIST_OF_BAMS} > /path/all_SNPs_pt${SLURM_ARRAY_TASK_ID}.vcf

But it gives several vcf files as per RUN and gives unknown in place of sample name in the header. I also want to include parameters:

--max-complex-gap         1 \
          --theta                   0.01 \
          --ploidy                  2 \
          --min-alternate-fraction  0.2 \
          --min-alternate-count     2 \
          --min-coverage            3 \
          --min-mapping-quality     1 \
          --min-base-quality        20 \
          --report-genotype-likelihood-max \
          --no-complex \
          --no-mnps \
          --no-indels \

But these seems to be outdated as it gives error if I include these parameters in the command line. Thank you!

[Ram]

You're running into a lot of problems combining SLURM variables and freebayes. Can you try not mixing them?

As a first step, instead of running the freebayes command, echo it:

RUN=${SLURM_ARRAY_TASK_ID}
LIST_OF_BAMS=/path/list_of_bam.txt
echo "freebayes  -f  /path/genome.fasta -L ${LIST_OF_BAMS} > /path/all_SNPs_pt${SLURM_ARRAY_TASK_ID}.vcf"
echo "freebayes  -f  /path/genome.fasta -L ${LIST_OF_BAMS} --max-complex-gap         1 \
          --theta                   0.01 \
          --ploidy                  2 \
          --min-alternate-fraction  0.2 \
          --min-alternate-count     2 \
          --min-coverage            3 \
          --min-mapping-quality     1 \
          --min-base-quality        20 \
          --report-genotype-likelihood-max \
          --no-complex \
          --no-mnps \
          --no-indels > /path/all_SNPs_pt${SLURM_ARRAY_TASK_ID}.vcf"

If the above echos the commands fine, then try remove the echos, add a set -x after your SLURM header lines and look at the exact freebayes commands being run - that should show you if where things fail.

[evelyn]

It says:

freebayes: unrecognized option '--no-complex'
did you mean --samples ?

[Ram]

Did the echo work fine? Is this what you got after removing the echo? Did you use set -x?

[evelyn]

Yes, echo worked fine and I got it after removing the echo.

[GenoMax]

Looks like version of freebayes you are using does not have --no-complex option. Run freebayes -h and take a look at possible options.

[Ram]

--no-complex is not a valid parameter. Are you looking for --throw-away-complex-obs? Similarly, --no-mnps and --no-indels are not valid. The closest I see are --throw-away-mnp-obs and --throw-away-indels-obs.

I found these by looking at freebayes -h (well, actually, I read through the CPP source code that prints that usage) - you could have read the manual yourself and learned this on your own.

Answer 1

[swbarnes2]

I don't think FreeBayes takes mutliple bams. Try making one giant bam with each sample having its own readgroup.

[h.mon]

And if you didn't add read groups when mapping, there is a small utility (bamaddrg) by the same main FreeBayes author which allows adding read groups to several bam files and streaming the output into freebayes.

[evelyn]

Can samtools merge be used to make one bam file and use that for freebayes?

[swbarnes2]

I'm pretty sure that what I used once, but it's been a while. Picard can add readgroups too.

[evelyn]

Thank you, I will use picard to add readgroups and then merge all the bam files to use for freebayes.

[Ram]

What does the samtools manual say about merging BAMs?