-
Notifications
You must be signed in to change notification settings - Fork 1
/
Fig4B_part2.py
156 lines (149 loc) · 7.66 KB
/
Fig4B_part2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 11 21:03:37 2022
@author: hromada
"""
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from usefulFunctions import COM1_colors, numID
import seaborn as sns
from matplotlib.colors import LinearSegmentedColormap
def gLV_Model1(x_vector, t=0, u_vector=None, a_matrix=None, antib_param_vector=None, antib_conc=None):
new_x_vector = np.zeros(len(x_vector))
for i in range(len(x_vector)):
new_x_vector[i] = x_vector[i] * (u_vector[i] + np.sum([a_matrix[i][j] * x_vector[j] for j in range(len(x_vector))])+(antib_param_vector[i] * antib_conc))
return new_x_vector
def simulate_gLV_Model1(x0_vector,Ts,u_vector, a_matrix, antib_param_vector, antib_conc):
return odeint(gLV_Model1, x0_vector, Ts, args=(u_vector, a_matrix, antib_param_vector, antib_conc))
def to_matrix(l, n):
return [l[i:i+n] for i in range(0, len(l), n)]
# load best fit parameters
speciesID={'CD':1,'ER':2,'DP':3,'FP':4,'BH':5,'CA':6,'PC':7,'EL':8,'CH':9,'BO':10,'BT':11,'BU':12,'BV':13,'CS':14}
gLV_params_df = pd.read_csv('gLV_parameters_MSB_2021.csv', header=None, index_col=None)
u_gLV_vector = gLV_params_df[0].tolist()[0::15]
aijs = gLV_params_df[0].tolist()
del aijs[0::15]
a_matrix = to_matrix(aijs, 14)
antib_param_df = pd.read_csv('Model1_scaled_parameter.csv', index_col=[0])
scaling_factor_dict = {'Metronidazole':24,'Vancomycin':12}
# load part2 data (communities part2)
E81_comms = {'SEHE81_1':[1],'SEHE81_2':[1,3,4,8],'SEHE81_3':[1,3,4,6,8],'SEHE81_4':[1,3,4,8,9],'SEHE81_5':[1,3,4,8,13],'SEHE81_6':[1,3,4,7,8,14]}
E81_df = pd.read_csv('Abundance_data_SuppFig3AB_part2.csv', index_col=[0]) #this is the same data as 'Abundance_data_SuppFig3AB' part 2 deposited at https://doi.org/10.5281/zenodo.7626486 , just formatted slightly differently
#duplicate 0 antibiotic condition to be Metr and Vanco
zero_df = E81_df[E81_df['Concentration']==0]
for ind in zero_df.index:
zero_df.at[ind,'Antibiotic'] = 'Metronidazole'
E81_df = pd.concat([E81_df, zero_df])
modelType_fraction_correct = {'Model1':[-1],'gLV_no_antib':[-1],'gLV_no_int':[-1],'Model1_num':[-1],'gLV_num':[-1],'gLV_no_int_num':[-1]}
#%% Simulate conditions from part2 data with expanded gLV model (Model1) and two null models (gLV_no_antib, gLV_no_int)
# output1: heatmap of prediction (accurate or inaccurate) for each concentration for each condition as .pdf file for Model1
# output2: save the fraction of correct prediction for Model1 and each null model in .csv file
tmpt=47
#modelType
for modelTypenum in [0,1,2]:
modelType = ['Model1','gLV_no_antib','gLV_no_int'][modelTypenum]
if modelType in ['Model1','gLV_no_antib','Model4']:
my_amatrix = a_matrix
if modelType in ['gLV_no_int']:
my_amatrix = [[] for i in range(len(a_matrix))]
for r in range(len(a_matrix)):
for v in range(len(a_matrix)):
if r==v:
my_amatrix[r].append(a_matrix[r][v]) #copy aiis
else:
my_amatrix[r].append(0) #remove all aijs
count_match = 0
count_unmatch = 0
count_nodata = 0
for antib in ['Metronidazole','Vancomycin']:
if modelType in ['Model1', 'gLV_no_int']:
b_list = []
for i in range(0,14):
sp = numID[i+1]
b_val = float(antib_param_df[(antib_param_df['Antibiotic']==antib)&(antib_param_df['Species']==sp)].iloc[0]['B'])
b_list.append(b_val)
if modelType in ['gLV_no_antib']:
b_list = [0 for i in range(0,14)]
conclist = {'Metronidazole':[0,0.375,0.75,1.5,3,6,12,24,48],'Vancomycin':[0,0.09375,0.1875,0.375,0.75,1.5,3,6]}[antib]
d_heatmap = []
m_heatmap = []
yticklabels = []
for s in range(6):
#Get general information
comm_name = 'SEHE81_'+str(s+1)
splist = E81_comms[comm_name]
yticklabels.append(comm_name)
initial_condition = [0 for i in range(14)]
#Data
heatmap_row = []
zerodf = E81_df[(E81_df['Condition']==comm_name)&(E81_df['Time']==tmpt)&(E81_df['Antibiotic']==antib)&(E81_df['Concentration']==0)]
CDODzero = np.mean(zerodf['sp1OD'])
for c in range(len(conclist)):
conc = conclist[c]
concdf = E81_df[(E81_df['Condition']==comm_name)&(E81_df['Time']==tmpt)&(E81_df['Antibiotic']==antib)&(E81_df['Concentration']==conc)]
if len(concdf)>0:
CDOD = np.mean(concdf['sp1OD'])
subMICfc = CDOD/CDODzero
if CDOD < 0.05:
heatmap_row.append(-1) #no growth
elif subMICfc > {'Metronidazole':1.21,'Vancomycin':1.05}[antib]:
heatmap_row.append(1)
else:
heatmap_row.append(0)
else:
heatmap_row.append(0.1)
d_heatmap.append(heatmap_row)
#Simulation
heatmap_row = []
initial_OD = 0.0022 #This is different than other experiments, but E81 and E90 kept initial OD constant btw species
for spnum in splist:
initial_condition[int(spnum)-1]=initial_OD
sim_zero_df = pd.DataFrame(simulate_gLV_Model1(initial_condition,[0,48],u_gLV_vector, my_amatrix,b_list,0))
CDODzero = sim_zero_df.loc[1][0]
for c in range(len(conclist)):
conc = conclist[c]
spnum = 1
scaling_factor = scaling_factor_dict[antib]
scaled_conc = conc/scaling_factor
#do simulation
my_sim_data = simulate_gLV_Model1(initial_condition,[0,48],u_gLV_vector, my_amatrix,b_list,scaled_conc)
sim_df = pd.DataFrame(my_sim_data)
CDOD = sim_df.loc[1][0]
subMICfc = CDOD/CDODzero
if CDOD < 0.05:
heatmap_row.append(-1) #no growth
elif subMICfc > {'Metronidazole':1.21,'Vancomycin':1.05}[antib]:
heatmap_row.append(1)
else:
heatmap_row.append(0)
m_heatmap.append(heatmap_row)
df1 = pd.DataFrame(m_heatmap)
df2 = pd.DataFrame(d_heatmap)
df3 = df1==df2
for r in range(len(df2)):
for v in range(len(df2.iloc[r])):
if df2.iloc[r][v] == 0.1:
count_nodata += 1
else:
if df3.iloc[r][v] == True:
count_match += 1
else:
count_unmatch += 1
if modelType == 'Model1':
fig3,ax = plt.subplots()
ng = np.ma.masked_where(df2.values == 0.1, df3)
myColors = ['gray','dodgerblue']
cmap = LinearSegmentedColormap.from_list('Custom', myColors, len(myColors))
sns.heatmap(df3, cmap=cmap, alpha=0.1, cbar=False, square=True,xticklabels = conclist, yticklabels = yticklabels)
ax.pcolor(np.arange(len(df3.columns)+1),np.arange(len(df3.index)+1), ng,cmap=cmap,linewidth=1,edgecolor='black', alpha=0.9)
plt.title(antib)
plt.show()
plt.close()
total = count_unmatch+count_match
frac = count_match/total
modelType_fraction_correct[modelType] = [frac]
modelType_fraction_correct[modelType+'_num'] = [count_match]
modelType_df = pd.DataFrame.from_dict(modelType_fraction_correct)
modelType_df.to_csv('part2_model_prediction_accuracy.csv')