bcftools mpileup no output (empty calls.bcf file)

Question

[RobertUt]

Hi everyone,

I am trying to make a consensus from a aligned bam file. This is my workflow:

bcftools mpileup -Ou -f reference_HG19.fa aligned.bam | bcftools call -mv -Ob -o calls.bcf

bcftools index calls.bcf

cat reference_HG19.fa | bcftools consensus calls.bcf > consensus.fa

However, the created consensus.fa file turned out to be identical with the reference_HG19.fa file. I looked at calls.bcf file with bcftools view, and it seems the following columns are all empty:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT aligned.bam

When I look at calls.bcf with samtools view I get the error "Aborted".

Does anyone have an idea how to solve this issues? Thank you very much!

Robert

[finswimmer]

Hello and welcome to biostars RobertUt ,

Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time.

Thank you.

Are you sure that you will have any variants in your aligned.bam?

When I look at calls.bcf with samtools view I get the error "Aborted".

AFAIK samtools view cannot read bcf files. bcftools view would be what you are looking for.

fin swimmer

[RobertUt]

Dear fin,

Thank's for your suggestions.

I just checked the bam file with samtools tview and found plenty of variants.

By the way, the only warning/error I get is about ploidy and a chromosome, which is missing in my bam file but present in the reference fasta.

Any other suggestions?

[GenoMax]

Please use ADD COMMENT/ADD REPLY when responding to existing posts to keep threads logically organized.

[finswimmer]

How did you produce your bam file? Please show us the exact command used.

[RobertUt]

The bam file wasn't generated in our on lab but was downloaded from:

https://www.ebi.ac.uk/ena/data/view/PRJEB3371

It's a whole genome sequencing file, generated with Illumina HiSeq 2000. I am also suspecting the bam file to be the problem, because the fasta worked fine in previous commands. Chromosome notation is consistent between bam and fasta.

[Kevin Blighe]

To which genome was the BAM aligned (check header)? Which genome ref are you using?

Just to humour me, could you try:

bcftools mpileup -Oz -f reference_HG19.fa aligned.bam | bcftools call -mv -Ob -o calls.bcf

...and:

bcftools mpileup -f reference_HG19.fa aligned.bam | head

[finswimmer]

Hello,

you have to make sure that the bam files are coordinate sorted before doing variant calling. At least the first sample isn't.

fin swimmer

[RobertUt]

Hi fin,

the bam file I'm using is coordinate sorted.

Robert

Answer 1

[RobertUt]

We finally found the solution: The downloaded bam was generated by paired end sequencing. For some unknown reason bcftools considers these read pairs to be "anomalous". By using the -A input option one prevents bcftools from skipping anomalous read pairs. Now it's working as it is supposed to.

[untitpoi]

Thanks for that answer, option -A solved my issue (void output) too!

[pwchan]

Thanks! I have a similar workflow (below), and using the -A input flag also solved my void output issue on my .vcf.gz file. However, my bcftools consensus command is still outputting a .fasta file that is identical to my input reference sequence, it appears no variants are being applied.

bcftools mpileup --max-depth 500 -f -A $referenceSeq $sample.plastome.rmdup.bam | bcftools call -mv -Oz -o $sample.calls.vcf.gz bcftools index $sample.calls.vcf.gz bcftools consensus -f $referenceSeq $sample.calls.vcf.gz -o $sample.plastome.fasta