Carbon Balance Calculations and Gas Transfer Rates

import os                                # navigating the operation system
import numpy as np                       # doing basic math
import pandas as pd                      # handling tables
from datetime import date                # getting the date for file saving
from scipy import stats                  # doing statistical analysis
from scipy.optimize import curve_fit     # data fitting
import matplotlib.pyplot as plt          # visualization
from iambcodes.rates import *            # improved visualization and fitting

Root_Dir = os.getcwd()
Target_Dir = Root_Dir  # Adjust if needed
File = 'Example_GrowthExp.xlsx' # Add Excel file with growth data
# All files are attached with the current date
Today = date.today().strftime('%y%m%d')
# Define image properties
ImageType = 'svg'      
fontsize = 15
plt.rc('font', size=fontsize)
kwargs = {'mew':5, 'lw':4}

# Substrate mol to g conversion parameter
mol2g = {'glc':1/.18, 'glyc':1/.092, 'ace':1/.059, 'eg':1/.0621}
# plus minus sign
pm = '\u00B1'

FileAddress = os.path.join(Target_Dir, File)
Xls_df = pd.ExcelFile(FileAddress)
print('Excel sheets:\n{}'.format(Xls_df.sheet_names))

# Choosing Excel sheet
Sheet = '10mM' 
# Define conversion factor, e.g. yeast: 0.62 (BNID 111182 in Bionumbers)

# Loading data from sheet into python, ignoring Metadata columns
df_orig = pd.read_excel(FileAddress, sheet_name=Sheet, usecols=lambda x: 'Metadata' not in x,)
df_orig = df_orig.dropna(how = 'all')
DelHeader = df_orig.columns[['Unnamed' in header for header in df_orig.columns]].values
df_orig = df_orig.drop(DelHeader, axis=1)

myMeta = pd.read_excel(FileAddress, sheet_name=Sheet, usecols=lambda x: 'Metadata' in x,)
TimeUnit = ''.join(myMeta.loc[myMeta['Metadata: Key']=='Time unit']['Metadata: Value'].values)
# add a biomass correction factor, e.g. from OD to gCDW, if none set to 1
BiomassUnit = ''.join(myMeta.loc[myMeta['Metadata: Key']=='Biomass unit']['Metadata: Value'].values)
if BiomassUnit == 'OD':
    ODConvert = True
    # read or define existing conversion factor for OD to gCDW/L
    BiomassCorr = myMeta.loc[myMeta['Metadata: Key']=='OD2CDW']['Metadata: Value'].values
    OD2DW = BiomassCorr if BiomassCorr.size > 0 else OD2DW
    BiomassUnit = 'gDW/L'
elif BiomassUnit == 'gDW/L':
    OD2DW = 1.
else:
    print('Use standard abbreviation for biomass unit: gDW/L, OD')
SubstrateUnit = ''.join(myMeta.loc[myMeta['Metadata: Key']=='Substrate unit']['Metadata: Value'].values)# 'mM'
# add a substrate correction factor, e.g. to convert from g/l to mM
SubstrateID = ''.join(myMeta.loc[myMeta['Metadata: Key']=='Substrate']['Metadata: Value'].values).lower()
# conversion of g/l concentration to mmol/l
SubstrateCorr = mol2g[SubstrateID] if SubstrateUnit=='g/L' else 1.
YieldCorr = 1/mol2g[SubstrateID]# if SubstrateUnit=='mM' else 1.

# Extract column ids from metadata
SubstrateColumns = myMeta['Metadata: Value'][myMeta['Metadata: Key']=='Substrate'].astype(str).values[0]
BiomassColumns = myMeta['Metadata: Value'][myMeta['Metadata: Key']=='Biomass Column ID'].astype(str).values[0]
TimeID = myMeta['Metadata: Value'][myMeta['Metadata: Key']=='Time Column ID'].astype(str).values[0]
BiomassID = [col for col in df_orig.columns if BiomassColumns in col]
SubstrateID = [col for col in df_orig.columns if SubstrateColumns in col]

df=df_orig.dropna(how = 'all').copy()

df['Biomass'] = np.mean(df[BiomassID].to_numpy(), axis=1)*OD2DW
df['Substrate'] = np.mean(df[SubstrateID].to_numpy(), axis=1)*SubstrateCorr

# Calculating substrate uptake rate
DeltaBm = np.array(df['Biomass']-df['Biomass'][0], dtype=float)
Sub_sl, Sub_ic, Sub_r, _, Sub_sl_sterr = stats.linregress(DeltaBm, df['Substrate'])
SubUp = round(-Sub_sl*Bm_sl,1)
SubUp_sterr = round(SubUp*np.sqrt( (Sub_sl_sterr/Sub_sl)**2 + (Bm_sl_sterr/Bm_sl)**2 ), 2)

# r2 calculation
r2_Sub, _ = FitR2(df['Substrate'], DeltaBm, [Sub_sl, Sub_ic], law='lin', inlog=False)

# yield calculation
Yield = round(Bm_sl,2) / (SubUp * YieldCorr) 
# calculating the standard error of the yield based on gaussian error propagation
Yield_sterr = round(Yield*np.sqrt( (SubUp_sterr/SubUp)**2 + (Bm_sl_sterr/Bm_sl)**2 ), 2)

print(f'Growth rate: {round(Bm_sl,2)}({pm}{round(Bm_sl_sterr,2)}) /{TimeUnit}, with R2: {round(r2_BM,2)}')
print(f'Substrate uptake: {SubUp}({pm}{SubUp_sterr}) mmol/gCDW/h (check units), with R2: {round(abs(Sub_r),2)}')
print(f'Biomass yield: {round(Yield,2)}{pm}{Yield_sterr} gDW/g')

# Substrate plots
LinReg = lambda x, m, b: m*x + b
plt.plot(DeltaBm, LinReg(DeltaBm, Sub_sl, Sub_ic), **kwargs)
plt.plot(DeltaBm, df['Substrate'],'go', **kwargs)
plt.legend([f'Slope={Sub_sl:.2f}mmol/gDW,\nIntercept={Sub_ic:.2f}mM'])
plt.xlabel('$\\Delta$ Biomass, gDW/L')
plt.ylabel('Substrate, mM')
plt.show()

plt.plot(df['time'], df['Substrate'], 'og', label='Substrate', **kwargs)
plt.plot(df['time'], Sub_ic+Sub_sl*DeltaBm, '-b', label=f'fit, rate={SubUp}{pm}{SubUp_sterr}mmol/gCDW/h', **kwargs)
plt.xlabel(f'Time, {TimeUnit}')
plt.ylabel('Substrate, mM')
plt.show()

prod = myMeta.loc[myMeta['Metadata: Key']=='Product', 'Metadata: Value'].squeeze()
Products = [p.strip() for p in prod.split(',')] if isinstance(prod, str) and prod.strip().lower() not in ('', 'nan') else []
if Products:
    Prod_df = pd.DataFrame()
    Prod_dict = {p: {} for p in Products}   # initialize keys

    for myProd in Products:
        myCols = [col for col in df_orig.columns if myProd in col]
        if len(myCols) == 0:
            print(f'Warning: no columns found for product {myProd}')
            continue

        myDat = df[myCols].to_numpy()
        Prod_df[f'{myProd}_mean'] = np.mean(myDat, axis=1)

        # linear fit of product vs DeltaBm
        Prod_sl, Prod_ic, Prod_r, _, Prod_sl_sterr = stats.linregress(DeltaBm, Prod_df[f'{myProd}_mean'])

        # compute rate and stderr safely
        rate = round(Prod_sl*Bm_sl, 1)
        if Prod_sl == 0 or Bm_sl == 0:
            rate_sterr = np.nan
        else:
            rate_sterr = round(rate * np.sqrt((Prod_sl_sterr / Prod_sl) ** 2 + (Bm_sl_sterr / Bm_sl) ** 2), 2)

        Prod_dict[myProd]['rate'] = rate
        Prod_dict[myProd]['rate_sterr'] = rate_sterr
        Prod_dict[myProd]['r2'] = round(abs(Prod_r), 2)

        print(f"Product {myProd} formation rate: {Prod_dict[myProd]['rate']}({pm}{Prod_dict[myProd]['rate_sterr']}) mmol/gCDW/h, R2: {Prod_dict[myProd]['r2']}")

    # create a consistent color for each product using a colormap
    cmap = plt.get_cmap('tab10', len(Products))
    colors = {p: cmap(i) for i, p in enumerate(Products)}

    for myProd in Products:
        col = colors[myProd]
        plt.plot(DeltaBm, Prod_df[f'{myProd}_mean'], 'o', color=col, label='_nolegend_', **kwargs)
        plt.plot(DeltaBm, Prod_dict[myProd]['rate']/Bm_sl * DeltaBm, '-', color=col, **kwargs)
    plt.xlabel('$\\Delta$ Biomass, gDW/L')
    plt.ylabel('Product concentration, mM')
    plt.legend(Products, loc='best')
    plt.show()

    plt.plot(df['time'], df['Substrate'], 'ok', label='_nolegend_', **kwargs)
    plt.plot(df['time'], Sub_ic+Sub_sl*DeltaBm, '-k', label=f'Glc:{SubUp}mmol/gCDW/h', **kwargs)
    for myProd in Products:
        col = colors[myProd]
        plt.plot(df['time'], Prod_df[f'{myProd}_mean'], 'o', color=col, label='_nolegend_', **kwargs)
        plt.plot(df['time'], Prod_dict[myProd]['rate']/Bm_sl * DeltaBm, '-', color=col, label=f'{myProd}:{Prod_dict[myProd]["rate"]}mmol/gDW/h', **kwargs)
    plt.xlabel('Time, h')
    plt.ylabel('Concentration, mM')
    plt.legend(loc='best')

    plt.tight_layout()
    plt.show()