<jupyter_start><jupyter_text>Requirements* download graph2vec code from  https://github.com/benedekrozemberczki/graph2vec* NCI1 and NCI109 datasets     Also set the path(s) in the cell below. On Windows, there could be issues with the graph2vec code. One way to fix this is to change “/“ to “\\\\” on lines 63 and 97 of graph2vec.py<jupyter_code>## pick one of the two datasets:

#datadir = '../Datasets/NCI1/'
datadir = '../Datasets/NCI109/'

## location of graph2vec python code
g2v = '~/graph2vec/'

## random state for reproducibility with train/test splits
RS = 99
import pandas as pd
import numpy as np
import igraph as ig
import json
import os
import matplotlib.pyplot as plt
import umap, umap.plot
%matplotlib inline
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier as knn 
from sklearn import metrics
from sklearn.metrics import roc_auc_score
from sklearn.ensemble import RandomForestClassifier as rfc 
from collections import Counter<jupyter_output><empty_output><jupyter_text>Prepare datasets for graph2vec* save each graph in a json file, where the filenames are 0-based integers* each file contains a list of edges and dictionary of node features<jupyter_code>## read data from NCI folder
edges = datadir+'DS_A.txt'
graph_id = datadir+'DS_graph_indicator.txt'
graph_label = datadir+'DS_graph_labels.txt'
node_label = datadir+'DS_node_labels.txt'

## read edges, build overall graph
X = np.array(pd.read_csv(edges,header=None))
E = [list(x) for x in X]
G = ig.Graph.TupleList(E, directed=True)

## add missing vertices (isolates)
vertices = set(G.vs['name'])

## read subgraph membership (1-based)
m = [int(x) for x in np.array(pd.read_csv(graph_id,header=None))]

## add difference
v = set(np.arange(1,len(m)+1))
diff = v.difference(vertices)
G.add_vertices(list(diff))

## for plotting
G.vs['size'] = 10
G.vs['color'] = 'darkgrey'
G.es['color'] = 'grey'
G.es['arrow_size'] = .33
G.vs['label_size'] = 8
## mapping nodes in order of name
idx = list(np.argsort(G.vs['name']))

## assign subgraph
for i in range(len(m)):
    G.vs[idx[i]]['graph'] = m[i]

## verify with graph label list
l = [G.vs[i]['graph'] for i in idx]
l == m
## read node labels
l = [int(x) for x in np.array(pd.read_csv(node_label,header=None))]
for i in range(len(l)):
    G.vs[idx[i]]['label'] = l[i]
G = G.as_undirected()

## read graph labels
gl = [int(x) for x in np.array(pd.read_csv(graph_label,header=None))]

## build the subgraphs and save json files in the data directory
## this can take several minutes

## we also count nodes and edges for EDA
vc = []
ec = []

for gp in np.arange(1,np.max(G.vs['graph'])+1):
    v = [v for v in G.vs if v['graph']==gp]
    sg = G.subgraph(v)
    vc.append(sg.vcount())
    ec.append(sg.ecount())
    sg_edges = [list(e.tuple) for e in sg.es]
    sg_features = {str(v.index):str(v['label']) for v in sg.vs}
    sg_json = {"edges":sg_edges,"features":sg_features}
    fn = datadir+str(gp-1)+'.json'
    with open(fn,'w') as fp:
        json.dump(sg_json,fp)
a = [vc[i] for i in range(len(vc)) if gl[i]==0]
b = [vc[i] for i in range(len(vc)) if gl[i]==1]
c = [ec[i] for i in range(len(ec)) if gl[i]==0]
d = [ec[i] for i in range(len(ec)) if gl[i]==1]

plt.subplots(1,2,figsize=(9,4))
plt.subplot(121)
plt.boxplot([a,b],labels=['0','1'],widths=.6, 
            flierprops = dict(marker='.', markerfacecolor='black', markersize=3,linestyle='none'),
            medianprops = dict(linestyle='-', linewidth=1.5, color='black'))
plt.ylabel('Count per graph',fontsize=14);
plt.xlabel('Label',fontsize=14)
plt.title('Number of nodes',fontsize=14);

plt.subplot(122)
plt.boxplot([c,d],labels=['0','1'],widths=.6, 
            flierprops = dict(marker='.', markerfacecolor='black', markersize=3,linestyle='none'),
            medianprops = dict(linestyle='-', linewidth=1.5, color='black'))
#plt.ylabel('Count per graph',fontsize=11)
plt.xlabel('Label',fontsize=14);
plt.title('Number of edges',fontsize=14);
#plt.savefig('nci_counts.eps');<jupyter_output><empty_output><jupyter_text>Run graph2vec -- 64 dimensions<jupyter_code>## run graph2vec 
cmd = 'python '+g2v+'graph2vec.py --input-path '+datadir+' --output-path '+datadir+\
'NCI64.csv --dimensions 64 --workers 1' 
x = os.system(cmd+' >/dev/null 2>&1')<jupyter_output><empty_output><jupyter_text>Classification in embedded space* we use 80% for training, 20% for testing* we use random forest classifiers<jupyter_code>## divide data into training and testing
D = pd.read_csv(datadir+'NCI64.csv')
D = np.array(D.drop(columns=['type']))
X_train, X_test, y_train, y_test = train_test_split(D, gl, test_size=0.2, random_state=RS)
## random forest classifier -- accuracy
rfc_mdl = rfc(n_estimators=100,criterion='entropy')
rfc_mdl.fit(X_train,y_train)
y_pred = rfc_mdl.predict(X_test)
print("Accuracy:",metrics.accuracy_score(y_test, y_pred))<jupyter_output><empty_output><jupyter_text>Re-run graph2vec -- 1024 dimensions<jupyter_code>## run graph2vec 
cmd = 'python '+g2v+'graph2vec.py --input-path '+datadir+' --output-path '+datadir+\
'NCI1024.csv --dimensions 1024 --workers 1' 
x = os.system(cmd+' >/dev/null 2>&1')<jupyter_output><empty_output><jupyter_text>Visualize embeddings * color w.r.t. graph label* we see a few outliers (what are those?), but in general, no clear class separation (but this is only 2D!)<jupyter_code>D = pd.read_csv(datadir+'NCI1024.csv')
D = np.array(D.drop(columns=['type']))
U = umap.UMAP(metric='cosine').fit(D)
umap.plot.points(U,labels=np.array(gl),theme='viridis',width=800,height=500);<jupyter_output><empty_output><jupyter_text>Classification in embedded space* we use 80% for training, 20% for testing* we use random forest classifiers* we also try mapping from 1024 to 64 dimensions* Project:  - divide training set into training and validation  - use this to select best model by trying several hyperparameters<jupyter_code>## random forest classifier
X_train, X_test, y_train, y_test = train_test_split(D, gl, test_size=0.2, random_state=RS)
rfc_mdl = rfc(n_estimators=100,criterion='entropy')
rfc_mdl.fit(X_train,y_train)
y_pred = rfc_mdl.predict(X_test)
print("Accuracy:",metrics.accuracy_score(y_test, y_pred))
## try in lower dimension -- same random_seed to get same split
U = umap.UMAP(n_components=64,metric='cosine').fit(D).embedding_
X_train, X_test, y_train, y_test = train_test_split(U, gl, test_size=0.2, random_state=RS)
rfc_mdl = rfc(n_estimators=100,criterion='entropy')
rfc_mdl.fit(X_train,y_train)
y_pred = rfc_mdl.predict(X_test)
print("Accuracy:",metrics.accuracy_score(y_test, y_pred))
## ROC curve
y_probs = rfc_mdl.predict_proba(X_test)[:, 1]
fpr, tpr, _ = metrics.roc_curve(y_test,  y_probs)
auc = metrics.roc_auc_score(y_test, y_probs)
plt.plot(fpr,tpr,label="ROC, auc="+str('%.3f' % auc),color='black')
plt.plot([0,1],[0,1],'--',label='Random',color='black')
plt.legend(loc=4,fontsize=14)
plt.xlabel('False Positive Rate',fontsize=14)
plt.ylabel('True Positive Rate',fontsize=14);
#plt.savefig('nci_64d_roc.eps');

## bootstrap CI
def bootstrap(y_test,y_probs,n_boot=1000):
    y_pred = np.array(y_probs)
    y_true = np.array(y_test)
    scores = []
    rng = np.random.RandomState(RS)
    for i in range(n_boot):
        indices = rng.randint(0, len(y_pred), len(y_pred))
        sc = roc_auc_score(y_true[indices], y_pred[indices])
        scores.append(sc)
    sorted_scores = np.array(scores)
    sorted_scores.sort()
    conf_lo = sorted_scores[int(0.05 * len(sorted_scores))]
    conf_up = sorted_scores[int(0.95 * len(sorted_scores))]
    return (conf_lo, conf_up)
bootstrap(y_test,y_probs)<jupyter_output><empty_output><jupyter_text>Subgraph features- we compute several statistics on the graphs - number of nodes, edges, density  - degree distribution  - number of components, size of the giant component - transitivity - degree assortativity - coreness distribution - node labels distribution - we compare those for the two classes of graphs- we perform binary classification using those features<jupyter_code>## example of a graph with label 0
gp = 800
v = [v for v in G.vs if v['graph']==gp]
sg = G.subgraph(v)
print('label:',gl[gp-1])
ig.plot(sg,bbox=(0,0,300,300))
#ig.plot(sg,target='nci_0.eps',bbox=(0,0,300,300))

## example of a graph with label 1
gp = 2300
v = [v for v in G.vs if v['graph']==gp]
sg = G.subgraph(v)
print('label:',gl[gp-1])
ig.plot(sg,bbox=(0,0,300,300))
#ig.plot(sg,target='nci_1.eps',bbox=(0,0,300,300))<jupyter_output><empty_output><jupyter_text>Build graph-based features fo each compound<jupyter_code>def degrees(sg,md):
    ctr = Counter(sg.degree())
    return [ctr[i] for i in range(md+1)]

def core(sg,mc):
    ctr = Counter(sg.coreness())
    return [ctr[i] for i in range(mc+1)]

def labels(sg,ml):
    ctr = Counter(sg.vs['label'])
    return [ctr[i+1] for i in range(ml)]

## Compute features for each graph
## this takes a few minutes as there are many graphs

L = []
md = np.max(G.degree())
mc = np.max(G.coreness())
ml = np.max(G.vs['label'])

for gp in np.arange(1,np.max(G.vs['graph'])+1):
    v = [v for v in G.vs if v['graph']==gp]
    sg = G.subgraph(v)
    ## node and edge counts, density
    x = [sg.vcount(),sg.ecount(),sg.ecount()/sg.vcount()]
    ## number of components, relative size of giant component
    x.extend([np.max(sg.clusters().membership)+1,sg.clusters().giant().vcount()/sg.vcount()])
    ## transitivity, assortativity
    x.extend([sg.transitivity_undirected(),sg.transitivity_avglocal_undirected()])
    ## assortativity
    a = sg.assortativity_degree()
    if np.isnan(a):
        a=0
    x.extend([a])
    ## degree distribution
    x.extend(degrees(sg,md))
    ## coreness distribution
    x.extend(core(sg,mc))
    ## node labels distribution
    x.extend(labels(sg,ml))
    L.append(x)
## store all features in a dataframe
col = ['nodes','edges','density','components','giant','trans','local trans','assort']
col.extend(['deg'+str(i) for i in np.arange(md+1)])
col.extend(['core'+str(i) for i in np.arange(mc+1)])
col.extend(['label'+str(i+1) for i in np.arange(ml)])
F = pd.DataFrame(L,columns=col)
F['label'] = gl

## EDA of the graphs
print(F.shape[0],'graphs')
print(sum(F['label']),'have label 1')
print('avg nodes:',np.mean(F['nodes']))
print('avg edges:',np.mean(F['edges']))<jupyter_output><empty_output><jupyter_text>Add those features are re-do binary classification<jupyter_code>## binary classification with all features
Fp = pd.concat([pd.DataFrame(U,columns=['emb_'+str(i) for i in np.arange(0,U.shape[1])]),F],axis=1)
## divide data into training and testing
A = np.array(Fp.drop(columns=['label']))
X_train, X_test, y_train, y_test = train_test_split(A, gl, test_size=0.2, random_state=RS)
features = list(Fp.columns) ## names
## random forest classifier
rfc_mdl = rfc(n_estimators=100,criterion='entropy')
rfc_mdl.fit(X_train,y_train)
y_pred = rfc_mdl.predict(X_test)
print("Accuracy:",metrics.accuracy_score(y_test, y_pred))
## ROC curve
y_probs = rfc_mdl.predict_proba(X_test)[:, 1]
fpr, tpr, _ = metrics.roc_curve(y_test,  y_probs)
auc = metrics.roc_auc_score(y_test, y_probs)
plt.plot(fpr,tpr,label="ROC, auc="+str('%.3f' % auc),color='black')
plt.plot([0,1],[0,1],'--',label='Random',color='black')
plt.legend(loc=4,fontsize=14)
plt.xlabel('False Positive Rate',fontsize=14)
plt.ylabel('True Positive Rate',fontsize=14);
#plt.savefig('nci_120d_roc.eps');

## bootstrap CI
bootstrap(y_test,y_probs)
## most important features and plot a few
x = np.argsort(rfc_mdl.feature_importances_)
top = [x[i] for i in np.arange(len(x)-1,-1,-1)]
[features[i] for i in top[:5]]

## manually pick a few features to show
f = ['deg3','emb_32']

plt.subplots(1,2,figsize=(9,4))

plt.subplot(121)
a = Fp[Fp['label']==0][f[0]]
b = Fp[Fp['label']==1][f[0]]
plt.boxplot([a,b],labels=['0','1'],widths=.6, 
            flierprops = dict(marker='.', markerfacecolor='black', markersize=3,linestyle='none'),
            medianprops = dict(linestyle='-', linewidth=1.5, color='black'))
plt.ylabel(f[0],fontsize=11)
plt.xlabel('Label',fontsize=11);

plt.subplot(122)
a = Fp[Fp['label']==0][f[1]]
b = Fp[Fp['label']==1][f[1]]
plt.boxplot([a,b],labels=['0','1'],widths=.6, 
            flierprops = dict(marker='.', markerfacecolor='black', markersize=3,linestyle='none'),
            medianprops = dict(linestyle='-', linewidth=1.5, color='black'))
plt.ylabel(f[1],fontsize=11)
plt.xlabel('Label',fontsize=11);
#plt.savefig('nci_features.eps');<jupyter_output><empty_output><jupyter_text>Unsupervised learning<jupyter_code>NCL = 10
from sklearn.cluster import KMeans
cl = KMeans(n_clusters=NCL).fit(U).labels_
K = pd.DataFrame(np.array([gl,cl,np.repeat(1,len(gl))]).transpose(),columns=['label=1','cluster','total'])
C = K.groupby(by='cluster').sum()
C['ratio'] = C['label=1']/C['total']
plt.figure(figsize=(4,3))
plt.plot(C['total'],C['ratio'],'.',color='black')
plt.xlabel('cluster size')
plt.ylabel('proportion with label 1');
#plt.savefig('kmeans_1.eps');
df = C.sort_values(by='ratio',ascending=False).round(decimals=3)
df
#print(df.to_latex(index=True))
## we typicalls find a small pure cluster of size 10
SmallClust = df.index[df['total'] == 10][0]
Um = umap.UMAP(metric='cosine').fit(D)
umap.plot.points(Um,labels=cl,theme='viridis',width=800,height=500);
## small, pure cluster:
CL = SmallClust
df = F.loc[[i for i in range(len(cl)) if cl[i]==CL]][['nodes','edges','assort','deg1','deg2','deg3','label']]
df
#print(df.to_latex(index=True))
## plot some of those -- should be similar
gp = df.index[0] 
v = [v for v in G.vs if v['graph']==gp+1]
sg = G.subgraph(v)
print('label:',gl[gp])
ig.plot(sg,bbox=(0,0,300,300))
## another one -- 
gp = df.index[1] 
v = [v for v in G.vs if v['graph']==gp+1]
sg = G.subgraph(v)
print('label:',gl[gp])
ig.plot(sg,bbox=(0,0,300,300))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17573646785662458, 'fraction_numerical': 0.02510520969380351, 'mean_word_length': 5.658454106280193, 'pattern_counts': {'":': 2, '<': 45, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '528037', 'license': 'MIT', 'n_tokens_mistral': 5543, 'n_tokens_neox': 5279, 'n_words': 1136, 'path': 'Python_Notebooks/Complementary_B.ipynb', 'repo_name': 'aidarahim/GraphMiningNotebooks', 'size': '13782'}