BCFtools consensus sequence

Question

[Gene_MMP8]

I am trying to incorporate variants in my ref.fa file using bcftools consensus. But the output consensus.fa has only one chromosome (the first one from the ref.fa file). Why are the other chromosomes not showing up in consensus .fa?
The command used was :

bcftools consensus -f hg38.fa compressed.vcf.gz > consensus.fa

where compressed.vcf.gz is my variants file.

Answer 1

[Kevin Blighe]

Have you indexed the hg38.fa file with samtools faidx?; is your VCF indexed with tabix -p vcf?; double-check that the chromosomes follow the same ID system in both files, e.g., chr1, chr2 or 1, 2, etc. In the VCF header, also double check that the chromosome contigs are listed in the header itself.

[davidenoma]

Hi Kevin, I tried this command:

samtools faidx GRCh38_latest_genomic.fna | bcftools consensus hg00096.vcf.gz > out.fa

And it gave the error that it applied 0 variants.

[Kevin Blighe]

Hi, I am not sure that the output of samtools faidx can be piped into bcftools consensus. What are you aiming to do?

What you probably want is:

samtools faidx GRCh38_latest_genomic.fna ;
tabix -p vcf hg00096.vcf.gz ;
cat GRCh38_latest_genomic.fna | bcftools consensus hg00096.vcf.gz > out.fa ;

[davidenoma]

Yes, it gave applied 0 variants in results again. However, I received a warning that
Sequence "NT_187681.1" is not in the vcf file for instance. As well as other sequences.

Your last line should have bcftools consensus though.

What happens when the contig cannot be found in the vcf file?

[Kevin Blighe]

Thanks, I have added 'consensus' Are the main contigs present?, i.e., chr1, chr2, et cetera? Basically, the contigs in the fna file and your VCF will have to have the same names, including any prefixes. You can check by simply viewing, e.g., the first few records of the VCF, and just the first line of the fna file.