text: user1813: 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?
  
  user73: This is quite a similar question to this one, 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.
  
  
original_text: 1813: <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>
   73: <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>
  