# Mechanics behind deviation from random distribution

The system we are working on is biological, more specifically the distribution of programmed DNA damage events across a chromosome. This can be thought of as 1D array (the chromosome) across which points can be chosen (the sites of intentional damage). We have mapped the positions of these events experimentally and initially asked whether or not they fit a random distribution - that is, damage can occur at any point along the chromosome with equal chance and any given sites of damage are independent of one another. By generating random distributions on MATLAB (randi), this turned out not to be the case.

By analysing the inter-point distances (IPDs) from both the real- and modelled-datas, the real-data is seen to deviate from a random distribution only below a certain IPD size, before rejoining the random distribution above it i.e. there are fewer shorter IPDs than would be expected by chance in the real data. Example IPD results: Red = random modelled distribution
Blue = real data
Y-axis = IPD size (log-scale)
X-axis = IPD number (IPDs are just plotted in numerical order)


The IPDs are plotted here on a log Y-axis and simply in increasing order as if it were a histogram. As you can see below a certain IPD size (Y-axis), the blue line deviates from the red line.

The hypothesis that we are testing (which has a sound biological basis) is that the position of one event depends on the ones already formed. Specifically, as soon as a site is chosen, it invokes a zone of repression around it, making the surrounding region less likely to be chosen as the next site. This effectively spaces out the events and explains the absence of shorter IPDs. This zone gradually reduces in intensity the further you move away from a chosen point - explaining the return to independence above a certain IPD distance. Question: Is there a mathematical method by which we could derive the shape of this zone from the random and real datasets alone? For example, by calculating its strength (ability to deviate from randomness) at each given point until it's effects are no longer seen?

The shape and scale of the triangle in the above diagram is the main thing I am trying to obtain (it is not necessarily a triangle).

We have a second model which simulates this hypothesis - and which delivers promising results however we need guidance on the shape, scale etc. of the repression zone otherwise it is rather down to trial and error and multiple different windows + parameters can fit.

I have seen something similar done before by binning the IPDs into a histogram, fitting a gamma probability function and then converting this to a hazard function but I am not a mathematician and I do not know if this is the correct method nor how to go about it.

I largely work in MATLAB so if somebody could provide some help in the form of MATLAB, that would be great - but any help at all would be most appreciated.

Data used in the plot:

Real IPDs:

7126.5
11311.5
12582.25
21499
25429.25
28876.5
29178.5
35545.25
37498.75
37881.5
38152
45464
47372.5
48047.5
52397
55563
57100.75
59372
61640.5
63822.5
66672.25
67010
68969
69071.5
69680.75
70136
70228.25
75124
75487.5
76186.5
80091.5
80279
80727.75
83397.25
84412.25
84481
85453.5
85483.25
88821
88862.25
89089.5
90453.25
92416.25
96658
97369.75
98573.25
104459.5
105307.25
107716.5
113079.5
113357.75
113750.25
113848
114834.25
114871
114919.25
116882
116899.75
117400.75
113384.191
116714.9387
119898.1004
123046.5264
126504.6261
130069.3977
133819.0782
137747.762
141858.6185
146088.6625
150264.6261
154671.6308
159430.2967
164407.1167
169531.1443
174883.6052
180484.1524
186826.807
193794.4646
201090.8222
209380.867
218202.6614
228206.8165
239754.5876
252495.3356
267223.6972
285275.7581
308050.18
335997.8885
393927.4475
431000.091

Modelled IPDs:

6309.250317
7485.019638
8691.132742
9875.024811
11093.9262
12328.9784
13540.43008
14760.67732
16018.67552
17243.509
18560.20364
19830.60355
21235.71334
22592.75188
23931.62058
25240.54551
26572.1846
27899.31413
29311.17773
30765.96211
32251.92515
33713.78512
35191.37822
36695.70116
38301.07903
39893.27382
41474.13555
43128.17872
44764.51525
46449.33501
48116.12259
49799.81561
51567.24913
53351.51996
55228.92877
57039.44196
58826.45323
60615.27354
62437.5259
64364.0891
66308.25836
68317.33777
70389.35974
72571.9451
74659.85927
76782.19429
79186.51912
81427.22249
83761.00059
86187.90023
88672.44356
91239.82722
93885.18499
96423.67933
99062.67598
101676.3844
104409.6901
107253.7768
110233.3544
113384.191
116714.9387
119898.1004
123046.5264
126504.6261
130069.3977
133819.0782
137747.762
141858.6185
146088.6625
150264.6261
154671.6308
159430.2967
164407.1167
169531.1443
174883.6052
180484.1524
186826.807
193794.4646
201090.8222
209380.867
218202.6614
228206.8165
239754.5876
252495.3356
267223.6972
285275.7581
308050.18
335997.8885
393927.4475
431000.091

• The method you suggested is a standard way for fitting distributions. It is not clear to me whether time is important to you or not, though it seems it is, in which case you may be dealing with a non-homogeneous process. This will be trickier. – mandata May 21 '15 at 13:33
• I'm not sure I fully understand what you are asking regarding time. Could you elaborate? – AnnaSchumann May 21 '15 at 13:50
• Can you describe a little your phenomenon? What are you measuring? Also, it appears that by "random" you mean a certain distribution which you have in mind. The variable can be random but from a different distribution, which may produce thinner tails than you expect. – Aksakal May 21 '15 at 14:05
• @mandata My apologies - this was poorly worded. I have updated it via an edit. I meant that the intensity of the zone gradually reduces the further you go away from a chosen point - not that it dissipates over time. – AnnaSchumann May 21 '15 at 14:14

The problem is that you assumed a certain random distribution of IPD and it's not fitting the empirical distribution. So, the formulation of your question is a little confusing given the explanation you given so far. The "deviation" is not from randomness, but of the empirical distribution from the assumed theoretical one.

You generate locations $x_i\sim U(0,1000)$, where 0 and 1000 are bounds. Hence, the IPD is $\Delta x_i=|x_i-x_{i-1}|$.

We can find the unconditional probability of a small IPD
$$P(\Delta x_i)<\varepsilon$$ for any given small $\varepsilon>0$ as follows:

$$P(\Delta x_i)<\varepsilon=\frac{\varepsilon}{500}-\frac{\varepsilon^2}{1,000,000}$$

This is a peculiar distribution. Here's its cumulative and density functions:  The x-axis is IPD, and y-axis is cumulative (left) and density (right) probability functions.

As you can see your choice of model (i.e. randi function), implies that the probability of a small distance is quite high, much higher than of a large IPD. Your biological phenomenon is probably not fitting into this model. You've got try some other model.

• We initially tested to see whether or not our experimentally determined distributions matched a random distribution. They do not - and we know this. We are now trying to formulate a new model using the hypothesis stated in the OP (which has a specific biological basis to it). The main problem we have is that multiple combinations of shapes/scales for the repression zone fit the data - and we need to know which is correct hence i'm asking if it's possible to derive the shape/scale from the datasets alone. – AnnaSchumann May 22 '15 at 11:56
• @AnnaSchumann, look for a distribution which has a density lower for left tail than that of uniform distribution induced. For instance, if you noticed there's a mode in the distribution of IPDs, then you could start with Poisson, Neg Binomial, Lognormal or even normal just to see if the fit gets better. – Aksakal May 22 '15 at 12:25
• Thanks! I will take a look now. What are the axis labels on the graphs in your above answer? I am not sure I understand them fully. – AnnaSchumann May 25 '15 at 14:52
• @AnnaSchumann, updated the answer – Aksakal May 25 '15 at 15:01
• Wouldn't the probabilities of certain IPD sizes depend upon how many numbers are chosen by randi in each iteration? In our system we typically have quite large boundaries but perhaps will only pick 3-5 event sites. Is there a way to model the probability distributions for different amounts of sites chosen? – AnnaSchumann May 25 '15 at 15:08