<jupyter_start><jupyter_text>Imports<jupyter_code># Import standard libraries
import csv
import OmicsIntegrator as oi
import pandas as pd
import numpy as np
import scipy.stats as ss
import networkx as nx
import matplotlib.pyplot as plt
from matplotlib_venn import venn2, venn2_circles, venn3, venn3_circles
import pickle
import tqdm
import time

# Import custom libraries
import drugs2
import prizes
import sensitivity_analysis as sensitivity
import sensitivity_sarspartners as svp
import neighborhoods as nbh<jupyter_output><empty_output><jupyter_text>Create prized list  Select IAV differentially expressed genes (with aging information)<jupyter_code># Load IAV DE genes data
iav_de_genes_file_name = '../Data/iav_genes.csv'
iav_genes_df = pd.read_csv(iav_de_genes_file_name)
# Minor reformatting
iav_genes_df.columns = ['name','log2FC_blanco','log2FC_ageing']
iav_genes_df = iav_genes_df.replace([np.inf, -np.inf], np.nan)
iav_genes_df = iav_genes_df.dropna(subset = ['log2FC_blanco','log2FC_ageing'], how='any')
iav_genes_df['name'] = iav_genes_df['name'].str.upper().str.strip()
# Keep only protein coding genes
protein_coding_genes_file_name = '../Data/protein_coding_ensembl_gene_id_hgnc_hg19.txt'
coding_genes = prizes.load_protein_coding_genes(protein_coding_genes_file_name)
iav_genes_df = iav_genes_df.merge(coding_genes, on = 'name', how = 'inner')
iav_genes_df.head()<jupyter_output><empty_output><jupyter_text>Create final prized list<jupyter_code>terminal_df = iav_genes_df
terminal_df = terminal_df.loc[terminal_df['log2FC_blanco']*terminal_df['log2FC_ageing']>0]
terminal_df.insert(1,'prize', np.abs(terminal_df['log2FC_blanco']))
print(len(terminal_df))
terminal_df.head()
# Save terminal df to tsv
terminal_df.to_csv(r'../Save_iav/terminals_ppi_analysis.tsv', header=True, index=None, sep='\t', quoting = csv.QUOTE_NONE, escapechar = '\t')<jupyter_output><empty_output><jupyter_text>Prepare the sentivity analysis for Stein tree parameters  W range<jupyter_code># Load prizes data
prizes_data = terminal_df
terminals = list(prizes_data['name'])
n_terminals = len(terminals)
# Load IREF interactome
interactome_file = "../Data/iRefIndex_v14_MIScore_interactome_C9.costs.allcaps.txt"
graph = oi.Graph(interactome_file)
# Distribution of cheapest path between any two terminals without penalty (g=-\infty)
network = graph.interactome_graph
shortest_dist_mat = np.zeros(shape = (n_terminals,n_terminals))
for ix_prot1 in tqdm.tqdm(np.arange(n_terminals)):
    time.sleep(0.01)
    for ix_prot2 in np.arange(ix_prot1+1, n_terminals, 1):
        shortest_dist_mat[ix_prot1,ix_prot2] = nx.dijkstra_path_length(network,
                                                                       source = terminals[ix_prot1],
                                                                       target = terminals[ix_prot2],
                                                                       weight = 'cost')
        shortest_dist_mat[ix_prot2,ix_prot2] = shortest_dist_mat[ix_prot1,ix_prot2]
# Plot heatmap of cheapest path
plt.matshow(np.transpose(shortest_dist_mat))
plt.xlabel('terminals')
plt.ylabel('terminals')
plt.colorbar()
# Plot histogram of cheapest path
plt.subplot(1,2,1)
flat_mat = shortest_dist_mat[np.triu_indices(n_terminals, k=1)]
plt.hist(flat_mat,40)
plt.xlabel('Cost of shortest path between two terminals')
plt.ylabel('Number of pairs of terminals')<jupyter_output><empty_output><jupyter_text>G range<jupyter_code># Create a dictionary containing edge costs under several choices of g
g_range = [0,1,2,3,4,5]
edge_costs_dict={'g = -Inf': graph.edge_costs}
for g in g_range:
    params = {"w": 0, "b": 1, "g": g, "edge_noise": 0, "dummy_mode": "terminals", "seed": 0, "skip_checks": False}
    graph._reset_hyperparameters(params)
    edge_costs_dict['g='+str(g)] = graph.costs
edge_costs_df = pd.DataFrame(data = edge_costs_dict)
# Boxplot of edge costs under several choices of g
plt.figure()
edge_costs_df.boxplot()
plt.yscale('log')
plt.ylabel('penalized edge cost')
plt.title('Boxplots of penalized edge costs in IREF for different values of g')
plt.show()<jupyter_output><empty_output><jupyter_text>B range<jupyter_code># For a range of g, compute the maximum penalized edge cost
g_range = [0,1,2,3,4,5]
max_penalized_edge_cost = pd.DataFrame(edge_costs_df.max())
max_penalized_edge_cost.columns=['max_penalized_edge_cost']
max_penalized_edge_cost<jupyter_output><empty_output><jupyter_text>Sensitivity analysis for Steiner tree parameters  Run sensitivity analysis<jupyter_code># Load set of virus interacting genes from Gordon et al.
virus_partners_file_name = "../Data/SARSCov_targets_df.tsv"
virus_interacting_genes = sensitivity.import_virus_partners(virus_partners_file_name)
# Parameters for sensitivity analysis
interactome_file_name = "../Data/iRefIndex_v14_MIScore_interactome_C9.costs.allcaps.txt"
prize_file_name = "../Save_iav/terminals_ppi_analysis.tsv"
# Graph hyperparameters
graph_params =  {
                "noise": 0.0, 
                "dummy_mode": "terminals", 
                "exclude_terminals": False, 
                "seed": 1,
                "pruning": 'strong',
                "verbosity_level": 0
                }
# Set sweeping parameters, i.e. configurations (w,b,g)
W_list = np.linspace(start = 0.2, stop = 2, num = 10)
B_list = np.array([5., 10., 15., 20., 25., 30., 35., 40., 45., 50.])
# Run sensitivity analysis
networks_dict = sensitivity.run_sensitivity_analysis(interactome_file_name, 
                                                     prize_file_name, 
                                                     graph_params, 
                                                     W_list, 
                                                     B_list, 
                                                     G=0)
# Save as pickle
with open("../Save_iav/networks_dict.pkl", "wb") as f: 
    pickle.dump(networks_dict, f)
# Add metatadata
networks_dict = pickle.load(open("../Save_iav/networks_dict.pkl", "rb"))
networks_dict = sensitivity.add_metadata(networks_dict, virus_interacting_genes)
# Make summary
networks_summary_df = sensitivity.make_summary(networks_dict, n_terminals, g=0)
networks_summary_df.head()<jupyter_output><empty_output><jupyter_text>Plot node stability heatmaps<jupyter_code># Stability of selected nodes
mat_allnodes = sensitivity.create_matrix_gene_overlap_between_networks(networks_summary_df, networks_dict)
plt.figure()
plt.matshow(mat_allnodes)
plt.xlabel('networks')
plt.ylabel('networks')
plt.title('Nodes stability')
plt.colorbar()
plt.show()
# Stability of selected terminals
mat_terminals = sensitivity.create_matrix_terminal_overlap_between_networks(networks_summary_df, networks_dict)
plt.figure()
plt.matshow(mat_terminals)
plt.xlabel('networks')
plt.ylabel('networks')
plt.title('Terminals stability')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Select robust parameters<jupyter_code># Select network corresponding to g=0, w=1.4 and b=40
index_selected = 67
paramstring_selected = networks_summary_df[networks_summary_df['index']==index_selected].index[0]
network_selected = networks_dict[paramstring_selected]
# Save selected network to file
oi.output_networkx_graph_as_interactive_html(network_selected, filename="../Save_iav/network_selected.html")
oi.output_networkx_graph_as_pickle(network_selected, filename= '../Save_iav/network_selected.pickle')
#oi.output_networkx_graph_as_graphml_for_cytoscape(robust_network, filename= '../Save_iav/network_selected.graphml')
# Characteristics of selected network
networks_summary_df[networks_summary_df['index']==index_selected]<jupyter_output><empty_output><jupyter_text>Add drug targets to selected network  Construct drug/target data<jupyter_code># Compute degree centrality for all nodes in the interactome
graph = oi.Graph(interactome_file_name)
centrality_dic = nx.degree_centrality(graph.interactome_graph)
# Load drug/target data from DrugCentral
drug_target_file_name = '../Data/drug.target.interaction.tsv'
drugcentral_df = drugs2.load_drug_target_data(drug_target_file_name, aff_cst_thresh=5)
drugcentral_df['degree_centrality'] = [centrality_dic[gene] if (gene in list(centrality_dic.keys())) else None for gene in drugcentral_df['gene']]
drugcentral_df.head()
# Construct table that gives the number of targets per drug
num_targets_df0 = drugcentral_df.groupby('drug', as_index=False)[['gene']].agg({
    ('num_targets','count'),
    ('num_terminal_targets',lambda gs: len(set(gs).intersection(set(terminals))))
})
max_centrality_df0 = drugcentral_df.groupby('drug', as_index=False)[['degree_centrality']].agg({
    ('max_degree_centrality',max)
})
num_targets_df = pd.DataFrame({'drug': num_targets_df0.index,
                               'num_targets': num_targets_df0['gene']['num_targets'],
                               'num_terminal_targets': num_targets_df0['gene']['num_terminal_targets'],
                               'max_target_centrality': max_centrality_df0['degree_centrality']['max_degree_centrality']}).reset_index(drop=True)
num_targets_df = num_targets_df.sort_values(by='num_targets', ascending=False, inplace=False)
num_targets_df.head(10)
# Histogram of number of targets per drug
plt.figure()
plt.hist(num_targets_df['num_targets'],50)
plt.yscale('log')
plt.xlabel('Number of targets')
plt.ylabel('Count')
plt.plot()
# Load L1000 drugs with correlations
embedded_drugs_file_name = '../Data/iav_correlations_autoencoder_space.txt'
bestdrugs_df = pd.read_csv(embedded_drugs_file_name, header=None)
bestdrugs_df.columns = ['drug','corr']
bestdrugs_df['drug'] = bestdrugs_df['drug'].str.strip("()' ").str.lower()
bestdrugs_df['corr'] = pd.to_numeric(bestdrugs_df['corr'].str.strip("() "))
bestdrugs_df.sort_values(by='corr', axis=0, ascending=False, inplace=True)

# Plot histogram of correlations
plt.figure()
plt.hist(bestdrugs_df['corr'],50)
plt.xlabel('anticorrelation')
plt.ylabel('count')
plt.show()

# Select top drugs (most anticorrelated with IAV signature)
bestdrugs_df = bestdrugs_df.iloc[:142] # so that we have the same number (142) of selected drugs as in the A549-SARS-Cov-2 analysis
print(bestdrugs_df.head())
print(bestdrugs_df.tail())
# Merge L1000 drugs with DrugCentral drug/target dataset
targets_and_drugs_df = drugcentral_df.merge(bestdrugs_df, on = 'drug', how = 'inner')
targets_and_drugs_df.head()<jupyter_output><empty_output><jupyter_text>Add drug target information to selected network<jupyter_code>network_selected = pickle.load(open('../Save_iav/network_selected.pickle', "rb"))
network_selected = drugs2.add_drug_info_to_selected_network(network_selected, targets_and_drugs_df)
# Save enriched network as pickle
oi.output_networkx_graph_as_pickle(network_selected, filename= '../Save_iav/network_selected_with_drug_info.pickle')
# Construct table of drug targets in the network
drug_targets_df = drugs2.drug_targets_in_selected_network(network_selected)
drug_targets_df.to_csv(r'../Save_iav/drug_targets_in_network.tsv', header=True, index=None, sep='\t', quoting = csv.QUOTE_NONE, escapechar = '\t')
drug_targets_df['affinity'] = pd.to_numeric(drug_targets_df['affinity'], errors='coerce')
drug_targets_df.dropna(subset=['affinity'], inplace=True)
drug_targets_df['terminal'] = drug_targets_df['name'].isin(terminals)
drug_targets_df
# Save drug/target dataframe to csv
drug_targets_df.to_csv(r'../Save_iav/final_drug_target_table.tsv', header=True, index=None, sep='\t', quoting = csv.QUOTE_NONE, escapechar = '\t')<jupyter_output><empty_output><jupyter_text>Compare to A549 cell type<jupyter_code># Drug/target dataframe for A549
drug_targets_df_a549_file_name = '../Save/final_drug_target_table.tsv'
drug_targets_df_a549 = pd.read_csv(drug_targets_df_a549_file_name, sep = '\t')
drug_targets_df_a549.head()
# Common gene targets
targets_iav = set(drug_targets_df['name'])
targets_a549 = set(drug_targets_df_a549['name'])

# plot Venn diagram IAV vs SARS-Cov-2
plt.figure()
venn2(subsets = [targets_iav,targets_a549],
      set_labels = ('IAV gene targets','SARS-Cov-2 gene targets'),
      set_colors=('r', 'g'), 
      alpha = 0.5)
#venn2_circles(subsets = [targets_iav,targets_a549])
plt.savefig('../Save_iav/venn2_iav_sars_genes.pdf', format='pdf')
plt.show()

print(sorted(targets_iav.intersection(targets_a549)))
print(sorted(targets_iav.difference(targets_a549)))
# Common drugs
drugs_iav = set(drug_targets_df['drug'])
drugs_a549 = set(drug_targets_df_a549['drug'])

# plot Venn diagram IAV vs SARS-Cov-2
plt.figure()
venn2(subsets = [drugs_iav,drugs_a549],
      set_labels = ('IAV drugs','SARS-Cov-2 drugs'),
      set_colors=('r', 'g'), 
      alpha = 0.5)
#venn2_circles(subsets = [targets_iav,targets_a549])
plt.savefig('../Save_iav/venn2_iav_sars_drugs.pdf', format='pdf')
plt.show()

a549_minus_iav = num_targets_df.loc[num_targets_df['drug'].isin(drugs_a549.difference(drugs_iav))]
iav_minus_a549 = num_targets_df.loc[num_targets_df['drug'].isin(drugs_iav.difference(drugs_a549))]
a549_inter_iav = num_targets_df.loc[num_targets_df['drug'].isin(drugs_a549.intersection(drugs_iav))]

print(sorted(set(a549_minus_iav['drug'])))
print(sorted(set(iav_minus_a549['drug'])))
print(sorted(set(a549_inter_iav['drug'])))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11902239709443099, 'fraction_numerical': 0.01679782082324455, 'mean_word_length': 5.4725759059745345, 'pattern_counts': {'":': 13, '<': 49, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '101569', 'license': 'MIT', 'n_tokens_mistral': 4953, 'n_tokens_neox': 4534, 'n_words': 795, 'path': 'Code_ppi/Code/SteinerTree_notebook-IAV.ipynb', 'repo_name': 'adamzabek/Drug-Repurposing', 'size': '13216'}