SNP Frequency

[kvarala]

I am trying to estimate the frequency of SNP occurence in my library as reported by the maq aligner. The quick and dirty way of doing this, I believe, is to calculate the number of bases of the reference genome that are covered by at least 3 uniquely mapped reads. That should match the cutoff specified in SNPFilter command of maq. Does anyone know if this information is available somewhere in the maq output files?
Also, please let me know if you think my logic for estimating the frequency is blatantly wrong. I understand that its not perfect but this is just my first hack at it.

[apfejes]

Hi Kvarla,

I've written several SNP callers myself. In principle, your method is good. I don't trust 3 independent fragments, however - I think it's not nearly enough. Generally, I'll use a minimum of 8x, and the base qualities of each aligned base must be over a given threshold (20 for scores out of 30, 27 for scores out of 40). Some people also require that you see fragments in both directions. YMMV.

I don't think MAQ is much different than that.

The only thing I'm not sure of is why you're doing an estimate - it's not hard to get the actual value of SNPs in your library.

Anthony

[kvarala]

I can easily find the number of SNPs. I'm trying to get the number of bases that are aligned at the same quality as the SNP calls but are identical to the reference. I don't think its accurate to estimate the SNP frequency as no. of SNPs called / total bases aligned since all the bases are not available for making that call.

[apfejes]

That makes sense. I typically just write out the number of bases observed with that level of coverage, since I use my own SNP callers. I'm not sure how one would get Maq to do that on it's own. I can point you to one of my SNP callers, since it's open source and produces similar results to maq, but I don't think that's what you're asking.

As far as I'm concerned, there really is no way to estimate the number of bases Maq will see without doing the actual calculation.

Anthony

[bioinfosm]

will maq's pileup be of help?
use the verbose command to output qualities as well, and then parse out the numbers

[kvarala]

ah.. of course. Thanks bioinfosm for the pointer. Serves me right for not reading the manual more carefully. Next time before I ask I'll make sure I RTFM.

[dvh]

Note there are some default filters on maq pileup output:
-m INT Maximum number of mismatches allowed for a read to be used [7]
-Q INT Maximum allowed number of quality values of mismatches [60]

We find these subtly change results - dont expect maq pileup to generate exactly what is in the .map file.

david

[sam253]

SNP frequency

Hi,

Could you please tell me how you use pileup file for SNP frequency. I am also trying to find out variant proportion in each position. Do you know how can I calculate that? e.g if If reference base is A and consensus base is T/G how can I know at what proportion each one of these bases are?

Thanks