%load_ext rpy2.ipython

%%R
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(rootSolve)
library(gridExtra)

"""
Explanation: 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

blog on adsorbtion
5-25% DNA adsorbed to polypropylene
Belotserkovkii
whitepaper
Gaillard and Strauss
http://www.sciencedirect.com/science/article/pii/S0003269703000514
http://ac.els-cdn.com/S0003269785712717/1-s2.0-S0003269785712717-main.pdf?_tid=5920ba7a-9acf-11e5-94e8-00000aacb362&acdnat=1449265285_3576aa92de9eb1df4e7f299ef05389ff
Beckman ultra-cfg tubes:
http://www.centrifugebybeckman.com/?page_id=788/?pi_ad_id=62362322836&gclid=CLHP0c3gxckCFYQfHwodfp8CBg
surface properties to 'prevent adherence of nucleic acids and proteins'

Diffusion of DNA

Literature resourse
themofisher cfg guide

PNaS paper on DNA diffusion

$D \sim L^{-v}$
$v_L$ = 0.571 +/- 0.014
linear DNA


L in microns
D = uM^2/sec



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
End of explanation
"""


%%R
SIPdbFile = '/home/nick/db/SIPdb/DOE_SIP_db.xlsx'
sheet.id = 'Fractions'

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

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

%%R
x = r_max__mm - r_min__mm 
hyp = tube_height__mm
asin(x / hyp)

"""
Explanation: Setting parameters
End of explanation
"""


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

"""
Explanation: Isoconcentration point

Formula 6.7 in Birnine and Rickwood 1978
End of explanation
"""


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

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

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

"""
Explanation: Fraction size range (based on previous run data)
End of explanation
"""


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

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

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

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

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

"""
Explanation: band width distribution

tube height of fraction = volume of fraction 

$100 ul = 0.1 cm^3$
$v = pir^2h$
$h = v / (pi * r^2)$
End of explanation
"""


%%R
# in cm
frac_size = 0.01
tube_diam = 1.3
DBL_size = 0.01

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

"""
Explanation: 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}{\pirb}$
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
End of explanation
"""


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

"""
Explanation: Stability of gradient
$t = 0.3(r_b - r_t)^2$

r unit: cm
t unit: hours
End of explanation
"""


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

"""
Explanation: 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} $
End of explanation
"""


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

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

"""
Explanation: Geometry of a angled cross sections in the ultracentrifugation tube
End of explanation
"""


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

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

"""
Explanation: 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!
End of explanation
"""


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

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

"""
Explanation: 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
      ^^
~~~
End of explanation
"""


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

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

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

"""
Explanation: Misc
End of explanation
"""


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

%%R
# converting cylinder volume to cylinder height (cm)
cyl_vol2height = function(v, r){v / (pi * r**2)}

cyl_vol2height(0.1, 0.65)

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

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

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

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

"""
Explanation: 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)$
End of explanation
"""


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

%%R
# BD range
df.SIPdb = df.SIPdb %>%
    group_by(sample_id, fractionation_date) %>%
    arrange(fraction_id) %>%
    mutate(BD_range = BD - lead(BD)) 

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

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

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

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

"""
Explanation: fraction volume ~ BD

determining from empirical data
End of explanation
"""


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

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

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

"""
Explanation: tube height ~ BD
End of explanation
"""


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

%%R
height2cylinderVol = function(h, r){
    # h = cylinder height
    # r = cylinder radius
    pi * r**2 * h
}

# test
r = 0.65
height2cylinderVol(r, r)

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

%%R

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

"""
Explanation: Converting tube height to vertcal tube BD

height ~ volume (up to that height)
BD ~ volume (theory?)
End of explanation
"""


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

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

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

"""
Explanation: DBL overlap in fractions

max delta BD between fraction in rotor orientation vs when the orientation set to verticle for fractionation
End of explanation
"""


%%R

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

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

"""
Explanation: 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
End of explanation
"""


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

"""
Explanation: 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)
End of explanation
"""


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

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

"""
Explanation: 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.9Dt}$
End of explanation
"""


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


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

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

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)

"""
Explanation: 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



Sandbox
End of explanation
"""


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

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

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

"""
Explanation: Converting tube volume to vertical tube height (python)
End of explanation
"""


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

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

"""
Explanation: finding continuous function for  BD ~ vertical_tube_height
End of explanation
"""


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

"""
Explanation: Plotting DBL index
End of explanation
"""


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

"""
Explanation: Notes

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


