<p>This is quite a similar question to <a href="https://bioinformatics.stackexchange.com/questions/10845/rnaseq-normalization-stabilization-gene-length-and-rlog">this one</a>, but I notice that you're specifically asking about a <em>probability distribution</em> rather than differential expression. That sounds like gene dispersion and/or shot noise; a good reference for me on that has been the documentation for DESeq2:</p>
<p><a href="http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html" rel="nofollow noreferrer">http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html</a></p>
<p>The gene dispersion issue is discussed in more detail in their 2014 paper:</p>
<p><a href="https://doi.org/10.1186/s13059-014-0550-8" rel="nofollow noreferrer">https://doi.org/10.1186/s13059-014-0550-8</a></p>
<blockquote>
<p>We model read counts K ij as following a negative binomial distribution (sometimes also called a gamma-Poisson distribution) with mean μ ij and dispersion α i . The mean is taken as a quantity q ij, proportional to the concentration of cDNA fragments from the gene in the sample, scaled by a normalization factor s ij , i.e., μ ij =s ij q ij . For many applications, the same constant s j can be used for all genes in a sample, which then accounts for differences in sequencing depth between samples. To estimate these size factors, the DESeq2 package offers the median-of-ratios method already used in DESeq. However, it can be advantageous to calculate gene-specific normalization factors s ij to account for further sources of technical biases such as differing dependence on GC content, gene length or the like, using published methods, and these can be supplied instead.</p>
</blockquote>

CC BY-SA 4.0
2021-04-28T23:28:19.220Z
2021-04-28T23:28:19.220Z
gringer
<p>I am looking for a probability distribution of a number of reads mapped to a particular gene in metagenomic sequencing (NGS, shotgun, likely illumina).</p>
<p>Naively one could model it via a binomial (or Poisson) distribution as:
<span class="math-container">$$
P_n = {N\choose n}\pi^n(1-\pi)^{N-n},
$$</span>
where</p>
<ul>
<li><span class="math-container">$n$</span> is the number of reads mapped to the gene</li>
<li><span class="math-container">$N$</span> is the sequencing depth (total number of reads)</li>
<li><span class="math-container">$\pi\propto l$</span> is the probability of read mapping to the particular gene, proportioal to the gene length (possibly adjusted for the read length)</li>
</ul>
<p>This is however complicated by such intricacies of sequencing as:</p>
<ul>
<li>finite length of the fragments, which can be longer or shorter than the gene length</li>
<li>paired-end sequencing, which means one has to be careful when counting evens where only one or both reads map to a gene</li>
<li>something else?</li>
</ul>
<p>Is there an established model for this case? Could you recommend reading on the subject?</p>

<ngs><read-mapping><reads><illumina><genome-sequencing>
Modeling number of reads mapped to a gene
# Modeling number of reads mapped to a gene

(Asked by: username_0 on 2021-04-28)

I am looking for a probability distribution of a number of reads mapped to a particular gene in metagenomic sequencing (NGS, shotgun, likely illumina).

Naively one could model it via a binomial (or Poisson) distribution as:
$$
P\_n = {N\choose n}\pi^n(1-\pi)^{N-n},
$$
where

* $n$ is the number of reads mapped to the gene
* $N$ is the sequencing depth (total number of reads)
* $\pi\propto l$ is the probability of read mapping to the particular gene, proportioal to the gene length (possibly adjusted for the read length)

This is however complicated by such intricacies of sequencing as:

* finite length of the fragments, which can be longer or shorter than the gene length
* paired-end sequencing, which means one has to be careful when counting evens where only one or both reads map to a gene
* something else?

Is there an established model for this case? Could you recommend reading on the subject?

**Tags:** ngs, read-mapping, reads, illumina, genome-sequencing

---

## 1 answers

### Answer by username_1 on 2021-04-28. Score: 2

This is quite a similar question to [this one](https://bioinformatics.stackexchange.com/questions/10845/rnaseq-normalization-stabilization-gene-length-and-rlog), but I notice that you're specifically asking about a *probability distribution* rather than differential expression. That sounds like gene dispersion and/or shot noise; a good reference for me on that has been the documentation for DESeq2:

<http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html>

The gene dispersion issue is discussed in more detail in their 2014 paper:

<https://doi.org/10.1186/s13059-014-0550-8>

> We model read counts K ij as following a negative binomial distribution (sometimes also called a gamma-Poisson distribution) with mean μ ij and dispersion α i . The mean is taken as a quantity q ij, proportional to the concentration of cDNA fragments from the gene in the sample, scaled by a normalization factor s ij , i.e., μ ij =s ij q ij . For many applications, the same constant s j can be used for all genes in a sample, which then accounts for differences in sequencing depth between samples. To estimate these size factors, the DESeq2 package offers the median-of-ratios method already used in DESeq. However, it can be advantageous to calculate gene-specific normalization factors s ij to account for further sources of technical biases such as differing dependence on GC content, gene length or the like, using published methods, and these can be supplied instead.


---

