Thanks for that answer, option -A solved my issue (void output) too!
• 0 views
•
link
written
by
untitpoi •
bcftools mpileup no output (empty calls.bcf file)
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
• 7,330 views
•
link
updated
by
pwchan •
•
written
by
RobertUt •
•
written
by
RobertUt •
1 answer
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.
• 1 views
•
link
written
by
RobertUt •
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
Log in to answer this question.
More posts like this
-
Extract SNPs present in mutants
written by pablo
-
How to generate a consensus sequence from BAM file with bcftools?
written by marongiu.luigi
-
Reality check for consensus fasta generated with bcftools
written by marongiu.luigi
-
Consensus sequence calling with normalisation of indels
written by michael.flower.14 •
-
bcftools mpileup before bcftools call
written by Egelbets •
-
Problems in Variant Calling Code
written by arjun.murthy779 •
-
getting consensus sequence with bcftools failed
written by Heinrich_M •
-
bcftools mpileup --max depth : What would be a reasonable setting for this option?
written by VBer
-
The fasta sequence does not match the REF allele at chr1:1298836 ?
written by RobertUt •
-
bcftools consensus does not produce an output
written by ieie •
Hello and welcome to biostars RobertUt ,
Please use the formatting bar (especially the
codeoption) 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?AFAIK
samtools viewcannot readbcffiles.bcftools viewwould be what you are looking for.fin swimmer
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?
Please use
ADD COMMENT/ADD REPLYwhen responding to existing posts to keep threads logically organized.How did you produce your bam file? Please show us the exact command used.
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.
To which genome was the BAM aligned (check header)? Which genome ref are you using?
Just to humour me, could you try:
...and:
Hello,
you have to make sure that the
bamfiles are coordinate sorted before doing variant calling. At least the first sample isn't.fin swimmer
Hi fin,
the bam file I'm using is coordinate sorted.
Robert