<jupyter_start><jupyter_text>In this notebook we will test neural-network classification models for classifying biological traits based on gene expression data. We will start by exploring different neural network architectures for a given trait, and then move to investigating the use of deep learning visualization techniques (saliency maps and activation maximization) in order to identify input patterns important for each classification model. --------------------------------  Before we start - handling input data:  1. Choosing the biological trait of focus in this analysisChoose one of the following traits:<jupyter_code>myTrait = 'stimulation'  # IFNb, LPS, dNS1, unstim
# myTrait = 'gender'      # F, M
# myTrait = 'ethnicity'   # African-American, Caucasian, East Asian, MULTI-RACIAL<jupyter_output><empty_output><jupyter_text>2. Loading the gene expression data:The gene expression data is provided as two seperated files: (a) The gene expression table, consisting of 414 genes measured across 2441 samples (b) The features table, including info such as stimultaion, gender and ethnicity origin of the 2441 samples<jupyter_code>features_table = pd.read_table('Input/features.txt')
features = features_table.set_index('geo_accession')
features = np.transpose(features)
GE_table = pd.read_table('Input/GE.txt')
GE = GE_table.set_index('ID_REF')
GE = np.transpose(GE)<jupyter_output><empty_output><jupyter_text>Normalizaition  - Choose whether or not to normalize the gene expression using Z-scores:<jupyter_code>normalize = True
# normalize = False

if normalize:
    GE = (GE - GE.mean())/(GE.std()) # normalizing GE by z-score for each gene
    print("Normalized GE")
else:
    print("Original GE")
num_genes = GE.shape[1]
num_samples = GE.shape[0]

print('Number of genes: %i' %num_genes)
print('Number of samples: %i' %num_samples)<jupyter_output>Number of genes: 414
Number of samples: 2441<jupyter_text>3. Processing the data for classification of the selected trait: Input X  - a vector of 414 genes.  Output Y  - multiclass labeling (or the corresponding one-hot codded representation, dummyY) of the selected trait (stimulation, gender or ethnicity).<jupyter_code># Creating the input matrix X and output matrix Y:
X = GE.values.astype(float)  
Y = features[myTrait]

classes, counts = np.unique(Y, return_counts=True)
num_classes = classes.size
print(dict(zip(classes, counts))) # Check the number of samples in each lable

# Encode class values as integers
encoder = LabelEncoder()
encoder.fit(Y)
encodedY = encoder.transform(Y)

# Convert integers to dummy variables
dummyY = np_utils.to_categorical(encodedY)<jupyter_output>{'IFNb': 432, 'LPS': 806, 'dNS1': 469, 'unstim': 734}<jupyter_text>--------------------------------  Part 1 - Testing neural networks architectures for the classification: Figure 1:  Three neural network architectures will be tested along this section:  ![Net%20architectures%20-%20new.png](attachment:Net%20architectures%20-%20new.png) [1] Shallow neural network - a classification model with one-hidden layer. The number of neurons within the hidden layer is given as a parameter.<jupyter_code>def build_shallow_network_model(X, dummyY, hidden_dim):
    input_dim = X.shape[1]  # 414 genes measures for each sample
    output_dim = dummyY.shape[1]  # 4 classes in case of stimulations & ethnicity, 2 classes in the case of gender
    
    inputs = Input(shape=(input_dim,))
    hidden = Dense(hidden_dim, activation='relu', kernel_regularizer=regularizers.l2(0.01))(inputs)
    outputs = Dense(output_dim, activation='softmax')(hidden)
    myModel = Model(inputs=inputs, outputs=outputs)
    myModel.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])  
    
    return myModel<jupyter_output><empty_output><jupyter_text>[2] Two-hidden-layers neural network - a classification model with two hidden layers. The number of neurons within each hidden layer are given as parameters.<jupyter_code>def build_two_layer_network_model(X, dummyY, hidden_dim, hidden_dim_2):
    input_dim = X.shape[1]  # 414 genes measured for each sample
    output_dim = dummyY.shape[1]  # 4 classes in case of stimulations & ethnicity, 2 classes in the case of gender
    
    inputs = Input(shape=(input_dim,)) 
    hidden1 = Dense(hidden_dim, activation='relu', kernel_regularizer=regularizers.l2(0.01))(inputs)
    hidden2 = Dense(hidden_dim_2, activation='relu', kernel_regularizer=regularizers.l2(0.01))(hidden1)
    predictions = Dense(output_dim, activation='softmax')(hidden2)
    myModel = Model(inputs=inputs, outputs=predictions)
    myModel.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
    
    return myModel<jupyter_output><empty_output><jupyter_text>[3] Logistic regression network - a multiclass model without hidden layers.<jupyter_code>def build_logistic_regression_model(X, dummyY):
    input_dim = X.shape[1]  # 414 genes measures for each sample
    output_dim = dummyY.shape[1]  # 4 classes in case of stimulations & ethnicity, 2 classes in the case of gender
    
    inputs = Input(shape=(input_dim,))
    outputs = Dense(output_dim, activation='softmax')(inputs)
    
    myModel = Model(inputs=inputs, outputs=outputs)
    myModel.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])   
    
    return myModel<jupyter_output><empty_output><jupyter_text>Analysis is conducted by spliting the dataset into 10-folds, training and testing a model on each of the splits. Mean and standard deviation of the accuracy, precision, recall and F1 scores across the 10 splits are reported. Stratified cross validation is used to preserve the proportion of each class in each split.<jupyter_code>def cross_validation_prediction(X, encodedY, dummyY, seed, hidden_dim, hidden_dim_2, debug):
    # Create 10-fold cross validation split:
    k = 10
    cv_scores = pd.DataFrame(np.zeros((k, 5)), columns = ['accuracy','precision','recall','f1','time'])
    kfold = StratifiedKFold(n_splits=k, shuffle=True, random_state=seed)

    i=0
    for train, test in kfold.split(X,encodedY):
        # Choose the neural network architecture of the model:
        myModel = None
        if hidden_dim > 0:
            if hidden_dim_2 == 0:
                if i==0:
                    print("shallow neural network with " + str(hidden_dim) + " dim hidden layer:")
                myModel = build_shallow_network_model(X, dummyY, hidden_dim)  
            else:
                if i==0:
                    print("two-layer neural network with " + str(hidden_dim) + " and " + str(hidden_dim_2) + " dim hidden layer:")
                myModel = build_two_layer_network_model(X, dummyY, hidden_dim, hidden_dim_2)
        else:
            if i==0:
                print("logistic regression model:")
            myModel = build_logistic_regression_model(X, dummyY)

        # Fit model according to train set:    
        print("   split %i" %i)
        start_time = time.time()
        myModel.fit(X[train], dummyY[train], epochs=10, batch_size=5, verbose=0)
        cv_scores['time'] = time.time() - start_time

        # Evaluate accuracy on test set:
        scores = myModel.evaluate(X[test], dummyY[test], verbose=0)
        cv_scores['accuracy'][i] = scores[1] * 100
        if debug==True:
            print("%s: %.2f%%" % (myModel.metrics_names[1], scores[1]*100))
        
        # Eveluate precision, recall and f1 on test set:
        probs = myModel.predict(X[test], batch_size=5, verbose=0)
        predicted = np.argmax(probs, axis=1)
        true = encodedY[test]

        cv_scores['precision'][i] = 100*precision_score(true, predicted, average='weighted')
        cv_scores['recall'][i] = 100*recall_score(true, predicted, average='weighted')
        cv_scores['f1'][i] = 100*f1_score(true, predicted, average='weighted')
        i = i+1

        if debug == True:
            classes, counts = np.unique(predicted, return_counts=True)
            print("predicted:" + str(dict(zip(classes, counts))))
            classes, counts = np.unique(true, return_counts=True)
            print("true:" + str(dict(zip(classes, counts))))
            print("")        
    return cv_scores<jupyter_output><empty_output><jupyter_text>Testing different neural network architectures for classification of biological traits:Here we test different neural network architectures to classify the chosen trait based on an input panel of the 414 genes expression levels for each of the 2441 sammple. The input is a matrix of shape (2441,414), the output is a dummy variable representaion of a label vector of shape (2441,1). Performance is measured based on 10-fold cross validation accuracy scores. Reported are the mean and standard deviation calculated across the splits.    * For shallow network: hidden_dim > 0 + hidden_dim_2 = 0    * For 2-layer network: hidden_dim > 0 + hidden_dim_2 > 0    * For logistic regression: hidden_dim = 0<jupyter_code>import warnings
warnings.filterwarnings("ignore")

architectures = [[2,0],[8,0],[16,0],[32,0],[64,0],[8,2],[16,8],[32,16],[0,0]]
mean_results = pd.DataFrame(np.zeros((len(architectures), 5)), columns = ['accuracy','precision','recall','f1','time'])
std_results = pd.DataFrame(np.zeros((len(architectures), 5)), columns = ['accuracy','precision','recall','f1','time'])

j=0
print("Testing " + str(len(architectures)) + " architectures for " + myTrait + ":")
print("----------------------------------------")
for architecture in architectures:
    cvscores = cross_validation_prediction(X, encodedY, dummyY, seed, architecture[0], architecture[1], debug=False)
    mean_results.iloc[j,:] = cvscores.mean()
    std_results.iloc[j,:] = cvscores.std()
    j=j+1
mean_results_copy = mean_results.copy()
new_results = pd.DataFrame(index=mean_results.index,columns = mean_results.columns)
for row in range(mean_results.shape[0]):
    for col in range(mean_results.shape[1]-1):
        new_results.iloc[row,col] = ("%.2f (+/- %.2f)" % (mean_results.iloc[row,col], std_results.iloc[row,col]))
    new_results.iloc[row,col+1] = ("%.2f" % mean_results.iloc[row,col+1])
    
new_results
 new_results.to_excel('Output/performance_analysis' + myTrait + '.xlsx', myTrait)<jupyter_output><empty_output><jupyter_text>---------------  Part 2 - Utilizing deep learning visualization methods to study input patterns For this part of the analysis we will use a trained model of shallow network with a hidden layer of 8 neurons:<jupyter_code>myModel = None
myModel = build_shallow_network_model(X, dummyY, 8)  
kfold = StratifiedKFold(n_splits=10, shuffle=True, random_state=seed)
for train, test in kfold.split(X,encodedY):
    continue

myModel.fit(X[train], dummyY[train], epochs=10, batch_size=5, verbose=0)
scores = myModel.evaluate(X[test], dummyY[test], verbose=0)
print("Model accuracy: " + str(scores[1]*100))<jupyter_output>Model accuracy: 98.76033038147224<jupyter_text>Convert the softmax activation output to a linear output (needed for both visualizations):<jupyter_code>from vis.utils import utils

myModel.layers[-1].activation = activations.linear
myModel = utils.apply_modifications(myModel)<jupyter_output><empty_output><jupyter_text>Saliency Maps:Saliency Maps are commonly used for revealing a trained neural netowrk's attention given to specific patterns in the input. This is conducted by quantifying the impact that a small change in an input component has on correctly classifying the input vector.Here, we conduct saliency maps analysis to quantify the amount of attention a trained model is giving to each of the genes in each of the sample. Samples are then grouped based on their trait class.Code is based on https://blog.keras.io/how-convolutional-neural-networks-see-the-world.html and https://raghakot.github.io/keras-vis/visualizations/saliency/, adapted for the case of an input vector (instead of an input image).<jupyter_code>def my_saliency_map(myModel, class_idx, sample_vec):    
    input_vec = myModel.input
    output_vec = myModel.output[:, class_idx]    
    
    # Loss is defined according to the class we want to maximize:
    loss = K.mean(myModel.output[:, class_idx])  
 
    # Compute the gradient of the output vector w.r.t the input vector + normalize to avoid very small/large gradients:
    grads = K.gradients(output_vec, input_vec)[0]
    grads /= (K.sqrt(K.mean(K.square(grads))) + 1e-5)
    
    # Calculate the loss and grads given the input vector:
    iterate = K.function([input_vec], [loss,grads])

    # Gradient descent on the given sample_vec:
    alpha = 0.1
    num_steps = 1
    
    sample_vec = sample_vec.reshape(1,num_genes)
    for i in range(num_steps):
        loss_value, grads_value = iterate([sample_vec])
        grads_value -= grads_value*alpha
    
    return grads_value[0]<jupyter_output><empty_output><jupyter_text>In the following cell we aim to visualize the saliency maps obtained for each of the samples on the test set, and compare them to the original gene expression, while grouping the sampls according to their labeling. We further calculate the averaged saliency maps for each class by averaging the attention scores of all samples within each class.<jupyter_code># Create dataframes of gene expression and saliency maps of size num_genes X num_samples:
genes = GE.columns.values.tolist()
expression_all_samples =  pd.DataFrame(genes, columns=['Gene'])
saliency_all_samples =  pd.DataFrame(genes, columns=['Gene'])

# Create a dataframe of averaged saliency maps of size num_genes X num_classes:
saliency_average_samples = pd.DataFrame(genes, columns=['Gene'])

# Saving the attention scores to an output file:
filepath = 'Output/' + myTrait + '_saliency_maps_v2.xlsx'
writer = pd.ExcelWriter(filepath)

# For every class - 
# (i)  Go through each relevant sample to retrieve gene expression and calculate saliency map 
# (ii) Caclulate the averaged saliency across all relevant samples
class_idx=0
for class_label in np.unique(Y): 
    print(class_label)
 
    i=0
    class_samples_indices = np.where(encodedY[test]==class_idx)[0]
    np.random.shuffle(class_samples_indices)   
    num_class_samples = len(class_samples_indices)
    
    class_expression = pd.DataFrame(np.zeros((num_genes, num_class_samples)),columns=[class_label + '_' + str(x) for x in range(num_class_samples)])
    class_saliency = pd.DataFrame(np.zeros((num_genes, num_class_samples)),columns=[class_label + '_' + str(x) for x in range(num_class_samples)])
    
    for sample_idx in class_samples_indices:       
        sample_vec = X[test][sample_idx]
        class_expression.iloc[:,i] = sample_vec
        
        saliency_vec = my_saliency_map(myModel, class_idx, sample_vec)
        class_saliency.iloc[:,i] = saliency_vec
        i+=1
    
    expression_all_samples = expression_all_samples.join(class_expression)
    saliency_all_samples = saliency_all_samples.join(class_saliency)
    saliency_average_samples[class_label] = np.mean(class_saliency,axis=1)

    class_saliency = saliency_all_samples['Gene'].to_frame().join(class_saliency)   
    class_saliency.to_excel(writer,class_label)
    writer.save()
        
    class_idx+=1<jupyter_output>IFNb
LPS
dNS1
unstim<jupyter_text>Figure 2: Plotting the gene expression measurements (top) vs. the calculated saliency maps (bottom) for each sample in test set. Samples are grouped by left to right according to class labels.<jupyter_code>sns.set(font_scale=0.3)

f,(ax1,ax2) = plt.subplots(2,1, figsize=(1.5,2.25), sharex=True)
g1 = sns.heatmap(expression_all_samples.iloc[:,1:],ax=ax1, cmap="BuPu", vmin=-1, vmax=1)
g2 = sns.heatmap(saliency_all_samples.iloc[:,1:],ax=ax2, cmap="BuPu", vmin=-1, vmax=1)
plt.savefig('Output/' + myTrait + '_GE_and_SM_fixed_scale_BuPu.jpg', dpi = 1200)<jupyter_output><empty_output><jupyter_text>Activation Maximization:Activation Maximization is a common method that can be used to generate an input instance that maximizes the activation of a specific class. Here we use the activation maximiaztion to identify gene expression patterns within the input that are important for the classification of each trait. We compare the activation maximiaztion scores calculated for each class to the averaged saliency maps calculate above.Code is based on https://blog.keras.io/how-convolutional-neural-networks-see-the-world.html and https://raghakot.github.io/keras-vis/visualizations/activation_maximization/, adapted for the case of an input vector (instead of an input image).<jupyter_code>def my_activation_maximization(myModel, class_idx):    
    input_vec = myModel.input
    output_vec = myModel.output[:, class_idx] 
    
    # Loss is defined according to the class we want to maximize:
    loss = K.mean(myModel.output[:, class_idx]) 
        
    # Compute the gradient of the loss w.r.t the input vector + normalize to avoid very small/large gradients:
    grads = K.gradients(loss, input_vec)[0]
    grads /= (K.sqrt(K.mean(K.square(grads))) + 1e-5)
    
    # Calculate the loss and grads given the input vector:
    iterate = K.function([input_vec], [loss,grads])

    # Start gradient ascent from a random input vector:
    rand_vec = np.random.random((1, input_vec.shape[1]))
    
    # Gradient ascent:
    alpha = 0.1
    num_steps = 5
    for i in range(num_steps):
        loss_value, grads_value = iterate([rand_vec])
        rand_vec += grads_value*1.
        
    rand_vec = np.transpose(rand_vec)
    
    return rand_vec
# Create a dataframe of activation maximization scores of size num_genes X num_classes:
genes = GE.columns.values.tolist()
act_max =  pd.DataFrame(genes, columns=['Gene'])

# Computing activation maximization for each class:
class_idx=0
for class_label in np.unique(Y): 
    print(class_label)
    
    act_max[class_label] = my_activation_maximization(myModel, class_idx)        
    class_idx+=1
    
filepath = 'Output/' + myTrait + '_ASM_and_AM_scores_v2.xlsx'
writer = pd.ExcelWriter(filepath)
saliency_average_samples.to_excel(writer,'averaged_saliency')
act_max.to_excel(writer,'activation_maximization')
writer.save()<jupyter_output>IFNb
LPS
dNS1
unstim<jupyter_text>Figure 3: Plotting the averaged saliency maps (left) vs. the calculated activation maximization scores (left) for each class label.<jupyter_code>sns.set(font_scale=0.3)

if myTrait=='gender':
    fig, (ax0,ax1) = plt.subplots(1, 2, figsize=(2.2,0.5), sharey=True)
else:
    fig, (ax0,ax1) = plt.subplots(1, 2, figsize=(2.2,1), sharey=True)
sns.heatmap(np.transpose(saliency_average_samples.iloc[:,1:]),cmap="YlGnBu", ax=ax0, vmin=-1, vmax=1, yticklabels=False)
sns.heatmap(np.transpose(act_max.iloc[:,1:]),cmap="YlGnBu", ax=ax1, vmin=-6, vmax=6, yticklabels=False)
plt.subplots_adjust(wspace = 0.2)
plt.savefig('Output/' + myTrait + '_ASM_and_AM_fixed_scale.jpg', dpi = 1200)
print("Correlation between activation maximization and averaged saliency maps:")
act_max.corrwith(saliency_average_samples)<jupyter_output>Correlation between activation maximization and averaged saliency maps:<jupyter_text>--------------------------------  Part 3 - Investigating classification based only on trait-specific genes In this section we show that the performance of a shallow network classification model is not compromised when using only the top scoring genes of each class as input (as opposed to using all 414 genes). Given the high correlation between saliency maps and activation maximization scores, we will focus on the top activation maximization genes.We define the  top-K genes for a given biological trait as the genes that obtained the highest scores when calculating activation maximization with respect to this trait (Part 2). Since activation maximization scores were calculated for each gene  per class, we average the scores of a gene across the classes to obtain a  trait-specific gene score. (Note that average is taken across classes, and not maximum, to prioritize genes that obtained high scores in more than one class; absolute value is taken for highly positive or highly negative scores both contribute to the classification).<jupyter_code>genes = GE.columns.values.tolist()
act_max_scores = act_max.copy()
act_max_scores['avg_abs_activation'] = (act_max[classes]).abs().mean(axis=1)
act_max_scores = act_max_scores.sort_values(by=['avg_abs_activation'], ascending=False)
sorted_genes_act_max = act_max_scores['Gene'].tolist()

act_max_scores.head(10)<jupyter_output><empty_output><jupyter_text>We further focus on exploring the similarity between the activation maximization scores and regression coefficients, which are traditionally used to measure effect sized in logistic regressions. We similarly define the top-K genes as those obtaining the strongest averaged regression coefficients.<jupyter_code>from sklearn.linear_model import LogisticRegression
clf = LogisticRegression(random_state=0, solver='lbfgs',multi_class='multinomial').fit(X, encodedY)

genes = GE.columns.values.tolist()
beta_df =  pd.DataFrame(genes, columns=['Gene'])
beta_df = beta_df.join(pd.DataFrame(clf.coef_.transpose()))

beta_df['max_abs_coef'] = (beta_df.iloc[:, 1:]).abs().max(axis=1)

filepath = 'Output/' + myTrait + '_betas_scores.xlsx'
writer = pd.ExcelWriter(filepath)
beta_df.to_excel(writer,'beta')
writer.save()

beta_df = beta_df.sort_values(by=['max_abs_coef'], ascending=False)
sorted_genes_standard_form = beta_df['Gene'].tolist()

beta_df.head(10)<jupyter_output><empty_output><jupyter_text>Choose the number of genes (the top-K threshold) which will be used as input to the classification model and calculate the accuracy of a shallow neural network with this reduced input:<jupyter_code># k = 1
# k = 2
# k = 3
# k = 5
k = 10
# k = 20
# k = 50
# k = 100
# k = 200
# Top-K genes based on activaiton mazimization:
top_k_act_max_GE = GE[sorted_genes_act_max[0:k]]
X_top_act_max = top_k_act_max_GE.values.astype(float)  

act_max_cvscores = cross_validation_prediction(X_top_act_max, encodedY, dummyY, seed, 8, 0, debug=False)
print("\nClassification based on top %i activation maximization genes: %.2f%% (+/- %.2f%%)" % (k, np.mean(act_max_cvscores['accuracy']), np.std(act_max_cvscores['accuracy'])))
# Top-K genes based on standard logistic regression:
standard_form_GE = GE[sorted_genes_standard_form[0:k]]
X_standard_form = standard_form_GE.values.astype(float)  

standard_form_cvscores = cross_validation_prediction(X_standard_form, encodedY, dummyY, seed, 8, 0, debug=False)
print("\nClassification based on top %i closed-form genes: %.2f%% (+/- %.2f%%)" % (k, np.mean(standard_form_cvscores['accuracy']), np.std(standard_form_cvscores['accuracy'])))<jupyter_output>shallow neural network with 8 dim hidden layer:
   split 0
   split 1
   split 2
   split 3
   split 4
   split 5
   split 6
   split 7
   split 8
   split 9

Classification based on top 10 closed-form genes: 98.56% (+/- 1.03%)<jupyter_text>Print the sorted list of genes based on either activation maximization scores or regression coefficients:<jupyter_code>sorted_genes_act_max
sorted_genes_standard_form<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10833441797908804, 'fraction_numerical': 0.017094017094017096, 'mean_word_length': 4.652280529671407, 'pattern_counts': {'":': 1, '<': 74, '<?xml version=': 0, '>': 78, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '266400', 'license': 'MIT', 'n_tokens_mistral': 7578, 'n_tokens_neox': 7058, 'n_words': 2061, 'path': 'analysis-notebook.ipynb', 'repo_name': 'roniwile/neural-network-visualization-methods-on-gene-expression-data', 'size': '23049'}