<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Setup1.1&nbsp;&nbsp;Load Dataset2&nbsp;&nbsp;Regulon enrichments3&nbsp;&nbsp;Functional Enrichments4&nbsp;&nbsp;Combine Annotations<jupyter_code>import sys, os
sys.path.append('../../scripts/')
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
import pandas as pd
from sklearn.metrics import r2_score
from itertools import combinations
import numpy as np
from scipy import stats
from enrichments import *
from core import *<jupyter_output><empty_output><jupyter_text>Setup  Load Dataset<jupyter_code>datasets = sorted([x for x in os.listdir(os.path.join(DATA_DIR,'iModulons/'))
            if '.' not in x])
# Thresholds were obtained from 2_identify_thresholds
cutoffs = {'MA-1': 550,
 'MA-2': 600,
 'MA-3': 350,
 'RNAseq-1': 700,
 'RNAseq-2': 300,
 'combined': 400}
class IcaData():
    def __init__(self,M,A,X,metadata,cutoff):
        self.M = pd.read_csv(M,index_col=0)
        self.M.columns = self.M.columns.astype(int)
        self.A = pd.read_csv(A,index_col=0)
        self.A.index = self.A.index.astype(int)
        self.X = pd.read_csv(X,index_col=0)
        self.metadata = pd.read_csv(metadata,index_col=0)
        self.thresholds = {k:self.get_threshold(k,cutoff) for k in self.M.columns}
    
    def show_enriched(self,k):
        gene_table = gene_info.copy()
        gene_table = gene_table.reindex(self.M.index)
        in_imod = abs(self.M[k]) > self.thresholds[k]
        weights = self.M.loc[in_imod,k]
        weights.name = 'weight'
        rows = gene_table.loc[in_imod]
        final_rows = pd.concat([weights,rows],axis=1)
        return final_rows.sort_values('weight')
    
    def get_threshold(self,k,cutoff):
        i=0
        genes = abs(self.M[k]).sort_values()
        k2,p = stats.normaltest(self.M[k])
        while k2 > cutoff:
            i -= 1
            k2,p = stats.normaltest(self.M.loc[genes.index[:i],k])
        imod_genes = genes.iloc[i:]
        return np.mean([genes.iloc[i], genes.iloc[i - 1]])

def load(dataset):
    # Define directories
    ds_dir = os.path.join(DATA_DIR,'iModulons',dataset)
    
    # Define files
    X_file = os.path.join(DATA_DIR,'processed_data',dataset+'.csv')
    M_file = os.path.join(ds_dir,'M.csv')
    A_file = os.path.join(ds_dir,'A.csv')
    metadata_file = os.path.join(DATA_DIR,'metadata',dataset+'_metadata.csv')
    
    return IcaData(M_file,A_file,X_file,metadata_file,cutoffs[dataset])
# Load datasets
objs = {}
for ds in tqdm(datasets):
    objs[ds] = load(ds)
summary_table = pd.DataFrame(dtype=int)
for ds in datasets:
    summary_table.loc[ds,'Samples'] = len(objs[ds].metadata)
    summary_table.loc[ds,'Conditions'] = len(objs[ds].metadata[['project_id','condition_id']].drop_duplicates())
    summary_table.loc[ds,'I-modulons'] = len(objs[ds].A.index)
    
    r2_list = []
    for name,group in objs[ds].metadata.groupby(['project_id','condition_id']):
        for i1,i2 in combinations(group.index,2):
            r2_list.append(r2_score(objs[ds].X[i1],objs[ds].X[i2]))
    
    summary_table.loc[ds,'Mean R2'] = np.mean(r2_list)
summary_table<jupyter_output><empty_output><jupyter_text>Regulon enrichments<jupyter_code>import warnings
warnings.simplefilter('ignore',category=FutureWarning)
def get_enrichments(ica_data):
    all_genes = set(ica_data.M.index)
    list2struct = []
    
    # Get TF enrichments for each component
    for k in tqdm(ica_data.M.columns):
        genes = set(ica_data.show_enriched(k).index)
        df = compute_enrichments(genes,all_genes,trn,fdr=1e-5)
        df['component'] = k
        df.index.name = 'description'
        df = df.reset_index(drop=False)
        list2struct.append(df)
        
    DF_enrich = pd.concat(list2struct).reset_index(drop=True)
    
    return DF_enrich
list2struct = []
for ds in datasets:
    df_enrich = get_enrichments(objs[ds])
    df_enrich['dataset'] = ds
    list2struct.append(df_enrich)

DF_all_enrich = pd.concat(list2struct).reset_index(drop=True)
DF_reg_enrich = DF_all_enrich.sort_values(['dataset','component',
                                            'pvalue','precision']).reset_index(drop=True)
DF_reg_enrich['type'] = 'regulatory'
col_order = ['dataset','component','type','description','pvalue','qvalue','precision','recall','f1score','TP']
DF_reg_enrich = DF_reg_enrich[col_order]
DF_reg_enrich.head()
DF_reg_enrich.drop_duplicates(['dataset','component']).to_csv(os.path.join(DATA_DIR,'iModulons','1TF_enrichments.csv'))
DF_reg_enrich.to_csv(os.path.join(DATA_DIR,'iModulons','all_reg_enrichments.csv'))<jupyter_output><empty_output><jupyter_text>Functional Enrichments<jupyter_code>DF_GO = pd.read_csv(os.path.join(GENE_DIR,'DF_GO.csv'),index_col=0)
go_dict = {}
for name,group in DF_GO.groupby('go_name'):
    genes = set(group.bnumber)
    go_dict[name] = genes
def get_go_enrichments(ica_data):
    all_genes = set(ica_data.M.index)
    enrich_list = []
    for k in tqdm(ica_data.M.columns):
        ic_genes = set(ica_data.show_enriched(k).index)
        
        list2struct = []
        for go,go_genes in go_dict.items():
            ((tp,fp),(fn,tn)) = contingency(go_genes,ic_genes,all_genes)
            odds,pval = stats.fisher_exact(((tp,fp),(fn,tn)),alternative='greater')
            if len(ic_genes) > 0:
                list2struct.append([go,k,tp,tp/len(ic_genes),tp/len(go_genes),pval])
            else:
                list2struct.append([go,k,tp,0,0,pval])
        df_pvals = pd.DataFrame(list2struct,columns = ['description','component','TP',
                                                       'precision','recall','pvalue'])
        enrich_list.append(FDR(df_pvals,.01))
    go_enrichments = pd.concat(enrich_list).reset_index(drop=True)
    go_enrichments.component = go_enrichments.component.astype(int)
    return go_enrichments
list2struct = []

for ds in datasets:
    df = get_go_enrichments(objs[ds])
    df['dataset'] = ds
    list2struct.append(df)
    
DF_all_go = pd.concat(list2struct).reset_index(drop=True)
DF_go_enrich = DF_all_go.sort_values(['dataset','component',
                                      'pvalue'])
DF_go_enrich['f1score'] = 2*DF_go_enrich.precision*DF_go_enrich.recall/(DF_go_enrich.precision+DF_go_enrich.recall)
DF_go_enrich['type'] = 'functional'
col_order = ['dataset','component','type','description','pvalue','qvalue','precision','recall','f1score','TP']
DF_go_enrich = DF_go_enrich[col_order]
DF_go_enrich.head()
DF_go_enrich.drop_duplicates(['dataset','component']).to_csv(os.path.join(DATA_DIR,'iModulons','GO_enrichments.csv'))
DF_go_enrich.to_csv(os.path.join(DATA_DIR,'iModulons','all_GO_enrichments.csv'))<jupyter_output><empty_output><jupyter_text>Combine Annotations<jupyter_code>DF_reg_enrich = pd.read_csv(os.path.join(DATA_DIR,'iModulons','1TF_enrichments.csv'),index_col=0)
DF_go_enrich = pd.read_csv(os.path.join(DATA_DIR,'iModulons','GO_enrichments.csv'),index_col=0)
# Combine regulatory and GO enrichments
DF_categories = pd.concat([DF_reg_enrich,DF_go_enrich]).sort_values(['dataset','component','qvalue'])

# Add names
DF_categories['name'] = DF_categories.dataset.str.cat(DF_categories.component.astype(str),sep='_')

# Filter out GO enrichments when regulatory enrichment exists
DF_categories = DF_categories.sort_values(['dataset','component','type','qvalue'],ascending=[1,1,0,1]).drop_duplicates(['dataset','component'])

# Reorganize columns
col_order = ['dataset','component','type','description','pvalue','qvalue','precision','recall','f1score','TP']
DF_categories = DF_categories.set_index('name')[col_order]

# Add empty rows for uncharacterized comps
other_comps = []
idx_list = []
for ds in datasets:
    for k in objs[ds].M.columns:
        idx = ds+'_'+str(k)
        if idx not in DF_categories.index:
            idx_list.append(idx)
            other_comps.append([ds,k,'uncharacterized']+[None]*7)
DF_other = pd.DataFrame(other_comps,columns=DF_categories.columns,index=idx_list)

DF_categories = pd.concat([DF_categories,DF_other]).sort_values(['dataset','component'])
DF_categories.to_csv(os.path.join(DATA_DIR,'iModulons','categories.csv'))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15029469548133595, 'fraction_numerical': 0.01068271119842829, 'mean_word_length': 4.426382411725516, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '381450', 'license': 'MIT', 'n_tokens_mistral': 3026, 'n_tokens_neox': 2933, 'n_words': 380, 'path': 'notebooks/iModulon_annotation/2_characterize_components.ipynb', 'repo_name': 'SBRG/xplatform_ica_paper', 'size': '8144'}