<jupyter_start><jupyter_text>Add some additional layers of information to the hierarchyRequires DDOT for ontology parsing https://github.com/michaelkyu/ddot/blob/master/examples/Tutorial.ipynb<jupyter_code># load required packages

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
import pandas as pd
import random

import ndex2
import getpass

from IPython.display import display

# latex rendering of text in graphs
import matplotlib as mpl
mpl.rc('text', usetex = False)
mpl.rc('font', family = 'serif')

from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Arial']

sns.set(font_scale=1.4)

sns.set_style('white')

sns.set_style("ticks", {"xtick.major.size": 15, "ytick.major.size": 15})
plt.rcParams['svg.fonttype'] = 'none'

import sys

# % matplotlib inline
# need ddot to parse the ontology
import ddot
from ddot import Ontology
# ddot only compatible with nx 1.11....
nx.__version__<jupyter_output><empty_output><jupyter_text>Load interactome (needs to be the same one used in previous analysis)<jupyter_code>interactome_uuid='275bd84e-3d18-11e8-a935-0ac135e8bacf' # for STRING high confidence
ndex_server='public.ndexbio.org'
ndex_user=None
ndex_password=<PASSWORD>
G_PC = ndex2.create_nice_cx_from_server(
            ndex_server, 
            username=ndex_user, 
            password=<PASSWORD>, 
            uuid=interactome_uuid
        ).to_networkx()

# relabel nodes.... this is a nx 1.11 thing...
# G_PC_data = pd.DataFrame.from_dict(dict(G_PC.nodes(data=True)), orient='index')
# G_PC = nx.relabel_nodes(G_PC,dict(zip(G_PC_data.index.tolist(),G_PC_data['name'].tolist())))
# print(G_PC.nodes()[0:5])
# nodes = G_PC.nodes()

# # print out interactome num nodes and edges for diagnostic purposes
# print('number of nodes:')
# print(len(G_PC.nodes()))
# print('\nnumber of edges:')
# print(len(G_PC.edges()))

# ------ nx 2.1 ------
nodes = G_PC.nodes()

# print out interactome num nodes and edges for diagnostic purposes
print('number of nodes:')
print(len(G_PC.nodes()))
print('\nnumber of edges:')
print(len(G_PC.edges()))<jupyter_output>number of nodes:
17185

number of edges:
420534<jupyter_text>Load ASD rare-common hierarchy from ndex<jupyter_code>ASD_uuid='d27ca9ef-bce5-11eb-9e72-0ac135e8bacf'
ndex_server='public.ndexbio.org'
ndex_user='<EMAIL>'
ndex_password=<PASSWORD>('NDEx password: ')
G_ASD = ndex2.create_nice_cx_from_server(
            ndex_server, 
            username=ndex_user, 
            password=<PASSWORD>, 
            uuid=ASD_uuid
        ).to_networkx()
nodes = G_ASD.nodes()

# print out interactome num nodes and edges for diagnostic purposes
print('number of nodes:')
print(len(G_ASD.nodes()))
print('\nnumber of edges:')
print(len(G_ASD.edges()))
G_ASD.nodes()
# add node attributes to dataframe for easier access
ASD_clust_df = pd.DataFrame.from_dict(dict(G_ASD.nodes(data=True)), orient='index')
# relabel nodes in G_ASD
# G_ASD=nx.relabel_nodes(G_ASD,dict(ASD_clust_df['CD_CommunityName']))
ASD_clust_df['system_ID']=ASD_clust_df.index.tolist()
# some columns are not the right type
ASD_clust_df['CD_MemberList_Size']=[int(x) for x in ASD_clust_df['CD_MemberList_Size'].tolist()]
ASD_clust_df['HiDeF_persistence']=[int(x) for x in ASD_clust_df['HiDeF_persistence'].tolist()]
ASD_clust_df.head()
G_ASD.nodes()<jupyter_output><empty_output><jupyter_text>Add some annotation the the systems - Fraction rare/common/neither seeds Do we need to pull down the other network for this?<jupyter_code>ASD_clust_df.head()<jupyter_output><empty_output><jupyter_text>System validation with mammalian phenotype ontology + mouse KO dataParse the ontology, data from http://www.informatics.jax.org/vocab/mp_ontology<jupyter_code>from scipy.stats import hypergeom
from statsmodels.stats import contingency_tables

import requests

# find human orthologs of mouse genes
import mygene
mg = mygene.MyGeneInfo()
# download MGI phenotype data
url = 'http://www.informatics.jax.org/downloads/reports/MGI_PhenoGenoMP.rpt'
r = requests.get(url,allow_redirects=True)
open('MGI_PhenoGenoMP.rpt','wb').write(r.content)
mgi_df = pd.read_csv('MGI_PhenoGenoMP.rpt',sep='\t',
                    names=['MGI_Allele_Accession_ID','Allele symbol','involves','MP','PMID','MGI_marker_accession_ID'])
# extract gene names
gene_name = [a.split('<')[0] for a in mgi_df['Allele symbol'].tolist()]
mgi_df['gene_name']=gene_name
mgi_df.index=mgi_df['gene_name']
mgi_df.head()
mouse_genes = list(np.unique(mgi_df['gene_name']))
# mouse_genes = [g for g in mouse_genes if (not g.startswith('+'))]
mg_mapped = mg.querymany(mouse_genes,as_dataframe=True,species=['mouse','human'],scopes='symbol',fields='symbol')
mg_mapped.head()
# drop genes with no human ortholog
print(len(mg_mapped))
mg_mapped = mg_mapped.dropna(subset=['symbol'])
print(len(mg_mapped))
# drop duplicates
mg_mapped = mg_mapped[~mg_mapped.index.duplicated(keep='first')]
print(len(mg_mapped))
mg_mapped.head()

mgi_df['human_ortholog']=mgi_df['gene_name'].map(dict(mg_mapped['symbol']))
mgi_df.head()
url = 'http://www.informatics.jax.org/downloads/reports/MPheno_OBO.ontology'
r = requests.get(url,allow_redirects=True)
open('MPheno_OBO.ontology','wb').write(r.content)
ddot.parse_obo('MPheno_OBO.ontology',
               'parsed_mp.txt',
              'id2name_mp.txt',
              'id2namespace_mp.txt',
              'altID_mp.txt')
MP2desc = pd.read_csv('id2name_mp.txt',sep='\t',
                      names=['MP','description'],index_col='MP')

MP2desc=MP2desc.loc[MP2desc.index.dropna()] # drop NAN from index
print(len(MP2desc))


MP2desc.head()

hierarchy = pd.read_table('parsed_mp.txt',
                          sep='\t',
                          header=None,
                          names=['Parent', 'Child', 'Relation', 'Namespace'])

hierarchy.head()
MPO = Ontology.from_table(
    table=hierarchy,
    parent='Parent',
    child='Child',
    add_root_name='MP:00SUPER',
    ignore_orphan_terms=True)
#MPO.clear_node_attr()
#MPO.clear_edge_attr()

# add description to node attribute
terms_keep = list(np.unique(hierarchy['Parent'].tolist()+hierarchy['Child'].tolist()))
MPO.node_attr=MP2desc.loc[terms_keep]

MPO
print(len(MPO.node_attr))
# find terms related to brain
MP_focal_list = []
for t in MPO.node_attr.index.tolist():
    descr_temp = MPO.node_attr.loc[t]['description']
    if descr_temp.find('nervous')>-1:
        print(descr_temp)
        MP_focal_list.append(t)
    elif descr_temp.find('neuron')>-1:
        print(descr_temp)
        MP_focal_list.append(t)
    elif descr_temp.find('synapt')>-1:
        print(descr_temp)
        MP_focal_list.append(t)
        
print(len(MP_focal_list))
# test for enrichment in root node
OR_p_list,OR_CI_list,log_OR_list = [],[],[]
num_genes_in_term_list=[]

MP_keep_list = []

# root node is the largest node
root_node = ASD_clust_df['CD_MemberList_Size'].sort_values(ascending=False).head(1).index.tolist()[0]
# add a negative control phenotype: abnormal skeleton physiology: MP:0001533
for MP_focal in MP_focal_list+['MP:0001533']:
    MP_desc_focal = dict(MP2desc['description'])[MP_focal]

    # focus the hierarchy on one branch, and look up all terms within that branch
    if len(MPO.parent_2_child[MP_focal])>0:
        MPO_focal = MPO.focus(MP_focal,verbose=False)
        focal_terms = MPO_focal.terms
    else: # if the term has no children, just look at that term
        focal_terms=[MP_focal]


    # check enrichment in root node
    focal_genes = ASD_clust_df['CD_MemberList'].loc[root_node].split(' ')
    

    mgi_temp = mgi_df[mgi_df['MP'].isin(focal_terms)]
    mgi_temp = mgi_temp.dropna(subset=['human_ortholog'])
    mgi_genes = list(np.unique(mgi_temp['human_ortholog']))
    mgi_genes = [g.upper() for g in mgi_genes]
    mgi_genes = list(np.intersect1d(mgi_genes,G_PC.nodes()))
    
    if (len(mgi_genes)>10)&(len(mgi_genes)<2000): # only test if there are at least 10 genes, and fewer than 2000 genes
        print('\n'+MP_desc_focal)
        print('number of genes in root node = '+str(len(focal_genes)))
        print('number of genes in focal MPO term = '+str(len(mgi_genes)))

        q00 = len(np.intersect1d(mgi_genes,focal_genes))
        print('number overlapping genes = '+str(q00))
        q01 = len(mgi_genes)-q00

        q10 = len(focal_genes)-q00
        q11 = len(G_PC.nodes())-q00-q01-q10

        table_temp = [[q00,q01],[q10,q11]]
        print(table_temp)

        CT= contingency_tables.Table2x2(table_temp)
        OR_p_temp = CT.log_oddsratio_pvalue()
        OR_CI_temp = CT.log_oddsratio_confint()
        log_OR_temp = CT.log_oddsratio
        print(OR_p_temp)
        print(OR_CI_temp)
        print(log_OR_temp)
        print(len(mgi_genes))

        OR_p_list.append(OR_p_temp)
        OR_CI_list.append(OR_CI_temp)
        log_OR_list.append(log_OR_temp)
        num_genes_in_term_list.append(len(mgi_genes))
        
        MP_keep_list.append(MP_focal)

        print(hypergeom.sf(q00,len(G_PC.nodes()),len(focal_genes),len(mgi_genes)))
    
OR_CI_lower, OR_CI_upper = zip(*OR_CI_list)

root_KO_df = pd.DataFrame({'OR_p':OR_p_list,'log_OR':log_OR_list,
                           'log_OR_CI_lower':OR_CI_lower,'log_OR_CI_upper':OR_CI_upper,
                          'num_genes_in_term':num_genes_in_term_list},
                          index=MP_keep_list)
root_KO_df['MP_description']=root_KO_df.index.map(dict(MP2desc['description']))

root_KO_df=root_KO_df.sort_values('OR_p')
# root_KO_df=root_KO_df.sort_values('log_OR',ascending=False)
root_KO_df.head(10)
plt.figure(figsize=(3,6))
terms_plot = root_KO_df.head(10).index.tolist()+['MP:0001533'] # top 10 + positive controls

plt.errorbar(root_KO_df.loc[terms_plot]['log_OR'],np.arange(len(terms_plot)),
            xerr=[np.subtract(root_KO_df.loc[terms_plot]['log_OR'],root_KO_df.loc[terms_plot]['log_OR_CI_lower']),
                   np.subtract(root_KO_df.loc[terms_plot]['log_OR_CI_upper'],root_KO_df.loc[terms_plot]['log_OR'])],color='k',fmt='o')
# plt.errorbar(np.arange(len(MP_focal_list)),np.exp(log_OR_list),
#              yerr=[np.subtract(np.exp(log_OR_list),np.exp(OR_CI_lower)),
#                    np.subtract(np.exp(OR_CI_upper),np.exp(log_OR_list))],color='k',fmt='o')
# plt.plot(np.arange(len(MP_focal_list)),np.exp(log_OR_list),'o',color='k')

plt.scatter(root_KO_df.loc[terms_plot]['log_OR'],np.arange(len(terms_plot)),s=root_KO_df.loc[terms_plot]['num_genes_in_term'],alpha=.5)

plt.yticks(np.arange(len(terms_plot)),root_KO_df.loc[terms_plot]['MP_description'])
plt.xticks([0,1,2])
plt.xlabel('log(OR) +- 95% CI')

plt.plot([0,0],[-.5,10.5],'--',color='gray')
plt.ylim([-.5,10.5])

plt.gca().invert_yaxis()
# tmp=plt.xticks(np.arange(len(MP_focal_list)),MP2desc['description'].loc[MP_focal_list].tolist(),rotation='vertical')

# plt.ylabel('OR +- 95% CI')

plt.savefig('figures/ASD_rare_common_mouse_KO_STRING_210608.png',dpi=300,bbox_inches='tight')

np.exp(np.subtract(np.exp(root_KO_df.loc[terms_plot]['log_OR']),np.exp(root_KO_df.loc[terms_plot]['log_OR_CI_lower'])))<jupyter_output><empty_output><jupyter_text>Test for enrichment in all systems<jupyter_code>MP_focal_top = root_KO_df.head(10).index.tolist() # record the top 10

MP_full_results_df=pd.DataFrame(index=ASD_clust_df.index.tolist())

for MP_focal in MP_focal_top:
    MP_desc_focal = dict(MP2desc['description'])[MP_focal]
    print(MP_desc_focal)

    # focus the hierarchy on one branch, and look up all terms within that branch
    if len(MPO.parent_2_child[MP_focal])>0:
        MPO_focal = MPO.focus(MP_focal)
        focal_terms = MPO_focal.terms
    else: # if the term has no children, just look at that term
        focal_terms=MP_focal


    hyper_p_list = []
    num_genes_list = []
    genes_id_list = []

    OR_p_list,OR_CI_list,log_OR_list=[],[],[]
    for focal_cluster in ASD_clust_df.index.tolist():
        #print(focal_cluster)
        mFocal_genes = ASD_clust_df['CD_MemberList'].loc[focal_cluster].split(' ')


        from scipy.stats import hypergeom
        M=len(list(G_PC.nodes())) # only keep genes in PCnet
        # Look up all entries matching focal_terms, and mFocal_genes
        mgi_temp = mgi_df[mgi_df['MP'].isin(focal_terms)]
        mgi_temp = mgi_temp.dropna(subset=['human_ortholog'])
        mgi_genes = list(np.unique(mgi_temp['human_ortholog']))
        new_index=[g.upper() for g in mgi_temp.index.tolist()] 
        mgi_temp.index=new_index

        N=len(np.intersect1d(list(np.unique(mgi_temp.index.tolist())),list(G_PC.nodes()))) # only keep genes in PCnet

        if len(np.intersect1d(mFocal_genes,mgi_temp.index.tolist()))>0:
            mgi_genes = mgi_temp.index.tolist()
            mgi_genes = list(np.intersect1d(mgi_genes,list(G_PC.nodes())))
            x = len(np.intersect1d(mFocal_genes,mgi_genes))
            #print(x)
            n=len(mFocal_genes)

            #mgi_temp['MP_description']=mgi_temp['MP'].map(dict(MP2desc['description']))



            #print(hypergeom.sf(x,M,n,N))
            hyper_p_list.append(hypergeom.sf(x,M,n,N))
            num_genes_list.append(x)
            genes_id_list.append(' '.join(list(np.intersect1d(mFocal_genes,mgi_genes))))

            q00 = len(np.intersect1d(mgi_genes,mFocal_genes))
            q01 = len(mgi_genes)-q00

            q10 = len(mFocal_genes)-q00
            q11 = len(list(G_PC.nodes()))-q00-q01-q10

            table_temp = [[q00,q01],[q10,q11]]

            CT= contingency_tables.Table2x2(table_temp)
            OR_p_temp = CT.log_oddsratio_pvalue()
            OR_CI_temp = CT.log_oddsratio_confint()
            log_OR_temp = CT.log_oddsratio
            #print(OR_p_temp)


            OR_p_list.append(OR_p_temp)
            OR_CI_list.append(OR_CI_temp)
            log_OR_list.append(log_OR_temp)


        else:
            hyper_p_list.append(1)
            num_genes_list.append(0)
            genes_id_list.append('')

            OR_p_list.append(1)
            OR_CI_list.append(0)
            log_OR_list.append(0)



    MP_focal_df = pd.DataFrame({MP_desc_focal+':-log(hyper_p)':-np.log10(hyper_p_list),
                                MP_desc_focal+':-log(OR_p)':-np.log10(OR_p_list),
                                MP_desc_focal+':log_OR':log_OR_list,
                                MP_desc_focal+':num_genes':num_genes_list,
                                MP_desc_focal+':gene_ids':genes_id_list},index=ASD_clust_df.index.tolist())

    if MP_desc_focal+':-log(hyper_p)' not in MP_full_results_df.columns.tolist():
        MP_full_results_df=MP_full_results_df.join(MP_focal_df)

MP_full_results_df.sort_values('abnormal neuron differentiation:num_genes',ascending=False).head()

MP_full_results_df.to_csv('tables/ASD_MP_full_results_210608.txt',sep='\t')
np.power(10,.05)
np.exp(3)<jupyter_output><empty_output><jupyter_text>Update hierarchy on ndexNEED SOPHIE'S HELP WITH THIS<jupyter_code># add new fields as node attributes to hierarchy
MP_full_results_dict = MP_full_results_df.to_dict()

node['n']
MP_full_results_df.head()
G_ASD_cx = ndex2.create_nice_cx_from_networkx(G_ASD)
# add node attributes to the nice_cx
# for key in MP_full_results_dict.keys():
#     print(key)
#     # set the datatypes
#     if key.find('gene_ids')>-1:
#         data_type='string'
#     else:
#         data_type='double'
#     nx.set_node_attributes(G_ASD_cx,key,MP_full_results_dict[key],type=data_type)
    

for node_id, node in G_ASD_cx.get_nodes():
    data = MP_full_results_df.loc[node['n']]
    print(data)
    for row, value in data.items():
        if row.find('gene_ids')>-1:
            data_type = "string"
            value=str(value)
        else:
            data_type = "double"
            value = str(value) # nice cx can only accept strings as values...
            if value=='inf': # check if inf, set to -1 if so
                value='-1'
            
        #print(value)
        #print(type(value))
        G_ASD_cx.set_node_attribute(node_id, row, value, type=data_type)
G_ASD_cx.nodeAttributes
G_ASD_cx.set_opaque_aspect('cartesianLayout',{})
#Upload to NDEx
SERVER = input('NDEx server (probably ndexbio.org): ')
USERNAME = input('NDEx user name: ')
PASSWORD = <PASSWORD>('NDEx password: ')
network_uuid = G_ASD_cx.upload_to(SERVER, USERNAME, PASSWORD)
network_uuid
# update the network on ndex
ndex_password=<PASSWORD>('NDEx password: ') # enter password again (not sure if necessary)
G_ASD_cx.update_to(ASD_uuid,ndex_server,ndex_user,ndex_password)
G_ASD_cx.update_to?<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13683382735466698, 'fraction_numerical': 0.015486056499909261, 'mean_word_length': 4.077395577395578, 'pattern_counts': {'":': 2, '<': 37, '<?xml version=': 0, '>': 44, 'https://': 1, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '589760', 'license': 'MIT', 'n_tokens_mistral': 6562, 'n_tokens_neox': 6199, 'n_words': 926, 'path': 'example_notebooks/.ipynb_checkpoints/ASD_rare_common_hierarchy_annotation-checkpoint.ipynb', 'repo_name': 'Ceofy/NetColoc', 'size': '16566'}