# Genomics Statistics Problem

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.

• It sounds like you're fishing for a statistic to prove a point. I wouldn't try to over-analyze your data. You might just consider stating that 18% of the sequences are in non-coding regions which is greater than the expected null hypothesis of 13% of the sequences being in non-coding regions. – conchoecia Jan 8 '19 at 6:10
• One could ask how big are those sequences of interest (SOI) and if there are some of them that span between coding and non-coding regions. Also the distribution of SOI could depend on TADs or other higher structures of your sequence (if it is a crhomosome size sequence). To model the distribution of SOI would require knowing a lot more about the SOI and the sequence you are working with and without a null model it cannot be tested if it is randomly following this model or not – llrs Jan 8 '19 at 7:55
• Hi again @conchoecia. Thanks for the feedback. With this point, however, we find ourselves in disagreement. There is a chance that if my results are due to random chance there may well be an uneven distribution. Say if my non-coding hits were 14% it would be obvious this would not be reliable evidence in support of my hypothesis. Furthermore, I have no problem with obtaining a negative result but however I interpret my results it needs to be statistically reliable, there is a chance that 18% is not significantly different to 13% and therefore cannot be used in favour of my hypothesis. – redoxio_ad_absurdum Jan 8 '19 at 15:58
• @llrs Hi. The sequences of interest are 29 bp (promotors) and I am searching chromosome FASTA. There are a some overlaps, but I have not counted them as non-coding for simplicity. The non-coding regions were defined by a different study (doi.org/10.1186/s13068-017-0742-z) Table 1. – redoxio_ad_absurdum Jan 8 '19 at 16:30
• Thanks everyone, question answered on cross-validated. stats.stackexchange.com/questions/386173/genomics-stats-problem/… – redoxio_ad_absurdum Jan 8 '19 at 16:59

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).

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.

• Hi, Michael. Thank you. Yes, I am trying binomial analysis now, don't worry about whether this results in a null hypothesis. I am simply trying to rule out a possible explanation for some results I obtained in the wet-lab. Whatever result I get is fine by me I am just happy to show I have properly investigated possible explanations. – redoxio_ad_absurdum Jan 8 '19 at 17:28
• Good to know. I'm definitely advocating DL for molecular solutions, they are very, very cool. – M__ Jan 8 '19 at 17:33
• Haha, I'm interested in DL however my data science, coding and statistics skills need a lot of work before I can use such a tool. But I read a lot about the applications and I'm keen to get to the point where I can use such tools. – redoxio_ad_absurdum Jan 8 '19 at 18:05