<jupyter_start><jupyter_text>Parse the results of new clixo, output relevant files for cytoscape analysis<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
import pandas as pd
import random


# latex rendering of text in graphs
import matplotlib as mpl

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

sns.set_style('white')

import sys
from IPython.display import display
% matplotlib inline
from ddot import Ontology<jupyter_output><empty_output><jupyter_text>Load the network genes<jupyter_code>ASD_CHD_df = pd.read_excel('../data/supplemental_tables_cell_systems_210416.xlsx',sheet_name='Table S4',skiprows=1)

print(len(ASD_CHD_df))
ASD_CHD_df.index=ASD_CHD_df['gene']
ASD_CHD_df.head()<jupyter_output>844<jupyter_text>Load PCnetDownload PCnet from https://ndexbio.org//network/f93f402c-86d4-11e7-a10d-0ac135e8bacf (UUID: f93f402c-86d4-11e7-a10d-0ac135e8bacf)<jupyter_code>G_pcnet = nx.read_gpickle('/Users/brinrosenthal/Documents/CCBB_tickets_data/PCnet/G_PCnet.gpickle')

print(len(G_pcnet.nodes()))
print(len(G_pcnet.edges()))
# make a subgraph from ASD_CHD_df genes
G_ASD_CHD = nx.subgraph(G_pcnet,ASD_CHD_df.index.tolist())
print(len(G_ASD_CHD.nodes()))
print(len(G_ASD_CHD.edges()))<jupyter_output>844
28559<jupyter_text>Load similarity matrix<jupyter_code>sim_rank_EL = pd.read_csv('new_clixo_files/sim_rank_cosine_EL.tsv',sep='\t',names=['node1','node2','sim'])
sim_rank_EL.head()<jupyter_output><empty_output><jupyter_text>Read in the ontology<jupyter_code>alpha='01'
beta='45'
filepath='new_clixo_files/ASDCHD.cosine..'+alpha+'..'+beta+'.ont'
ont = Ontology.from_table('new_clixo_files/ASDCHD_parsed_alpha.cosine.'+alpha+'_beta_.'+beta+'.txt',parent=0,child=1)
ont<jupyter_output><empty_output><jupyter_text>Add the fraction validated genes to ontology terms, write out for input to cytoscape<jupyter_code>term_sizes = pd.Series(ont.term_sizes)
term_sizes.index = pd.Series(ont.terms_index).index.tolist()
term_sizes=term_sizes.sort_values(ascending=False)
term_sizes.head()
# need to use ont_connected to get actual genes in term (trace ontology)

ont_connected = pd.DataFrame(ont.connected())
ont_connected.index=pd.Series(ont.genes_index).index.tolist()+pd.Series(ont.terms_index).index.tolist()
ont_connected.columns=pd.Series(ont.genes_index).index.tolist()+pd.Series(ont.terms_index).index.tolist()
ont_connected.head()
gene2term_full = ont_connected.iloc[0:len(ont.genes_index)].T.iloc[len(ont.genes_index):].T
gene2term_full.head()
# save out the pathways which have at least one validated xenopus gene
xenopus_validated_genes = ['SCN2A','KMT2A','PTPN11','KMT2D','KANSL1','KAT6A','MAPT']
xenopus_pathways = gene2term_full.loc[xenopus_validated_genes]
xenopus_pathways.sum().sort_values(ascending=False).to_csv('new_clixo_files/xenopus_validated_pathways_210510.tsv',sep='\t')
xenopus_pathways.sum().head()
# write out the gene2term matrix
# gene2term_full.to_csv('new_clixo_files/gene2term_full_'+alpha+'_'+beta+'.csv')<jupyter_output><empty_output><jupyter_text>Run enrichment on each cluster, after discussion with TreyEnrich each term separately, using gprofilerAlso run GO alignment, but this happens in a separate notebook (clixo_GO_alignment/ASD_CHD_new_clixo_cosine_GO_alignment_200428.ipynb)<jupyter_code>from gprofiler import GProfiler
gp = GProfiler("MyToolName/0.1")
# loop over terms in gene2term_full
gene2term_full.head()
import os

savedir_profile = 'gprofiler_clixo_ASD_CHD_cosine_alpha_'+alpha+'_beta_'+beta
if not os.path.isdir(savedir_profile):
    os.mkdir(savedir_profile)

cluster2top_pathway = pd.DataFrame(columns=['cluster','pathway','term_ID',
                                            'term_p-value','term_recall','term_precision','num_genes_in_cluster'])

for focal_cluster in gene2term_full.columns.tolist():
    print(focal_cluster)
    focal_genes = gene2term_full[gene2term_full[focal_cluster]].index.tolist()
    print(len(focal_genes))
    if len(focal_genes)>5: # only run enrichment if > 5 genes in cluster
        # ---- USE THIS FOR BACKGROUND OF all PCnet genes ----
        gp_results = pd.DataFrame(gp.gprofile(focal_genes,custom_bg = G_pcnet.nodes(),correction_method=gp.THR_FDR))
        
        # ---- USE THIS FOR BACKGROUND OF 844 ASDCHD genes ----
        # gp_results = pd.DataFrame(gp.gprofile(focal_genes,custom_bg = gene2term_full.index.tolist(),correction_method=gp.THR_FDR))
        
        if len(gp_results)>0: # make sure there are some significant enrichments

            gp_results.columns = ["query.number", "significant", "p.value", "term.size",
                                  "query.size", "overlap.size", "recall", "precision",
                                  "term.id", "domain", "subgraph.number", "term.name",
                                  "relative.depth", "intersection"]
            print(gp_results[['p.value','precision','recall','term.size','overlap.size','term.id','term.name']].head())

            # uncomment if we want to save the results
            writer = pd.ExcelWriter(savedir_profile+'/cluster_'+str(focal_cluster)+'.xlsx')
            gp_results.to_excel(writer)
            writer.save()

            cluster2top_pathway = cluster2top_pathway.append(pd.DataFrame({'cluster':[focal_cluster],
                                                                          'pathway':gp_results.loc[0]['term.name'],
                                                                           'term_ID':gp_results.loc[0]['term.id'],
                                                                          'term_p-value':gp_results.loc[0]['p.value'],
                                                                           'term_recall':gp_results.loc[0]['recall'],
                                                                           'term_precision':gp_results.loc[0]['precision'],
                                                                          'num_genes_in_cluster':[len(focal_genes)]}),
                                                            ignore_index=True)
            
        else:
            
            cluster2top_pathway = cluster2top_pathway.append(pd.DataFrame({'cluster':[focal_cluster],
                                                                      'pathway':'None',
                                                                       'term_ID':'None',
                                                                      'term_p-value':1.0,
                                                                        'term_recall':'None',
                                                                        'term_precision':'None',
                                                                      'num_genes_in_cluster':[len(focal_genes)]}),
                                                        ignore_index=True)
            
    else:
        cluster2top_pathway = cluster2top_pathway.append(pd.DataFrame({'cluster':[focal_cluster],
                                                                      'pathway':'None',
                                                                       'term_ID':'None',
                                                                      'p-value':1.0,
                                                                        'term_recall':'None',
                                                                        'term_precision':'None',
                                                                      'num_genes_in_cluster':[len(focal_genes)]}),
                                                        ignore_index=True)

#gp_results = pd.DataFrame(gp.gprofile(focal_genes,correction_method=gp.THR_FDR))
cluster2top_pathway.head()
cluster2top_pathway.index=cluster2top_pathway['cluster']
cluster2top_pathway.head()<jupyter_output><empty_output><jupyter_text>add top gprofiler term per cluster for further analysis/ visualization in cytoscape<jupyter_code># add top pathway to ont2numVal
ont2numVal=ont2numVal.join(cluster2top_pathway[['num_genes_in_cluster','term_p-value',
                                                'term_precision','term_recall','pathway','term_ID']])
ont2numVal.to_csv('new_clixo_files/ont2numVal_cosine_alpha'+alpha+'_beta'+beta+'.txt',sep='\t')
# load precomputed ont2numVal so we don't have to enrich each pathway again
# ont2numVal = pd.read_csv('new_clixo_files/ont2numVal_cosine_alpha'+alpha+'_beta'+beta+'.txt',sep='\t',index_col='term')
# ont2numVal.head()
ont2numVal.head()<jupyter_output><empty_output><jupyter_text>write out ontology to cytoscape compatible format<jupyter_code># used to use the full ontology, but now we use the tree form for plotting
# ont.to_table('../clixo_hierarchy/new_clixo_no_alignment/ASD_CHD_cosine_'+alpha+'_'+beta+'.txt')
# write out ontology in tree form (this is what is used in making the figures)
tree = ont.get_tree()
tree_df = pd.DataFrame(zip(*tree)).T
tree_df.columns=['child','parent']
tree_df.to_csv('new_clixo_files/tree_ASDCHD_cosine_'+alpha+'_'+beta+'.csv')
tree_df.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11545195885314073, 'fraction_numerical': 0.014664040271394178, 'mean_word_length': 2.5272095715939793, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 38, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 68, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '301726', 'license': 'MIT', 'n_tokens_mistral': 2936, 'n_tokens_neox': 2786, 'n_words': 402, 'path': 'NB7_clixo_analysis/NB7B_parse_new_clixo_results_210419.ipynb', 'repo_name': 'ucsd-ccbb/ASD_combined_network_analysis', 'size': '9138'}