<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Initial investigation of drug-gene networks
 Brin Rosenthal (sbrosenthal@ucsd.edu) 
 April 15, 2016
Prototype for tool to be added to Search
Goals
Step1:  Test drug_gene_heatprop module
This section runs the inferred drug heat propagation module from a list of seed genes, and returns a list of genes ranked by their 'heat', or proximity to the seed gene set.  These are the genes which we think will be most related to the seed genes.
For a more detailed, step by step description of the process, continue reading past this section.
Step2:  More detailed description of methods below...
Load the drug bank database, and create a network out of it
Network is bipartite with types
Step3:  What is this drug-gene graph like?
how sparse is it?
Are there genes/drugs that have many connections?
Step4:  But we probably want to focus on within-cluster interactions, instead of the whole graph
Download a sample cluster from geneli.st (Adrenocortical carcinoma cluster 250250)
Extract a subnetwork from total drug-gene network containing only the genes from this cluster, and associated drugs
Plot this subnetwork
Step5:  Above, we plot the drug-gene interaction network for our sample cluster
We're showing only the genes that have associated drugs, in the cluster
This is one option for exploring drug-gene interaction space, as the sparseness of drugs/genes per cluster allows for easy visualization
Another option... heat propagation
Run heat propagation from seed nodes on a sample cluster, to prioritize genes (and their associated drugs) similar to seed node set
Some questions to resolve
Step6:  First let's plot the focal cluster of interest (Adrenocortical carcinoma cluster 250250)
Step7:  Now we will convert the cluster correlation matrix back to network form
Step8:  Now let's look up the drugs associated with these genes to see if there are any good candidates
<ASSISTANT_TASK:>
Python Code:
# import some useful packages
import numpy as np
import matplotlib.pyplot as plt
import seaborn
import networkx as nx
import pandas as pd
import random
import json
# latex rendering of text in graphs
import matplotlib as mpl
mpl.rc('text', usetex = False)
mpl.rc('font', family = 'serif')
% matplotlib inline
Explanation: Initial investigation of drug-gene networks
 Brin Rosenthal (sbrosenthal@ucsd.edu) 
 April 15, 2016
Prototype for tool to be added to Search
Goals:
Input a gene list
Use the DrugBank database to suggest drugs related to genes in input list
Note: data files and code for this notebook may be found in the 'data' and 'source' directories
End of explanation
# load the module
import sys
sys.path.append('../source')
import drug_gene_heatprop
import imp
imp.reload(drug_gene_heatprop)
path_to_DB_file = '../drugbank.0.json.new' # set path to drug bank file
path_to_cluster_file = 'sample_matrix.csv' # set path to cluster file
seed_genes = ['LETM1','RPL3','GRK4','RWDD4A'] # set seed genes (must be in cluster)
gene_drug_df = drug_gene_heatprop.drug_gene_heatprop(seed_genes,path_to_DB_file,path_to_cluster_file,
                                                    plot_flag=True)
gene_drug_df.head(25)
Explanation: Test drug_gene_heatprop module
This section runs the inferred drug heat propagation module from a list of seed genes, and returns a list of genes ranked by their 'heat', or proximity to the seed gene set.  These are the genes which we think will be most related to the seed genes.
For a more detailed, step by step description of the process, continue reading past this section.
End of explanation
def load_DB_data(fname):
    '''
    Load and process the drug bank data
    
    '''
    with open(fname, 'r') as f:
        read_data = f.read()
    f.closed
    si = read_data.find('\'\n{\n\t"source":')
    sf = read_data.find('\ncurl')
    DBdict = dict()
    # fill in DBdict
    while si > 0:
        
        db_temp = json.loads(read_data[si+2:sf-2])
        DBdict[db_temp['drugbank_id']]=db_temp
        # update read_data
        read_data = read_data[sf+10:]
        
        si = read_data.find('\'\n{\n\t"source":')
        sf = read_data.find('\ncurl')
        
    return DBdict
DBdict = load_DB_data('/Users/brin/Documents/DrugBank/drugbank.0.json.new')
# make a network out of drug-gene interactions
DB_el = []
for d in DBdict.keys():
    node_list = DBdict[d]['node_list']
    for n in node_list:
        DB_el.append((DBdict[d]['drugbank_id'],n['name']))
    
    
G_DB = nx.Graph()
G_DB.add_edges_from(DB_el)
gene_nodes,drug_nodes = nx.bipartite.sets(G_DB)
gene_nodes = list(gene_nodes)
drug_nodes = list(drug_nodes)
Explanation: More detailed description of methods below...
Load the drug bank database, and create a network out of it
Network is bipartite with types:
Drugs
Genes which are acted on by each drug
End of explanation
print('--> there are '+str(len(gene_nodes)) + ' genes with ' + str(len(drug_nodes)) + ' corresponding drugs')
DB_degree = pd.Series(G_DB.degree())
DB_degree.sort(ascending=False)
plt.figure(figsize=(18,5))
plt.bar(np.arange(70),DB_degree.head(70),width=.5)
tmp = plt.xticks(np.arange(70)+.4,list(DB_degree.head(70).index),rotation=90,fontsize=11)
plt.xlim(1,71)
plt.ylim(0,200)
plt.grid('off')
plt.ylabel('number of connections (degree)',fontsize=16)
Explanation: What is this drug-gene graph like?
how sparse is it?
Are there genes/drugs that have many connections?
End of explanation
# load a sample cluster for network visualization
sample_genes = pd.read_csv('/Users/brin/Documents/DrugBank/sample_cluster.csv',header=None)
sample_genes = list(sample_genes[0])
# also include neighbor genes
neighbor_genes = [nx.neighbors(G_DB,x) for x in sample_genes if x in G_DB.nodes()]
neighbor_genes = [val for sublist in neighbor_genes for val in sublist]
sub_genes = []
sub_genes.extend(sample_genes)
sub_genes.extend(neighbor_genes)
G_DB_sample = nx.subgraph(G_DB,sub_genes)
drug_nodes = list(np.intersect1d(neighbor_genes,G_DB.nodes()))
gene_nodes = list(np.intersect1d(sample_genes,G_DB.nodes()))
# return label positions offset by dx
def calc_pos_labels(pos,dx=.03):
    # input node positions from nx.spring_layout()
    
    pos_labels = dict()
    for key in pos.keys():
        pos_labels[key] = np.array([pos[key][0]+dx,pos[key][1]+dx])
    
    return pos_labels
pos = nx.spring_layout(G_DB_sample,k=.27)
pos_labels = calc_pos_labels(pos)
plt.figure(figsize=(14,14))
nx.draw_networkx_nodes(G_DB_sample,pos=pos,nodelist = drug_nodes,node_shape='s',node_size=80,alpha=.7,label='drugs')
nx.draw_networkx_nodes(G_DB_sample,pos=pos,nodelist = gene_nodes,node_shape='o',node_size=80,node_color='blue',alpha=.7,label='genes')
nx.draw_networkx_edges(G_DB_sample,pos=pos,alpha=.5)
nx.draw_networkx_labels(G_DB_sample,pos=pos_labels,font_size=10)
plt.grid('off')
plt.legend(fontsize=12)
plt.title('Adrenocortical carcinoma cluster 250250',fontsize=16)
Explanation: But we probably want to focus on within-cluster interactions, instead of the whole graph
Download a sample cluster from geneli.st (Adrenocortical carcinoma cluster 250250)
Extract a subnetwork from total drug-gene network containing only the genes from this cluster, and associated drugs
Plot this subnetwork
End of explanation
sample_mat = pd.read_csv('/Users/brin/Documents/DrugBank/sample_matrix.csv',index_col=0)
print(sample_mat.head())
idx_to_node = dict(zip(range(len(sample_mat)),list(sample_mat.index)))
sample_mat = np.array(sample_mat)
sample_mat = sample_mat[::-1,0:-1] # reverse the indices for use in graph creation
Explanation: Above, we plot the drug-gene interaction network for our sample cluster
We're showing only the genes that have associated drugs, in the cluster
This is one option for exploring drug-gene interaction space, as the sparseness of drugs/genes per cluster allows for easy visualization
Another option... heat propagation
Run heat propagation from seed nodes on a sample cluster, to prioritize genes (and their associated drugs) similar to seed node set
Some questions to resolve:
- ### How should we handle negative edge weights?
- ### Should we return drugs associated with individual genes, or drugs most associated with total input gene list?
End of explanation
plt.figure(figsize=(7,7))
plt.matshow(sample_mat,cmap='bwr',vmin=-1,vmax=1,fignum=False)
plt.grid('off')
plt.title('Adrenocortical carcinoma cluster 250250',fontsize='16')
Explanation: First let's plot the focal cluster of interest (Adrenocortical carcinoma cluster 250250)
End of explanation
G_cluster = nx.Graph()
G_cluster = nx.from_numpy_matrix(np.abs(sample_mat))
G_cluster = nx.relabel_nodes(G_cluster,idx_to_node)
pos = nx.spring_layout(G_cluster,k=.4)
seed_genes = ['STIM2','USP46','FRYL','COQ2'] #['STIM2','USP46'] # input gene list here
plt.figure(figsize=(10,10))
nx.draw_networkx_nodes(G_cluster,pos=pos,node_size=20,alpha=.5,node_color='blue')
nx.draw_networkx_nodes(G_cluster,pos=pos,nodelist=seed_genes,node_size=50,alpha=.7,node_color='red',linewidths=2)
nx.draw_networkx_edges(G_cluster,pos=pos,alpha=.03)
plt.grid('off')
plt.title('Sample subnetwork: pre-heat propagation',fontsize=16)
Wprime = network_prop.normalized_adj_matrix(G_cluster,weighted=True)
Fnew = network_prop.network_propagation(G_cluster,Wprime,seed_genes)
plt.figure(figsize=(10,10))
nx.draw_networkx_edges(G_cluster,pos=pos,alpha=.03)
nx.draw_networkx_nodes(G_cluster,pos=pos,node_size=20,alpha=.8,node_color=Fnew[G_cluster.nodes()],cmap='jet',
                       vmin=0,vmax=.005)
nx.draw_networkx_nodes(G_cluster,pos=pos,nodelist=seed_genes,node_size=50,alpha=.7,node_color='red',linewidths=2)
plt.grid('off')
plt.title('Sample subnetwork: post-heat propagation',fontsize=16)
N = 50
Fnew.sort(ascending=False)
print('Top N hot genes: ')
Fnew.head(N)
# plot the hot subgraph in gene-gene space
G_cluster_sub = nx.subgraph(G_cluster,list(Fnew.head(N).index))
pos = nx.spring_layout(G_cluster_sub,k=.5)
plt.figure(figsize=(10,10))
nx.draw_networkx_nodes(G_cluster_sub,pos=pos,node_size=100,node_color=Fnew[G_cluster_sub.nodes()],cmap='jet',
                      vmin=0,vmax=.005)
nx.draw_networkx_edges(G_cluster_sub,pos=pos,alpha=.05)
pos_labels = calc_pos_labels(pos,dx=.05)
nx.draw_networkx_labels(G_cluster_sub,pos=pos_labels)
plt.grid('off')
plt.title('Sample cluster: hot subnetwork \n (genes most related to input list)', fontsize=16)
Explanation: Now we will convert the cluster correlation matrix back to network form
End of explanation
top_N_genes = list(Fnew.head(N).index)
top_N_genes = list(np.setdiff1d(top_N_genes,seed_genes)) # only keep non-seed genes
top_N_genes = Fnew[top_N_genes]
top_N_genes.sort(ascending=False)
top_N_genes = list(top_N_genes.index)
drug_candidates_list = seed_genes # build up a list of genes and drugs that may be related to input list
for g in top_N_genes:
    if g in G_DB.nodes():  # check if g is in drugbank graph
        drug_candidates_list.append(g)
        drug_neighs_temp = list(nx.neighbors(G_DB,g))
        drug_candidates_list.extend(drug_neighs_temp)
    
# make a subgraph of these drug/gene candidates
G_DB_sub = nx.subgraph(G_DB,drug_candidates_list)
# define drug_nodes and gene_nodes from the subgraph
drug_nodes = list(np.intersect1d(neighbor_genes,G_DB_sub.nodes()))
gene_nodes = list(np.intersect1d(sample_genes,G_DB_sub.nodes()))
plt.figure(figsize=(12,12))
pos = nx.spring_layout(G_DB_sub)
pos_labels = calc_pos_labels(pos,dx=.05)
nx.draw_networkx_nodes(G_DB_sub,pos=pos,nodelist=gene_nodes,node_size=100,alpha=.7,node_color='blue',label='genes')
nx.draw_networkx_nodes(G_DB_sub,pos=pos,nodelist=drug_nodes,node_size=100,alpha=.7,node_color='red',node_shape='s',label='drugs')
nx.draw_networkx_edges(G_DB_sub,pos=pos,alpha=.5)
nx.draw_networkx_labels(G_DB_sub,pos=pos_labels,font_color='black')
plt.grid('off')
ax_min = np.min(pos.values())-.3
ax_max = np.max(pos.values())+.3
plt.xlim(ax_min,ax_max)
plt.ylim(ax_min,ax_max)
plt.legend(fontsize=14)
#plt.axes().set_aspect('equal')
plt.title('Genes in hot subnetwork with associated drugs', fontsize=16)
Explanation: Now let's look up the drugs associated with these genes to see if there are any good candidates
End of explanation
<END_TASK>