Modeling number of reads mapped to a gene

Question

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?

Answer 1

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.