Bioassays mit Python auswerten:

Wir beschreiben auf dieser Seite eine Anwendung mit Python für die Auswertung von Bioassay-Daten. Wir gehen davon aus, dass unsere Daten anhand einer 4-parametrigen logistischen Regression modelliert werden und via der gewichteten Kleinste-Quadrate Schätzmethode ausgewertet werden. Um die gewichtete Kleinste-Quadrate Schätzmethode anzuwenden, muss man im Voraus die Gewichte spezifizieren. Diese werden in einem a-priori Schritt anhand von historischen Daten bestimmt. Die Gewichte stellen die Inverse der Varianz dar. Dementsprechend modellieren wir in dem a-priori Schritt die Varianz als eine Funktion der Dosis. Die Modellierung der Varianz aus historischen Daten wird als eine eigenständige Python-Applikation entwickelt. Wir präsentieren eine kurze Dokumentation der Applikation und einführende Python-Befehle:

  1. Collection and import of all files representing the historical assay data.
  2. For each assay:
    • Computation of the log10 of the dose.
    • Computation of the empirical standard deviation for each treatment and dose.
    • Computation of the natural logarithm of the standard deviation +1 for each treatment and dose e.g. log(std+1).
  3. Assemblage of all assays to one dataset.
  4. Computation of the simple linear regression of the log(std+1) vs log10 dose.
  5. Computation of the quadratic regression of the log(std+1) vs log10 dose (resp. the quadrat of log10 dose).
  6. Comparison of the achieved R² between linear regression and quadratic regression
  7. Outputs and graphical representations from the fitted model are delivered.
  8. Final decision table containing the model with the highest R² and the estimated parameters of the linear regression or of the quadratic regression is printed to the output.

Python Code

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
import scipy.stats
from sklearn import linear_model
import statsmodels.formula.api as smf

#read all files from local folder
files = os.listdir('dir')
derived={}

for i in range(len(files)):

    #read sas dataset
    filepath='dir'+files[i]
    sasdata=pd.read_sas(filepath)
    sasdata=sasdata.rename(columns={'F1': 'ST', 'F2': 'Dose','F3':'Response1','F4':'Response2'})

    #data preparations

    #derived dataset
    derived[i]=sasdata

final_df=pd.concat([derived[i] for i in range(6)], ignore_index=True)

#Linear Regression
# short via scipy
scipy.stats.linregress(np.log(final_df['Dose']),np.log(final_df['std_response']+1))

# via sklearn
x=final_df['Dose'].values
x = x.reshape(144, 1)
x = np.log(x)
y=final_df['std_response'].values
y=y.reshape(144,1)
y = np.log(y+1)
regr = linear_model.LinearRegression()
regr.fit(x, y)
plt.scatter(x, y, color='black')
plt.plot(x, regr.predict(x), color='blue', linewidth=3)
plt.xticks(())
plt.yticks(())
plt.show()

#via statsmodels
final_df=final_df.assign(log_std = np.log(final_df['std_response']+1),log_dose = np.log(final_df['Dose']))
sm_model=smf.ols('log_std~log_dose',data=final_df).fit()
print(sm_model.summary())

#coeffcient determination for linear regression
lr_r2=regr.score(x,y)