Goal:

  • Basic assessment on whether the diffusive boundary layer of the ultracentrifugation tubes may be causing the broad OTU abundance distributions (most span the entire gradient)

Literative on diffusive boundary layers (DBL)

DNA binding to plastics

Diffusion of DNA

  • Literature resourse

  • sedimentation coefficient of DNA (S)

    • $S = 2.8 + (0.00834 * M^{0.479})$
      • where
        • M = molecular weight of DNA
    • OR $S = 2.8 + (0.00834 * (L*666)^{0.479})$
      • where
        • L = length of DNA
  • Svedberg's equation
    • $s/D = \frac{M(1-\bar{V}p)}{RT}$
    • where
      • s = sedimentation coefficient
      • D = diffusion coefficient
      • M = molecular weight
      • $\bar{V} = 1/\rho_p$
        • $\rho_p$ = density of the sphere
      • p = density of the liquid
      • R = universal gas constant
      • T = absolute temperature
  • Finding diffusion coefficient of DNA in CsCl
    • $D = \frac{RT}{M(1-\bar{V}p)}*s$
      • where
        • R = 8.3144598 (J mol^-1 K^-1)
        • T = 293.15 (K)
        • p = 1.7 (Buckley lab gradients)
        • $\bar{V} = 1/\rho_p$
          • $\rho_p$ = 1.99
        • $s = 2.8 + (0.00834 * (L*666)^{0.479})$
          • L = DNA length (bp)
  • Einstein-Smoluchowski relation
    • $t = \frac{z^2}{0.9 * D}$
      • where
        • t = time (sec)
        • z = mean deviation of molecules from starting position
        • D = diffusion coefficient (cm^2 s^-1)

Init


In [1]:
%load_ext rpy2.ipython

In [135]:
%%R
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(rootSolve)
library(gridExtra)

Setting parameters


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


Tube angle from axis of rotation: 27.95319 

In [83]:
%%R
x = r_max__mm - r_min__mm 
hyp = tube_height__mm
asin(x / hyp)


[1] 0.4878751

Isoconcentration point

  • Formula 6.7 in Birnine and Rickwood 1978

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


Isoconcentration point: 3.781204 (cm)

Fraction size range (based on previous run data)


In [8]:
%%R
# loading SIPdb info
df.SIPdb = read_excel(SIPdbFile, sheet=sheet.id)
df.SIPdb %>% head(n=3) %>% as.data.frame


      sample_id gradient_date fractionation_date fraction_id well_id     RI
1 13C-Ami.D3.R1    2015-01-27         2015-01-30           1      A1 1.4075
2 13C-Ami.D3.R1    2015-01-27         2015-01-30           2      B1 1.4071
3 13C-Ami.D3.R1    2015-01-27         2015-01-30           3      C1 1.4067
  RI_corrected       BD volume experiment
1       1.4059 1.770113     NA    fullCyc
2       1.4055 1.765742     NA    fullCyc
3       1.4051 1.761371     NA    fullCyc

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


    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
-0.00219  0.00328  0.00437  0.00440  0.00546  0.01748      120 

In [10]:
%%R
# fraction size range (1st & 3rd quartiles)
frac_size_range__BD = c(0.003, 0.005)

band width distribution

  • tube height of fraction = volume of fraction

$100 ul = 0.1 cm^3$

$v = pi*r^2*h$

$h = v / (pi * r^2)$


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


The height of each fraction in vertically oriented tube: 0.07533962 (cm)
(Just applies to cylinder)

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


Tube bottom volume: 0.5751733 (cm^3)
Number of fractions in tube bottom: 5.751733 

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


Source: local data frame [6 x 3]
Groups: sample_id [3]

       sample_id fractionation_date BD_heaviest
           (chr)             (time)       (dbl)
1 12C-Con.D14.R1         2015-02-23    1.769020
2 12C-Con.D14.R1         2015-07-17    1.771206
3  12C-Con.D1.R2         2015-02-16    1.769020
4  12C-Con.D1.R2         2015-06-26    1.773391
5 12C-Con.D30.R1         2015-03-09    1.772298
6 12C-Con.D30.R1         2015-05-04    1.766835

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


   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.07534 0.37670 0.37670 0.39240 0.45200 0.52740 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   5.000   5.000   5.208   6.000   7.000 

Figure description:

  • rug plot is showing the heaviest gradient fraction actually used further in SIP pipeline

Notes:

  • similar to Flamm et al., 1972 on fixed angle rotors
  • near the bottom of the tube, the density changes most rapidly
    • smaller bandwidths at bottom of the tube

Converting verticle bandwidth to angled bandwidth

  • Flamm et al., 1966

Width of band in angled rotor (at speed):

$w = \frac{V}{\pi*r*b}$

Width of band in vertical position:

$w' = \frac{V}{\pi*r^2}$

Where:

  • V = common volume
  • r = tube radius
  • b = band ellipse major axis (at speed)

surface to volume ratio of each fraction

Rough approximation


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='')


Percent of fraction that is DBL: 3.053254%

Stability of gradient

$t = 0.3(r_b - r_t)^2$

  • r unit: cm
  • t unit: hours

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


Time span of stability: 1.51875 (h)

Tube position of fraction (during spin; at equilibrium)

  • angled tube: this sets up the diffusive boundary layer
  • Determining tube location of each density band:
    • For density range of interest:
      • determine location range in tube
  • Determining where each density band touched on the tube during ultracentrifugation:
    • gradient is basically a set of ellipsoids (angled cross sections) of the tube (cylinder)
      • Need to determine the ellipsoid main axis diameter for each position in the gradient
    • For each location range in tube (see above):
      • Determine the ellipsoid major axis length
      • Determine the tube side wall locations (bottom & top) for the ellipsoid
        • This is the location of the diffusive boundary layer for the ellispoid

Tube location of each density band

Density at a given radius from axis of rotation:

$ x = \big(\frac{w^2}{2*\beta^o} * (r_d^2 - I^2)\big) + D$

  • where:
    • $x$ = density at a given radius
    • $w^2$ = angular velocity
    • $\beta^o$ = beta coef
    • $r_d$ = radius of desired density (cm)
    • $I$ = isocencentration point (cm)
    • $D$ = average density of gradient

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


isoconcentration point: 3.781204 
radius range for BD-min to BD-max:  3.289207 to 4.895445 

Geometry of a angled cross sections in the ultracentrifugation tube


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


ellipsoid major axis length in tube (cylinder) cross section: 2.664616 (cm)

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


radius position where round bottom comes into play: 3.55 

Notes:

  • Need to take into account round bottom because it falls into the radius range for BD min-max

Tube position of band (angled tube top & bottom) as a function distance from axis of rotation

Thanks Sam!


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


[1] 1.767843 3.434764
[1] 1.66692

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


[1] 1.767843 3.434764
[1] 1.66692

tube_band_position ~ BD

  • For tube in rotor with gradient at equilibrium, where on the tube height does each gradient band touch

    • Both bottom of tube and top
           /  /
    top-->/  /
         /| /
        / |/
       (  /<--bottom
        ^^

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


[1] 0.65
    tube_low__cm tube_high__cm      BD particle_pos
97     0.2029038     0.5099015 1.76923     4.826341
98     0.2487061     0.4514961 1.77023     4.839784
99           NaN           NaN 1.77123     4.853190
100          NaN           NaN 1.77223     4.866559
101          NaN           NaN 1.77323     4.879892
102          NaN           NaN 1.77423     4.893188

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


Misc


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


The height of each fraction in vertically oriented tube: 0.07533962 (cm)
(Just applies to cylinder)

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


Tube bottom volume: 0.5751733 (cm^3)
Number of fractions in tube bottom: 5.751733 

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


Source: local data frame [6 x 3]
Groups: sample_id [3]

       sample_id fractionation_date BD_heaviest
           (chr)             (time)       (dbl)
1 12C-Con.D14.R1         2015-02-23    1.769020
2 12C-Con.D14.R1         2015-07-17    1.771206
3  12C-Con.D1.R2         2015-02-16    1.769020
4  12C-Con.D1.R2         2015-06-26    1.773391
5 12C-Con.D30.R1         2015-03-09    1.772298
6 12C-Con.D30.R1         2015-05-04    1.766835

Tube position of fraction in vertical position

  • For each density band:
    • Determine the height in the tube
    • This will rely on empirical data to determine the relationship between BD and tube height (vertical position)

Relationship between density and tube height

  • volume ~ BD
  • tube height ~ volume

tube height ~ volume

Cylinder volume:

$v = \pi r^2$ h

Sphere volume:

$v = \frac{4}{3} \pi r^3$

Spherical cap volume:

$V = \frac{\pi * h^2}{3} (3r-h)$


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)


[1] 0.2360648

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)


[1] 0.07533962

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


[1] 0.2360648
[1] 0.5932149
[1] 0.65

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)


 [1] 0.0000000 0.2360648 0.3449723 0.4348016 0.5161206 0.5932149 0.6687044
 [8] 0.7440440 0.8193836 0.8947232 0.9700628 1.0454025 1.1207421 1.1960817
[15] 1.2714213 1.3467609 1.4221006 1.4974402 1.5727798 1.6481194 1.7234590
[22] 1.7987986 1.8741383 1.9494779 2.0248175 2.1001571 2.1754967 2.2508364
[29] 2.3261760 2.4015156 2.4768552 2.5521948 2.6275344 2.7028741 2.7782137
[36] 2.8535533 2.9288929 3.0042325 3.0795722 3.1549118 3.2302514

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


fraction volume ~ BD

  • determining from empirical data

In [35]:
%%R
# loading SIPdb info
df.SIPdb = read_excel(SIPdbFile, sheet=sheet.id)
df.SIPdb %>% head(n=3) %>% as.data.frame


      sample_id gradient_date fractionation_date fraction_id well_id     RI
1 13C-Ami.D3.R1    2015-01-27         2015-01-30           1      A1 1.4075
2 13C-Ami.D3.R1    2015-01-27         2015-01-30           2      B1 1.4071
3 13C-Ami.D3.R1    2015-01-27         2015-01-30           3      C1 1.4067
  RI_corrected       BD volume experiment
1       1.4059 1.770113     NA    fullCyc
2       1.4055 1.765742     NA    fullCyc
3       1.4051 1.761371     NA    fullCyc

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


Tube bottom volume: 0.5751733 (cm^3)
Number of fractions in tube bottom: 5.751733 

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


/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.

  res = super(Function, self).__call__(*new_args, **new_kwargs)

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


Call:
lm(formula = y ~ x + I(x^2) + I(x^3), data = dataset)

Coefficients:
(Intercept)            x       I(x^2)       I(x^3)  
    -2761.6       4943.8      -2931.8        576.5  


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


[1] 3.338287
[1] 2.0457
[1] 0.8140625
[1] 0.0528
WARNING: BD value 1.9 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0
[1] 0

tube height ~ BD


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


WARNING: BD value 1.811638 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0
WARNING: BD value 1.810545 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0
WARNING: BD value 1.809452 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0
WARNING: BD value 1.811638 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0
WARNING: BD value 1.810545 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0
WARNING: BD value 1.813823 is > 1.809, where tube volume = 0 Setting BD to 1.809
WARNING: tube volume < 0 for BD:  1.809 
 Setting volume to 0

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


Converting tube height to vertcal tube BD

  • height ~ volume (up to that height)
  • BD ~ volume (theory?)

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)


[1] 0.8627599

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)


[1] 0.5751733

In [ ]:
%%R

tubeVolume2BD = function(){
    # function to convert tube volume to BD
    
}

DBL overlap in fractions

  • max delta BD between fraction in rotor orientation vs when the orientation set to verticle for fractionation

In [44]:
%%R
## df of tube position for gradient in rotor (assuming equilibrium)
angle.tube.pos %>% head(n=4)


  tube_low__cm tube_high__cm      BD particle_pos
1     1.368140      3.817911 1.67323     3.289207
2     1.326125      3.775896 1.67423     3.308901
3     1.284360      3.734131 1.67523     3.328479
4     1.242838      3.692609 1.67623     3.347942

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


Notes

  • For BD bands in vertical tube (blue dots), the DBL is determined by the angled tube bands that hit that same tube position (green and red points).
  • For instance, the heaviest vertical tube band (could be thought of as a fraction with BD range of 0.001) has a DBL consisting of DNA from angled tube gradient bands with BD of ~1.69 & ~1.76.
    • Thus, both very light and heavy contaminating DNA should be diffusing into the fraction containing DNA of the corrected measured BD (measured via refractive index during fractionation).

Conclusions

  • For most of the fractions, the DBL should be contributing both very light and very heavy DNA.
    • The difference in BD of contaminating vs 'native' DNA is greatest for heaviest fractions (getting very light contaminating DNA) vs the lightest fractions (getting only 'medium' DNA).

Modeling DBL diffusion

  • Assumptions/params
    • Certain amount of DNA absorbed into tube wall
    • Thickness of DBL
    • Diffusion coef of DNA in CsCl
  • Modeling
    • Diffusion from a tube surface

In [107]:
%%R

(0.0362)**2 / (0.9 * 2e-5)

#F = function(z, D) z**2/(0.9*D)
#curve(F, )


[1] 72.80222

DNA diffusion

  • sedimentation coefficient of DNA (S)
    • $S = 2.8 + (0.00834 * M^{0.479})$
      • where
        • M = molecular weight of DNA
    • OR $S = 2.8 + (0.00834 * (L*666)^{0.479})$
      • where
        • L = length of DNA
  • Svedberg's equation
    • $s/D = \frac{M(1-\bar{V}p)}{RT}$
    • where
      • s = sedimentation coefficient
      • D = diffusion coefficient
      • M = molecular weight
      • $\bar{V} = 1/\rho_p$
        • $\rho_p$ = density of the sphere
      • p = density of the liquid
      • R = universal gas constant
      • T = absolute temperature
  • Finding diffusion coefficient of DNA in CsCl ($\mu m^2 / s$)
    • $D = \frac{RT}{M(1-\bar{V}p)}*s$
      • where
        • R = 8.3144598 (J mol^-1 K^-1)
        • T = 293.15 (K)
        • p = 1.7 (Buckley lab gradients)
        • $\bar{V} = 1/\rho_p$
          • $\rho_p$ = 1.99
        • $s = 2.8 + (0.00834 * (L*666)^{0.479})$
          • L = DNA length (bp)

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


Notes

  • Smaller DNA fragment length = height diffusion coefficient, with exponential-like increase at small lengths
    • very similar function to DNA_length ~ time_to_equilibrium

Calculating diffusion from DBL

  • Einstein-Smoluchowski relation
    • $t = \frac{z^2}{0.9 * D}$
      • where
        • t = time (sec)
        • z = mean deviation of molecules from starting position
        • D = diffusion coefficient (cm^2 s^-1)
    • rewritten: $z = \sqrt{0.9*D*t}$

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


Notes

  • according to the above figures, the essentially all of the DNA fragments should be diffused into the fractions within a couple of minutes at the most

In [ ]:



Sandbox


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


[1] 1.079949
[1] 1.28833

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)


[1] 1.079949 1.288330

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)


          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
[1,]  9.610939 8.544161 7.477382 6.410604 5.343825 4.277047 3.210268
[2,] 10.835989 9.769211 8.702432 7.635654 6.568875 5.502097 4.435318

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]:
0.0

Converting tube volume to vertical tube height (python)


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]:
0.07533961803166643

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]:
array([ 0.        ,  0.23603985,  0.34495023,  0.43482548,  0.51613246])

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]:
array([ 0.        ,  0.23603985,  0.34495023,  0.43482548,  0.51613246,
        0.59233273,  0.66767235,  0.74301197,  0.81835158,  0.8936912 ,
        0.96903082,  1.04437044,  1.11971006,  1.19504967,  1.27038929,
        1.34572891,  1.42106853,  1.49640815,  1.57174776,  1.64708738,
        1.722427  ,  1.79776662,  1.87310624,  1.94844585,  2.02378547,
        2.09912509,  2.17446471,  2.24980433,  2.32514394,  2.40048356,
        2.47582318,  2.5511628 ,  2.62650242,  2.70184203,  2.77718165,
        2.85252127,  2.92786089,  3.00320051,  3.07854012,  3.15387974,
        3.22921936])

finding continuous function for BD ~ vertical_tube_height


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


Call:
lm(formula = BD ~ vert_tube_pos, data = df)

Coefficients:
  (Intercept)  vert_tube_pos  
       1.6520        -0.0217  


Call:
lm(formula = BD ~ vert_tube_pos + I(vert_tube_pos^2), data = df)

Coefficients:
       (Intercept)       vert_tube_pos  I(vert_tube_pos^2)  
          1.669765           -0.039001            0.002332  


Call:
lm(formula = BD ~ vert_tube_pos + I(vert_tube_pos^2) + I(vert_tube_pos^3), 
    data = df)

Coefficients:
       (Intercept)       vert_tube_pos  I(vert_tube_pos^2)  I(vert_tube_pos^3)  
         1.6740702          -0.0468151           0.0048450          -0.0002039  


Call:
lm(formula = BD ~ vert_tube_pos + I(vert_tube_pos^2) + I(vert_tube_pos^3) + 
    I(vert_tube_pos^4), data = df)

Coefficients:
       (Intercept)       vert_tube_pos  I(vert_tube_pos^2)  I(vert_tube_pos^3)  
         1.677e+00          -5.616e-02           1.009e-02          -1.174e-03  
I(vert_tube_pos^4)  
         5.606e-05  

Plotting DBL index


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


Notes

  • The dashed line shows the BD values of fragments in the DBL for the vertical tube gradient BD band of 1.7

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


[1] 7.142857
[1] 40.81633
[1] 74.4898

In [ ]: