I am doing amplicon sequencing of a virus across many different regions. Lets say I have 20k unique variants of the same virus that I put into my pcr assay and after sequencing and amplifications I am left with, 19k variants that appear. But many of the variants appear at really low read count and I don't know if that read is real or just noise. I have a threshold T that I say if a variant appears above T I'm counting it as real. Out of the 20k variants about only 15% of variants appear to be "real". I have some data on different input species counts (2k, 10k, 20k, 40k), and their corresponding fraction of "real".

My question is how do I determine which input variants count is best? I only looked at 4 different values and maybe the correct value isn't sampled, I obviously can't try everything between 2k and 40k. I think I want a balance between largest fraction of "real" variants and total number of variants.

Is there any research (or better yet code) that answers/discusses this problem?

Thanks in advanced.

  • $\begingroup$ @M__ Thanks for the comment adjusted my question. $\endgroup$ Jul 15, 2023 at 23:22
  • 1
    $\begingroup$ Added doing amplicon sequencing across different regions $\endgroup$ Jul 15, 2023 at 23:59
  • $\begingroup$ Thanks thats slightly clearer. There is an issue here on the nomenclature. 15+K viral species is 6 times more viruses than designated by the ICTV (International Committee on the Taxonomy of Viruses), i.e. more viruses than the total number of species ever declared. Could you kindly define the term species? What you have stated is not taxonomically possible. $\endgroup$
    – M__
    Jul 16, 2023 at 1:04
  • $\begingroup$ I see, using species not exactly clearly. They are all variants of the same virus. Will adjust $\endgroup$ Jul 16, 2023 at 1:26

1 Answer 1


The underlying answer to the question is have you experimentally controlled for in vitro artefacts? Thats taking a cloned virus and performing the same RT-PCR. If not thats makes things very difficult. If there isn't a cloned virus thats a problem. You might need to ask your wet-lab colleagues what old school "cloning a virus" means, or Biology SE will explain this. The reason this is important is given below.

The first step to the question is simply to plot a histogram against copy number. What you will find is a heavily skewed distribution where the copy number of one variant will represent >40% (at a wild guess) of all the viral copy number. Moreover, a handful (less than 100) of variants will comprise 90% of all the copy number. Once that is established it is much easier to assess what is a true variant. What I'm saying is the difference is one log minimum between what you think you are looking at and what is there.

The sources of artefact:

  • RT-PCR artefacts, particularly if there is an RT step: Taq polymerase alone has an error rate of 1 in 1000 and RT has a very high error rate, worse it's the first step.

What DNA polymerase was used?

The key issue is separating the RT-PCR artefacts from naturally occurring RNA polymerase errors in situ is not possible without a control. If a frequency analysis cut-off strategy is used it needs to be an ultra-conservative threshold, i.e. I mean ridiculously high and this means lot of variants are going to be dumped.

Why? Again it's back to the problem whether the variants were generated in vitro, i.e. artefacts, or in situ, i.e. resulting from RNA polymerase. Generating a cut-off within a skewed distribution is not trivial per se but extremely difficult when the ratio between artefact and genuine variant has not been controlled.

There appears some sort of control. Maximising the threshold: thats a function of the control, so the control needs to accommodate a model of the different numbers of input variables; it will not be a linear function. Basically what input will give me the best representation of non-artifactual variants, that is what the control has to find out. Thus, for every input, you want a control.

Again, if the control is lacking, in any context, then the rule of thumb is to be ultra-conservative, not trying to maximise the number of variants, but minimise them, particularly if the control isn't perfect. Most virological reviewers raise concern with the numbers been reported here on the information given (obviously we don't know the virus).

  • NGS sequencing is another source of artefact - I should have mentioned there are some famous disasters here. I'm sure you NGS is really good, but a small number of investigators have deposited a disproportionate number of artefact ridden sequences - that everyone knows about.

from comments (again)

A control of 4 is tiny in data science especially a GLM approach. However, other studies would have been performed - possibly lots of them if this was e.g. HIV - and controls would have been run. There are two major strategies (really 3) for biologically cloning a virus. Really you'll need 100+ controls, that may not be feasible but as many within 4 to 100 as are available. GLM would be a good approach IMO. It would need stating that this does not truly control the rate of artefacts, but it's much closer.

Final comments ... don't want to assess thresholds.

I gave this question some thought and I think it would work via machine learning (what I specialise in). However, it may be difficult to find appropriate published controls and/or augmentation. Personally, I'm satisfied this would work if done correctly, but implementation is fairly advance ML. Given the amount of time I've spent on thinking this question through I wouldn't use an alternative approach, because I think it would be valid. Augmentation BTW is a generic term in machine learning and deep learning to increase the size of the training set.

  • $\begingroup$ Yes, they have controlled for some in vitro artefacts. I have plotted histograms of frequencies and each of the input variants. I got the expected picture you described, one variant taking a large portion of the reads a small number with a higher read count and a lot of low read count ones. I don't think I want to find a cutoff for the distribution, let's say I have one. I'm looking for a way to choose the number of input variants that maximizes observations over that threahold $\endgroup$ Jul 16, 2023 at 3:22
  • $\begingroup$ Following up on your edit. For each input I have a control but I only have 4 inputs. It is too expensive/time consuming to sequence too many more. Do I just use the controls on the 4 inputs I have? Or is there some way I can simulate behavior for unknown inputs? $\endgroup$ Jul 16, 2023 at 5:43
  • 1
    $\begingroup$ ok I think, I'm starting to pick up what you are saying. I have a few different options for what is considered a good observation. A hard threshold, item using some calculus (knee/inflection), as well as some methods based on spike ins (these are ideal probably). But lets say for each input we can compute the fraction of observations above the threshold. Then fit a GLM to these 4 points? 4 seems a bit small for a GLM right? $\endgroup$ Jul 16, 2023 at 5:54
  • $\begingroup$ Ok I think we are slightly talking about different things. I don't want to evaluate the controls or use them to select a threshold. I want to decide how many unique variants of my virus to sequence at once. Is that what your comment is describing and I'm just missing something? $\endgroup$ Jul 19, 2023 at 5:41
  • $\begingroup$ Ok @TheNumber23 answered in comments above. Upvotes are encouraged. $\endgroup$
    – M__
    Jul 19, 2023 at 7:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.