This is a test version of Biostars. For the public version, visit https://www.biostars.org.
getting consensus sequence with bcftools failed

I want to extract the mapped regions out of bam files. So I tried bcftools. I have the current version of bcftools (V.19). I tried the following steps

bcftools mpileup -Ou -f reference.fa alignments.bam | bcftools call -mv -Oz -o calls.vcf.gz

after that a note appeared

Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid

but this file has the regions which I want, so I go on with:

bcftools index calls.vcf.gz

cat reference.fa | bcftools consensus calls.vcf.gz > consensus.fa

Here also the note appear, and the whole fasta regions are displayed in consensus.fa

Note: the --sample option not given, applying all records

But when doing this:

cat reference.fa | bcftools consensus -s calls.vcf.gz > consensus.fa

nothing works. I tried a lot of stuff and have no idea where the mistake is.

Hope you can help me.

genome bcftools samtools rna-seq vcf

The -s option of bcftools consensus is meant for specifying the sample name (i.e. if the input is a multi-sample VCF file, this tells bcftools which sample to pick). Can you see what the sample name is in the VCF file?

Okay thanks, than -s is not necessary. The sample has the following header:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  mapped_male_replicate.bam 
W003_0-161190     154995  .       G       A7.30814 .       DP=2;VDB=0.18;SGB=-0.453602;MQ0F=0;AC=2;AN=2;DP4=0,0,2,0;MQ=20  GT:PL  1/1:36,6,0 
W003_0-161190     155048  .       T       G       19.4636 . DP=3;VDB=0.698947;SGB=-0.511536;MQ0F=0;AC=2;AN=2;DP4=0,0,3,0;MQ=20     GT:PL   1/1:49,9,0

Is there something wrong with the vcf file?

0 answers

No answers yet.

Log in to answer this question.