Question: Generate consensus sequence from BAM without ambiguity codes
0
gravatar for Tonor
9 months ago by
Tonor420
UK
Tonor420 wrote:

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

samtools mpileup -uf REF BAM | bcftools call -c | vcfutils.pl 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
1

Past post that may be useful: How to create consensus from bam file without IUPAC code?
https://github.com/cbg-ethz/ConsensusFixer

ADD REPLYlink written 9 months ago by genomax59k

Consensus Fixer looks great - thanks very much

ADD REPLYlink written 9 months ago by Tonor420
1
gravatar for swbarnes2
9 months ago by
swbarnes24.5k
United States
swbarnes24.5k wrote:

You can edit vcfutils.pl 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);

to

 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.

Help
Access

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