small fix data (SHK0073) and other analysis files

Author: PHvanLent

figures/optimization_results/sample_size_exp0.png

@@ -0,0 +1,179 @@
+"""Runs specific pathway model combinatorial experiment analysis (linear modelling)"""
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from source.visualization import plot_dbtl_cycles
+from matplotlib.colors import ListedColormap
+from sklearn.linear_model import LinearRegression
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.linear_model import LinearRegression
+from sklearn.pipeline import make_pipeline
+import statsmodels.formula.api as smf
+import pandas as pd
+import statsmodels.formula.api as smf
+
+name = "pathwayA"
+target = 'Y_product_A'
+base_dir = f"results/optimization_results/stratified_sampling/{name}/experiment1"
+output_dir = f"results/statistics_combinatorial_exp/stratified_sampling/{name}/"
+dirs = os.listdir(base_dir)
+
+
+#pathway A settings
+Number_of_Samples = [50, 100, 200]
+Screening_Ratio = [1, 2, 4, 6]
+Positions = [4, 6, 8, 10]
+Features = [6, 12, 19]
+
+# # # pathway B settings
+# Number_of_Samples = [50, 100,200]
+# Screening_Ratio = [2, 4, 6]
+# Positions = [4, 6, 8, 10]
+# Features = [6, 10,12]
+#
+# # #pathway C settings
+# Number_of_Samples = [50, 100,200]
+# Screening_Ratio = [2, 4, 6]
+# Positions = [4, 6, 8, 10]
+# Features = [6, 10, 12, 17]
+# #
+# # #pathway D settings
+# Number_of_Samples = [50, 100,200,300]
+# Screening_Ratio = [2, 4, 6]
+# Positions = [4, 6, 8, 10]
+# Features = [5, 10, 15, 21]
+# #
+# # pathway E settings
+# Number_of_Samples = [50, 100, 200, 300]
+# Screening_Ratio = [1, 2, 4, 6]
+# Positions = [4, 6, 8, 10]
+# Features = [5, 10, 15, 20, 25]
+
+
+option_combinations = []
+list_of_list_combinations = []
+for F in Features:
+    for P in Positions:
+        for N in Number_of_Samples:
+            for R in Screening_Ratio:
+                grid_str = f"S{R*N}X4N{N}F{F}P{P}"
+                option_combinations.append(grid_str)
+                list_of_list_combinations.append([R*N, N, F, P])
+
+print("number of combinations", len(list_of_list_combinations))
+
+results = {}
+for k, option in enumerate(option_combinations):
+    option_values = []
+    for dir in dirs:
+        working_dir = os.path.join(base_dir, dir)
+        files = os.listdir(working_dir)
+        for file in files:
+
+            if option in file:
+                file_path = os.path.join(working_dir, file)
+                process = pd.read_csv(file_path, index_col=0)
+
+                max_val = process[process.index.str.startswith('cycle4')][target].median() # calculates the maximum after 5 rounds
+
+                #calculates an auc
+                test = process.copy()
+                test.index = test.index.astype(str)
+                test["cycle"] = test.index.str.extract(r"(cycle\d+)",expand=False)
+                cycle_auc = test.groupby("cycle")[target].mean().sum()
+                max_val = cycle_auc
+
+                option_values.append(max_val)
+    results[option]=option_values
+
+gridsearch_results = pd.DataFrame.from_dict(results,
+                                            orient='index').iloc[:,:10] #only takes 10 repeats (sometimes more experiments were done)
+
+
+print("number of missing experiments", sum(gridsearch_results.isna().sum(axis=1)))
+
+
+
+
+# First, make a DataFrame with column names
+mean_performance = gridsearch_results.mean(axis=1)
+scenarios = pd.DataFrame(np.array(list_of_list_combinations), columns = ['S','N','F','P'])
+scenarios['mean_performance'] = np.array(mean_performance)
+scenarios= scenarios.sort_values('mean_performance', ascending=False)
+scenarios['SNratio'] = scenarios['S']/scenarios['N']
+
+
+df = scenarios.copy()
+
+# Build formula with interactions
+# Example: 'target ~ A + B + C + A:B + A:C + B:C'
+features = ['SNratio','F','P']
+target = 'mean_performance'
+# formula = f"{target} ~ " + " + ".join(features) + " + " + " + ".join([f"{a}:{b}" for i, a in enumerate(features) for b in features[i+1:]])
+formula = 'mean_performance ~ SNratio + F + P'
+
+
+model = smf.ols(formula, data=df).fit()
+print(model.summary())
+with open(f"{output_dir}/{name}_lm_SNratio_interaction.txt", "w") as f:
+    f.write(model.summary().as_text())
+
+model.summary()
+from statsmodels.stats.multitest import multipletests
+
+# Assume you already have a fitted model
+pvals = model.pvalues
+
+# Apply Benjamini-Hochberg (FDR) or Bonferroni correction
+corrected = multipletests(pvals, alpha=0.05, method='fdr_bh')  # or method='bonferroni'
+
+# Output adjusted p-values
+
+adjusted_pvals = pd.DataFrame({
+    'coef': model.params,
+    'raw_pval': pvals,
+    'adj_pval': corrected[1],
+    'significant (FDR)': corrected[0]
+})
+
+adjusted_pvals.to_csv(f"{output_dir}/mht_{name}_lm_with_interaction.csv")
+
+
+import statsmodels.formula.api as smf
+
+# First, make a DataFrame with column names
+df = scenarios.copy()
+
+# Build formula with interactions
+# Example: 'target ~ A + B + C + A:B + A:C + B:C'
+features = df.columns[:-1]
+target = df.columns[-1]
+formula = 'mean_performance ~ S + N + F + P'
+
+
+model = smf.ols(formula, data=df).fit()
+print(model.summary())
+with open(f"{output_dir}/{name}_lm_SandN_interaction.txt", "w") as f:
+    f.write(model.summary().as_text())
+
+model.summary()
+from statsmodels.stats.multitest import multipletests
+
+# Assume you already have a fitted model
+pvals = model.pvalues
+
+# Apply Benjamini-Hochberg (FDR) or Bonferroni correction
+corrected = multipletests(pvals, alpha=0.05, method='fdr_bh')  # or method='bonferroni'
+
+# Output adjusted p-values
+import pandas as pd
+adjusted_pvals = pd.DataFrame({
+    'coef': model.params,
+    'raw_pval': pvals,
+    'adj_pval': corrected[1],
+    'significant (FDR)': corrected[0]
+})
+
+
+adjusted_pvals.to_csv(f"{output_dir}/mht_{name}_lm_without_interaction.csv")