I am a bioinformatician and work on RNA-Seq data. The data contains a lot of reads (of length 80 bp in my case). These reads are fragments of those genes that were expressed. I map them back to my reference genome which helps answer many questions (gene expression, snp variants, alternative splicing etc...).
I work on alternative splicing, where I am interested in characterizing events where same genes express in more than 1 ways leading to different transcripts. Let's say I my data consists of 4 libraries (or biological replicates) for each of the 2 plant subspecies (closely related ones belonging to same species). And, I would like to look at exon skipping events (where certain exons from coding genes sometimes gets excluded).
My question is about the statistical analysis towards characterizing these exon skipping events between the two subspecies over those replicates. The basis of my question (or line of thinking and therefore the complexity), however, comes from statistical analysis on RNA-Seq data for gene expression. So, I'll frame the problem from there. Please bear with me.
For every gene, $g_{i}$, its expression is = number of reads mapped to that gene (or count data). For one gene, $g_{1}$, say, the data would look like this.
Lib1 Lib2 Lib3 Lib4
sA 400 420 600 250
sB 180 229 60 125
Since its count data, given the concentration $q_{1}$ of gene $g_{1}$, biostatisticians usually model the number of reads, $r_{1}$ that map to a $g_{1}$ as a poisson model. It can be shown that $r_{1} \propto q_{i}$, or $r_{1} = sq_{1}$ (poisson parameter). That is, $p(Y = r_{1} | sq_{1}) \sim poisson$. One question here: This means, that if I were to measure for the same gene, the read count and note down the concentration, then, I would end up in a poisson model, right?
From here, since we can not know the concentration $q_{1}$, they make use of the biological replicates, to estimate the variation. If $Q_{1}$ represents the concentrations over (all) the biological replicates, then under the assumption that this follows a gamma distribution (for the sake of mathematical convenience), the read counts $R_{1}$ for gene $g_{1}$ would then follow a negative binomial distribution. This is an overdispersed poisson model, it seems. Then a statistical test is performed to estimate if the difference in expression is statistically significant under a given $\alpha$. There are R-packages available that implement this.
Now, coming back to my problem, as opposed to looking at every gene, for every exon, $E_{i}$, I find out the number of times this exon is spliced out (= not included in the transcript) and the number of times it is included, in all libraries, in both subspecies.
So, for an exon, $E_{1}$, the data would look like this, for example: The first number in each entry is the number of exon skipping events and the second is the normal events.
Lib1 Lib2 Lib3 Lib4
sA 2, 80 1, 65 0, 40 2, 66
sB 10, 120 0, 22 8, 90 4, 90
The difference is that I have two values for each entry. I would like to find out, for each exon, if there is a difference between the two subspecies. My immediate thought was to calculate the ratio, for each entry, and then use existing model, however, they are on based on count data. So, my questions are these:
1) In general, is it known/ possible to estimate what distribution the ratio of two counts would follow?
2) Are there any other approaches ( for example, generalized linear models) which could from these data (including dispersion from replicates) help me calculate if the occurrence of exon skipping events between the two subspecies are statistically significant?
PS: In case somethings are not clear enough, I'd love to clarify them (I am by no means a statistician). I would appreciate any thoughts you guys have on this problem. Thanks again!
