MiBiPreT example: Isotope Analysis with Amersfoort data¶

Diagnostic plots for data analysis on microbial biodegredation at the contaminant Amersfoort site. Author: Alraune Zech

Data based on the PhD thesis of Johan van Leeuwen, 2021 'Biodegredation of mono- and polyaromatic hydrocarbons in a contaminated aquifer originating from a former Pintsch gas factory site' which is equivalent to the manuscript of van Leeuwen et al., 2022 'Anaerobic degradation of benzene and other aromatic hydrocarbons in a tar-derived plume: Nitrate versus iron reducing conditions', J. of Cont. Hydrol. The data was provided by Johan van Leeuwen.

Background: Amersfoort contaminant site¶

Close to the train station in Amersfoort, the Netherlands, the subsurface is contaminated with organic hydrocarbons forming a NAPL. Contamination originates from decades of operating manufactured gas plant, which dumped tar by-products in waste lagoons. The tar is a DNAPL and has spread into the underlying shallow unconfined aquifer. Sample wells were installed to measure various characteristics of the subsurface. The raw data contains measurements on

• environmental conditions, such as pH, redox potential, concentrations of oxygen, nitrate, etc
• contaminant concentration such as BTEX, indene, indane, naphtalene and multiple other (typically cyclic) petroleum hydrocarbons
• metabolite concentration, i.e. byproducts of degredation processes of contaminant
• isotope measurments for specific contaminants and samples
• counts of genes (RNA/DNA) of mibrobiota that is know to perform biodegredation as well as functional enzymes know to be responsible for biodegredation

In [1]:

import mibiscreen as mbs

import mibiscreen as mbs

In [2]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt import numpy as np import pandas as pd

Data loading¶

In [3]:

file_path = './amersfoort.xlsx'

file_path = './amersfoort.xlsx'

Load and standardize data of isotopes:

In [4]:

isotopes_raw,units = mbs.load_excel(file_path,
                                   sheet_name = 'isotopes',
                                   verbose = False)

isotopes,units = mbs.standardize(isotopes_raw,
                                reduce = True,
                                verbose=False)

isotopes_raw,units = mbs.load_excel(file_path, sheet_name = 'isotopes', verbose = False) isotopes,units = mbs.standardize(isotopes_raw, reduce = True, verbose=False)

Lambda regression and Lambda regression plot for single molecule¶

In [5]:

C_data,H_data = mbs.extract_isotope_data(isotopes,'toluene')
results = mbs.Lambda_regression(C_data,
                                H_data,
                                validate_indices = True,
                                verbose = True,
                                )
mbs.Lambda_plot(**results,
                title = 'toluene',
               fit_color = 'k',
               marker_color = 'C1')

C_data,H_data = mbs.extract_isotope_data(isotopes,'toluene') results = mbs.Lambda_regression(C_data, H_data, validate_indices = True, verbose = True, ) mbs.Lambda_plot(**results, title = 'toluene', fit_color = 'k', marker_color = 'C1')

==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 24.92
______________________________________________________________

Out[5]:

(<Figure size 375x280 with 1 Axes>,
 <Axes: title={'center': 'toluene'}, xlabel='$\\delta^{{13}}$C', ylabel='$\\delta^2$H'>)

Lambda regression and Lambda regression plot for all molecules (separately)¶

In [6]:

molecules = ['benzene','toluene','ethylbenzene','pm_xylene','naphthalene','indene']
molecules_analysis = []

for j,molecule in enumerate(molecules):

    C_data,H_data = mbs.extract_isotope_data(isotopes,molecule)

    results = mbs.Lambda_regression(C_data,
                                H_data,
                                validate_indices = True,
                                verbose = True,
                                )

    molecules_analysis.append(results)

    mbs.Lambda_plot(**results,title = molecule)
    # Lambda_plot(**results,save_fig = 'Amersfoort_isotope_Lambda_{}.png'.format(molecule),dpi = 300)

molecules = ['benzene','toluene','ethylbenzene','pm_xylene','naphthalene','indene'] molecules_analysis = [] for j,molecule in enumerate(molecules): C_data,H_data = mbs.extract_isotope_data(isotopes,molecule) results = mbs.Lambda_regression(C_data, H_data, validate_indices = True, verbose = True, ) molecules_analysis.append(results) mbs.Lambda_plot(**results,title = molecule) # Lambda_plot(**results,save_fig = 'Amersfoort_isotope_Lambda_{}.png'.format(molecule),dpi = 300)

==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 18.37
______________________________________________________________
==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 24.92
______________________________________________________________
==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 26.26
______________________________________________________________
==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 20.16
______________________________________________________________
==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 6.38
______________________________________________________________
==============================================================
 Running function 'Lambda_regression()' on data
==============================================================
Lambda value, being the slope of the linear fit is 
 identified with 1.66
______________________________________________________________

Lambda regression and plot of results of all molecules in one plot¶

Produce plot with matplotlib for indiviual adaptions

In [7]:

molecules = ['benzene','toluene','naphthalene','indene']
molecules_analysis_set = [molecules_analysis[i] for i in [0,1,4,5]]
fig, axes = plt.subplots(figsize=[5.5,5],
                          ncols = 2,
                          nrows = 2)
ax = axes.flat
for j,molecule in enumerate(molecules):

    results =  molecules_analysis_set[j]

    x,y = results['delta_C'], results['delta_H']

    # Plot the scatter plot
    ax[j].scatter(x,y,color = 'C{}'.format(j),ec = '0.5',s=50, zorder = 3,label = 'data')
    ax[j].set_title(molecule,bbox=dict(facecolor="oldlace",boxstyle="round"))
    ax[j].grid(True,zorder = 0)

                    
    # plot trendlines
    polynomial = np.poly1d(results['coefficients'])
    trendline_x = np.linspace(np.min(x), np.max(x), 50)
    trendline_y = polynomial(trendline_x)
    ax[j].plot(trendline_x, trendline_y,c= 'k', label=r'$\Lambda = {:.0f}$'.format(results['coefficients'][0]))

    if j%2 == 0:
        ax[j].set_ylabel(r'$\delta^2$H')

    if j >= len(ax)-2:
        ax[j].set_xlabel(r'$\delta^{{13}}$C')

    ax[j].legend()

fig.tight_layout()
plt.savefig("Amersfoort_isotope_Lambda.png",dpi=300)

molecules = ['benzene','toluene','naphthalene','indene'] molecules_analysis_set = [molecules_analysis[i] for i in [0,1,4,5]] fig, axes = plt.subplots(figsize=[5.5,5], ncols = 2, nrows = 2) ax = axes.flat for j,molecule in enumerate(molecules): results = molecules_analysis_set[j] x,y = results['delta_C'], results['delta_H'] # Plot the scatter plot ax[j].scatter(x,y,color = 'C{}'.format(j),ec = '0.5',s=50, zorder = 3,label = 'data') ax[j].set_title(molecule,bbox=dict(facecolor="oldlace",boxstyle="round")) ax[j].grid(True,zorder = 0) # plot trendlines polynomial = np.poly1d(results['coefficients']) trendline_x = np.linspace(np.min(x), np.max(x), 50) trendline_y = polynomial(trendline_x) ax[j].plot(trendline_x, trendline_y,c= 'k', label=r'$\Lambda = {:.0f}$'.format(results['coefficients'][0])) if j%2 == 0: ax[j].set_ylabel(r'$\delta^2$H') if j >= len(ax)-2: ax[j].set_xlabel(r'$\delta^{{13}}$C') ax[j].legend() fig.tight_layout() plt.savefig("Amersfoort_isotope_Lambda.png",dpi=300)