Literature resourse
sedimentation coefficient of DNA (S)
In [1]:
%load_ext rpy2.ipython
In [135]:
%%R
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(rootSolve)
library(gridExtra)
In [80]:
%%R
SIPdbFile = '/home/nick/db/SIPdb/DOE_SIP_db.xlsx'
sheet.id = 'Fractions'
In [81]:
%%R
# tube characteristics
tube_diam__mm = 13
tube_height__mm = 48
tube_round_bottom_height__mm = 6.5
tube_capacity__ml = 4.7
tube_composition = 'polypropylene'
# rotor
rotor_id = 'TLA-110'
r_min__mm = 26
r_ave__mm = 37.2
r_max__mm = 48.5
frac_tube_angle = 90
# cfg run
## rpm of run
rpm = 55000
## angular velocity (w^2)
angular_velocity = 17545933.74
## average particle density
ave_gradient_density = 1.70
## beta^o
BetaO = 1.14e9 # CsCl at density of 1.70
## position of particle at equilibrium
particle_at_eq = 3.78
## max 13C shift
max_13C_shift_in_BD = 0.036
## min BD (that we care about)
min_GC = 13.5
min_BD = min_GC/100.0 * 0.098 + 1.66
## max BD (that we care about)
max_GC = 80
max_BD = max_GC / 100.0 * 0.098 + 1.66 # 80.0% G+C
max_BD = max_BD + max_13C_shift_in_BD
# diffusive boundary layer (DBL)
DBL_size_range__micron = c(10,100)
# misc
fraction_vol__cm3 = 0.1
In [82]:
%%R
# rotor angle
## sin(x) = opp / hypo
## x = sin**-1(opp/hypo)
rad2deg = function(rad) {
return((180 * rad) / pi)
}
deg2rad = function(deg){
return(deg * pi / 180)
}
x = r_max__mm - r_min__mm
hyp = tube_height__mm
rotor_tube_angle = rad2deg(asin(x / hyp))
cat("Tube angle from axis of rotation:", rotor_tube_angle, "\n")
In [83]:
%%R
x = r_max__mm - r_min__mm
hyp = tube_height__mm
asin(x / hyp)
In [7]:
%%R
# isoconcentration point
r_top = r_min__mm / 10
r_bottom = r_max__mm / 10
I__cm = sqrt((r_top**2 + r_top * r_bottom + r_bottom**2)/3)
cat('Isoconcentration point:', I__cm, '(cm)\n')
In [8]:
%%R
# loading SIPdb info
df.SIPdb = read_excel(SIPdbFile, sheet=sheet.id)
df.SIPdb %>% head(n=3) %>% as.data.frame
In [9]:
%%R -w 550 -h 300
# BD range
df.SIPdb = df.SIPdb %>%
group_by(sample_id, fractionation_date) %>%
arrange(fraction_id) %>%
mutate(BD_range = BD - lead(BD))
df.SIPdb$BD_range %>% summary %>% print
ggplot(df.SIPdb, aes(BD_range)) +
geom_histogram(binwidth=0.001) +
labs(x = 'gradient fraction size (delta BD)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [10]:
%%R
# fraction size range (1st & 3rd quartiles)
frac_size_range__BD = c(0.003, 0.005)
In [11]:
%%R
frac_size__ul = 100
frac_vol__cm3 = frac_size__ul / 1000
tube_radius__cm = tube_diam__mm / (2 * 10)
frac_height = frac_vol__cm3 / (pi * tube_radius__cm **2)
cat('The height of each fraction in vertically oriented tube:', frac_height, '(cm)\n')
cat('(Just applies to cylinder)\n')
In [12]:
%%R
# round bottom volume = (4/3 * pi * r**3) / 2
tube_radius__cm = tube_diam__mm / (2*10)
bottom_volume__cm3 = (4/3 * pi * tube_radius__cm**3) / 2
cat('Tube bottom volume:', bottom_volume__cm3, '(cm^3)\n')
cat('Number of fractions in tube bottom:', bottom_volume__cm3 / 0.1, '\n')
In [13]:
%%R
df.SIPdb.f = df.SIPdb %>%
filter(BD <= max_BD) %>%
arrange(fraction_id) %>%
group_by(sample_id, fractionation_date) %>%
summarize(BD_heaviest = first(BD)) %>% head %>% print
In [14]:
%%R
# plotting fraction vs fraction height in verticle tube
## calculating culmulative fraction height
df.SIPdb = df.SIPdb %>%
ungroup() %>%
mutate(frac_height = frac_height) %>%
group_by(sample_id, fractionation_date) %>%
arrange(fraction_id) %>%
mutate(cum_frac_height = cumsum(frac_height))
## first fraction kept (max BD fraction used for each gradient)
heavy.frac = df.SIPdb %>%
filter(BD <= max_BD) %>%
arrange(fraction_id) %>%
group_by(sample_id, fractionation_date) %>%
summarize(BD_heaviest = first(BD))
heavy.frac = inner_join(heavy.frac, df.SIPdb, c('sample_id' = 'sample_id',
'fractionation_date' = 'fractionation_date',
'BD_heaviest' = 'BD')) %>%
select(sample_id, fractionation_date, BD_heaviest, cum_frac_height, BD_range, fraction_id)
## plotting
ggplot(df.SIPdb, aes(cum_frac_height, BD_range)) +
geom_density2d() +
geom_rug(data=heavy.frac, sides='b', size=1, alpha=0.25, aes(color=fraction_id)) +
#geom_point(alpha=0.25) +
#stat_smooth(formula= y ~ poly(x, 2), method='lm') +
labs(x='Height of fraction\nfrom tube bottom (cm)', y='BD span of fraction') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [15]:
%%R
# summary of where in tube the heaviest fraction is used (cm from the bottom)
heavy.frac$cum_frac_height %>% summary %>% print
heavy.frac$fraction_id %>% summary %>% print
In [230]:
%%R
# in cm
frac_size = 0.01
tube_diam = 1.3
DBL_size = 0.01
In [231]:
%%R
# cylinder volume = pi * r**2 * h
tube_radius = tube_diam / 2
frac_vol = pi * tube_radius**2 * frac_size
nonDBL_vol = pi * (tube_radius - DBL_size)**2 * frac_size
DBL_vol = frac_vol - nonDBL_vol
DBL_to_frac = DBL_vol / frac_vol * 100
cat('Percent of fraction that is DBL: ', DBL_to_frac, '%\n', sep='')
In [18]:
%%R
r_t = r_min__mm / 10
r_b = r_max__mm / 10
t = 0.3 * (r_b - r_t)**2
cat('Time span of stability:', t, '(h)\n')
Density at a given radius from axis of rotation:
$ x = \big(\frac{w^2}{2*\beta^o} * (r_d^2 - I^2)\big) + D$
Radius from axis of rotation for a given density:
$ r_d = \sqrt{\frac{(x - D) * 2\beta^o}{w^2} + I^2} $
In [19]:
%%R
BD2radius = function(x, D, BetaO, w2, I){
sqrt( ((x-D)*2*BetaO/w2) + I^2)
}
min_BD_r = BD2radius(min_BD, ave_gradient_density, BetaO, angular_velocity, I__cm)
max_BD_r = BD2radius(max_BD, ave_gradient_density, BetaO, angular_velocity, I__cm)
cat('isoconcentration point:', I__cm, '\n')
cat('radius range for BD-min to BD-max: ', min_BD_r, 'to', max_BD_r, '\n')
In [20]:
%%R
tube_diam__cm = tube_diam__mm / 10
elps_main_axis__cm = tube_diam__cm / deg2rad(rotor_tube_angle)
cat('ellipsoid major axis length in tube (cylinder) cross section:', elps_main_axis__cm, '(cm)\n')
In [21]:
%%R
# radius position where rounded tube bottom affects ellipsoid major axis
## should just be round bottom diameter
round_bottom_diam__cm = tube_diam__mm / 10
r_round_bottom = r_max__mm / 10 - round_bottom_diam__cm
cat('radius position where round bottom comes into play:', r_round_bottom, '\n')
In [22]:
%%R
## Equation for lower point of band in the rounded part of tube
distFromAxis2angleTubePosOLD = function(x, r, D, A){
# x = a distance from the axis of rotation
# r = radius of cfg tube
# D = max tube distance from axis of rotation
# A = angle of tube to axis of rotation (degrees)
# Equation for finding the lower point of the band
if(x >= D-(r*aspace::cos_d(A))-r) {
if (x >= D-(r-r*aspace::sin_d(A))) {
# d = D - x
# hv = sqrt(d*(2*r-d))
# a = (90-A)-aspace::asin_d(hv/r)
d = x-(D-r)
a = A-aspace::asin_d(d/r)
}else if (x <= D-r){
# d = r-((D-x)-r)
# hv = sqrt(d*(2*r-d))
# a = (90+A)-aspace::asin_d(hv/r)
d = x-(D-r)
a = A-aspace::asin_d(d/r)
}else{
d = x-(D-r)
a = A-aspace::asin_d(d/r)
}
LowH = r-r*aspace::cos_d(a)
#print(LowH) ## This band will be in the rounded part
}else{
d = D-(r*aspace::cos_d(A))-r-x
hc = d/aspace::sin_d(A)
LowH = r+hc
# print(LowH) ## This band will be in the cylinder part
}
# Equation for finding the upper band
if(x > D-(r-r*aspace::cos_d(A))) {
d = D - x
hv = sqrt(d*(2*r-d))
a = (90-A)+aspace::asin_d(hv/r)
HighH = r-r*aspace::cos_d(a)
#print(HighH) ## This band will be in the rounded part
}else{
d = D-(r-r*aspace::cos_d(A))-x
hc = d/aspace::sin_d(A)
HighH = r+hc
#print(HighH) ## This band will be in the cylinder part
}
return(c(LowH, HighH))
}
# test
r = 0.65 # radius of tube
D = 4.85 # distance from axis of rotation to furthest part of tube
A = 37.95 # angle of tube to axis of rotation (degrees)
x = 3 # some distance from axis of rotation (from equation)
pos = distFromAxis2angleTubePosOLD(x, r, D, A)
pos %>% print
delta = pos[2] - pos[1]
delta %>% print
In [23]:
%%R
## Equation for lower point of band in the rounded part of tube
distFromAxis2angleTubePos = function(x, r, D, A){
# x = a distance from the axis of rotation
# r = radius of cfg tube
# D = max tube distance from axis of rotation
# A = angle of tube to axis of rotation (degrees)
# Equation for finding the lower point of the band
if(x >= D-(r*aspace::cos_d(A))-r) {
d = x-(D-r)
a = A-aspace::asin_d(d/r)
LowH = r-r*aspace::cos_d(a)
#print(LowH) ## This band will be in the rounded part
}else{
d = D-(r*aspace::cos_d(A))-r-x
hc = d/aspace::sin_d(A)
LowH = r+hc
# print(LowH) ## This band will be in the cylinder part
}
# Equation for finding the upper band
if(x > D-(r-r*aspace::cos_d(A))) {
d = x-(D-r)
a = (A)-(180-aspace::asin_d(d/r))
HighH = r-r*aspace::cos_d(a)
#print(HighH) ## This band will be in the rounded part
}else{
d = D-(r-r*aspace::cos_d(A))-x
hc = d/aspace::sin_xd(A)
HighH = r+hc
#print(HighH) ## This band will be in the cylinder part
}
return(c(LowH, HighH))
}
# test
r = 0.65 # radius of tube
D = 4.85 # distance from axis of rotation to furthest part of tube
A = 37.95 # angle of tube to axis of rotation (degrees)
x = 3 # some distance from axis of rotation (from equation)
pos = distFromAxis2angleTubePos(x, r, D, A)
pos %>% print
delta = pos[2] - pos[1]
delta %>% print
In [24]:
%%R
# params
tube_radius__cm = tube_diam__mm / (2*10)
tube_radius__cm %>% print
r_max__cm = r_max__mm / 10
## sequence of BD values in gradient
BDs = seq(min_BD, max_BD, 0.001)
## converting BD to position in gradient (assuming equilibrium)
rads = sapply(BDs, BD2radius,
D=ave_gradient_density,
BetaO=BetaO,
w2=angular_velocity,
I=I__cm)
# converting to band position in gradient to position in vertical tube height
angle.tube.pos = sapply(rads, distFromAxis2angleTubePos,
r=tube_radius__cm,
D=r_max__cm,
A=rotor_tube_angle)
angle.tube.pos = angle.tube.pos %>% t %>% as.data.frame
colnames(angle.tube.pos) = c('tube_low__cm', 'tube_high__cm')
angle.tube.pos$BD = BDs
angle.tube.pos$particle_pos = rads
angle.tube.pos %>% tail
In [25]:
%%R -w 600
# plotting
angle.tube.pos.g = angle.tube.pos %>%
gather(tube_pos, tube_pos__cm, tube_low__cm, tube_high__cm)
ggplot(angle.tube.pos.g, aes(particle_pos, tube_pos__cm, color=tube_pos, group=BD)) +
geom_point() +
geom_line(alpha=0.5) +
labs(x='Distance from axis (cm)', y='Band position on tube wall (cm)') +
theme_bw() +
theme(
text = element_text(size=18)
)
In [26]:
%%R
frac_size__ul = 100
frac_vol__cm3 = frac_size__ul / 1000
tube_radius__cm = tube_diam__mm / (2 * 10)
frac_height = frac_vol__cm3 / (pi * tube_radius__cm **2)
cat('The height of each fraction in vertically oriented tube:', frac_height, '(cm)\n')
cat('(Just applies to cylinder)\n')
In [27]:
%%R
# round bottom volume = (4/3 * pi * r**3) / 2
tube_radius__cm = tube_diam__mm / (2*10)
bottom_volume__cm3 = (4/3 * pi * tube_radius__cm**3) / 2
cat('Tube bottom volume:', bottom_volume__cm3, '(cm^3)\n')
cat('Number of fractions in tube bottom:', bottom_volume__cm3 / 0.1, '\n')
# need to determine the tube height of each of these first 6 fractions
In [28]:
%%R
df.SIPdb.f = df.SIPdb %>%
filter(BD <= max_BD) %>%
arrange(fraction_id) %>%
group_by(sample_id, fractionation_date) %>%
summarize(BD_heaviest = first(BD)) %>% head %>% print
In [29]:
%%R
# relationship: volume ~ height
r = 0.65
f = function(x){x**3 - 3*r*x**2 + 3*0.1/pi}
uniroot.all(f, c(0, r*2)) %>% print
curve((pi*x**2)/3 * (3*r - x), from=-1, to=4)
In [30]:
%%R
# converting cylinder volume to cylinder height (cm)
cyl_vol2height = function(v, r){v / (pi * r**2)}
cyl_vol2height(0.1, 0.65)
In [177]:
%%R
# converting sphere cap volume to sphere height (cm)
sphr_vol2height = function(v, r){
# h**3 - 3*r*h**2 + (3v / pi) = 0
f = function(x){x**3 - 3*r*x**2 + 3*v/pi}
roots = uniroot.all(f, c(0, r*2))
if(length(roots) > 1){
root_str = paste(roots, sep=',')
cat('ERROR: Number of roots > 1 for volume:', v, '\n roots:', root_str, '\n')
stop()
} else if (length(roots) == 0){
cat('WARNING: no roots for volume:', v, '\n')
roots = NA
}
return(roots[1])
}
# testing
sphere_half_vol = (4/3 * pi * 0.65**3)/2
sphr_vol2height(0.1, 0.65) %>% print
sphr_vol2height(0.5, 0.65) %>% print
sphr_vol2height(sphere_half_vol, 0.65) %>% print
In [32]:
%%R
# -- check on calculation of sphere hight from bottom --
# height of rounded tube bottom
h_bottom = tube_diam__mm / (2*10)
# volume of rounded portion of tube
v_bottom = (4/3 * pi * h_bottom**3)/2
stopifnot(sphr_vol2height(v_bottom, 0.65) == h_bottom)
In [33]:
%%R
# function to convert volume filling cfg tube to tube height
tubeVol2height = function(v, r=0.65){
# v = volume (cm^3)
# r = sphere radius
stopifnot(length(v) == 1)
sphr_cap_vol = (4/3 * pi * r**3)/2
if(v <= sphr_cap_vol){
# height does not extend to cylinder
h = sphr_vol2height(v, r)
} else {
# height = sphere_cap + cylinder
sphr_cap_height = sphr_vol2height(sphr_cap_vol, r)
h = sphr_cap_height + cyl_vol2height(v - sphr_cap_vol, r)
}
return(h)
}
# test
tube_radius__cm = tube_diam__mm / (2*10)
vol = 0.1 # 100 ul
vols = seq(0, 4, 0.1)
sapply(vols, tubeVol2height, r=tube_radius__cm)
In [34]:
%%R
# plotting
n_fracs = 32
x = seq(0, 0.1*(n_fracs+1), 0.1)
df = data.frame('x' = x, 'y' = sapply(x, tubeVol2height))
ggplot(df, aes(x, y)) +
geom_point(color='blue') +
geom_line(color='blue') +
stat_function(fun = function(x) x, linetype='dashed', alpha=0.7) +
labs(x='Culmulative fraction volume (ml)',
y='Tube height (cm)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [35]:
%%R
# loading SIPdb info
df.SIPdb = read_excel(SIPdbFile, sheet=sheet.id)
df.SIPdb %>% head(n=3) %>% as.data.frame
In [36]:
%%R
# BD range
df.SIPdb = df.SIPdb %>%
group_by(sample_id, fractionation_date) %>%
arrange(fraction_id) %>%
mutate(BD_range = BD - lead(BD))
In [37]:
%%R
# round bottom volume = (4/3 * pi * r**3) / 2
tube_radius__cm = tube_diam__mm / (2*10)
bottom_volume__cm3 = (4/3 * pi * tube_radius__cm**3) / 2
cat('Tube bottom volume:', bottom_volume__cm3, '(cm^3)\n')
cat('Number of fractions in tube bottom:', bottom_volume__cm3 / 0.1, '\n')
In [38]:
%%R
# plotting fraction vs fraction height in verticle tube
## calculating culmulative fraction height
df.SIPdb = df.SIPdb %>%
ungroup() %>%
mutate(frac_vol = fraction_vol__cm3) %>%
group_by(sample_id, fractionation_date) %>%
arrange(fraction_id) %>%
mutate(culm_frac_vol = cumsum(frac_vol))
## plotting
ggplot(df.SIPdb, aes(BD, culm_frac_vol)) +
geom_point() +
geom_smooth() +
labs(y='Culmulative fraction volume (cm^3)', x='fraction BD') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [39]:
%%R
# curve fit for: culm_frac_vol ~ fraction_BD
dataset = df.SIPdb %>%
select(BD, culm_frac_vol) %>%
rename('x' = BD, 'y' = culm_frac_vol) %>%
as.data.frame
fit = lm(y~x+I(x^2)+I(x^3), data=dataset)
fit %>% print
dataset$y_pred = predict(fit)
ggplot(dataset, aes(x)) +
geom_point(aes(y=y)) +
geom_line(aes(y=y_pred), color='red') +
labs(y='Culmulative fraction volume (cm^3)', x='fraction BD') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [40]:
%%R
BD2tubeVol = function(BD){
if(BD > 1.809){
cat('WARNING: BD value', BD, 'is > 1.809, where tube volume = 0',
'Setting BD to 1.809\n')
BD = 1.809
}
tube_vol = 545.9*BD**3 - 2771.9*BD**2 + 4664.9*BD - 2599.5
if(tube_vol < 0){
cat('WARNING: tube volume < 0 for BD: ', BD,
'\n Setting volume to 0\n')
tube_vol = 0
}
return(tube_vol)
}
# check
BD2tubeVol(1.65) %>% print
BD2tubeVol(1.7) %>% print
BD2tubeVol(1.75) %>% print
BD2tubeVol(1.8) %>% print
BD2tubeVol(1.9) %>% print
In [41]:
%%R
BD_gradient = rev(seq(min_BD, max_BD, 0.001))
# wrapper function for BD2tubeVol & tubeVol2height
BD2tubeHeight = function(BD, r){
stopifnot(length(BD) == 1)
v = BD2tubeVol(BD)
tubeVol2height(v, r)
}
h = sapply(BD_gradient, BD2tubeHeight, r=tube_radius__cm)
df = data.frame('BD' = BD_gradient, 'height' = h)
ggplot(df, aes(BD, h)) +
geom_point() +
labs(x='Buoyant density',
y='Tube height (cm)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [42]:
%%R -w 400 -h 350
# plotting relationship between BD and tube height for SIPdb data
df.SIPdb = df.SIPdb %>%
mutate(tube_height__cm = sapply(BD, BD2tubeHeight, r=tube_radius__cm))
## plotting
ggplot(df.SIPdb, aes(BD, tube_height__cm, color=-fraction_id)) +
#geom_point(alpha=0.5, shape='O') +
geom_point() +
geom_vline(xintercept=min_BD, linetype='dashed', alpha=0.5) +
geom_vline(xintercept=max_BD, linetype='dashed', alpha=0.5) +
labs(title='All SIPdb gradient data',
y='Height of fraction\nfrom tube bottom (cm)',
x='Fraction BD') +
theme_bw() +
theme(
text = element_text(size=16),
legend.position = 'none'
)
In [43]:
%%R -w 600
# do the BD bands expand further up in the gradient?
df.SIPdb = df.SIPdb %>%
group_by(sample_id, fractionation_date) %>%
arrange(fraction_id) %>%
mutate(delta_tube_height = lead(tube_height__cm) - tube_height__cm)
df.SIPdb.f = df.SIPdb %>%
filter(delta_tube_height > 0)
ggplot(df.SIPdb.f, aes(BD, delta_tube_height, color=tube_height__cm)) +
geom_point() +
scale_x_continuous(limits=c(1.65, 1.80)) +
scale_y_continuous(limits=c(0, 0.2)) +
labs(title='All SIPdb data',
x='Buoyant density',
y='Change in tube height between fractions (cm)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [69]:
%%R
height2sphereCapVolume = function(h, r){
# converting height of spherical cap to volume
# h = height in sphere
# r = radius of sphere
(pi*h**2)/3 * (3*r - h)
}
# test
r = 0.65
half_sphere = 4/3 * pi * r**3 / 2 # volume of half a sphere
half_sphere = round(half_sphere, 6)
sphere_cap = height2sphereCapVolume(r, r)
sphere_cap = round(sphere_cap, 6)
stopifnot(half_sphere == sphere_cap)
In [70]:
%%R
height2cylinderVol = function(h, r){
# h = cylinder height
# r = cylinder radius
pi * r**2 * h
}
# test
r = 0.65
height2cylinderVol(r, r)
In [73]:
%%R
tubeHeight2tubeVolume = function(h, r=0.65){
# function to convert vertical tube height to vertical tube volume (up to that height)
# h = tube height (cm)
# r = tube radius (cm) # assuming spherical tube bottom
stopifnot(length(h) == 1)
half_sphere_vol = 4/3 * pi * r**3 / 2
if(h <= r){
# height does not extend to cylinder; just calculating for spherical bottom
v = height2sphereCapVolume(h, r)
} else {
# volume = sphere_half_volume + cylinder_volume
v = half_sphere_vol + height2cylinderVol(h - r, r)
}
return(v)
}
# test
tubeHeight2tubeVolume(r)
In [ ]:
%%R
tubeVolume2BD = function(){
# function to convert tube volume to BD
}
In [44]:
%%R
## df of tube position for gradient in rotor (assuming equilibrium)
angle.tube.pos %>% head(n=4)
In [59]:
%%R -w 600
# converting BD to tube height for vertical gradient
angle.tube.pos$tube_vert__cm = sapply(angle.tube.pos$BD, BD2tubeHeight, r=tube_radius__cm)
# plotting
df.g = angle.tube.pos %>%
gather(tube_angle_pos, tube_angle_pos__cm, tube_low__cm, tube_high__cm)
ggplot(df.g, aes(tube_angle_pos__cm, tube_vert__cm, color=tube_angle_pos, group=BD)) +
geom_line(alpha=0.25, color='black') +
geom_point() +
labs(x='Band position on angled tube (cm)',
y='Band position in vertical tube (cm)') +
theme_bw() +
theme(
text = element_text(size=18)
)
In [94]:
%%R -w 600
# converting BD to tube height for vertical gradient
angle.tube.pos$tube_vert__cm = sapply(angle.tube.pos$BD, BD2tubeHeight, r=tube_radius__cm)
# plotting
df.g = angle.tube.pos %>%
gather(tube_pos, tube_pos__cm, tube_low__cm, tube_high__cm, tube_vert__cm)
# min/max of vertical tube height
df.min.max = angle.tube.pos %>%
select(BD, tube_vert__cm) %>%
group_by() %>%
mutate(min_tube_vert__cm = min(tube_vert__cm),
max_tube_vert__cm = max(tube_vert__cm)) %>%
distinct(min_tube_vert__cm, max_tube_vert__cm)
ggplot(df.g, aes(BD, tube_pos__cm, color=tube_pos, group=BD)) +
geom_line(alpha=0.25, color='black') +
geom_point() +
geom_hline(data=df.min.max, aes(yintercept=min_tube_vert__cm), alpha=0.5, linetype='dashed') +
geom_hline(data=df.min.max, aes(yintercept=max_tube_vert__cm), alpha=0.5, linetype='dashed') +
labs(y='Tube vertical height (cm)', x='Buoyant density') +
theme_bw() +
theme(
text = element_text(size=18)
)
Conclusions
In [107]:
%%R
(0.0362)**2 / (0.9 * 2e-5)
#F = function(z, D) z**2/(0.9*D)
#curve(F, )
In [13]:
%%R -h 300
length2MW = function(L){ L * 666 }
length2sedCoef = function(L){
2.8 + (0.00834 * (L*666)**0.479)
}
MW2diffuseCoef = function(L, p, R=8.3144598, T=293.15){
V = 1/1.99
M = length2MW(L)
s = length2sedCoef(L)
(R*T)/(M*(1-V*p)) * s
}
# test
L = seq(100, 50000, 100)
p = 1.7
D = sapply(L, MW2diffuseCoef, p=p)
df = data.frame('L' = L, 'D' = D)
# plotting
ggplot(df, aes(L, D)) +
geom_point() +
geom_line(alpha=0.5) +
theme_bw() +
theme(
text = element_text(size=16)
)
In [54]:
%%R
# converting D to cm^2/s
df$D_cm = df$D * 1e-5
# time periods (sec)
t = seq(1, 300, 10)
# calculating z (cm)
ES = function(D, t){
sqrt(0.9 * D * t)
}
df2 = expand.grid(df$D_cm, t)
colnames(df2) = c('D_cm', 't')
df2$z = mapply(ES, df2$D_cm, df2$t)
tmp = expand.grid(df$L, t)
# adding variable
df2$L = tmp$Var1
df2$t_min = df2$t / 60
df2$z_uM = df2$z / 1e-5
## plotting
ggplot(df2, aes(t_min, z_uM, color=L, group=L)) +
#geom_point(size=1.5) +
geom_line() +
labs(x='Time (minutes)',
y='mean deviation of molecules\nfrom starting position (uM)') +
scale_color_continuous('DNA fragment\nlength (bp)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [55]:
%%R -w 800
## plotting
ggplot(df2, aes(L, z_uM, color=t_min, group=t_min)) +
#geom_point(size=1.5) +
geom_line() +
labs(x='DNA fragment length (bp)',
y='mean deviation of molecules\nfrom starting position (uM)') +
scale_color_continuous('Time\n(minutes)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [ ]:
In [535]:
%%R
## Libraries
library(aspace)
## Given values
r = 0.65 # radius of tube
D = 4.85 # distance from axis of rotation to furthest part of tube
A = 37.95 # angle of tube to axis of rotation (degrees)
x = 4.32 # some distance from axis of rotation (from equation)
## Equation for lower point of band in the rounded part of tube
if(x >= D-(r*cos_d(A))-r) {
d = D - x
hv = sqrt(d*(2*r-d))
if (x >= D-(r-r*sin_d(A))) {
a = (90-A)-asin_d(hv/r)
}
else if (x < D-(r-r*sin_d(A))) {
a = -(90-A)-asin_d(hv/r)
}
LowH = r-r*cos_d(a)
print(LowH)
}
## Equation for lower point of band in column part of tube
if(x < D-(r*cos_d(A))-r) {
d = D-x-r
hc = d/sin_d(A)
LowH = r+hc
print(LowH)
}
## Equation for upper point of band in rounded part of tube
if(x > D-(r-r*cos_d(A))) {
d = D - x
hv = sqrt(d*(2*r-d))
a = (90-A)+asin_d(hv/r)
HighH = r-r*cos_d(a)
print(HighH)
}
## Equation for upper point of band in cylindrical part of tube
if(x <= D-(r-r*cos_d(A))) {
d = D-(r-r*cos_d(A))-x
hc = d/sin_d(A)
HighH = r+hc
print(HighH)
}
In [6]:
%%R
## Equation for lower point of band in the rounded part of tube
distFromAxis2angleTubePos = function(x, r, D, A){
# x = a distance from the axis of rotation
# r = radius of cfg tube
# D = max tube distance from axis of rotation
# A = angle of tube to axis of rotation (degrees)
if(x >= D-(r*aspace::cos_d(A))-r) {
d = D - x
hv = sqrt(d*(2*r-d))
if (x >= D-(r-r*aspace::sin_d(A))) {
a = (90-A)-aspace::asin_d(hv/r)
}
else if (x < D-(r-r*aspace::sin_d(A))) {
a = -(90-A)-aspace::asin_d(hv/r)
}
LowH = r-r*aspace::cos_d(a)
#print(LowH)
}
## Equation for lower point of band in column part of tube
if(x < D-(r*aspace::cos_d(A))-r) {
d = D-x-r
hc = d/aspace::sin_d(A)
LowH = r+hc
#print(LowH)
}
## Equation for upper point of band in rounded part of tube
if(x > D-(r-r*aspace::cos_d(A))) {
d = D - x
hv = sqrt(d*(2*r-d))
a = (90-A)+aspace::asin_d(hv/r)
HighH = r-r*aspace::cos_d(a)
#print(HighH)
}
## Equation for upper point of band in cylindrical part of tube
if(x <= D-(r-r*aspace::cos_d(A))) {
d = D-(r-r*aspace::cos_d(A))-x
hc = d/aspace::sin_d(A)
HighH = r+hc
#print(HighH)
}
return(c(LowH, HighH))
}
# test
r = 0.65 # radius of tube
D = 4.85 # distance from axis of rotation to furthest part of tube
A = 37.95 # angle of tube to axis of rotation (degrees)
x = 4.32 # some distance from axis of rotation (from equation)
distFromAxis2angleTubePos(x, r, D, A)
In [7]:
%%R
x = seq(0,3, 0.5)
r = 0.65 # radius of tube
D = 4.85 # distance from axis of rotation to furthest part of tube
A = 27.95 # angle of tube to axis of rotation (degrees)
sapply(x, distFromAxis2angleTubePos, r=r, D=D, A=A)
In [169]:
import numpy as np
from scipy.integrate import quad
# numpy-dependent code
def distFromAxis2angleTubeVol(x, r, A, D):
# x = distance from axis of rotation
# r = tube radius
# A = tube angle
# D = max distance from axis of rotation
a = np.deg2rad(A)
p1 = r-(r*np.cos(a))
p2 = r+(r*np.cos(a))
R12 = p2-p1
d = D-x
#print p1, p2
def SphVol(X, r, p2, R12):
return (X*((2*r)-X)/2) * ((2*np.pi*((p2-X)/R12)) - np.sin(2*np.pi*((p2-X)/R12)))
def CylWedVol(X, r, a, p1):
return 2*(X-r+((X-p1)/np.cos(a))) / np.tan(a) * np.sqrt(r**2-X**2)
if x < (D-p2):
h1 = (p2-x)/np.sin(a)
h2 = ((D-p1)-x)/np.sin(a)
area1 = (2/3)*np.pi*r**3
area2 = (1/2)*np.pi*r**2*(h1+h2)
area = area1+area2
#print h1, h2, area1, area2
elif x > (D-p2):
area1 = (1/3)*np.pi*p1**2*(3*r-p1)
area2 = quad(SphVol, p1, d, args=(r, p2, R12))
area3 = quad(CylWedVol, r, (r-((d-p1)/np.cos(a))), args=(r, a, p1))
area = area1+area2+area3
elif x > (D-p1):
area = (1/3)*np.pi*d**2*(3*r-d)
return area
# test
r = 0.65 # radius of tube (cm)
D = 4.85 # distance from axis of rotation to furthest part of tube
A = 27.95 # angle of tube to axis of rotation (degrees)
x = 3.0 # some distance from axis of rotation (from equation)
distFromAxis2angleTubeVol(x, r, D, A)
Out[169]:
In [170]:
# converting cylinder volume to height
def cylVol2height(v, r):
# v = volume (ml)
# r = tube radius (cm)
h = v / (np.pi * r**2)
return h
# test
cylVol2height(0.1, 0.65)
Out[170]:
In [173]:
# converting sphere cap volume to sphere height
from scipy import optimize
def sphereCapVol2height(v, r):
# v = volume (ml)
# r = tube radius (cm)
# h**3 - 3*r*h**2 + (3v / pi) = 0
f = lambda x : x**3 - 3*r*x**2 + 3*v/np.pi
try:
root = optimize.brentq(f, 0, r*2, maxiter=1000)
except ValueError:
msg = 'WARNING: not roots for volume {}\n'
sys.stderr.write(msg.format(v))
root = np.nan
return(root)
# test
sphereCapVol2heightV = np.vectorize(sphereCapVol2height)
heights = np.arange(0, r**2, 0.1)
sphereCapVol2heightV(heights, 0.65)
Out[173]:
In [121]:
# convert liquid volume in vertical cfg tube to tube height
def tubeVol2height(v, r):
# v = volume (ml)
# r = tube radius (cm)
sphere_cap_vol = (4/3 * np.pi * r**3)/2
if v <= sphere_cap_vol:
# height does not extend to cylinder
h = sphereCapVol2height(v, r)
else:
# height = sphere_cap + cylinder
sphere_cap_height = sphereCapVol2height(sphere_cap_vol, r)
h = sphere_cap_height + cylVol2height(v - sphere_cap_vol, r)
return(h)
# test
vol = 0.1 # 100 ul
vols = np.arange(0, 4+vol, vol)
tubeVol2heightV = np.vectorize(tubeVol2height)
tubeVol2heightV(vols, r=0.65)
Out[121]:
In [123]:
BD_vert = """BD,vert_tube_pos
1.5,9.48114816107
1.501,8.82939732893
1.502,8.5139232829
1.503,8.27086657273
1.504,8.06562489985
1.505,7.88464192336
1.506,7.72092856414
1.507,7.57032006486
1.508,7.43009721563
1.509,7.29836842716
1.51,7.17375489298
1.511,7.05521480775
1.512,6.94193821822
1.513,6.83328063528
1.514,6.72871925446
1.515,6.62782303112
1.516,6.53023160098
1.517,6.43564004601
1.518,6.34378763861
1.519,6.25444936373
1.52,6.16742942407
1.521,6.08255618938
1.522,5.9996782158
1.523,5.91866107041
1.524,5.83938477069
1.525,5.76174169905
1.526,5.68563488916
1.527,5.610976606
1.528,5.53768716048
1.529,5.46569391276
1.53,5.39493042902
1.531,5.32533576355
1.532,5.25685384425
1.533,5.18943294389
1.534,5.12302522289
1.535,5.0575863321
1.536,4.99307506643
1.537,4.92945306131
1.538,4.86668452591
1.539,4.80473600771
1.54,4.74357618404
1.541,4.68317567694
1.542,4.62350688823
1.543,4.5645438521
1.544,4.50626210304
1.545,4.44863855724
1.546,4.39165140575
1.547,4.335280018
1.548,4.27950485459
1.549,4.22430738813
1.55,4.16967003135
1.551,4.11557607165
1.552,4.06200961137
1.553,4.00895551324
1.554,3.9563993504
1.555,3.90432736057
1.556,3.85272640399
1.557,3.80158392469
1.558,3.75088791479
1.559,3.70062688162
1.56,3.65078981723
1.561,3.60136617031
1.562,3.55234582006
1.563,3.50371905201
1.564,3.45547653553
1.565,3.40760930293
1.566,3.36010873002
1.567,3.31296651792
1.568,3.26617467618
1.569,3.21972550695
1.57,3.17361159019
1.571,3.12782576985
1.572,3.08236114087
1.573,3.03721103706
1.574,2.9923690197
1.575,2.94782886677
1.576,2.90358456291
1.577,2.85963028993
1.578,2.81596041788
1.579,2.77256949663
1.58,2.72945224793
1.581,2.68660355792
1.582,2.64401847003
1.583,2.60169217832
1.584,2.55962002113
1.585,2.5177974751
1.586,2.47622014947
1.587,2.43488378073
1.588,2.39378422748
1.589,2.35291746564
1.59,2.3122795838
1.591,2.27186677887
1.592,2.23167535195
1.593,2.19170170432
1.594,2.15194233377
1.595,2.1123938309
1.596,2.07305287583
1.597,2.03391623486
1.598,1.9949807574
1.599,1.95624337303
1.6,1.91770108862
1.601,1.87935098569
1.602,1.84122425323
1.603,1.80344903231
1.604,1.76610234807
1.605,1.72923269162
1.606,1.69287464061
1.607,1.65705312048
1.608,1.62178546514
1.609,1.58708264414
1.61,1.5529501074
1.611,1.51938843709
1.612,1.48639389738
1.613,1.4539589297
1.614,1.42207261839
1.615,1.39072114057
1.616,1.35988820676
1.617,1.32955549481
1.618,1.29970307742
1.619,1.27030984166
1.62,1.24135389839
1.621,1.21281297865
1.622,1.18466481413
1.623,1.15688749867
1.624,1.12945982786
1.625,1.10236161426
1.626,1.07557397568
1.627,1.0490795947
1.628,1.02286294766
1.629,0.996910501947
1.63,0.971210880723
1.631,0.945754994636
1.632,0.920536140475
1.633,0.895550067055
1.634,0.870795009008
1.635,0.846271689462
1.636,0.821983292901
1.637,0.797935409751
1.638,0.774135954503
1.639,0.75059505936
1.64,0.727324945596
1.641,0.70433977494
1.642,0.681655483401
1.643,0.65928960005
1.644,0.637261053281
1.645,0.615589967122
1.646,0.594297450135
1.647,0.573405379394
1.648,0.55293618197
1.649,0.53267033275
1.65,0.512331039339
1.651,0.492227332309
1.652,0.472352659119
1.653,0.452702217388
1.654,0.433272723828
1.655,0.414062219375
1.656,0.395069905454
1.657,0.376296007692
1.658,0.357741664502
1.659,0.339408838938
1.66,0.321300253076
1.661,0.303419345041
1.662,0.28577024968
1.663,0.268357804786
1.664,0.251187585845
1.665,0.234265973321
1.666,0.217600257475
1.667,0.201198786109
1.668,0.185071159176
1.669,0.169228467638
1.67,0.153683553869
1.671,0.13845121426
1.672,0.123548105236
1.673,0.108991638882
1.674,0.0947955527856
1.675,0.0809527300524
1.676,0.0673390171276
1.677,0.0537732836508
1.678,0.0402458116705
1.679,0.0267562792572
1.68,0.0133043689708
"""
In [166]:
%%R -i BD_vert
con = textConnection(BD_vert)
df = read.csv(con)
close(con)
df = df %>%
select(vert_tube_pos, BD)
# curve fitting
fit = lm(BD~vert_tube_pos, data=df)
fit %>% print
df$y_pred_poly1 = predict(fit)
fit = lm(BD~vert_tube_pos+I(vert_tube_pos^2), data=df)
fit %>% print
df$y_pred_poly2 = predict(fit)
fit = lm(BD~vert_tube_pos+I(vert_tube_pos^2)+I(vert_tube_pos^3), data=df)
fit %>% print
df$y_pred_poly3 = predict(fit)
fit = lm(BD~vert_tube_pos+I(vert_tube_pos^2)+I(vert_tube_pos^3)+I(vert_tube_pos^4), data=df)
fit %>% print
df$y_pred_poly4 = predict(fit)
# plotting
ggplot(df, aes(vert_tube_pos)) +
geom_point(aes(y=BD), size=1.5) +
geom_line(aes(y=y_pred_poly1), color='green')+ #, size=2, alpha=0.3) +
geom_line(aes(y=y_pred_poly2), color='blue')+ #, size=2, alpha=0.3) +
geom_line(aes(y=y_pred_poly3), color='red')+ #, size=2, alpha=0.3) +
geom_line(aes(y=y_pred_poly4), color='purple')+ #, size=2, alpha=0.3) +
geom_hline(yintercept=1.74, linetype='dashed', alpha=0.5) +
theme_bw() +
theme(
text = element_text(size=16)
)
In [207]:
%%R -w 600
inFile = '/home/nick/notebook/SIPSim/t/genome100/DBL_index.txt'
df = read.delim(inFile, sep='\t') %>%
gather(pos, vert_grad_BD, vert_gradient_BD_low, vert_gradient_BD_high)
# example
df.ex = data.frame('DBL_BD' = c(1.667, 1.733), 'vert_grad_BD' = c(1.7, 1.7))
# plot
ggplot(df, aes(DBL_BD, vert_grad_BD, color=pos, group=DBL_BD)) +
#geom_point() +
geom_line(color='black', size=1) +
geom_point(data=df.ex, color='red', size=4) +
geom_line(data=df.ex, aes(group=vert_grad_BD), color='red', linetype='dashed', size=1.2) +
geom_vline(xintercept=1.774, linetype='dashed', alpha=0.5, color='blue') + # theoretical max fragment BD
scale_y_reverse() +
labs(x='BD of DBL', y='BD of vertical gradient\n(during fractionation)') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [205]:
%%R
# how big of a GC spread?
BD2GC = function(x) (x - 1.66) / 0.098 * 100
BD2GC(1.667) %>% print
BD2GC(1.7) %>% print
BD2GC(1.733) %>% print
In [ ]: