I am trying to reproduce this chart:

enter image description here

which depicts the probability of various a gene is insert in a mutant screen by a transposon. Apparently the graph is constructed from the following equation:

enter image description here

where Pi s the probability of finding one insert within a given gene, x is the average length of a gene transcript (1.7 kb), y is the total length of transcribed region in the gene space (113.1 Mb), Rtr is the probability of insertion into a transcription region within the gene space (0.573), Rgs is the probability of insertion into gene space in the whole genome (0.786), and N is the total number of Tnt1 inserts required.

However when I try and use this equation to generate the data, say for a 1.7Kbp gene I get this:

N <- seq(1, 20, 1) 
P <- 1 - (1 - ( 1.7 / 113.1 )) ^ (0.573 * 0.786 * N)
plot(N, P, type = "l")

enter image description here

Where the values are nothing like the chart. I have tried converting the kb and mb values to the same units but this did not help. I feel that I must be interpreting the equation wrong and must be super obvious and staring me in the face but I just cant see it.

To be clear, I am fine with recreating the aesthetic effects I just need help regenerating the data.

  • $\begingroup$ The notation is too small for me to tell, but I think N goes from 1 to 9*10^3, not from 1 to 20. Which if you plot that for the x-axis range the chart looks closer to the one you've posted. $\endgroup$ – Andy W Jun 26 '15 at 4:57

Andy W is correct.

I misunderstood the x-axis scale. Its actually 10^5, but I thought the exponential only applied to the last value on the x-axis to represent a saturation of transposons converging on prob = 1. To kind of demonstrate the asymptote.

Here is a demonstration that it is in fact correct.

N <- seq(1, 9e+05, 1) 
P <- 1 - (1 - ( 0.5 / 113.1 )) ^ (0.573 * 0.786 * N)
q <- 1 - (1 - ( 5 / 113.1 )) ^ (0.573 * 0.786 * N)
r <- 1 - (1 - ( 1.5 / 113.1 )) ^ (0.573 * 0.786 * N)
plot(N[1:5000], P[1:5000], type = "l")
lines(N[1:5000], q[1:5000], col = "red")
lines(N[1:5000], r[1:5000], col = "blue")

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.