Question: Generate consensus sequence from BAM without ambiguity codes
gravatar for Tonor
9 months ago by
Tonor420 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 • 703 views
ADD COMMENTlink modified 9 months ago by swbarnes24.5k • written 9 months ago by Tonor420

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

ADD REPLYlink written 9 months ago by genomax59k

Consensus Fixer looks great - thanks very much

ADD REPLYlink written 9 months ago by Tonor420
gravatar for swbarnes2
9 months ago by
United States
swbarnes24.5k 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 9 months ago by swbarnes24.5k

Brilliant thanks for that

ADD REPLYlink written 9 months ago by Tonor420
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: 2044 users visited in the last hour