mpileup consensus behaviour

[rwness]

I am wondering if anyone has any experience calling a consensus genotype from mpileup data. I am curious how the suggested command from samtools documentation treats indels and what it does with 'heterozygous' indels.

samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

Thanks

[rwness]

If anyone is curious - the answer to my question is that this consensus will include NO indels. Any lines in the VCF/BCF that contain indels -- heterozygous or homozygous -- will be ignored.

I wrote a parser for the vcf created by this command:

$samtools mpileup -uf ref.fa aln.bam | bcftools view -cg -

that will generate a consensus sequence where heterozygous indels are ignored and heterozygous sites are represented with ambiguity codes. It will also output 2 additional consensus sequences in which heterozygous SNPs and indels are included with gaps to keep the sequences aligned.

[abeeaton]

I've been sticking with the old pileup -c, so I can get the heterozygous sites with ambiguity codes and genotype quality scores. How were you able to get this out of the mpileup/vcf output?

[swbarnes2]

The new way posted above, using the all-points vcf and vcfutils.pl will get make ambiguous letters in the consensus. Here's the relevant part of the code:

my %het = (AC=>'M', AG=>'R', AT=>'W', CA=>'M', CG=>'S', CT=>'Y',
GA=>'R', GC=>'S', GT=>'K', TA=>'W', TC=>'Y', TG=>'K');


$q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
if ($q < 0) {
$_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
$b = ($_ < .5 || $alt eq '.')? $ref : $alt;
$q = -$q;
} else {
$b = $het{"$ref$alt"};
$b ||= 'N';
}

In English, the first line means "check the FQ"; the FQ is negative when the SNP is homozygous, and positive when it's mixed, and the bigger the absolute value, the more confident the SNP.

So if it's < 0, it does the first part of code (it checks against the AF1, and if the AF1 is > 0.5, which it should be for a homozygous SNP, it sets $b as the alternate letter, if for some reason the AF1 is < .5, it sets $b as the old reference letter.)

If the FQ is positive, then the SNP should be mixed, and it concatenates the two letters, and checks the hash above to know what single letter to set $b to. Then it adds $b to the growing sequence.

$q, which is derived from the FQ, ends up being the quality score, though it gets tweaked a little; it adds 33.449 to the figure, then converts it to a letter, capping it at a quality of 126.

$q = int($q + 33 + .499);
$q = chr($q <= 126? $q : 126);

Gaps are handled as they were in the old program, where they are NOT added in, there is just a window of lowercase letters around them. Personally, I made a little perl script, and I feed it the genome, and a conservatively filtered list of SNPs, and I put the changes in that way.

[lh3]

Thanks for the nice explanation, swbarnes2. Just want to add that the key purpose of generating consensus is to distinguish regions that cannot be called and regions that are called confidently as homozygous reference. Heterozygotes found in the consensus sequence actually have lower accuracy because fewer filters were applied.

[nupurgupta]

I am trying to understand how to set $q. In the line

$q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);

Does it mean get $q from command line (1st paramter) in case there is an FQ value?

Thanks very much