Question: Generate consensus sequence from BAM without ambiguity codes
gravatar for Tonor
16 days ago by
Tonor410 wrote:

To generate a consensus sequence from a BAM file I use the command:

samtools mpileup -uf REF BAM | bcftools call -c | vcf2fq > SEQ

samtools/vcfutils gives ambiguity codes when there is a mixed base population at a site, even if the secondary base is quite a low frequency (say 20%) - is there anyway to avoid ambiguity codes? or to increase the frequency at which they would be used? i.e. if I have 60% A and 40% G at a position - i want my consensus to be A.

Also, anyway to get samtools to give ma fasta file instead of fastq - I have a script that converts it - but was wondering if vcfutils can give fasta automatically.

Happy to use alternative to samtools to get a consenus sequence.

sequence • 105 views
ADD COMMENTlink modified 16 days ago by swbarnes23.3k • written 16 days ago by Tonor410

Past post that may be useful: How to create consensus from bam file without IUPAC code?

ADD REPLYlink written 16 days ago by genomax44k

Consensus Fixer looks great - thanks very much

ADD REPLYlink written 16 days ago by Tonor410
gravatar for swbarnes2
16 days ago by
United States
swbarnes23.3k wrote:

You can edit to have it output fastas instead of fastqs. You want to change these lines

 print "\@$chr\n"; &v2q_print_str($seq);
 print "+\n"; &v2q_print_str($qual);


 print "\>$chr\n"; &v2q_print_str($seq);
 #print "+\n"; &v2q_print_str($qual);

Probably best to make a copy of the whole vcf2fq part, which starts with

sub vcf2fq {
  my %opts = (d=>3, D=>100000, Q=>10, l=>5);

and ends with

  &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l});

And change 'sub vcf2fq' to 'sub vcf2fa' Then you'll have both functionalities.

Unfortunately, putting in a filter for when to not output the het letter is trickier.

ADD COMMENTlink written 16 days ago by swbarnes23.3k

Brilliant thanks for that

ADD REPLYlink written 15 days ago by Tonor410
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 799 users visited in the last hour