Thanks Michael,
This is an interesting idea. Usually, we resolve PCR/optical duplicates with the Picard MarkDuplicates command, which just chooses the read with the highest average quality. This approximates what you want, I think.
Is the Picard MarkDuplicates command run by default? Or is this a separate command from a different package?
Also, I find it hard to believe that sequence errors would be causing you trouble at 10%, as long as you're filtering for quality and have decent coverage. PCR might in limited cases, and Picard effectively takes care of that.
10% is fine, but we have problems at 1% and lower. I think even the defaults in tallyVariants() use 1% as a cutoff if I'm not mistaken.

What you want is to specify the cycleBreaks argument to VariantTallyParam. It allows you to define bins by cycle. So if you made bins for all 100bp, you could do things like generate a consensus by read family. A family corresponding to position X would be cycle bin #1 at X, cycle bin #2 at X+1, etc. Then just pick the alt (or ref) with the highest count. Might be tricky to implement efficiently for the whole genome. Of course, this assumes you've got single end reads, otherwise this won't work, because the mate is not considered.
Cool, I will try that out, that may be just what I was looking for. It might be tricky for the whole genome, but we are restricting ourselves to just the mitochondrial genome. Still gives us several thousand start sites to work from but should be easier than the whole genome. As for single end reads, can I just filter on the strand (i.e. '+' or '-')? Or will I have to perform separate alignments for each of the paired reads?


	[[alternative HTML version deleted]]

