# Modeling number of reads mapped to a gene

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?