BCFtools consensus sequence
1
1
Entering edit mode
3.4 years ago
Gene_MMP8 ▴ 210

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.

next-gen • 2.6k views
ADD COMMENT
3
Entering edit mode
3.3 years ago

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.

ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

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?

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

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