<jupyter_start><jupyter_text>Creating Data Report<jupyter_code>#profile = ProfileReport(data, \
#                        title="Data Profile")
#profile.to_file("./data_report.html")<jupyter_output><empty_output><jupyter_text>Finding and counting missing values<jupyter_code>data.isna().sum()<jupyter_output><empty_output><jupyter_text>Plotting two important columns in data<jupyter_code>fig, ax = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(12,10))
#ax = ax.ravel() 
mask = data['simplified.significance'] == 'benign'
ax[0,0].hist(data[mask]['foldx.ddg'].values)
ax[0,0].set_ylabel('benign')
ax[0,0].set_title('foldx.ddg')
ax[0,1].hist(data[mask]['covar.2d'].values)
ax[0,1].set_title('covar.2d')
mask = data['simplified.significance'] == 'pathogenic'
ax[1,0].hist(data[mask]['foldx.ddg'].values)
#ax[1,0].set_title('Histogram: foldx.ddg (pathogenic)')
ax[1,0].set_ylabel('pathogenic')
ax[1,1].hist(data[mask]['covar.2d'].values)
ax[0,1].set_title('covar.2d')
plt.setp([a.get_xticklabels() for a in ax[0, :]], visible=False)
plt.setp([a.get_yticklabels() for a in ax[:, 1]], visible=False)
fig.tight_layout()
fig.subplots_adjust(top=0.88)
#plt.tight_layout()
#plt.savefig('./histograms.png')<jupyter_output>/home/dispel/anaconda3/lib/python3.7/site-packages/numpy/lib/histograms.py:829: RuntimeWarning: invalid value encountered in greater_equal
  keep = (tmp_a >= first_edge)
/home/dispel/anaconda3/lib/python3.7/site-packages/numpy/lib/histograms.py:830: RuntimeWarning: invalid value encountered in less_equal
  keep &= (tmp_a <= last_edge)<jupyter_text>Extracting features from amino acids using NLF Encode<jupyter_code>nlf = pd.read_csv('https://raw.githubusercontent.com/dmnfarrell/epitopepredict/master/epitopepredict/mhcdata/NLF.csv',index_col=0)

def nlf_encode(seq):    
    '''
    '''
    
    x = pd.DataFrame([nlf[i] for i in seq]).reset_index(drop=True)  
    #show_matrix(x)
    e = x.values.flatten()
    return e
nlf_wt = []
for aa in data['AA.wt'].values:
    nlf_wt.append(nlf_encode(aa))
nlf_wt = np.array(nlf_wt)
nlf_mut = []
for aa in data['AA.mut'].values:
    nlf_mut.append(nlf_encode(aa))
nlf_mut = np.array(nlf_mut)
std = StandardScaler()
nlf_wt = std.fit_transform(nlf_wt)
nlf_mut = std.fit_transform(nlf_mut)<jupyter_output><empty_output><jupyter_text>One hot encoding and scaling foldx.ddg and covar.2d features<jupyter_code>def to_hot(df, col):
    '''
    '''
    na_mask = df[col].notna()
    temp = np.zeros((len(df), 2))
    mu = df[col].mean()
    sd = df[col].std()
    df[col][na_mask] = (df[col][na_mask] - mu)/sd
    temp[np.where(na_mask == True),0] = 1
    temp[np.where(na_mask == True),1] = df[na_mask][col]
    
    return temp
onehot_foldx = to_hot(data, 'foldx.ddg')
onehot_covar = to_hot(data, 'covar.2d')<jupyter_output>/home/dispel/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:8: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy<jupyter_text>Scaling Allele Frequency feature<jupyter_code>data_allele = data['Allele.Frequency']
mu = data_allele.mean()
sd = data_allele.std()
data_allele = (data_allele - mu)/sd 
def to_onehot(data, col):
    '''
    '''
    from sklearn.preprocessing import OneHotEncoder
    from sklearn.preprocessing import LabelEncoder
    
    ohe = OneHotEncoder()
    le = LabelEncoder()
    ohe_fit = ohe.fit_transform(data[col].values.reshape(-1,1)).toarray()
    
    return ohe_fit
out_inh = to_onehot(data, 'inheritance')
out_disease = to_onehot(data, 'disease.cat')
cus = np.hstack((out_inh, out_disease))<jupyter_output><empty_output><jupyter_text>Combining all features into a new matrix<jupyter_code>features = np.c_[nlf_wt, nlf_mut, onehot_foldx, onehot_covar, cus, data_allele.values]<jupyter_output><empty_output><jupyter_text>Preparing data for machine learning<jupyter_code>le = LabelEncoder()

mask1 = data['simplified.significance'] == 'pathogenic'
mask2 = data['simplified.significance'] == 'benign'
labels = data['simplified.significance'][(mask1 == True) | (mask2 == True)]

y_le = le.fit_transform(labels)
new_data = features[np.where((mask1 == True) | (mask2 == True))]
X_train, X_test, y_train, y_test = train_test_split(new_data, y_le, \
                                                   test_size=0.2, random_state=1)<jupyter_output><empty_output><jupyter_text>GP Optimization to find best parameters<jupyter_code>from skopt.space import Real, Integer, Categorical
from skopt.utils import use_named_args
from sklearn.model_selection import cross_val_score
from skopt import gp_minimize
from sklearn.svm import SVC
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Real(10**-3, 10**3, "uniform", name='gamma'),\
          Categorical(['poly', 'rbf', 'linear', 'sigmoid'], name="kernel")]
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Real(10**-3, 10**3, "uniform", name='gamma'),\
          Categorical(['rbf', 'linear', 'sigmoid', 'poly'], name='kernel')]

clf = SVC()

@use_named_args(space)
def objective(**params):
    clf.set_params(**params)

    return -np.mean(cross_val_score(clf, X_train, y_train, cv=5, n_jobs=-1,
                                    scoring="neg_mean_absolute_error"))
res = gp_minimize(objective, space, n_calls=20, acq_func='EI',\
                  n_random_starts=2, random_state=1)
res.x
clf = SVC(C=res.x[0], gamma=res.x[1], kernel=res.x[2])
clf.fit(X_train, y_train)
clf.score(X_train, y_train), clf.score(X_test, y_test)
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import auc
from sklearn.metrics import plot_roc_curve

cv = StratifiedKFold(n_splits=6)
classifier = SVC(C=res.x[0], gamma=res.x[1], \
                     kernel=res.x[2],\
                 random_state=1)
tprs = []
aucs = []
mean_fpr = np.linspace(0, 1, 100)
fig, ax = plt.subplots(figsize=(12,8))
for i, (train, test) in enumerate(cv.split(new_data, y_le)):
    classifier.fit(new_data[train], y_le[train])
    viz = plot_roc_curve(classifier, new_data[test], y_le[test],
                         name='ROC fold {}'.format(i),
                         alpha=0.3, lw=1, ax=ax)
    interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)
    interp_tpr[0] = 0.0
    tprs.append(interp_tpr)
    aucs.append(viz.roc_auc)
    
ax.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',
        label='Chance', alpha=.8)

mean_tpr = np.mean(tprs, axis=0)
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
std_auc = np.std(aucs)
ax.plot(mean_fpr, mean_tpr, color='b',\
        label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),\
        lw=2, alpha=.8)
std_tpr = np.std(tprs, axis=0)
tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
ax.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,\
                label=r'$\pm$ 1 std. dev.')

ax.set(xlim=[-0.05, 1.05], ylim=[-0.05, 1.05],\
       title="Receiver operating characteristic SVC")
ax.legend(loc="lower right")
plt.savefig('./roc_auc_svc.png')
plt.show()
from skopt.space import Real, Integer, Categorical
from skopt.utils import use_named_args
from sklearn.model_selection import cross_val_score
from skopt import gp_minimize
from sklearn.linear_model import LogisticRegression
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Real(10**-3, 10**3, "uniform", name='gamma'),\
          Categorical(['poly', 'rbf', 'linear', 'sigmoid'], name="kernel")]
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Categorical(['l1', 'l2'], name='penalty'),\
         Categorical([ 'liblinear', 'saga'], name='solver')]

lr = LogisticRegression()

@use_named_args(space)
def objective(**params):
    lr.set_params(**params)

    return -np.mean(cross_val_score(lr, X_train, y_train, cv=5, n_jobs=-1,
                                    scoring="neg_mean_absolute_error"))
res_lr = gp_minimize(objective, space, n_calls=20, acq_func='EI',\
                  n_random_starts=2, random_state=1)
res_lr.x
lr = LogisticRegression(C=res_lr.x[0], penalty=res_lr.x[1],
                       solver=res_lr.x[2])
lr.fit(X_train, y_train)
lr.score(X_train, y_train), lr.score(X_test, y_test)
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import auc
from sklearn.metrics import plot_roc_curve

cv = StratifiedKFold(n_splits=6)
classifier = LogisticRegression(C=res_lr.x[0], penalty=res_lr.x[1], \
                     solver=res_lr.x[2],\
                 random_state=1)
tprs = []
aucs = []
mean_fpr = np.linspace(0, 1, 100)
fig, ax = plt.subplots(figsize=(12,8))
for i, (train, test) in enumerate(cv.split(new_data, y_le)):
    classifier.fit(new_data[train], y_le[train])
    viz = plot_roc_curve(classifier, new_data[test], y_le[test],
                         name='ROC fold {}'.format(i),
                         alpha=0.3, lw=1, ax=ax)
    interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)
    interp_tpr[0] = 0.0
    tprs.append(interp_tpr)
    aucs.append(viz.roc_auc)
    
ax.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',
        label='Chance', alpha=.8)

mean_tpr = np.mean(tprs, axis=0)
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
std_auc = np.std(aucs)
ax.plot(mean_fpr, mean_tpr, color='b',\
        label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),\
        lw=2, alpha=.8)
std_tpr = np.std(tprs, axis=0)
tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
ax.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,\
                label=r'$\pm$ 1 std. dev.')

ax.set(xlim=[-0.05, 1.05], ylim=[-0.05, 1.05],\
       title="Receiver operating characteristic Logistic Regression")
ax.legend(loc="lower right")
plt.savefig('./roc_auc_lr.png')
plt.show()<jupyter_output><empty_output><jupyter_text>Base<jupyter_code>data.head()
mask1 = data['foldx.ddg'].notna()
mask2 = data['covar.2d'].notna()
data = data[mask1 & mask2]
X = data[['foldx.ddg', 'covar.2d']].values
y = data['simplified.significance'].values
X = X[np.where((y == 'benign') | (y == 'pathogenic'))]
y = y[np.where((y == 'benign') | (y == 'pathogenic'))]
le = LabelEncoder()
y_le_base = le.fit_transform(y)
X_train, X_test, y_train, y_test = train_test_split(X, y_le_base, test_size=0.2,\
                                                   random_state=1)
lr = LogisticRegression()
lr.fit(X_train, y_train)
lr.score(X_train, y_train), lr.score(X_test, y_test)
from skopt.space import Real, Integer, Categorical
from skopt.utils import use_named_args
from sklearn.model_selection import cross_val_score
from skopt import gp_minimize
from sklearn.linear_model import LogisticRegression
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Real(10**-3, 10**3, "uniform", name='gamma'),\
          Categorical(['poly', 'rbf', 'linear', 'sigmoid'], name="kernel")]
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Categorical(['l1', 'l2'], name='penalty'),\
         Categorical([ 'liblinear', 'saga'], name='solver')]

lr = LogisticRegression()

@use_named_args(space)
def objective(**params):
    lr.set_params(**params)

    return -np.mean(cross_val_score(lr, X_train, y_train, cv=5, n_jobs=-1,
                                    scoring="neg_mean_absolute_error"))
res_lr_base = gp_minimize(objective, space, n_calls=20, acq_func='EI',\
                  n_random_starts=2, random_state=1)
res_lr_base.x
lr = LogisticRegression(C=res_lr_base.x[0], penalty=res_lr_base.x[1],
                       solver=res_lr_base.x[2])
lr.fit(X_train, y_train)
lr.score(X_train, y_train), lr.score(X_test, y_test)
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import auc
from sklearn.metrics import plot_roc_curve

cv = StratifiedKFold(n_splits=6)
classifier = LogisticRegression(random_state=1)
tprs = []
aucs = []
mean_fpr = np.linspace(0, 1, 100)
fig, ax = plt.subplots(figsize=(12,8))
for i, (train, test) in enumerate(cv.split(X, y_le_base)):
    classifier.fit(X[train], y_le_base[train])
    viz = plot_roc_curve(classifier, X[test], y_le_base[test],
                         name='ROC fold {}'.format(i),
                         alpha=0.3, lw=1, ax=ax)
    interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)
    interp_tpr[0] = 0.0
    tprs.append(interp_tpr)
    aucs.append(viz.roc_auc)
    
ax.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',
        label='Chance', alpha=.8)

mean_tpr = np.mean(tprs, axis=0)
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
std_auc = np.std(aucs)
ax.plot(mean_fpr, mean_tpr, color='b',\
        label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),\
        lw=2, alpha=.8)
std_tpr = np.std(tprs, axis=0)
tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
ax.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,\
                label=r'$\pm$ 1 std. dev.')

ax.set(xlim=[-0.05, 1.05], ylim=[-0.05, 1.05],\
       title="Receiver operating characteristic Logistic Regression Base")
ax.legend(loc="lower right")
plt.savefig('./roc_auc_lr_base.png')
plt.show()
svc = SVC()
svc.fit(X_train, y_train)
svc.score(X_train, y_train), svc.score(X_test, y_test)
from skopt.space import Real, Integer, Categorical
from skopt.utils import use_named_args
from sklearn.model_selection import cross_val_score
from skopt import gp_minimize
from sklearn.svm import SVC
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Real(10**-3, 10**3, "uniform", name='gamma'),\
          Categorical(['poly', 'rbf', 'linear', 'sigmoid'], name="kernel")]
'''
space  = [Real(10**-3, 10**3, "uniform", name='C'),\
          Real(10**-3, 10**3, "uniform", name='gamma'),\
          Categorical(['rbf', 'linear', 'sigmoid', 'poly'], name='kernel')]

clf = SVC()

@use_named_args(space)
def objective(**params):
    clf.set_params(**params)

    return -np.mean(cross_val_score(clf, X_train, y_train, cv=5, n_jobs=-1,
                                    scoring="neg_mean_absolute_error"))
res_svc_base = gp_minimize(objective, space, n_calls=20, acq_func='EI',\
                  n_random_starts=2, random_state=1, verbose=True)
res_svc_base.x
clf = SVC(C=res_svc_base.x[0], gamma=res_svc_base.x[1], \
          kernel=res_svc_base.x[2])
clf.fit(X_train, y_train)
clf.score(X_train, y_train), clf.score(X_test, y_test)
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import auc
from sklearn.metrics import plot_roc_curve

cv = StratifiedKFold(n_splits=6)
classifier = SVC(C=res_svc_base.x[0], gamma=res_svc_base.x[1], \
                     kernel=res_svc_base.x[2],\
                 random_state=1)
tprs = []
aucs = []
mean_fpr = np.linspace(0, 1, 100)
fig, ax = plt.subplots(figsize=(12,8))
for i, (train, test) in enumerate(cv.split(X, y_le_base)):
    classifier.fit(X[train], y_le_base[train])
    viz = plot_roc_curve(classifier, X[test], y_le_base[test],
                         name='ROC fold {}'.format(i),
                         alpha=0.3, lw=1, ax=ax)
    interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)
    interp_tpr[0] = 0.0
    tprs.append(interp_tpr)
    aucs.append(viz.roc_auc)
    
ax.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',
        label='Chance', alpha=.8)

mean_tpr = np.mean(tprs, axis=0)
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
std_auc = np.std(aucs)
ax.plot(mean_fpr, mean_tpr, color='b',\
        label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),\
        lw=2, alpha=.8)
std_tpr = np.std(tprs, axis=0)
tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
ax.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,\
                label=r'$\pm$ 1 std. dev.')

ax.set(xlim=[-0.05, 1.05], ylim=[-0.05, 1.05],\
       title="Receiver operating characteristic SVC base")
ax.legend(loc="lower right")
plt.savefig('./roc_auc_svc_base.png')
plt.show()
data = data.dropna()
AC = data[['covar.2d', 'Allele.Frequency']].values
patho = data['simplified.significance'].values
AC = AC[np.where((patho == 'pathogenic') | (patho == 'benign'))]
patho = patho[np.where((patho == 'pathogenic') | (patho == 'benign'))]
le = LabelEncoder()
patho_le = le.fit_transform(patho)
std = StandardScaler()
AC_std = std.fit_transform(AC)
X_train, X_test, y_train, y_test = train_test_split(AC_std, patho_le, test_size=0.2,\
                                                   random_state=1)
lr = LogisticRegression()
lr.fit(X_train, y_train)
lr.score(X_train, y_train), lr.score(X_test, y_test)
clf = SVC()
clf.fit(X_train, y_train)
clf.score(X_test, y_test)
color = np.where(patho == 'pathogenic', 'red', 'blue')
plt.figure(figsize=(12,8))
plt.scatter(AC[:,1][patho == 'pathogenic'], AC[:,0][patho == 'pathogenic'], \
            c=color[patho == 'pathogenic'], alpha=0.4, label='pathogenic')
plt.scatter(AC[:,1][patho == 'benign'], AC[:,0][patho == 'benign'], \
            c=color[patho == 'benign'], alpha=0.4, label='benign')
plt.xlabel('Allele.Frequency')
plt.ylabel('covar.2d')
plt.title('Alle.Frequency and covar.2d Scatter Plot')
plt.legend()
plt.savefig('./scatterplot.png')
import statsmodels.api as sm
data = sm.datasets.scotland.load(as_pandas=False)
X_train = sm.add_constant(X_train)
mod = sm.GLM(y_train, X_train, \
             family=sm.families.Binomial(link=sm.genmod.families.links.logit))
mod_results = mod.fit()
print(mod_results.summary())
X_test = sm.add_constant(X_test)
preds = (mod.predict(mod_results.params, X_test) > 0.5) * 1
from sklearn.metrics import accuracy_score
accuracy_score(y_test, preds)
lr = LogisticRegression()
lr.fit(X_train, y_train)
lr.score(X_test, y_test)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17068520802960638, 'fraction_numerical': 0.023943030167096556, 'mean_word_length': 4.085543199315654, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 41, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '140572', 'license': 'Apache-2.0', 'n_tokens_mistral': 7417, 'n_tokens_neox': 7062, 'n_words': 1227, 'path': 'Notebook.ipynb', 'repo_name': 'anilk991/BioHackathon2020', 'size': '17834'}