0
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.