-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmalaria_regresion.py
195 lines (144 loc) · 7.28 KB
/
malaria_regresion.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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
# Eda Ehler, molekulární informatika
# qsar analýza léků proti malárii
# (predikce účinnosti látky - regrese vypočítaných prediktorů na IC50 (aktivitu))
############################################################
# !!!důležité - spouštím přes anacondu (nainstalovan rdkit)#
# - mám udělané virt.env "molinf" #
# - spustit z příkazové řádky příkazem "activate molinf" #
# -> "python malaria_regresion.py" #
############################################################
# IMPORTY - pandas, numpy, složky rdkitu a scikit-learnu
#--------
import pandas as pd
import numpy as np
from rdkit.Chem import AllChem
from rdkit.Chem import Descriptors
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVR # vypočet je moc pomalý, nepoužívám
from sklearn.cross_validation import KFold, cross_val_score, train_test_split, cross_val_predict
from sklearn.linear_model import LinearRegression, BayesianRidge, SGDRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score, explained_variance_score
# PROMĚNNÉ - input/output
#------------------------
bioactivity = pd.read_csv("bioactivity.txt", sep="\t")
malaria = pd.read_csv("malaria.txt", sep="\t")
out_file = "malaria_output.txt"
# FUNKCE
#-------
def convert_to_nM(radek):
# pokud má molekula aktivitu v ug.mL-1, převede ji na nM
if radek["ACTIVITY_UNITS"] == "ug.mL-1":
act = int(radek["ACTIVITY"])
molwt = AllChem.CalcExactMolWt(radek["MOL_OBJECT"])
radek["ACTIVITY"] = (act/molwt)*1000000
radek["ACTIVITY_UNITS"] = "nM"
return radek
else:
return radek
# funkce na cross-validaci, upraveno dle Garreta, Moncecchi (2013) Learning scikit-learn. Chapter 2, p.29
def evaluate_cross_validation(clf, X, y, K):
# create a k-fold cross-validation iterator
cv = KFold(len(y), K, shuffle=False, random_state=0)
# by default the score used is the one returned by score method of the estimator (accuracy)
scores = []
for train_index, test_index in cv:
X_train, X_test = X[train_index], X[test_index]
y_train, y_test = y[train_index], y[test_index]
scores.append(train_and_evaluate(clf, X_train, X_test, y_train, y_test))
#print(scores)
mean_abs_err = sum([x[0] for x in scores])/K
mean_sqr_err = sum([x[1] for x in scores])/K
avg_r2 = sum([x[2] for x in scores if (x[2] >= -1 and x[2] <= 1)])/K # pouze od -1 do 1
explained_var = sum([x[3] for x in scores if (x[3] >= -1 and x[3] <= 1)])/K # pouze od -1 do 1
print("Mean abs error:", mean_abs_err)
print("Mean sqr error:", mean_sqr_err)
print("Coeficient of determination (r2):", avg_r2)
print("Explained variance:", explained_var)
return (mean_abs_err, mean_sqr_err, avg_r2, explained_var)
# funkce na trénování a evaluaci, upraveno dle Garreta, Moncecchi (2013) Learning scikit-learn. Chapter 2, p.29
def train_and_evaluate(clf, X_train, X_test, y_train, y_test):
clf.fit(X_train, y_train)
#print("Accuracy on training set:")
#print(clf.score(X_train, y_train))
#print("Accuracy on testing set:")
#print(clf.score(X_test, y_test))
y_predicted = clf.predict(X_test)
vysledek = (mean_absolute_error(y_test, y_predicted),
mean_squared_error(y_test, y_predicted),
r2_score(y_test, y_predicted),
explained_variance_score(y_test, y_predicted))
#print("mean_absolute_error:")
#print(vysledek[0])
#print("mean_squared_error:")
#print(vysledek[1])
#print("r2_score:")
#print(vysledek[2])
#print("explained_variance_score:")
#print(vysledek[3])
return vysledek
# VÝPOČTY
#--------
# vyberu jen relevantní sloupce (jméno, smiles na vypočet fingerprintů, aktivitu, jednotky)
bioactivity = bioactivity[["CMPD_CHEMBLID", "CANONICAL_SMILES", "STANDARD_VALUE", "STANDARD_UNITS"]]
malaria = malaria[["CHEMBLID", "CANONICAL_SMILES", "A1_STANDARD_VALUE", "A1_STANDARD_UNITS"]]
# pro rychlé testy vezmu jen 30 prvních řádků##
#bioactivity = bioactivity[:30] #
#malaria = malaria[:30] #
###############################################
# přejmenovat sloupce stejně
jmena_sloupcu = ["CHEMBLID", "SMILES", "ACTIVITY", "ACTIVITY_UNITS"]
bioactivity.columns = jmena_sloupcu
malaria.columns = jmena_sloupcu
# spojit bioactivity a malaria do jedné tabulky, a vyhodit stejná jména
bio_mal = pd.concat([bioactivity, malaria])
bio_mal = bio_mal.drop_duplicates(subset="CHEMBLID")
# dropne všechny řádky, kde je nějaký "NaN" (DataFrame.dropna(axis=0, how='any', thresh=None, subset=None, inplace=False))
bio_mal = bio_mal.dropna()
#print(bio_mal.head)
# add RDKit mol objects
bio_mal["MOL_OBJECT"] = bio_mal.SMILES.apply(AllChem.MolFromSmiles)
jmena_sloupcu.append("MOL_OBJECT")
# přehodí všechny řádky na stejné jednotky: ug/ml -> g/dm3 -> nM
bio_mal = bio_mal.apply(convert_to_nM, axis=1)
# vypočítá deskriptory z MOL_OBJECT proměnné a přidá je do tabulky
for jmeno, funkce in Descriptors.descList: # slovnik nazev:funkce na vypocet deskriptoru z mol objektu
bio_mal[jmeno] = bio_mal.MOL_OBJECT.apply(funkce)
#print(bio_mal.head)
# dropne sloupce s malým standard deviation (uniformní nebo skoro uniformní deskriptory)
# http://stackoverflow.com/questions/31799187/drop-columns-with-low-standard-deviation-in-pandas-dataframe
threshold = 0.05
data_ready = bio_mal.drop(bio_mal.std()[bio_mal.std() < threshold].index.values, axis=1)
#print(data_ready.shape)
# pro další použití vezmu jen numpy matice prediktorů (variabilní deskriptory)
# a target qsar analýzy -> to je aktivita (de facto hodnota IC50)
prediktory = np.asarray(data_ready.drop(jmena_sloupcu, axis=1))
targety = np.asarray(data_ready.ACTIVITY)
# standardizace
prediktory = StandardScaler().fit_transform(prediktory)
"""
print(prediktory)
input("prediktory")
print(targety)
input("targety")
"""
# REGRESE
#--------
regressors_dict = { "BayesianRidge": BayesianRidge(),
"LinearRegression": LinearRegression(n_jobs=-1),
"SGDRegressor": SGDRegressor(loss='squared_loss', penalty=None, random_state=42),
#"SVR_linear": SVR(kernel='linear'), # linear function
#"SVR_poly": SVR(kernel='poly'), # polynomial funtion
#"SVR_rbf": SVR(kernel='rbf')} #radial basis function
}
# výpočet a tisk výsledků
with open(out_file, mode="w", encoding="utf-8") as output:
output.write("Regressor".ljust(20) + "\tmean_abs_err\tmean_sqr_err\tavg_r2\texplained_var\n")
for name, regressor in regressors_dict.items():
print(name)
scores = evaluate_cross_validation(regressor, prediktory, targety, 5)
output.write(name.ljust(20) + "\t")
output.write(str(round(scores[0], 3)) + "\t")
output.write(str(round(scores[1], 3)) + "\t")
output.write(str(round(scores[2], 3)) + "\t")
output.write(str(round(scores[3], 3)) + "\n")
print("---")