standard conditions from: Clay et al., 2003. Eur Biophys J
((2 * 3.14159 * R)/60)^2
In [18]:
import numpy as np
In [3]:
%load_ext rpy2.ipython
In [26]:
%%R
library(ggplot2)
library(dplyr)
In [16]:
# angular velocity: our setup
angular_vel_f = lambda R: (2 * 3.14159 * R / 60)
print angular_vel_f(55000)
# angular velocity: De Sario et al., 1995
print angular_vel_f(35000)
radius_max - radius_min = width_of_tube
In [26]:
%%R -w 14 -h 6 -u in
library(ggplot2)
library(reshape)
library(grid)
# radius top,bottom (cm)
r.top = 57.9 / 10
r.bottom = 71.1 / 10
# isoconcentration point
I = sqrt((r.top^2 + r.top * r.bottom + r.bottom^2)/3)
# rpm
R = 35000
# particle density
D = 1.70
# beta^o
B = 1.14e9
# dna in bp from 0.1kb - 100kb
L = seq(100, 100000, 100)
# angular velocity
## 2*pi*rpm / 60
w = ((2 * 3.14159 * R)/60)^2
# DNA GC content
G.C = seq(0.1, 0.9, 0.05)
# Molecular weight
# M.W in relation to GC content (dsDNA)
A = 313.2
T = 304.2
C = 289.2
G = 329.2
GC = G + C
AT = A + T
#GC2MW = function(x){ x*GC + (1-x)*AT + 157 } # assuming 5' monophosphate on end of molecules
GC2MW = function(x){ x*GC + (1-x)*AT }
M.W = sapply(G.C, GC2MW)
# buoyant density
## GC_fraction = (p - 1.66) / 0.098
GC2buoyant.density = function(x){ (x * 0.098) + 1.66 }
B.D = GC2buoyant.density(G.C)
# radius of the isoconcentration point from cfg center (AKA: r.p)
## position of the particle at equilibrium
buoyant.density2radius = function(x){ sqrt( ((x-D)*2*B/w) + I^2 ) }
P = buoyant.density2radius(B.D)
# calculating S
S.fun = function(L){ 2.8 + (0.00834 * (L*M.W)^0.479) }
S = t(sapply( L, S.fun ))
# calculating T
T = matrix(ncol=17, nrow=length(L))
for(i in 1:ncol(S)){
T[,i] = 1.13e14 * B * (D-1) / (R^4 * P[i]^2 * S[,i])
}
## formating
T = as.data.frame(T)
colnames(T) = G.C
T$dna_size__kb = L / 1000
T.m = melt(T, id.vars=c('dna_size__kb'))
colnames(T.m) = c('dna_size__kb', 'GC_content', 'time__h')
#T.m$GC_content = as.numeric(as.character(T.m$GC_content))
## plotting
p = ggplot(T.m, aes(dna_size__kb, time__h, color=GC_content, group=GC_content)) +
geom_line() +
scale_y_continuous(limits=c(0,175)) +
labs(x='DNA length (kb)', y='Time (hr)') +
scale_color_discrete(name='GC content') +
#geom_hline(yintercept=66, linetype='dashed', alpha=0.5) +
theme( text = element_text(size=18) )
#print(p)
# plotting at small scale
p.sub = ggplot(T.m, aes(dna_size__kb, time__h, color=GC_content, group=GC_content)) +
geom_line() +
scale_x_continuous(limits=c(0,5)) +
scale_y_continuous(limits=c(0,175)) +
labs(x='DNA length (kb)', y='Time (hr)') +
scale_color_discrete(name='GC content') +
#geom_hline(yintercept=66, linetype='dashed', alpha=0.5) +
theme(
text = element_text(size=14),
legend.position = 'none'
)
vp = viewport(width=0.43, height=0.52, x = 0.65, y = 0.68)
print(p)
print(p.sub, vp=vp)
In [196]:
%%R
# gas constant
R = 8.3144621e7 #J / mol*K
# temp
T = 273.15 + 23 # 23oC
# rotor speed (rpm)
S = 55000
# beta^o
beta = 1.14 * 10^-9
#beta = 1.195 * 10^-10
# G
G = 7.87 * 10^10 #cgs
# angular velocity
## 2*pi*rpm / 60
omega = 2 * pi * S /60
# GC
GC = seq(0,1,0.1)
# lengths
lens = seq(1000, 100000, 10000)
# molecular weight
GC2MW.dry = function(x){
A = 313.2
T = 304.2
C = 289.2
G = 329.2
GC = G + C
AT = A + T
x*GC + (1-x)*AT
}
M.dry = sapply(GC, GC2MW)
GC2MW.dryCS = function(n){
#n = number of bases
#base pair = 665 daltons
#base pair per dry cesium DNA = 665 * 4/3 ~= 882
return(n * 882)
}
M.dryCS = sapply(lens, GC2MW.dryCS)
# BD
GC2BD = function(x){
(x * 0.098) + 1.66
}
rho = sapply(GC, GC2BD)
# sd
calc_s.d = function(p=1.72, L=50000, T=298, B=1.195e9, G=7.87e-10, M=882){
R = 8.3145e7
x = (100 / 0.098)^2 * ((p*R*T)/(B^2*G*L*M))
return(x)
}
# run
p=seq(1.7, 1.75, 0.01)
L=seq(1000, 50000, 1000)
m = outer(p, L, calc_s.d)
rownames(m) = p
colnames(m) = L
In [85]:
%%R
# gas constant
R = 8.3144621e7 #J / mol*K
# temp
T = 273.15 + 23 # 23oC
# rotor speed (rpm)
S = 55000
# beta^o
beta = 1.14 * 10^-9
#beta = 1.195 * 10^-10
# G
G = 7.87 * 10^10 #cgs
# angular velocity
## 2*pi*rpm / 60
omega = 2 * pi * S /60
# GC
GC = seq(0,1,0.1)
# lengths
lens = seq(1000, 100000, 10000)
# molecular weight
GC2MW.dry = function(x){
A = 313.2
T = 304.2
C = 289.2
G = 329.2
GC = G + C
AT = A + T
x*GC + (1-x)*AT
}
#GC2MW = function(x){ x*GC + (1-x)*AT }
M.dry = sapply(GC, GC2MW.dry)
GC2MW.dryCS = function(n){
#n = number of bases
#base pair = 665 daltons
#base pair per dry cesium DNA = 665 * 4/3 ~= 882
return(n * 882)
}
M.dryCS = sapply(lens, GC2MW.dryCS)
# BD
GC2BD = function(x){
(x * 0.098) + 1.66
}
rho = sapply(GC, GC2BD)
# sd
calc_s.d = function(p=1.72, L=50000, T=298, B=1.195e9, G=7.87e-10, M=882){
R = 8.3145e7
x = (100 / 0.098)^2 * ((p*R*T)/(B^2*G*L*M))
return(sqrt(x))
}
# run
p=seq(1.7, 1.75, 0.01)
L=seq(500, 50000, 500)
m = outer(p, L, calc_s.d)
rownames(m) = p
colnames(m) = L
In [86]:
%%R
heatmap(m, Rowv=NA, Colv=NA)
In [87]:
%%R -w 500 -h 350
df = as.data.frame(list('fragment_length'=as.numeric(colnames(m)), 'GC_sd'=m[1,]))
#df$GC_sd = sqrt(df$GC_var)
ggplot(df, aes(fragment_length, GC_sd)) +
geom_line() +
geom_vline(xintercept=4000, linetype='dashed', alpha=0.6) +
labs(x='fragment length (bp)', y='G+C s.d.') +
theme(
text = element_text(size=16)
)
Notes:
In [70]:
%%R
calc_s.d = function(p=1.72, L=50000, T=298, B=1.195e9, G=7.87e-10, M=882){
R = 8.3145e7
sigma_sq = (p*R*T)/(B^2*G*L*M)
return(sqrt(sigma_sq))
}
# run
p=seq(1.7, 1.75, 0.01)
L=seq(500, 50000, 500)
m = outer(p, L, calc_s.d)
rownames(m) = p
colnames(m) = L
head(m)
In [66]:
%%R
heatmap(m, Rowv=NA, Colv=NA)
In [79]:
%%R -w 500 -h 350
BD50 = 0.098 * 0.5 + 1.66
df = as.data.frame(list('fragment_length'=as.numeric(colnames(m)), 'BD_sd'=m[BD50,]))
ggplot(df, aes(fragment_length, BD_sd)) +
geom_line() +
geom_vline(xintercept=4000, linetype='dashed', alpha=0.6) +
labs(x='fragment length (bp)', y='G+C s.d.') +
theme(
text = element_text(size=16)
)
In [ ]: