Fasta file not matching VCF
1
0
Entering edit mode
4.8 years ago

Hi,

I've used bcftools consensus to create fasta files from individual vcf files. However, some sites in the fasta file do not correspond to the vcf. Using the --sample flag resolved some of these problems, but not all. For example this site

chr22   2600340 .   C   T   1638.22 PASS    .   GT:AD:DP:GQ:PL  1/1:0,9:9:27:405,27,0

should have T, but

samtools faidx 8L19766_sam_sample.fasta chr22:2600340-2600340

yield a C.

The vcf files have been called using GATKs HaplotypeCaller. I have also attempted using GATKs FastaAlternateReferenceMaker to create fasta files, but this made the discrepancies far worse.

What could be causing these issues?

Thanks so much for any help!

fasta BCFtools SAMtools • 1.3k views
ADD COMMENT
0
Entering edit mode
4.8 years ago
JC 13k

Your VCF it indicates reference is "C", and that's what you see in your Fasta.

To visualize the alternative call, you need to view your BAM file: samtools tview FILE.bam chr22:2600340-2600340

ADD COMMENT
0
Entering edit mode

No, if the VCF indicates that the individual has the alternative allele at a certain site, as is the case here, bcftools consensus should place the alternative allele at that position in the fasta file. This is how it works for most of my sites, but not all.

ADD REPLY
1
Entering edit mode

Do you have insertions or deletions in your vcf? If so, then chr22:260034 in the consensus fasta will not represent the same position in the original fasta.

ADD REPLY
0
Entering edit mode

Can you please give us the exact bcftools consensus command you're using?

ADD REPLY
0
Entering edit mode

so, you create a consensus "patching" your reference genome? how do you generate that?

ADD REPLY
0
Entering edit mode

No, there are no insertions or deletions .

This is the command I use:

cat house_sparrow_genome_assembly-18-11-14_masked.fa \
| bcftools consensus --sample 8L64869 8L64869_SNPs.vcf.gz \
> 8L64869.fasta
ADD REPLY

Login before adding your answer.

Traffic: 820 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