Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-28T23:28:19.220Z
Answers.LastActivityDate: 2021-04-28T23:28:19.220Z
Answers.UserDisplayName: gringer
Body: <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>
  
Tags: <ngs><read-mapping><reads><illumina><genome-sequencing>
Title: Modeling number of reads mapped to a gene
ThreadText: # 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.
  
  
  ---
  
  