Genomics Statistics Problem

Question

I have used a python script to identify target sequences in a DNA sequence file. There are two classes of sequence: coding and non-coding.

I have identified 728 sequences of interest. 597 of these fall into the coding regions and 131 of these fall into the non-coding regions. This is the equivalent of 18% non-coding. The total non-coding region in the sequence file is 13%.

Is there a statistical tool to demonstrate the python script identified target sequences in a non-random fashion?
If the script identified sequences that were randomly distributed then 13% of them would have been found in the non-coding region, from a total of 728 this seems like it should be reliable.

Answer 1

If your null hypothesis is that the event of being selected by your program is statistically independent from the event of being non-coding, you can use Fisher's exact test to reject the null hypothesis.

Suppose you have 10,000 (I'm making that up but presumably you have the real number) total sequences, 8,700 coding and 1,300 noncoding. Given that your program selects 728 sequences, you would expect around 13%*728 ~= 95 noncoding sequences. But your program picked 131 noncoding sequences and you want to know if 131 is far enough away from 95 to conclude that the difference is not just due to chance.

The probability of selecting 131 or more of these noncoding sequences (given that we are choosing 728 total sequences) is given by the so-called hypergeometric distribution. The two-sided test will also incorporate the probably of getting an equally unlikely low number (in this case, somewhere around 57) or less.

In R for example, you can calculate this using the 2-by-2 contingency table as follows:

fisher.test(rbind(c(597, 131), c(8103, 1169)))

which gives:

p-value = 7.399e-05

So in this situation (and assuming the probability model described), there's around a 0.0074% chance of observing a ratio as "extreme" as the one you observed, which indicates it's unlikely that the two events are independent.

Note: As others have mentioned, a binomial test could be appropriate too. If your number of total sequences is very large, then the two results will converge. But if it isn't, then Fisher's exact test will be more appropriate (essentially since you're sampling without replacement).

Answer 2

A Poisson or binomial would suffice.

I'm not at all sure the difference will significantly deviate from randomness.

Perl: 95% confidence interval of a non-parametric bootstrap for e.g. 1000 replicates is an alternative approach. I'd be more comfortable coding that solution in Perl than I would in Python, so is unlikely to be of much use here. Anyway, I would use a 2-dimensional array for your total data-set, the second dimension being coding/non-coding, I would then use the inbuilt random-number generator and multiply by the size of the array to randomly pluck out each element. 95% CI of the mean works by ranking the data and removing the top and bottom 2.5% and taking the interval for a single replicate (which is repeated 999 times more).

There will be a deep learning approach to this data, which would be easier to implement in Python, but I'd need to think about this because it is not simple classification.

Again my warning is that given the sample size, the result at best might be borderline. There are likely to be lots of interesting in your data, but coding/non-coding may not be one of them.