bcftools consensus doesn't allow filtering of phased vs. unphased heterozygous sites

[elmoremh]

I want unphased heterozygous sites to be entered into my consensus sequences as either Ns or maybe the appropriate IUPAC ambiguity code. It's not clear that there's a way to communicate this to bcftools consensus.

It's possible there's a way to do this with expressions, but if so I haven't figured it out. (Can I parse text within a column? None of the examples in the documentation lead me to think that.)

Here's test data (please use individuals 10171 or 10004):
phase_filtering_test.tar.gz

[pd3]

This should work, have you tried to run as follows?

bcftools consensus -s <sample> -I file.vcf.gz -f ref.fa

[elmoremh]

I was using bcftools 1.5, so I hadn't seen the -I option, so thank you! But I want to separate the two haplotypes in my partially-phased diploid data as much as possible. I want to be able to separate the two haplotypes where variants are phased and put in Ns or IUPAC where variants are unphased.

Below, it looks like -H 1 overrules the -I option. I want some way to tell the program to distinguish between phased and unphased variants, only giving an N or IUPAC code for unphased heterozygous sites.

Thank you so much for you help with this matter already.

[elmoremh@holylogin03 ~]$ /n/home15/elmoremh/sw/bcftools-1.9/bcftools consensus -H 1 -I -s 10171 ~/AmanitaBASE_all_individuals_phased_HD1.recode.vcf.gz -f ~/10511_HD1.fasta  > 10171_HD1_H1_bcftools-consensus.fasta
[elmoremh@holylogin03 ~]$ more 10171_HD1_H1_bcftools-consensus.fasta
>Contig87:398638-399360
TGCATTCCACTATCCGCCAACGATTATTGCGGGTCGAAGAGGATCTACTCGATGC**TAA
TCGGTGAAGAATCGGCAATGAATGTTCTCATTGATTCCTGGGCATTATTGTGCACTGACC
TTCAAACGGCACTCGACGCTCAGCTATCGGATGAGGACACTAAACACCTTGCGCATGTCA
CTGCTTCTCGCGTGCAGATAATTGCGGAAAGCCTTATGGAGGTCCGCTCACAGGCTGATG
CAATGGCTTATTCACTGGAGGTCCTAGAA*CAATGATGTCGTCGTTGACCATATCAAAGG
ATAACGAGGGTACAACCACCCGCGGCTCCCTATCACCTACTCCAGGTCC*TCTCCACGTT
GCGCTATCACGTCGTACATCAGGACCGCTTATGAGTGGTTGGTACGTAACCTCCATAATC
CTTACCCACCGAAGCATGTCAAGGAGGCGATCGCGCGCGACACGGGGTCCGCTATCAAAG
ACCTGGATAACTGGTTTATCGATGCACGAAGAAGGATTGGCTGGAATCAGCTACGTCGCA
AGCGTTTTTCAAATAAAAGGCAAGAAATCGTGCGAGCTGCTACAAAGTTTTTCGCTGGTG
CAGGATGTTCGGAATTATTGGATCCTGATACCGAGATGGAATTTGTGGCGATCGAAGCCT
CTGCCTTTCAACTCTATCCGAATCAGATTACATCGCGACATGAGCTGGCCAAACAAACAA
AC
[elmoremh@holylogin03 ~]$ /n/home15/elmoremh/sw/bcftools-1.9/bcftools consensus -I -s 10171 ~/AmanitaBASE_all_individuals_phased_HD1.recode.vcf.gz -f ~/10511_HD1.fasta  > 10171_HD1_H1_bcftools-consensus.fasta
[elmoremh@holylogin03 ~]$ more 10171_HD1_H1_bcftools-consensus.fasta
>Contig87:398638-399360
TGCATTCCACTATCCGCCAACGATTATTGCGGGTCGAAGAGGATCTACTCGATGC**TAA
TCGGTGAAGAATCGGCAATGAATGTTCTCATTGATTCCTGGGCATTATTGTGCACTGACC
TTCAAACGGCACTCGACGCTCAGCTATCGGATGAGGACACTAAACACCTTGCGCATGTCA
CTGCTTCTCGCGTGCAGATAATTGCGGAAAGCCTTATGGAGGTCCGCTCACAGGCTGATG
CAATGGCTTATTCACTGGAGGTCCTAGAA*CAATGATGTCGTCGTTGACCATATCAAAGG
ATAACGARGNTACAACYACCCRYGRCTCCCTATCWCCTACTCCAGGTCC*TCTCCRCGTT
GCGCTATCACGTCGTACATCAGGACCGCTTATGAGTGGYTGKTRCGTAACCTCCAYAATC
CTTACCCACCGAARCATSTCAAGGAGKYGATCGCRCGSGACACGGGKTCCKCTAYCAAAG
ACCTGGATAACTGGTTTATCGATGCACGWAGAAGGATYGGCTGGAATCAGCTACGTCGCA
AGCGTTTTTCAAATAAAMGGCAAGAAATCGTGCGAGCTGCTACMAAGTTTTTCGCTGGTG
CAGGATGYYCKGAWTYATTGGATCCTGATACCGAGATGGAATTTGYRGCGATCGAAGCCT
CYRCCTTTCAACTCTATCCGAAYMRGATTACMTCRCGACWYGAGNNNGCCMAACARACAA
RC

[pd3]

This should be now possible with the --haplotype 1pIu or --haplotype 2pIu option:

-H, --haplotype 1pIu,2pIu: first/second allele for phased and IUPAC code for unphased GTs

Let me know in case of problems.

[elmoremh]

It works! Thank you so much!!! You've saved me a lot of time and frustration.