I have a number of tumor/normal pairs from a cancer sequencing project (exomes sequenced with Illumina HiSeq). They are not necessarily sequenced to the same depth. I am interested in assessing the number of somatic mutations in each sample, and comparing these counts between samples.
One way to do this would be to run a somatic variant caller on the samples, and count the number of variants identified in each. However, at any fixed parameterization of the caller, samples sequenced to less depth will erroneously tend to give fewer mutations, since the variant caller will have less confidence in any particular call.
What I'd like instead is something that gives me an unbiased estimate of the number of mutations in a sample, as well as an error bar on that number. That way, as the depth increases, I'd expect to get a tighter confidence interval, but not an ever-increasing number of mutations.
Does anyone know of a method that does this or have a suggestion on the correct theoretical way to think about it?
