Generating vcf files per individual instead of one big file from bam files

Question

[salman_96]

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

[Pierre Lindenbaum]

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.

[salman_96]

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

[Pierre Lindenbaum]

and I cannot find any way to loop through

what have you tried ?

[salman_96]

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

[Pierre Lindenbaum]

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

[salman_96]

Hi,

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

[salman_96]

No vcf files generated

[Pierre Lindenbaum]

because a for loop should end with done

https://www.cyberciti.biz/faq/bash-for-loop/

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

[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

[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

Answer 1

[salman_96]

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