Notes

Panel A is a diagram.

To see how we got from the raw electrochemical scans to the datasets used here, please see the following processing notebooks:

Then see how those data were analyzed in these notebooks for supplemental figures S6 and S7:

These supplemental figure notebooks produced model coefficients that are used in this notebook.

Setup packages and plotting for the notebook:

Fig. 5B PYO soak - transfer

Fig. 5C PCA soak - transfer

Fig. 5D - SWV signal

Let’s read in the SWV data and plot it:

Fig. 5E - GC signal

Let’s read in the GC data and plot it:

Fig. 5F - GC vs. SWV signal

Now we will plot the SWV and GC peak currents against each other

Fig. 5G - Blank vs. biofilm PYO retention and \(D_{loss}\)

To plot the SWV signals over time with the \(D_{loss}\) model fits, we will read in the data, model coefficients, and predictions for the blank and ∆phz biofilm datasets. Let’s read in those 6 files:

# Model coefficients
df_dphz_nls_1 <- read_csv("../supplement/Fig_S7/phz2019_dPHZ_Dphys_nls_coefs.csv") %>% 
  filter(exp == 1 & run == 1 & term == 'a')

df_blank_nls_75 <- read_csv("../supplement/Fig_S7/phz2019_blank_Dphys_nls_coefs.csv") %>% 
  filter(PHZadded == '75uM' & term == 'a')

# The signals are the peak current data points
# we are subtracting the constant background signal from model fit
df_dphz_swv_decay <- read_csv("../../processing/processed_data/phz_eDNA_2019_signals_long.csv") %>% 
  filter(echem == 'SWV' & exp == 1 & run == 1 & reactor == 'transfer' & electrode == 'i1') %>% 
  mutate(IDA = 'biofilm', signal = signal - df_dphz_nls_1$estimate) %>% 
  select(time, signal, IDA)

df_blank_decay <- read_csv("../../processing/processed_data/phz_eDNA_2019_swv_blank_tran_time_signals.csv") %>% 
  filter(PHZadded == '75uM') %>% 
  mutate(IDA = 'blank', signal = signal - df_blank_nls_75$estimate) %>%
  select(time, signal, IDA)

# Predictions are best fit lines generated from the model
# we are subtracting the constant background signal from model fit
df_dphz_preds <- read_csv("../supplement/Fig_S7/phz2019_dPHZ_Dphys_preds.csv") %>% 
  filter(exp == 1 & run == 1) %>% 
  select(time, pred, pred_low, pred_high) %>% 
  mutate(pred = pred - df_dphz_nls_1$estimate,
         pred_low = pred_low - df_dphz_nls_1$estimate, 
         pred_high = pred_high - df_dphz_nls_1$estimate)%>% 
  mutate(IDA = 'biofilm')

df_blank_preds <- read_csv("../supplement/Fig_S7/phz2019_blank_Dphys_preds.csv") %>% 
  filter(PHZadded == '75uM') %>% 
  select(time, pred, pred_low, pred_high) %>% 
  mutate(pred = pred - df_blank_nls_75$estimate,
         pred_low = pred_low - df_blank_nls_75$estimate, 
         pred_high = pred_high - df_blank_nls_75$estimate) %>% 
  mutate(IDA = 'blank')

df_preds <- bind_rows(df_dphz_preds, df_blank_preds)

df_decays <- bind_rows(df_dphz_swv_decay, df_blank_decay)

We also are subtracting out the constant background signal as determined by the model intercept, so that both curves are easier to visually compare.

Now let’s make the final plot:

Fig. 5H - \(D_{ap}\) vs. \(D_{loss}\)

In order to construct figure 5F we must do a few analysis steps using data analysed for Figs S6 and S7. This is somewhat complicated, because to calculate \(D_{loss}\) we need the \(D_{ap}\) value and the SWV soak value (I0) for each run. Therefore we will do the following:

  1. Take biofilm linear models (SWV vs. GC) from fig. S6 and calculate \(D_{ap}\).
  2. Take blank linear model (SWV vs. GC) from Fig. S7 and calculate \(D_{ap}\).
  3. Take nonlinear models from fig. S6 and join with \(D_{ap}\) values and initial current I0, which is taken to be the peak current from the soak reactors before each run. Then calculate \(D_{loss}\) from nonlinear fit coefficient, \(D_{ap}\) and I0.

\(D_{ap}\)

First let’s read in the linear model coefficients from Figs S6 and S7. Recall that this linear model was used to fit the SWV vs. GC data plots and that the slope is proportional to \(\sqrt{D_{ap}}\).

exp_id run_id term estimate std.error statistic p.value conf.low conf.high IDA
Biofilm 1 Rep 1 signal_SWV 0.1513328 0.0039875 37.951641 0.0000000 0.1427183 0.1599473 biofilm
Biofilm 1 Rep 2 signal_SWV 0.1045085 0.0034865 29.975017 0.0000000 0.0969763 0.1120407 biofilm
Biofilm 1 Rep 3 signal_SWV 0.1203080 0.0079031 15.222956 0.0000000 0.1032345 0.1373816 biofilm
Biofilm 2 Rep 1 signal_SWV 0.2791415 0.0036863 75.723946 0.0000000 0.2711777 0.2871052 biofilm
Biofilm 2 Rep 2 signal_SWV 0.2125194 0.0077260 27.507105 0.0000000 0.1958284 0.2292104 biofilm
Biofilm 2 Rep 3 signal_SWV 0.2230086 0.0077181 28.894223 0.0000000 0.2063346 0.2396825 biofilm
Biofilm 1 Rep 1 (Intercept) 0.0000000 0.0000000 -14.499185 0.0000000 0.0000000 0.0000000 biofilm
Biofilm 1 Rep 2 (Intercept) 0.0000000 0.0000000 -10.551620 0.0000001 0.0000000 0.0000000 biofilm
Biofilm 1 Rep 3 (Intercept) -0.0000001 0.0000000 -9.616970 0.0000003 -0.0000001 0.0000000 biofilm
Biofilm 2 Rep 1 (Intercept) 0.0000000 0.0000000 -18.031183 0.0000000 0.0000000 0.0000000 biofilm
Biofilm 2 Rep 2 (Intercept) 0.0000000 0.0000000 -10.724706 0.0000001 0.0000000 0.0000000 biofilm
Biofilm 2 Rep 3 (Intercept) 0.0000000 0.0000000 -9.337288 0.0000004 0.0000000 0.0000000 biofilm
Blank 1 Rep 1 (Intercept) -0.0000002 0.0000000 -7.205283 0.0055105 -0.0000003 -0.0000001 blank
Blank 1 Rep 1 signal_SWV 0.1987799 0.0033905 58.628598 0.0000109 0.1879898 0.2095699 blank


So you can see we have estimated the slopes for each dataset (term ‘signal_SWV’ ‘estimate’), including confidence intervals. Now let’s redefine the function we used in Fig. S7 to calculate \(D_{ap}\) from the slope of the SWV vs. GC line.

Recall that for this type of data:

\[D_{ap} = \frac{1}{\pi t_p} \left( \frac{m A \psi}{S} \right) ^2\]

Now we’ll convert slope into \(D_{ap}\) estimates (with confidence intervals) using our function dap_from_swvGC(), and we will go ahead and plot the biofilm and blank values:

\(D_{loss}\)

Now we will calculate \(D_{loss}\), so we first need the model coefficients from the SWV decay signals from Fig S6. We also need the peak SWV current from each scan in the soak reactor before it was transferred - this will serve as our value for initial current I0. Lastly, we also need the \(D_{ap}\) values calculated above. Below we will read in those data and put them all together in one data frame that we can work with.

# Biofilm SWV decay fits and soak I0's

df_dphz_nls <- read_csv("../supplement/Fig_S7/phz2019_dPHZ_Dphys_nls_coefs.csv")
df_dphz_i0 <- read_csv("../../processing/processed_data/phz_eDNA_2019_swv_signals.csv") %>% 
  filter(reactor == 'soak' & electrode == 'i1') %>% 
  mutate(run = as.double(run)) %>% 
  select(exp, run, i0 = signal)

# Join all biofilm data together

df_dphz_nls_i0 <- left_join(df_dphz_nls, df_dphz_i0, by = c('exp','run')) %>% mutate(IDA = 'biofilm')

# Blank SWV decay fits and soak I0's

df_blank_nls <- read_csv("../supplement/Fig_S7/phz2019_blank_Dphys_nls_coefs.csv")
df_blank_i0 <- read_csv("../../processing/processed_data/phz_eDNA_2019_swv_blank_signals.csv") %>% 
  filter(reactor == 'soak') %>% 
  select(PHZadded, i0 = signal)

# Join all blank data together

df_blank_nls_i0 <- left_join(df_blank_nls, df_blank_i0, by = c('PHZadded')) %>% 
  mutate(exp = 3, IDA = 'blank') %>% 
# Add an ID variable for matching below
  mutate(run = case_when(
    PHZadded == '10uM' ~ 4,
    PHZadded == '25uM' ~ 4,
    PHZadded == '50uM' ~ 4,
    PHZadded == '75uM' ~ 4,
    PHZadded == '100uM' ~ 4
  ))

# Join the biofilm and blank data

df_all_nls_i0 <- bind_rows(df_dphz_nls_i0, df_blank_nls_i0) %>% filter(term == 'b') %>% 
  select(exp, run,b_estimate = estimate, b_low = conf.low, b_high = conf.high, i0, IDA, PHZadded)

# Add the Dap values for the biofilm and blank data

df_all_dap_i0 <- left_join(df_all_nls_i0, df_all_dap %>% select(exp, run, IDA, dap, dap_high, dap_low),
                           by = c('exp','run','IDA'))

df_all_dap_i0 %>% kable() %>% kable_styling() %>% scroll_box(height = '300px')
exp run b_estimate b_low b_high i0 IDA PHZadded dap dap_high dap_low
1 1 2.8e-06 2.8e-06 2.9e-06 3.40e-06 biofilm NA 4.00e-06 4.40e-06 3.50e-06
1 2 1.8e-06 1.8e-06 1.9e-06 3.50e-06 biofilm NA 1.90e-06 2.20e-06 1.60e-06
1 3 2.1e-06 2.0e-06 2.1e-06 3.80e-06 biofilm NA 2.50e-06 3.30e-06 1.80e-06
2 1 2.1e-06 2.1e-06 2.2e-06 6.60e-06 biofilm NA 1.35e-05 1.42e-05 1.27e-05
2 2 4.3e-06 4.1e-06 4.5e-06 6.60e-06 biofilm NA 7.80e-06 9.10e-06 6.60e-06
2 3 2.3e-06 2.2e-06 2.3e-06 5.00e-06 biofilm NA 8.60e-06 9.90e-06 7.40e-06
3 4 3.0e-07 2.0e-07 3.0e-07 3.90e-06 blank 25uM 6.80e-06 7.60e-06 6.10e-06
3 4 7.0e-07 7.0e-07 7.0e-07 6.70e-06 blank 50uM 6.80e-06 7.60e-06 6.10e-06
3 4 7.0e-07 7.0e-07 7.0e-07 9.30e-06 blank 75uM 6.80e-06 7.60e-06 6.10e-06
3 4 1.2e-06 1.1e-06 1.3e-06 1.21e-05 blank 100uM 6.80e-06 7.60e-06 6.10e-06


Recall that each of the decays with the expression:

\[y = b (x)^{-0.5} + a\]

So the variable b is the model coefficient we will use here.

As defined in the supplementary text, we can define \(D_{loss}\) as follows:

\[D_{loss} = \frac{I_0^2 D_{ap} t_s}{\pi b^2}\]

Now that we have all the required values, we can define a function to calculate \(D_{loss}\):

We will calculate \(D_{loss}\) using a reasonable guess for the scan time parameter \(t_s\), 6ms. The scan time is not a known parameter, so there is uncertainty in this guess.

Now we can compare \(D_{loss}\) and \(D_{ap}\) graphically:

Let’s go ahead and calculate the mean values for \(D_{loss}\) and \(D_{ap}\) in the biofilm and the blank:

## # A tibble: 2 x 3
##   IDA       mean_dap   mean_dloss
##   <chr>        <dbl>        <dbl>
## 1 biofilm 0.00000637 0.0000000680
## 2 blank   0.00000683 0.00000195

Overall, this seems reasonable, but you can see that the blank \(D_{loss}\) is lower than the blank \(D_{ap}\). We aimed to obtain a conservative estimate of biofilm \(D_{loss}\), so given the uncertainty in the scan time, we decided to choose a value that made the blank \(D_{loss}\) equal to the blank \(D_{ap}\).

Redoing the above analysis with a value of \(t_s = 21ms\) yields the following results:

And comparing \(D_{loss}\) and \(D_{ap}\) now yields the plot for figure 5F:

We can calculate the mean \(D_{ap}\) and \(D_{loss}\) values:

## # A tibble: 2 x 3
##   IDA       mean_dap  mean_dloss
##   <chr>        <dbl>       <dbl>
## 1 biofilm 0.00000637 0.000000238
## 2 blank   0.00000683 0.00000681

And we can see that there is a statistically significant difference between the biofilm \(D_{ap}\) and \(D_{loss}\) whether we pool all the technical replicates:

## # A tibble: 2 x 4
##   IDA     conf_int_low conf_int_high p_value
##   <chr>          <dbl>         <dbl>   <dbl>
## 1 biofilm   0.00000248           Inf 0.00969
## 2 blank    -0.00000338           Inf 0.495

Or by comparing \(D_{ap}\) and \(D_{loss}\) within each biological replicate:

## # A tibble: 3 x 5
## # Groups:   exp [3]
##     exp IDA     conf_int_low conf_int_high p_value
##   <dbl> <chr>          <dbl>         <dbl>   <dbl>
## 1     1 biofilm  0.000000943           Inf  0.0234
## 2     2 biofilm  0.00000443            Inf  0.0155
## 3     3 blank   -0.00000338            Inf  0.495

Create figure

Old Fig. 5:

Revised Fig. 5:

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.15.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.7.4   hms_0.5.3         modelr_0.1.5     
##  [4] broom_0.5.2       kableExtra_1.1.0  cowplot_0.9.4    
##  [7] viridis_0.5.1     viridisLite_0.3.0 knitr_1.23       
## [10] forcats_0.4.0     stringr_1.4.0     dplyr_0.8.3      
## [13] purrr_0.3.3       readr_1.3.1       tidyr_1.0.0      
## [16] tibble_2.1.3      ggplot2_3.3.0     tidyverse_1.3.0  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2       lattice_0.20-38  assertthat_0.2.1 digest_0.6.21   
##  [5] utf8_1.1.4       R6_2.4.0         cellranger_1.1.0 backports_1.1.4 
##  [9] reprex_0.3.0     evaluate_0.14    httr_1.4.1       highr_0.8       
## [13] pillar_1.4.2     rlang_0.4.6      readxl_1.3.1     rstudioapi_0.10 
## [17] Matrix_1.2-15    rmarkdown_1.13   labeling_0.3     splines_3.5.3   
## [21] webshot_0.5.1    munsell_0.5.0    compiler_3.5.3   xfun_0.7        
## [25] pkgconfig_2.0.3  mgcv_1.8-27      htmltools_0.4.0  tidyselect_0.2.5
## [29] gridExtra_2.3    fansi_0.4.0      crayon_1.3.4     dbplyr_1.4.2    
## [33] withr_2.1.2      grid_3.5.3       nlme_3.1-137     jsonlite_1.6    
## [37] gtable_0.3.0     lifecycle_0.1.0  DBI_1.0.0        magrittr_1.5    
## [41] scales_1.0.0     cli_1.1.0        stringi_1.4.3    fs_1.3.1        
## [45] xml2_1.2.2       ellipsis_0.3.0   generics_0.0.2   vctrs_0.3.1     
## [49] tools_3.5.3      glue_1.3.1       yaml_2.2.0       colorspace_1.4-1
## [53] rvest_0.3.5      haven_2.2.0