<jupyter_start><jupyter_text>Functions<jupyter_code>%matplotlib inline
import os
import tables
import glob
import sys
sys.path.insert(1, os.path.dirname(os.getcwd()) + "/utils/")
import scipy
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import tensorflow as tf
import tensorflow.keras as K
import sklearn.metrics as skm


tf_version = tf.__version__  # ToDo use packaging.version
if tf_version <= '1.13.1':
    from utils.LocallyDirectedConnected import LocallyDirected1D
elif tf_version >= '2.0':
    from utils.LocallyDirectedConnected_tf2 import LocallyDirected1D
else:
    print("unexpected tensorflow version")
    from utils.LocallyDirectedConnected_tf2 import LocallyDirected1D
    
# import LocallyDirectedConnected_tf2
# when using tensorflow 2:
# import tensorflow2.LocallyDirectedConnected as LocallyDirectedConnected


from scipy.sparse import coo_matrix
from scipy import stats
from sklearn.metrics import precision_recall_curve
from sklearn.utils.fixes import signature
from sklearn.metrics import average_precision_score
from utils import sensitivity, specificity, evaluate_performance
tf.keras.backend.set_epsilon(0.0000001)


def weighted_binary_crossentropy(y_true, y_pred): 
    y_true = K.backend.clip(y_true, 0.0001, 1)
    y_pred = K.backend.clip(y_pred, 0.0001, 1)

    return K.backend.mean(-y_true * K.backend.log(y_pred + 0.0001) * weight_positive_class - (1 - y_true) * K.backend.log(
        1 - y_pred + 0.0001) * weight_negative_class)


def get_valdata(datapath):
    yval = pd.read_csv(datapath + "yval_"+studyname+".csv")
    h5file = tables.open_file(datapath  + studyname + '_genotype_processed.h5', "r")
    ybatch = yval["labels"]
    xbatchid = np.array(yval["tot_index"].values, dtype=np.int64)
    xbatch = h5file.root.data[xbatchid, :]
    ybatch = np.reshape(np.array(ybatch), (-1, 1))
    h5file.close()
    return (xbatch, ybatch)


def get_testdata(datapath):
    ytest = pd.read_csv(datapath + "ytest_"+studyname+".csv")
    h5file = tables.open_file(datapath  + studyname + '_genotype_processed.h5', "r")
    ybatch = ytest["labels"]
    xbatchid = np.array(ytest["tot_index"].values, dtype=np.int64)
    xbatch = h5file.root.data[xbatchid, :]
    ybatch = np.reshape(np.array(ybatch), (-1, 1))
    h5file.close()
    return (xbatch, ybatch)


class train_data_generator(K.utils.Sequence):

    def __init__(self, datapath, batch_size, trainsize, startindex, stopindex, shuffle = True):
        self.datapath = datapath
        self.batch_size = batch_size
        self.ytrainsize = trainsize
        self.startindex = startindex
        self.stopindex = stopindex
        self.shuffledindexes = np.arange(trainsize)
        if shuffle:
            np.random.shuffle(self.shuffledindexes)

    
    def __len__(self):
        return int(np.ceil(self.ytrainsize / float(self.batch_size)))

    
    def __getitem__(self, idx):

            batchindexes = self.shuffledindexes[idx * self.batch_size:((idx + 1) * self.batch_size)]
            ytrain = pd.read_csv(datapath + "ytrain_"+studyname+".csv")
            h5file = tables.open_file(self.datapath + studyname + '_genotype_processed.h5', "r")
            ybatch = ytrain["labels"].iloc[batchindexes]
            xbatchid = np.array(ytrain["tot_index"].iloc[batchindexes], dtype=np.int64)
            xbatch = h5file.root.data[xbatchid, :]
            ybatch = np.reshape(np.array(ybatch), (-1, 1))
            h5file.close()
            return xbatch, ybatch
    

    
    def on_epoch_begin(self):
        'Updates indexes after each epoch'
        np.random.shuffle(self.shuffledindexes)
    
    def on_epoch_end(self):
        'Updates indexes after each epoch'
        np.random.shuffle(self.shuffledindexes)
        
        
class val_data_generator(K.utils.Sequence):

    def __init__(self, datapath, batch_size, valsize, startindex, stopindex):
        self.datapath = datapath
        self.batch_size = batch_size
        self.yvalsize = valsize
        self.startindex = startindex
        self.stopindex = stopindex

    def __len__(self):
        val_len = int(np.floor(self.yvalsize / float(self.batch_size)))
        return val_len 
        

    def __getitem__(self, idx):
        yval = pd.read_csv(self.datapath + "yval_"+studyname+".csv")
        h5file = tables.open_file(self.datapath + studyname + '_genotype_processed.h5', "r")
        ybatch = yval["labels"].iloc[idx * self.batch_size:((idx + 1) * self.batch_size)]
        xbatchid = np.array(yval["tot_index"].iloc[idx * self.batch_size:((idx + 1) * self.batch_size)], dtype=np.int64)
        xbatch = h5file.root.data[xbatchid,:]
        ybatch = np.reshape(np.array(ybatch), (-1, 1))
        h5file.close()
        return xbatch, ybatch


def Lasso(inputsize, l1_value):
    inputs = K.Input((inputsize,), name='inputs')
    x1 = K.layers.BatchNormalization(center=False, scale=False, name="inter_out")(inputs)
    x1 = K.layers.Dense(units=1, kernel_regularizer=K.regularizers.l1(l1value))(x1)
    x1 = K.layers.Activation("sigmoid")(x1)
    model = K.Model(inputs=inputs, outputs=x1)
    return model


def GenNet_gene_layer(inputsize):
    mask = scipy.sparse.load_npz( datapath + '/SNP_gene_mask.npz')
    
    input_ = K.Input((inputsize,), name='input_layer')
    input_layer = K.layers.Reshape(input_shape=(inputsize,), target_shape=(inputsize, 1))(input_)
    
    gene_layer = LocallyDirectedConnected.LocallyDirected1D(mask=mask, filters=1, input_shape=(inputsize, 1), name="gene_layer")(input_layer)
    gene_layer = K.layers.Flatten()(gene_layer)
    gene_layer = K.layers.Activation("tanh")(gene_layer)
    gene_layer = K.layers.BatchNormalization(center=False, scale=False, name="inter_out")(gene_layer)

    output_layer = K.layers.Dense(units=1, name="output_layer")(gene_layer)
    output_layer = K.layers.Activation("sigmoid")(output_layer)

    model = K.Model(inputs=input_, outputs=output_layer)
    return model

def GenNet_gene_layer_l1(inputsize, l1_value = 0.01):
    mask = scipy.sparse.load_npz( datapath + '/SNP_gene_mask.npz')
    
    input_ = K.Input((inputsize,), name='input_layer')
    input_layer = K.layers.Reshape(input_shape=(inputsize,), target_shape=(inputsize, 1))(input_)
    
    gene_layer = LocallyDirectedConnected.LocallyDirected1D(mask=mask, filters=1,
                                                            input_shape=(inputsize, 1),
                                                            name="gene_layer",
                                                            activity_regularizer=tf.keras.regularizers.l1(l=0.01))(input_layer)
    gene_layer = K.layers.Flatten()(gene_layer)
    gene_layer = K.layers.Activation("tanh")(gene_layer)
    gene_layer = K.layers.BatchNormalization(center=False, scale=False, name="inter_out")(gene_layer)

    output_layer = K.layers.Dense(units=1, name="output_layer",activity_regularizer=tf.keras.regularizers.l1(l=0.01),
                                  kernel_regularizer = tf.keras.regularizers.l1(l=l1_value) )(gene_layer)
    output_layer = K.layers.Activation("sigmoid")(output_layer)

    model = K.Model(inputs=input_, outputs=output_layer)
    return model<jupyter_output><empty_output><jupyter_text>Main<jupyter_code>jobid = 1
modeltype = "GenNet_gene_layer"

optimizer = "Adadelta" # good to start with
batch_size = 32

namescore = 'score' + str(jobid)

basepath = os.getcwd()[:-4]
datapath = basepath + "/processed_data/"
studyname = str(np.load(datapath + "studyname.npy"))

epochs = 200
l1_value = 0.01
weight_positive_class = 3 # adjust for imbalanced datasets
weight_negative_class = 1

print(studyname)
print(weight_positive_class)
print(weight_negative_class)

train_size = len(pd.read_csv(datapath + "ytrain_"+studyname+".csv"))
val_size = len(pd.read_csv(datapath + "yval_"+studyname+".csv"))


if optimizer == "Adam":
    lr_opt = 0.0006 # seems to work in most cases
    optimizer = tf.keras.optimizers.Adam(lr = lr_opt)
if optimizer == "Adadelta":
    optimizer = tf.keras.optimizers.Adadelta()
    

folder = (str(studyname)+ "__" +str(jobid) )

h5file = tables.open_file(datapath  + studyname + '_genotype_processed.h5', "r")
data_shape = h5file.root.data.shape
inputsize = h5file.root.data.shape[1]
startindex = 0
stopindex= -1
h5file.close()

rfrun_path = "//media/avanhilten/pHDD1TB/SCZ/results/" + folder + "/"
if not os.path.exists(rfrun_path):
    print("Runpath did not exist but is made now")
    os.mkdir(rfrun_path)

print("jobid =  " + str(jobid))
print("folder = " + str(folder))
print("batchsize = " + str(batch_size))
print("n_features " + str(inputsize))

if modeltype == "GenNet_gene_layer":
    model = GenNet_gene_layer(inputsize=int(inputsize))    
if modeltype == "GenNet_gene_layer_l1":
    model = GenNet_gene_layer_l1(inputsize=int(inputsize))    
if modeltype == "Lasso":
    model = Lasso(inputsize=int(inputsize), l1_value = l1_value)  
    
model.compile(loss=weighted_binary_crossentropy, optimizer=optimizer, metrics=["accuracy",sensitivity, specificity])

print(model.summary())
model_summary = str(model.to_json())

with open(rfrun_path + '/experiment_stats_results_.txt', 'a') as f:
    f.write('gtname = ' + str(studyname))
    f.write('\n jobid = ' + str(jobid))
    f.write('\n model = ' + str(modeltype))
    f.write('\n batchsize = ' + str(batch_size))
    f.write('\n weightnegative = ' + str(weight_negative_class))
    f.write('\n weightpositive = ' + str(weight_positive_class))

with open(rfrun_path + '/experiment_summary_model.txt', 'w') as fh:
    model.summary(print_fn=lambda x: fh.write(x + '\n'))

csv_logger = K.callbacks.CSVLogger(rfrun_path + 'log.csv', append=True, separator=';')
earlystop =K.callbacks.EarlyStopping(monitor='val_loss', min_delta=0.0001, patience=25, verbose=1, mode='auto')
saveBestModel = K.callbacks.ModelCheckpoint(rfrun_path + "bestweight_job.h5", monitor='val_loss',
                                          verbose=1, save_best_only=True, mode='auto')
reduce_lr = K.callbacks.ReduceLROnPlateau(monitor='val_loss', factor=0.2,
                              patience=5, min_lr=0.001)
# %%
if os.path.exists(rfrun_path + '/bestweight_job.h5'):
    print('loading weights')
    model.load_weights(rfrun_path + '/bestweight_job.h5')
    print("evaluate over " + str(data_shape[0]) + " patients")
    xval, yval = get_valdata(datapath)
    evaluate_val = model.evaluate(xval, yval)
else:
    history = model.fit_generator(generator=train_data_generator(datapath=datapath, batch_size=batch_size, trainsize=int(train_size),
                                                        startindex=startindex, stopindex=stopindex),
                        shuffle=True,
                        epochs=epochs,
                        verbose=1,
                        callbacks=[earlystop, saveBestModel, csv_logger, reduce_lr],
                        workers=5,
                        use_multiprocessing=True,
                        validation_data=val_data_generator(datapath=datapath, batch_size=batch_size, valsize=val_size,
                                                            startindex=startindex, stopindex=stopindex)
                        )
    plt.plot(history.history['loss'])
    plt.plot(history.history['val_loss'])
    plt.title('model loss')
    plt.ylabel('loss')
    plt.xlabel('epoch')
    plt.legend(['train', 'val'], loc='upper left')
    plt.show()

print("Finished")



print("evaluate over patients with the model from the latest epoch")
xval, yval = get_valdata(datapath)
pval = model.predict(xval)
evaluate_performance(yval, pval)

print('loading weights')
model.load_weights(rfrun_path + '/bestweight_job.h5')
print("evaluate over the patients with the best model")
pval = model.predict(xval)
evaluate_performance(yval, pval)
xtest, ytest = get_testdata(datapath)
ptest = model.predict(xtest)
evaluate_performance(ytest, ptest)
np.save(rfrun_path + "/ptest.npy", ptest)
colors = ['#7dcfe2','#4b78b5','darkgrey','dimgray','#7dcfe2','#4b78b5','darkgrey','dimgray','#7dcfe2','#4b78b5','darkgrey','dimgray','#7dcfe2','#4b78b5','darkgrey','dimgray','#7dcfe2','#4b78b5','darkgrey','dimgray','#7dcfe2','#4b78b5','darkgrey','dimgray']


plt.figure(figsize=(20, 10))
num_annotated = 10
gene_end  = np.load(datapath + "gene_end.npy")
gene_middle = []
for i in  range(len(gene_end) - 1):
    gene_middle.append((gene_end[i]+ gene_end[i+1])/2)



y = model.layers[6].get_weights()[0]

y = y / max(y)
x = np.arange(len(y))

plt.ylim(bottom=-1.2, top=1.19)
plt.xlim(0, len(x) + 100)

plt.xlabel("Genes (colored per chromosome)", size=24)
plt.ylabel("Weights", size=24)

gene_overview = pd.read_csv(datapath + "/gene_overview.csv", sep = "\t")
gene_overview['mean'] = y
gene_overview = gene_overview.sort_values(["CHR", "bp"], ascending = (True, True))
gene_overview["pos"] = x

    
x = gene_overview["pos"].values
y = gene_overview["mean"].values

for i in range(len(gene_end) - 1):
    plt.scatter(x[gene_end[i]:gene_end[i + 1]], y[gene_end[i]:gene_end[i + 1]],  c=colors[i])
    
gene_overview_annotate = gene_overview.sort_values("mean", ascending=False).head(num_annotated)
top_hits_genes = gene_overview.sort_values("mean", ascending=False).copy()


for i in range(num_annotated):
    plt.annotate(gene_overview_annotate["gene"].iloc[i],
                 (gene_overview_annotate["pos"].iloc[i], gene_overview_annotate["mean"].iloc[i]),
                 xytext=(gene_overview_annotate["pos"].iloc[i] + 100,
                         gene_overview_annotate["mean"].iloc[i] + np.random.randint(100) * 0.0006), size=16)


plt.ylim(bottom = 0)
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.spines.left'] = True
plt.yticks( size=16)
plt.xticks(gene_middle, np.arange(len(gene_middle))+1, size=16)
plt.savefig(datapath + 'Manhattan_genes_blue.png', dpi = 300)
plt.show()
import h5py
from tensorflow.keras.models import Model

ytrain = np.reshape(np.array(pd.read_csv(datapath + "ytrain_"+studyname+".csv")["labels"]), (-1, 1))
ptrain = np.squeeze(model.predict_generator(train_data_generator(datapath=datapath, batch_size=100, trainsize=int(train_size),
                                                startindex=startindex, stopindex=stopindex, shuffle = False), workers=8, use_multiprocessing=True ))

np.save(rfrun_path + "/ptrain.npy", ptrain)


intermediate_layer_model = Model(inputs=model.input,
                                 outputs=model.get_layer(name='inter_out').output)

intermediate_layer_model.compile(loss=weighted_binary_crossentropy, optimizer=optimizer,
                                 metrics=["accuracy", sensitivity, specificity])

intermediate_output = np.squeeze(intermediate_layer_model.predict_generator(train_data_generator(datapath=datapath, batch_size=100, trainsize=int(train_size),
                                               startindex=startindex, stopindex=stopindex, shuffle = False), workers=8, use_multiprocessing=True ))



np.save(rfrun_path + "/intermediate_output.npy", intermediate_output)
norm_mean = np.mean(intermediate_output, axis=0)
norm_std = np.std(intermediate_output, axis=0)

plt.figure()
plt.hist(norm_mean)
plt.xlabel("mean")
plt.ylabel("number")
plt.title("hist mean of activations")
plt.savefig(rfrun_path + "mean_act_hist")
plt.show()

plt.figure()
plt.hist(norm_std)
plt.xlabel("std")
plt.ylabel("number")
plt.title("hist std of activations")
plt.savefig(rfrun_path + "std_act_hist")
plt.show()

np.save(rfrun_path + "/ptrain.npy", ptrain)
print("mean ptrain = " + str(np.mean(ptrain)))
print("min ptrain = " + str(np.min(ptrain)))
print("max ptrain = " + str(np.max(ptrain)))

print("\n f1 = " + str(skm.f1_score(ytrain, ptrain.round())))
print("\n confusion matrix")
cm = skm.confusion_matrix(ytrain, ptrain.round())
print(cm)
ptrain.max()

fpr, tpr, thresholds = skm.roc_curve(ytrain, ptrain)
roc_auc = skm.auc(fpr, tpr)

optimal_idx = np.argmax(tpr - fpr)
optimal_threshold = thresholds[optimal_idx]
print(optimal_threshold)

lw = 2
Specificity1 = cm[0, 0] / (cm[0, 0] + cm[0, 1])
print('Specificity : ', Specificity1)

Sensitivity1 = cm[1, 1] / (cm[1, 0] + cm[1, 1])
print('Sensitivity : ', Sensitivity1)

plt.figure()
plt.plot(fpr, tpr, color='darkorange',
         lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.plot(1 - Specificity1, Sensitivity1, color='b', marker='o')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate / Recall')
plt.ylabel('True Positive Rate / Precision')
plt.title('AUCROC and Precision-Recall curve:')

print("\n train: loss, acc, sensitivity, specificity")
# print(evaluate_train)
print("\n AUC = " + str(roc_auc) + '\n')
print(skm.classification_report(y_pred=np.round(ptrain), y_true=ytrain))

with open(rfrun_path + '/experiment_stats_results_.txt', 'a') as f:
    f.write(str(cm))
    f.write("\n")
    f.write("\n mean ptrain = " + str(np.mean(ptrain)))
    f.write("\n f1 = " + str(skm.f1_score(ytrain, ptrain.round())))
    f.write("\n auc = " + str(roc_auc))

average_precision = average_precision_score(ytrain, ptrain)
precision, recall, thresholds = precision_recall_curve(ytrain, ptrain)

# In matplotlib < 1.5, plt.fill_between does not have a 'step' argument
step_kwargs = ({'step': 'post'}
               if 'step' in signature(plt.fill_between).parameters
               else {})
plt.step(recall, precision, color='b', alpha=0.2,
         where='post', label='Average Precision (area = {0:0.2f})'.format(
        average_precision))
plt.fill_between(recall, precision, alpha=0.2, color='b', **step_kwargs)

plt.savefig(rfrun_path + studyname + "_ROCtrain.png")
plt.legend(loc="lower right")
plt.show()    



data = h5py.File(rfrun_path + "/bestweight_job.h5", "r")
dname = [s for s in [key for key in data["model_weights"].keys()] if "output" in s]
weights_dense = np.array(data["model_weights"][dname[0]][dname[0]]["kernel:0"]).flatten()
data.close()


x= np.arange(len(weights_dense))

results_fold = pd.DataFrame(np.array([weights_dense, x]).T, columns=['weights_dense', "pos"])
results_fold["norm_mean"]  = norm_mean
results_fold["norm_std"] = norm_std
results_fold.to_csv(rfrun_path + "results_fold.csv")

compare_weight_stdcor = results_fold.weights_dense.values * results_fold.norm_std.values

x = np.arange(len(compare_weight_stdcor))
compare_weight_stdcor = np.abs(compare_weight_stdcor) / np.max(np.abs(compare_weight_stdcor))


plt.figure(figsize=(20, 10))
num_annotated = 10
gene_end = np.load(datapath + "gene_end.npy")

for i in range(len(gene_end) - 1):
    plt.scatter(x[gene_end[i]:gene_end[i + 1]], compare_weight_stdcor[gene_end[i]:gene_end[i + 1]], c=colors[i])
        
plt.ylim(bottom= 0 , top=1.2)
plt.xlim(0, len(compare_weight_stdcor )+ 100)
plt.title("Gene Importance", size=36)
plt.xlabel("Gene coordinate (colored per chromosome)", size=18)
plt.ylabel("Importance ", size=18)

gene_overview = pd.read_csv(datapath + "gene_overview.csv", sep = "\t")
gene_overview['mean'] = compare_weight_stdcor
gene_overview["pos"] = x
gene5_overview = gene_overview.sort_values("mean", ascending=False).head(num_annotated)
top_hits_genes = gene_overview.sort_values("mean", ascending=False).copy()


for i in range(num_annotated):
    plt.annotate(gene5_overview["gene"].iloc[i],
                 (gene5_overview["pos"].iloc[i], gene5_overview["mean"].iloc[i]),
                 xytext=(gene5_overview["pos"].iloc[i] + 100,
                         gene5_overview["mean"].iloc[i]), size=16)

gene5_overview = gene_overview.sort_values("mean", ascending=True).head(num_annotated)


plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')

plt.savefig(rfrun_path + "geneimportance_train_10f_c_wide.png", bbox_inches='tight', pad_inches=0)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14527333131984294, 'fraction_numerical': 0.018071076210611093, 'mean_word_length': 4.043665905052044, 'pattern_counts': {'":': 5, '<': 11, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '595257', 'license': 'Apache-2.0', 'n_tokens_mistral': 7449, 'n_tokens_neox': 7006, 'n_words': 1153, 'path': 'jupyter_notebooks/4_GenNet_run_network.ipynb', 'repo_name': 'reneevdw/GenNet', 'size': '19866'}