<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Binding Site Prediction
In this notebook we perform various machine learning methods and compare various aspects of machine learning paradigms
Step1:  Data Import and pre-processing
Step2:  Data Pre-processing
Step3:  Analysis I
Step4:  Some PSSM shows deviation from expected normal distribution -- which should be the case in neutral information setting due to Central Limit Theorem (CLT).
Analysis II
Step5:  With normalization
Step6:  Normalization over each feature space reduces the complexity of the problem which in turn improves the result.
SVM on PSSM + other putative useful features
We now test the performance on SVM using non-zero knowledge approach. We scourged through the complete list of features acquired by PredictProtein and include features that might have certain influence on DNA/RNA binding
Step7:  Analysis III
Step8:  Decision Tree
Decision Trees (DTs) are a non-parametric supervised learning method used for classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features.
Reference
L. Breiman, J. Friedman, R. Olshen, and C. Stone. Classification and Regression Trees. Wadsworth, Belmont, CA, 1984.
Step9:  Random Forest
In random forests, each tree in the ensemble is built from a sample drawn with replacement (i.e., a bootstrap sample) from the training set. In addition, when splitting a node during the construction of the tree, the split that is chosen is no longer the best split among all features. Instead, the split that is picked is the best split among a random subset of the features. As a result of this randomness, the bias of the forest usually slightly increases (with respect to the bias of a single non-random tree) but, due to averaging, its variance also decreases, usually more than compensating for the increase in bias, hence yielding an overall better model.
[SKL]
Reference
Step10:  Extremely Randomized Tree
In extremely randomized trees (see ExtraTreesClassifier and ExtraTreesRegressor classes), randomness goes one step further in the way splits are computed. As in random forests, a random subset of candidate features is used, but instead of looking for the most discriminative thresholds, thresholds are drawn at random for each candidate feature and the best of these randomly-generated thresholds is picked as the splitting rule. This usually allows to reduce the variance of the model a bit more, at the expense of a slightly greater increase in bias
Step11:  Random Forest on selected features
Step12:  While there is a significant accuracy improvement going from Decision Tree to Random Forest, the resulting prediction from Extremely Random Forest only improves the accuracy by the margin. Likewise, manually handpicking the features does not seem to improve the performance of the accuracy.
IV
<ASSISTANT_TASK:>
Python Code:
## matrix and vector tools
import pandas as pd
from pandas import DataFrame as df
from pandas import Series
import numpy as np
## sklearn
from sklearn.datasets import make_friedman1
from sklearn.feature_selection import RFE
from sklearn.svm import SVR
from sklearn.svm import SVC
from sklearn.model_selection import cross_val_score
from sklearn.datasets import make_blobs
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier
from sklearn.feature_selection import VarianceThreshold
# matplotlib et al.
from matplotlib import pyplot as plt
%matplotlib inline
Explanation: Binding Site Prediction
In this notebook we perform various machine learning methods and compare various aspects of machine learning paradigms:
Zero-knowledge vs. domain-knowledge based prediction
Single algorithms vs. ensemble methods
Prediction over normalized vs. non-normalized space
We also reviewed several machine learning algorithms such as Support Vector Machine including its variants (c-SVN, regressive SVN etc), tree based methods (decision tree, random forest, extremely random forest etc) and other ensembel methods (AdaBoost)
End of explanation
dna = df.from_csv('../../data/training_data_binding_site_prediction/dna_big.csv')
## embed class
dna = dna.reset_index(drop=False)
dna['class_bool'] = dna['class'] == '+'
dna['class_num'] = dna.class_bool.apply(lambda x: 1 if x else 0)
## added protein ID and corresponding position
dna['ID'] = dna.ID_pos.apply(lambda x: ''.join(x.split('_')[:-1]))
dna['pos'] = dna.ID_pos.apply(lambda x: x.split('_')[-1])
## data columns
dna.columns
## print available features
for feature in dna.columns[:-6]:
    
    print feature
dna
Explanation: Data Import and pre-processing
End of explanation
## create column-wise normalized data-set
dna_norm = dna.copy()
for col in dna_norm[dna_norm.columns[1:][:-6]].columns:
    dna_norm[col] = (dna_norm[col] - dna_norm[col].mean()) / (dna_norm[col].std() + .00001)
dna_norm
Explanation: Data Pre-processing: normalization
We apply Student's t normalization on each column:
$$X' = \frac{X - \hat{X}}{s}$$
This set would be used in parallel to normal dataset for comparison
End of explanation
# extract dataset and prediction
X = dna[dna.columns[1:][:-6]]
X = X[[x for x in X.columns.tolist() if 'pssm' in x]]
X = X.iloc[range(1000)]
y = dna['class_bool']
y = y[range(1000)]
# apply RFE on linear c-SVM
estimator = SVC(kernel="linear")
selector = RFE(estimator, 5, step=1)
selector = selector.fit(X, y)
print selector.ranking_
# redid previous routine on the whole data
pssm_rank = pd.DataFrame()
cat = dna['class']
for i in range(dna.index.size / 1000):
    
    this_cat = cat[range(i * 1000, (i + 1) * 1000)]
    if this_cat.unique().size > 1:
    
        X = dna[dna.columns[1:][:-6]]
        X = X[[c for c in X.columns.tolist() if 'pssm' in c]]
        X = X.iloc[range(i * 1000, (i + 1) * 1000)]
        y = dna['class_bool']
        y = y[range(i * 1000, (i + 1) * 1000)]
        estimator = SVC(kernel="linear")
        selector = RFE(estimator, 5, step=1)
        selector = selector.fit(X, y)
        print selector.ranking_
    
        pssm_rank[str(i)] = selector.ranking_
    
pssm_rank.index = [c for c in X.columns.tolist() if 'pssm' in c]
##sort PSSM features by its predictive power
rank_av = [np.mean(pssm_rank.ix[i]) for i in pssm_rank.index]
arg_rank_av = np.argsort(rank_av)
pssm_rank_sorted = pssm_rank.ix[pssm_rank.index[arg_rank_av]]
pssm_rank_sorted['RANK_AV'] = np.sort(rank_av)
pssm_rank_sorted
# plot average rank of all HSSP values
plt.hist([np.mean(pssm_rank.ix[i]) for i in pssm_rank.index], bins=60, alpha=.5)
plt.title("Histogram of Average HSSP Features Rank (RFE on linear SVM)")
fig = plt.gcf()
fig.set_size_inches(10, 6)
Explanation: Analysis I: Measuring the variation of significance among PSSM features
We want to see whether certain evolution patterns have any influence on DNA binding mechanism. We apply Recursive Feature Elimination (RFE) to rank our all PSSM-based features according to their predictive power using linear SVM (that is, SVM with linear kernel function).
Reference:
Guyon, I., Weston, J., Barnhill, S., & Vapnik, V., “Gene selection for cancer classification using support vector machines”, Mach. Learn., 46(1-3), 389–422, 2002.
End of explanation
X = dna[dna.columns[1:][:-6]]
y = dna.class_num
## train c-SVM
clf_svm1 = SVC(kernel='rbf', C=0.7)
clf_svm1.fit(X[dna.fold == 0], y[dna.fold == 0])
## predict class
pred = clf_svm1.predict(dna[dna.fold == 1][dna.columns[1:][:-6]])
truth = dna[dna.fold == 1]['class_num']
tp = pred[(np.array(pred) == 1) & (np.array(truth) == 1)].size
tn = pred[(np.array(pred) == 0) & (np.array(truth) == 0)].size
fp = pred[(np.array(pred) == 1) & (np.array(truth) == 0)].size
fn = pred[(np.array(pred) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
print "Size of (-)- and (+)-sets:\n\t(+)\t %d\n\t(-)\t%d" % (truth[truth == 1].index.size, truth[truth == 0].index.size)
Explanation: Some PSSM shows deviation from expected normal distribution -- which should be the case in neutral information setting due to Central Limit Theorem (CLT).
Analysis II: Support Vector Machine
We use c-SVM on two models :
all features model
handpicked features model
Reference
“Support-vector networks”, C. Cortes, V. Vapnik - Machine Learning, 20, 273-297 (1995).
“Automatic Capacity Tuning of Very Large VC-dimension Classifiers”, I. Guyon, B. Boser, V. Vapnik - Advances in neural information processing 1993.
SVM on all features
Without Normalization
End of explanation
X_norm = dna_norm[dna_norm.columns[1:][:-6]]
y = dna_norm.class_num
## train c-SVM
clf_svm2 = SVC(kernel='rbf', C=0.7)
clf_svm2.fit(X_norm[dna_norm.fold == 0], y[dna_norm.fold == 0])
## predict class
pred2 = clf_svm2.predict(dna_norm[dna_norm.fold == 1][dna_norm.columns[1:][:-6]])
truth = dna_norm[dna_norm.fold == 1]['class_num']
tp = pred2[(np.array(pred2) == 1) & (np.array(truth) == 1)].size
tn = pred2[(np.array(pred2) == 0) & (np.array(truth) == 0)].size
fp = pred2[(np.array(pred2) == 1) & (np.array(truth) == 0)].size
fn = pred2[(np.array(pred2) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: With normalization
End of explanation
## hand-pick features
features = [x for x in dna.columns[1:][:-6] if 'pssm' in x] +\
            [x for x in dna.columns[1:][:-6] if 'glbl_aa_comp' in x] +\
            [x for x in dna.columns[1:][:-6] if 'glbl_sec' in x] +\
            [x for x in dna.columns[1:][:-6] if 'glbl_acc' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_mass' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_hyd' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_cbeta' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_charge' in x] +\
            [x for x in dna.columns[1:][:-6] if 'inf_PP' in x] +\
            [x for x in dna.columns[1:][:-6] if 'isis_bin' in x] +\
            [x for x in dna.columns[1:][:-6] if 'isis_raw' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profbval_raw' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_sec_raw' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_sec_bin' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_acc_bin' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_normalize' in x] +\
            [x for x in dna.columns[1:][:-6] if 'pfam_within_domain' in x] +\
            [x for x in dna.columns[1:][:-6] if 'pfam_dom_cons' in x]
X_norm = dna_norm[features]
y = dna_norm.class_num
## train c-SVM
clf_svm3 = SVC(kernel='rbf', C=0.7)
clf_svm3.fit(X_norm[dna_norm.fold == 0], y[dna_norm.fold == 0])
## predict class
pred3 = clf_svm3.predict(X_norm[dna_norm.fold == 1])
truth = dna_norm[dna_norm.fold == 1]['class_num']
tp = pred3[(np.array(pred3) == 1) & (np.array(truth) == 1)].size
tn = pred3[(np.array(pred3) == 0) & (np.array(truth) == 0)].size
fp = pred3[(np.array(pred3) == 1) & (np.array(truth) == 0)].size
fn = pred3[(np.array(pred3) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: Normalization over each feature space reduces the complexity of the problem which in turn improves the result.
SVM on PSSM + other putative useful features
We now test the performance on SVM using non-zero knowledge approach. We scourged through the complete list of features acquired by PredictProtein and include features that might have certain influence on DNA/RNA binding
End of explanation
X = dna[dna.columns[1:][:-6]]
y = dna.class_num
Explanation: Analysis III: Tree-based Classificators
We picked three tree-based algorithms: Decision Tree (DT), Random Forest (RT) and Extremely Random Forest (ERT). From left to right, the algorithm allows more complexity into the models by introducing more randomness and biased than the previous algorithm.
End of explanation
# compute cross validated accuracy of the model
clf_t1 = DecisionTreeClassifier(max_depth=None, min_samples_split=2,
                             random_state=0)
scores = cross_val_score(clf_t1, X, y, cv=5)
print scores
print scores.mean()                             
clf_t1.fit(X[dna.fold == 0], y[dna.fold == 0])
pred_t1 = clf_t1.predict(X[dna.fold == 1])
truth = dna[dna.fold == 1]['class_num']
tp = pred_t1[(np.array(pred_t1) == 1) & (np.array(truth) == 1)].size
tn = pred_t1[(np.array(pred_t1) == 0) & (np.array(truth) == 0)].size
fp = pred_t1[(np.array(pred_t1) == 1) & (np.array(truth) == 0)].size
fn = pred_t1[(np.array(pred_t1) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: Decision Tree
Decision Trees (DTs) are a non-parametric supervised learning method used for classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features.
Reference
L. Breiman, J. Friedman, R. Olshen, and C. Stone. Classification and Regression Trees. Wadsworth, Belmont, CA, 1984.
End of explanation
# compute cross validated accuracy of the model
clf_t2 = RandomForestClassifier(n_estimators=10, max_depth=None,
                                 min_samples_split=2, random_state=0)
scores = cross_val_score(clf_t2, X, y, cv=5)
print scores
print scores.mean()
clf_t2.fit(X[dna.fold == 0], y[dna.fold == 0])
pred_t2 = clf_t2.predict(X[dna.fold == 1])
truth = dna[dna.fold == 1]['class_num']
tp = pred_t2[(np.array(pred_t2) == 1) & (np.array(truth) == 1)].size
tn = pred_t2[(np.array(pred_t2) == 0) & (np.array(truth) == 0)].size
fp = pred_t2[(np.array(pred_t2) == 1) & (np.array(truth) == 0)].size
fn = pred_t2[(np.array(pred_t2) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: Random Forest
In random forests, each tree in the ensemble is built from a sample drawn with replacement (i.e., a bootstrap sample) from the training set. In addition, when splitting a node during the construction of the tree, the split that is chosen is no longer the best split among all features. Instead, the split that is picked is the best split among a random subset of the features. As a result of this randomness, the bias of the forest usually slightly increases (with respect to the bias of a single non-random tree) but, due to averaging, its variance also decreases, usually more than compensating for the increase in bias, hence yielding an overall better model.
[SKL]
Reference:
Breiman, “Random Forests”, Machine Learning, 45(1), 5-32, 2001.
Breiman, “Arcing Classifiers”, Annals of Statistics 1998.
End of explanation
# compute cross validated accuracy of the model
clf_t3 = ExtraTreesClassifier(n_estimators=10, max_depth=None,
                              min_samples_split=2, random_state=0)
scores = cross_val_score(clf_t3, X, y, cv=5)
print scores
print scores.mean()
clf_t3.fit(X[dna.fold == 0], y[dna.fold == 0])
pred_t3 = clf_t3.predict(X[dna.fold == 1])
truth = dna[dna.fold == 1]['class_num']
tp = pred_t3[(np.array(pred_t3) == 1) & (np.array(truth) == 1)].size
tn = pred_t3[(np.array(pred_t3) == 0) & (np.array(truth) == 0)].size
fp = pred_t3[(np.array(pred_t3) == 1) & (np.array(truth) == 0)].size
fn = pred_t3[(np.array(pred_t3) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: Extremely Randomized Tree
In extremely randomized trees (see ExtraTreesClassifier and ExtraTreesRegressor classes), randomness goes one step further in the way splits are computed. As in random forests, a random subset of candidate features is used, but instead of looking for the most discriminative thresholds, thresholds are drawn at random for each candidate feature and the best of these randomly-generated thresholds is picked as the splitting rule. This usually allows to reduce the variance of the model a bit more, at the expense of a slightly greater increase in bias:
[SKL]
Reference:
P. Geurts, D. Ernst., and L. Wehenkel, “Extremely randomized trees”, Machine Learning, 63(1), 3-42, 2006.
End of explanation
features = [x for x in dna.columns[1:][:-6] if 'pssm' in x] +\
            [x for x in dna.columns[1:][:-6] if 'glbl_aa_comp' in x] +\
            [x for x in dna.columns[1:][:-6] if 'glbl_sec' in x] +\
            [x for x in dna.columns[1:][:-6] if 'glbl_acc' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_mass' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_hyd' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_cbeta' in x] +\
            [x for x in dna.columns[1:][:-6] if 'chemprop_charge' in x] +\
            [x for x in dna.columns[1:][:-6] if 'inf_PP' in x] +\
            [x for x in dna.columns[1:][:-6] if 'isis_bin' in x] +\
            [x for x in dna.columns[1:][:-6] if 'isis_raw' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profbval_raw' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_sec_raw' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_sec_bin' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_acc_bin' in x] +\
            [x for x in dna.columns[1:][:-6] if 'profphd_normalize' in x] +\
            [x for x in dna.columns[1:][:-6] if 'pfam_within_domain' in x] +\
            [x for x in dna.columns[1:][:-6] if 'pfam_dom_cons' in x]
X = dna[features]
y = dna.class_num
# compute cross validated accuracy of the model
clf_t4 = RandomForestClassifier(n_estimators=10, max_depth=None,
                                 min_samples_split=2, random_state=0)
scores = cross_val_score(clf_t4, X, y, cv=5)
print scores
print scores.mean()
clf_t4.fit(X[dna.fold == 0], y[dna.fold == 0])
pred_t4 = clf_t4.predict(X[dna.fold == 1])
truth = dna[dna.fold == 1]['class_num']
tp = pred_t4[(np.array(pred_t4) == 1) & (np.array(truth) == 1)].size
tn = pred_t4[(np.array(pred_t4) == 0) & (np.array(truth) == 0)].size
fp = pred_t4[(np.array(pred_t4) == 1) & (np.array(truth) == 0)].size
fn = pred_t4[(np.array(pred_t4) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: Random Forest on selected features
End of explanation
X = dna[dna.columns[1:][:-6]]
y = dna.class_num
# compute cross validated accuracy of the model
ada = AdaBoostClassifier(n_estimators=100)
scores = cross_val_score(ada, X, y, cv=5)
print scores
print scores.mean()  
ada.fit(X[dna.fold == 0], y[dna.fold == 0])
pred_ada = ada.predict(X[dna.fold == 1])
truth = dna[dna.fold == 1]['class_num']
tp = pred_ada[(np.array(pred_ada) == 1) & (np.array(truth) == 1)].size
tn = pred_ada[(np.array(pred_ada) == 0) & (np.array(truth) == 0)].size
fp = pred_ada[(np.array(pred_ada) == 1) & (np.array(truth) == 0)].size
fn = pred_ada[(np.array(pred_ada) == 0) & (np.array(truth) == 1)].size
cm = "Confusion Matrix:\n\tX\t\t(+)-pred\t(-)-pred\n" +\
            "\t(+)-truth\t%d\t\t%d\n" +\
            "\t(-)-truth\t%d\t\t%d"
    
print cm % (tp, fn, fp, tn)
Explanation: While there is a significant accuracy improvement going from Decision Tree to Random Forest, the resulting prediction from Extremely Random Forest only improves the accuracy by the margin. Likewise, manually handpicking the features does not seem to improve the performance of the accuracy.
IV: Other Ensemble Learning Methods
AdaBoost
The core principle of AdaBoost is to fit a sequence of weak learners (i.e., models that are only slightly better than random guessing, such as small decision trees) on repeatedly modified versions of the data. The predictions from all of them are then combined through a weighted majority vote (or sum) to produce the final prediction.
[SKL]
Reference:
Y. Freund, and R. Schapire, “A Decision-Theoretic Generalization of On-Line Learning and an Application to Boosting”, 1997.
J. Zhu, H. Zou, S. Rosset, T. Hastie. “Multi-class AdaBoost”, 2009.
End of explanation
<END_TASK>