This is a test version of Biostars. For the public version, visit https://www.biostars.org.
how to get one vcf file for freebayes from many bam files

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!

snp
• 5,626 views • link
updated by Ram Ram • written by evelyn evelyn

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.

• 0 views • link
written by Ram Ram

It says:

freebayes: unrecognized option '--no-complex'
did you mean --samples ?
• 0 views • link
written by evelyn evelyn

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

• 0 views • link
written by Ram Ram

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

• 0 views • link
written by evelyn evelyn

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

• 0 views • link
updated by Ram Ram • written by GenoMax GenoMax

--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.

• 0 views • link
written by Ram Ram

1 answer

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

• 0 views • link
written by swbarnes2 swbarnes2

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.

• 0 views • link
written by h.mon h.mon

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

• 0 views • link
written by evelyn evelyn

What does the samtools manual say about merging BAMs?

• 0 views • link
written by Ram Ram

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

• 0 views • link
written by swbarnes2 swbarnes2

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

• 0 views • link
written by evelyn evelyn

Log in to answer this question.

More posts like this
  • How to solve Freebayes memory issue
    written by anikcropscience • anikcropscience
    Hello, I am running `Freebayes` variant calling pipeline. I have 1440 bam files from genotype by sequencing data. The genome size is around 540 Mb. …
  • FreeBayes detection capacity
    written by quentinperriere • quentinperriere
    Hello, I'm currently utilizing FreeBayes for variant detection with reads generated through Oxford Nanopore technology. However, I've encountered an issue where not all known variants …
  • multiple reference sequences
    written by emilydolivo97 • emilydolivo97
    Hello , Can I use freebayes on multiple reference sequences . These reference sequences are merged in one FASTA file and each reference start with …
  • freebayes not reading the chunks of regions file
    written by evafinegan • evafinegan
    Hello, I am using freebayes for SNP calling using: RUN=${SLURM_ARRAY_TASK_ID} INPUT_DIR=/path/ OUTPUT_DIR=/path/ INPUT=$(ls $INPUT_DIR/*.txt | sed -n ${RUN}p) SAMPLE=$(basename "$INPUT" .txt) freebayes -f genome.fasta \ …
  • freebayes error variant calling
    written by evelyn evelyn
    Hello, I am using `freebayes` for variant calling using: list=list.txt freebayes -f transcripts.fasta -L $list --ploidy 2 \ --min-alternate-fraction 0.2 \ --min-alternate-count 2 \ --min-coverage …
  • freebayes problem in sample name in vcf
    written by evelyn evelyn
    Hello, I am using following command to create a `vcf` file with `freebayes`: freebayes -f transcripts.fasta duplicates.bam --max-complex-gap 1 \ --theta 0.01 \ --ploidy 2 …
  • freebayes SNP calling question
    written by evelyn evelyn
    Hello, I am using freebayes for making vcf files using: 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 this error: unable …
  • Error in freebayes SNP calling
    written by evelyn evelyn
    Hello, I am using `freebayes` to make a `vcf` from a list of `.bam` files (duplicates removed): RUN=${SLURM_ARRAY_TASK_ID} LIST_OF_BAMS=/path/list_of_bam.txt freebayes \ -f /path/genome.fasta \ -L …
  • SNP calling specific genes with GATK
    written by gubrins gubrins
    Good afternoon, I'm trying to do the SNP calling of some files from a target-capture sequencing experiment with GATK. I don't have any problem but …
  • SNPs mining using freebayes: problem in using multiple bam file of different ploidy
    written by ram.mainali • ram.mainali
    My objective is to identify the SNPs on resistance gene/haplotype. I have resistance specific haplotype read and mapped bam file (ploidy = 1), susceptible haplotype …