Answers.Body: <p>The underlying answer to the question is have you experimentally controlled for <em>in vitro</em> 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 &quot;cloning a virus&quot; means, or Biology SE will explain this. The reason this is important is given below.</p>
  <p>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 &gt;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.</p>
  <p>The sources of artefact:</p>
  <ul>
  <li>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.</li>
  </ul>
  <p>What DNA polymerase was used?</p>
  <p>The key issue is separating the RT-PCR artefacts from naturally occurring RNA polymerase errors <em>in situ</em> 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.</p>
  <p><strong>Why?</strong> Again it's back to the problem whether the variants were generated <em>in vitro</em>, i.e. artefacts, or <em>in situ</em>, i.e. resulting from RNA polymerase. Generating a cut-off within a skewed distribution is not trivial <em>per se</em> but extremely difficult when the ratio between artefact and genuine variant has not been controlled.</p>
  <hr />
  <p>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.</p>
  <p>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).</p>
  <ul>
  <li>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.</li>
  </ul>
  <hr />
  <p>from comments (again)</p>
  <p>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.</p>
  <hr />
  <p>Final comments ... don't want to assess thresholds.</p>
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-16T03:22:49.083Z
Answers.Text: 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
Answers.UserDisplayName: TheNumber23
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-16T05:43:32.953Z
Answers.Text: 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?
Answers.UserDisplayName: TheNumber23
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-16T05:54:57.970Z
Answers.Text: 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?
Answers.UserDisplayName: TheNumber23
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-19T05:41:22.090Z
Answers.Text: 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?
Answers.UserDisplayName: TheNumber23
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-19T07:38:49.227Z
Answers.Text: Ok @TheNumber23 answered in comments above. Upvotes are encouraged.
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-16T03:10:58.563Z
Answers.LastActivityDate: 2023-07-19T07:49:56.943Z
Answers.LastEditDate: 2023-07-19T07:49:56.943Z
Answers.UserDisplayName: M__
Body: <p>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 &quot;real&quot;. I have some data on different input species counts (2k, 10k, 20k, 40k), and their corresponding fraction of &quot;real&quot;.</p>
  <p>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 &quot;real&quot; variants and total number of variants.</p>
  <p>Is there any research (or better yet code) that answers/discusses this problem?</p>
  <p>Thanks in advanced.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-15T23:22:51.367Z
Comments.Text: @M__ Thanks for the comment adjusted my question.
Comments.UserDisplayName: TheNumber23
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-15T23:59:43.137Z
Comments.Text: Added doing amplicon sequencing across different regions
Comments.UserDisplayName: TheNumber23
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-16T01:04:56.793Z
Comments.Text: 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.
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-16T01:26:06.220Z
Comments.Text: I see, using species not exactly clearly. They are all variants of the same virus. Will adjust
Comments.UserDisplayName: TheNumber23
Tags: <ngs><variant-calling><virus><pcr><qpcr>
Title: Choosing the best number of species for assay
ThreadText: # Choosing the best number of species for assay
  
  (Asked by: username_0 on 2023-07-15)
  
  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.
  
  **Tags:** ngs, variant-calling, virus, pcr, qpcr
  
  > **Comments:**
  > **username_0**: @username_1 Thanks for the comment adjusted my question.  
  > **username_0**: Added doing amplicon sequencing across different regions  
  > **username_1**: 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.  
  > **username_0**: I see, using species not exactly clearly. They are all variants of the same virus. Will adjust  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2023-07-16. Score: 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.
  
  > **Comments:**
  > **username_0**: 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  
  > **username_0**: 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?  
  > **username_0**: 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?  
  > **username_0**: 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?  
  > **username_1**: Ok @username_0 answered in comments above. Upvotes are encouraged.  
  
  ---
  
  