This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Generating vcf files per individual instead of one big file from bam files

Hi,

I have multiple bam files that I used to generate vcf file for variant calling. I used the command below

bcftools mpileup -O v -f $REF  bam/*.bam > genotypes.vcf

Instead of getting one file each sample, I am getting one big file that contains information on all samples.

Is there any way to either split them into samples or re-run the command to get one file per sample?

Regards

vcf bam variant
• 3,003 views • link
updated by Ram Ram • written by salman_96 • salman_96 
bcftools mpileup -O v -f $REF  bam/sample1.bam > sample1.vcf
bcftools mpileup -O v -f $REF  bam/sample2.bam > sample2.vcf

use a loop or better, a workflow manager.

'bcftools call' should also be invoked after mpileup.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Yes, I did write only the first command that generates genotype vcf files. Unfortunately, there are multiple files and I cannot find any way to loop through

• 0 views • link
written by salman_96 • salman_96

and I cannot find any way to loop through

what have you tried ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I tried a one liner to avoid large genotype file. I get no results'

for bam in $(ls *.bam) ; do bcftools mpileup -O v -f GRCh38.primary_assembly.genome.fa $bam | bcftools call --ploidy 1 -vm -O v > VARIANTS.vcf

• 0 views • link
written by salman_96 • salman_96

and there is no error message ? the vcf are all empty ? no header ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Hi,

Yes no error message. Just the command returns ">"

• 0 views • link
written by salman_96 • salman_96

Yes, The following command did work when used done but it is generating same named file for every sample. So at the end the process resulted in only one file which was the last one in the samples because of loop.

• 0 views • link
written by salman_96 • salman_96

The bam files are named somehow like this

CL03_T1STARpairedAligned.sortedByCoord.out.bam CL05_T1STARpairedAligned.sortedByCoord.out.bam CL05_T1_2DSTARpairedAligned.sortedByCoord.out.bam CL09T1_2DSTARpairedAligned.sortedByCoord.out.bam

• 0 views • link
written by salman_96 • salman_96

This is my script

/bin/bash

cd /The_Variant_call_project/bam

for filename in *.bam

    do
    echo "working with file $filename"

base=$(basename $filename .bam) echo "base name is $base"

for bam in $(ls *.bam) ; do bcftools mpileup -O v -f GRCh38.primary_assembly.genome.fa $base | bcftools call --ploidy 1 -vm -O v > $final_variants

done

The error I run into is Failed to open file "CL48T3STARpairedAligned.sortedByCoord.out" : No such file or directory No I know that the file name should end with .bam but here I get CL48T3STARpairedAligned.sortedByCoord.out instead of CL48T3STARpairedAligned.sortedByCoord.out.bam.

I would appreciate if you would add any suggestion. I do believe there is some problem with the base name.

Regards

• 0 views • link
written by salman_96 • salman_96

1 answer

Hi I was able to solve the problem. Had to make some changes in the last command

bcftools mpileup -O v -f GRCh38.primary_assembly.genome.fa ${infile} | bcftools call --ploidy GRCh38 -vm -O v > ${final_variants}

Thank you for the help. I appreciate it

• 0 views • link
written by salman_96 • salman_96

Log in to answer this question.

More posts like this