【ChEMBL/STRING】化合物の標的Pathwayを取得【正確確率検定】

with open(f'{TEMP_DIR}/pathwayTerm_to_symbols.pkl', 'rb') as tf:
    ptm_syms = pickle.load(tf)

with open(f'{TEMP_DIR}/pathwayTerm_to_symbols.pkl', 'rb') as tf:
    ptm_syms = pickle.load(tf)

# n_proteins
proteins_string_chembl = list(set(pathways_proteins) | set(compounds_proteins))
n_proteins = len(proteins_string_chembl)
print(f'# of unique proteins found in either string dataset or chembl dataset: {n_proteins}')

# library
import pandas as pd

# contingency table: each compound x each pathway
def make_contingency_table(
    target_proteins,
    proteins_in_pathway,
    n_proteins
):
    _b = set(target_proteins) & set(proteins_in_pathway)
    _a = set(proteins_in_pathway) - _b
    _d = set(target_proteins) - _b
    
    a, b, d = len(list(_a)), len(list(_b)), len(list(_d))
    c = n_proteins - (a + b + d)
    
    df_table = pd.DataFrame([[a, b], [c, d]])
    df_table.columns = ['~is_target', 'is_target']
    df_table.index = ['isin_pathway', '~isin_pathway']
    df_table = df_table.T
    
    return df_table

# library
from scipy.stats import fisher_exact

# fisher's exact test
def calc_p(table):
    odds_ratio, p = fisher_exact(table, alternative='less')
    return p

# make table & tests for all inputs
def enrichment_analysis(
    inkList, ptmList, n_proteins,
    ink_syms, ptm_syms
):
    pvals = []
    
    for i, ik in enumerate(inkList):
        pvals.append([])

        for pt in ptmList:
            table = make_contingency_table(ink_syms[ik], ptm_syms[pt], n_proteins)
            pvals[-1].append(calc_p(table))

    return pvals

# library
import datetime
import pickle
import time
import os


# const
TEMP_DIR = 'temp'


# 結果の保存先のディレクトリ名に用いるタイムスタンプ出力関数
def now():
    return str(datetime.datetime.now().strftime('%Y-%m-%d %H-%M-%S'))

# make table & tests for all inputs
def enrichment_analysis(
    inkList, ptmList, n_proteins,
    ink_syms, ptm_syms,
    rep_interval=50
):
  	# 保存先の作成
	no = now()
    os.makedirs(f'{TEMP_DIR}/pvals/{no}/{len(inkList)}', exist_ok=True)
    
    # 計算時間測定用
    st0 = time.time()
    st = time.time()
    
    # enrichment analysis
    pvals = []    
    for i, ik in enumerate(inkList):
        pvals.append([])

        for pt in ptmList:
            table = make_contingency_table(ink_syms[ik], ptm_syms[pt], n_proteins)
            pvals[-1].append(calc_p(table))

        # save & stdout
        if (i + 1) % rep_interval == 0:
            # save temp array
            with open(f"{TEMP_DIR}/pvals/{len(inkList)}/pvals.{i + 1}.pkl", "wb") as tf:
                pickle.dump(pvals[(i + 1) - rep_interval:(i + 1)], tf)
            
            # stdout (time)
            print(f'{i + 1} / {len(inkList)}: Elapsed time was {time.time() - st:.1f} sec.')
            st = time.time()
    
    # save results
    with open(f"{TEMP_DIR}/pvals/{len(inkList)}/pvals.{i + 1}.pkl", "wb") as tf:
        pickle.dump(pvals[len(pvals) - (len(pvals) % rep_interval):], tf)        

    # stdout (time)
    print(f'Finished: {time.time() - st0:.1f} sec.')
    
    return pvals

# ChEMBLから取得した化合物のInChIKeyリスト
inks = list(ink_syms.keys())

# STRINGから取得したpathwayのterm(ID)のリスト
ptms = [pt for pt in ptm_syms.keys() if 3 <= len(ptm_syms[pt])]

# enrichment analysis
pvals = enrichment_analysis(
    inks, ptms, n_proteins,
    ink_syms, ptm_syms,
    rep_interval=100
)

# library
import numpy as np

# dataframe
df_pvals = pd.DataFrame(pvals)
df_pvals.index = inks
df_pvals.columns = ptms

# score
df_pvals = df_pvals.drop_duplicates()
df_score = df_pvals.map(lambda x: -np.log10(x))

# save
no = now()
os.makedirs(f'{TEMP_DIR}/scores/{no}/df_score', exist_ok=True)
with open(f'{TEMP_DIR}/scores/{no}/df_score.pkl', 'wb') as tf:
    pickle.dump(df_score, tf)