bcftools mpileup no output (empty calls.bcf file)
1
2
Entering edit mode
3.5 years ago
RobertUt ▴ 70

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

bcftools samtools mpileup consensus calls.bcf • 3.1k views
ADD COMMENT
0
Entering edit mode

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

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

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Hi fin,

the bam file I'm using is coordinate sorted.

Robert

ADD REPLY
4
Entering edit mode
3.4 years ago
RobertUt ▴ 70

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.

ADD COMMENT
2
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 2315 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6