0
$\begingroup$

What is the correct way to compute a standard error of the ratio between random variables following multinomial and negative binomial distribution?

I asked a question about computing variance in the next-generation sequencing experiment and now I would like to check the standard error of the ratio of molecules in one of my experimental replicates.

$N_i$ is the sum of all molecules derived from the gene $i$ from all available replicates and $Y_{1i}$ is the number of molecules in the replicate 1. I assume that for each gene molecules are uniformly distributed across replicates.

How to compute the standard error of: $\frac{Y_{1i}}{N_i}\ $ ?

$\endgroup$

1 Answer 1

1
$\begingroup$

I'm not sure that the multinomial modeling used in response to your other question is the best or simplest way to think about this. And I think that dealing with the ratio that you present might be an unnecessary and distracting complication.

This is much more simply considered as a Poisson distribution, in which individual events occur at a constant average rate and are independently distributed in time, space, or volume.

If you start with a properly mixed reservoir with $N_i$ molecules and deliver a volume that contains on average $\bar Y_i$ molecules of type $i$ per replicate, then $Y_i$ follows a Poisson distribution with mean value $\bar Y_i$. With a Poisson distribution the variance of the distribution is exactly equal to the mean, $\bar Y_i$. The standard deviation of $Y_i$ equals $\sqrt{\bar Y_i}$.

So if you have on the average 100 molecules per replicate, you have a standard deviation of 10, or a percentage standard deviation (100 * SD/mean) of 10%. With 10,000 molecules per replicate you have a standard deviation of 100, a percentage standard deviation of only 1%. With 4 molecules per replicate you have a standard deviation of 2, a percentage standard deviation of 50%.

You can use the Poisson distribution to provide more details about the how many replicates are expected to contain each specific number of molecules, which can be used to check that the assumptions of the Poisson distribution are met.

If you have multiple types of molecules and they are distributed independently, then each of them can be described as above. If they are not distributed independently (e.g., some bind to each other) then things are more complicated.

You can scale that by the total number of molecules, $N_i$, if you wish, but then the result will depend on how many replicates you have. I think what's most important for next-generation sequencing that might get down close to single-molecule levels is the variance per replicate regardless of how much sample you started with or how many replicates you chose to perform.

$\endgroup$
4
  • $\begingroup$ Thank you. Main main concern and motivation of this question was setting a threshold to discard genes with low number of mapped molecules. I was wondering, whether the number of replicates should decrease the threshold (more replicates, higher certainty in quantification of low-counts genes). But based of you wrote it should predominantly depend on the depth of the sequencing rather than number of replicates. Am I right? $\endgroup$
    – hibernicah
    Aug 8, 2018 at 22:49
  • $\begingroup$ @hibernicah the total read count/sequencing depth is the critical factor. It should be possible to pool read-count information from multiple replicates to get more information for low-expression genes, but you might need to adjust for systematic processing differences from replicate to replicate (maybe a bigger problem than having low numbers of copies per se). Poisson sampling SD even at 25 molecules per replicate is only 20%, similar to typical biologic variability. Look into the Bioconductor project or other well established systems for analyzing NGS data. $\endgroup$
    – EdM
    Aug 8, 2018 at 23:56
  • $\begingroup$ From what I explored so far, it seems that there is no consensus on how to filter low-counts genes. I'm using DESeq2 to for my differential expression analysis and I observe that pre-filtering results in about 5-16% differences in the number of recognised differentially expressed genes (regardless of using build-in independent filtering on low-counts genes), which I suppose may be caused by highly noisy properties of my dataset (expected from the experimental conditions). $\endgroup$
    – hibernicah
    Aug 9, 2018 at 0:39
  • $\begingroup$ @hibernicah how much do you really care about differential expression of very low count genes? Of what practical value might they be? With so many more genes than conditions in most experiments, you might be able to find more highly expressed genes that serve your purposes. $\endgroup$
    – EdM
    Aug 9, 2018 at 1:23

Your Answer

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