Hi,
I have a general question concerning genotype calling and thought that maybe some of you have encountered the same problem as me and could give opinion on this. I have a number of sequences, each with reads mapped to it. References come from an assembly of transcriptome from single individual. I would like to call heterozygous sites in each sequence, I don't want to get only a set of most likely SNPs but give the most probable genotype to each site so I could further estimate for example diversity. Unfortunately some parts of sequences have poor coverage. I used samtools and than I ran into some problems (although I think the software doesn’t really matter here, it could be a problem with any software which calls genotypes using probabilistic method). Since the qualities of heterozygotes are much lower than of homozygotes in the same low- covered regions I am hesitant to apply any quality threshold because I will underestimate the diversity. If I apply coverage threshold (take e.g. only sites covered with more than x reads) I may exclude low-expressed sequences in favor of highly expressed sequences and that could also bias the results (and also I will lose a lot of data). On the other hand if I choose the genotype with highest posterior probability I may have quite a lot of false positives especially in low-covered regions. I am not sure if there is the only right solution to this, but maybe I’m missing something that is already taken into account in the model. If anyone has already experienced similar dilemma and has maybe an idea how to properly deal with it I would like to hear it.
I have a general question concerning genotype calling and thought that maybe some of you have encountered the same problem as me and could give opinion on this. I have a number of sequences, each with reads mapped to it. References come from an assembly of transcriptome from single individual. I would like to call heterozygous sites in each sequence, I don't want to get only a set of most likely SNPs but give the most probable genotype to each site so I could further estimate for example diversity. Unfortunately some parts of sequences have poor coverage. I used samtools and than I ran into some problems (although I think the software doesn’t really matter here, it could be a problem with any software which calls genotypes using probabilistic method). Since the qualities of heterozygotes are much lower than of homozygotes in the same low- covered regions I am hesitant to apply any quality threshold because I will underestimate the diversity. If I apply coverage threshold (take e.g. only sites covered with more than x reads) I may exclude low-expressed sequences in favor of highly expressed sequences and that could also bias the results (and also I will lose a lot of data). On the other hand if I choose the genotype with highest posterior probability I may have quite a lot of false positives especially in low-covered regions. I am not sure if there is the only right solution to this, but maybe I’m missing something that is already taken into account in the model. If anyone has already experienced similar dilemma and has maybe an idea how to properly deal with it I would like to hear it.