<jupyter_start><jupyter_text>Load gene mappings<jupyter_code>with open(os.path.join(conf.GENES_METADATA_DIR, 'genes_mapping_simplified-0.pkl'), 'rb') as f:
    genes_mapping_0 = pickle.load(f)

with open(os.path.join(conf.GENES_METADATA_DIR, 'genes_mapping_simplified-1.pkl'), 'rb') as f:
    genes_mapping_1 = pickle.load(f)<jupyter_output><empty_output><jupyter_text>Load MultiXcan associations<jupyter_code>genes_associations_filename = os.path.join(conf.GENE_ASSOC_DIR, 'smultixcan-mashr-zscores.pkl.xz')
display(genes_associations_filename)

genes_associations = pd.read_pickle(genes_associations_filename)
display(genes_associations.shape)
display(genes_associations.head())<jupyter_output><empty_output><jupyter_text>Load PheWAS catalog<jupyter_code>phewas_catalog = pd.read_csv(os.path.join(conf.DATA_DIR, 'phewas-catalog.csv.gz'), dtype={'phewas code': str})
phewas_catalog.shape
phewas_catalog[phewas_catalog['phewas code'].isna()].head()
phewas_catalog[phewas_catalog['gene_name'].isna()].head()
phewas_catalog[phewas_catalog['gene_name'].isna()].shape
phewas_catalog = phewas_catalog.dropna(subset=['gene_name', 'phewas code'])
phewas_catalog.shape
phewas_catalog['gene_name'].unique().shape
phewas_catalog['phewas code'].unique().shape
phewas_catalog = phewas_catalog.assign(gene_id=phewas_catalog['gene_name'].apply(lambda x: genes_mapping_1[x] if x in genes_mapping_1 else None))
phewas_catalog = phewas_catalog.dropna(subset=['gene_name', 'gene_id', 'phewas code'])
phewas_catalog.shape
phewas_catalog.head()
phewas_catalog.sort_values('phewas phenotype').head()<jupyter_output><empty_output><jupyter_text>Genes in common<jupyter_code>shared_gene_ids = \
    set(phewas_catalog['gene_id'].values)\
    .intersection(genes_associations.index)
len(shared_gene_ids)<jupyter_output><empty_output><jupyter_text>HPO to MIM<jupyter_code>hpo_to_mim = pd.read_csv(os.path.join(conf.DATA_DIR, 'hpo-to-omim-and-phecode.csv.gz'), dtype={'phecode': str})
hpo_to_mim.shape
hpo_to_mim.head()<jupyter_output><empty_output><jupyter_text>Load silver standard to map from UKB to MIM<jupyter_code>omim_silver_standard = pd.read_csv(os.path.join(conf.DATA_DIR, 'omim_silver_standard.tsv'), sep='\t')
ukb_to_mim_map = omim_silver_standard[['trait', 'pheno_mim']].dropna()
ukb_to_mim_map.shape
ukb_to_mim_map.head()<jupyter_output><empty_output><jupyter_text>Read gwas2gene results<jupyter_code>from glob import glob

import rpy2.robjects as robjects
from rpy2.robjects import pandas2ri
pandas2ri.activate()
readRDS = robjects.r['readRDS']
f_files = glob(os.path.join(conf.OMIM_SILVER_STANDARD_GWAS_TO_GENE_DIR, '*.rds'))
display(len(f_files))

if len(f_files) != len(omim_silver_standard['trait'].unique()):
    print(f'WARNING: some files are not there. {len(omim_silver_standard.trait.unique())} expected, {len(f_files)} found.')
gwas2genes_results = {}

for f in f_files:
    f_base = os.path.basename(f)
    f_code = f_base.split('.')[0]
    
    #print(f_base)
    rds_contents = readRDS(f)
    
    if len(rds_contents[1]) > 0:
        f_gene_list = list(rds_contents[1][0].iter_labels())
    else:
        print(f'{f_code}: empty')
        f_gene_list = []
    
    gwas2genes_results[f_code] = f_gene_list
gwas2gene_all_genes = []

for k in gwas2genes_results.keys():
    gwas2gene_all_genes.extend(gwas2genes_results[k])

display(len(gwas2gene_all_genes))

gwas2gene_all_genes = set(gwas2gene_all_genes)
display(len(gwas2gene_all_genes))

gwas2gene_all_genes = shared_gene_ids.intersection(gwas2gene_all_genes)
display(len(gwas2gene_all_genes))
pd.Series(list(gwas2gene_all_genes)).head()<jupyter_output><empty_output><jupyter_text>Universe<jupyter_code>from entity import Trait
_ukb_traits = []
_ukb_traits_phecodes = []
_ukb_gene_available = []

for t in ukb_to_mim_map['trait'].unique():
    t_code = Trait(full_code=t).code
    if t_code not in gwas2genes_results:
        print(t_code)
        continue
    
    for g in gwas2genes_results[t_code]:
        _ukb_traits.append(t)
        _ukb_gene_available.append(g)
df = pd.DataFrame({'trait': _ukb_traits, 'gene': _ukb_gene_available})
df.shape
df.drop_duplicates().shape
df.head()<jupyter_output><empty_output><jupyter_text>Add MIM/Phecode<jupyter_code># add mim
_tmp = pd.merge(df, ukb_to_mim_map, on='trait', how='inner')
display(_tmp.shape)
display(_tmp.head())
_tmp[_tmp['pheno_mim'].isna()].shape
# mim to phecode
_tmp = pd.merge(_tmp, hpo_to_mim[['phecode', 'dID']].dropna(), left_on='pheno_mim', right_on='dID', how='inner').drop(columns=['dID'])
display(_tmp.shape)
display(_tmp.head())
_tmp[_tmp['phecode'].isna()].shape
_tmp.head()
# phecode to phewas catalog
_tmp = pd.merge(_tmp, phewas_catalog[['phewas code', 'gene_id']],
                left_on=['phecode', 'gene'], right_on=['phewas code', 'gene_id'],
                how='left').drop(columns=['phewas code'])
display(_tmp.shape)
_tmp.head()
_tmp[_tmp['gene_id'].isna()].shape
_tmp = _tmp.drop_duplicates(subset=['trait', 'gene', 'gene_id'])
_tmp.shape
_tmp.head(30)
_tmp[_tmp['gene_id'].isna()].shape
def _assign_true_class(x):
    tc = ~pd.isnull(x['gene_id'])
    idx = [0]
    if tc.shape[0] > 1 and tc.any():
        idx = np.where(tc)[0]
    return pd.Series({
        'pheno_mim': ', '.join(x.iloc[idx]['pheno_mim'].astype(str)),
        'phecode': ', '.join(x.iloc[idx]['phecode'].astype(str)),
        'true_class': int(tc.any()),
    })
_tmp2 = _tmp.groupby(['trait', 'gene']).apply(_assign_true_class)
_tmp2.shape
_tmp2.head()
assert not _tmp2.loc['M41-Diagnoses_main_ICD10_M41_Scoliosis', 'ENSG00000012504']['true_class']
assert not _tmp2.loc['M41-Diagnoses_main_ICD10_M41_Scoliosis', 'ENSG00000141665']['true_class']
assert _tmp2.loc['M41-Diagnoses_main_ICD10_M41_Scoliosis', 'ENSG00000112137']['true_class']
_tmp2['true_class'].value_counts()
18632 / 20837
2205 / 20837<jupyter_output><empty_output><jupyter_text>Add score<jupyter_code>_genes_unstacked = genes_associations.unstack()
_genes_unstacked.shape
_genes_unstacked.head()
classifier_table = _tmp2.assign(score=_genes_unstacked)
classifier_table.shape
classifier_table.head()
classifier_table[classifier_table['score'].isna()].shape
classifier_table = classifier_table.dropna(subset=['score'])
classifier_table.shape
N_TESTS = classifier_table.reset_index().drop_duplicates(subset=['trait', 'gene']).shape[0]
display(N_TESTS)

PVALUE_THRESHOLD = (0.05 / (N_TESTS))
display(PVALUE_THRESHOLD)

ZSCORE_THRESHOLD = np.abs(stats.norm.ppf(PVALUE_THRESHOLD / 2))
display(ZSCORE_THRESHOLD)
classifier_table = classifier_table.assign(predicted_class=(classifier_table['score'] > ZSCORE_THRESHOLD).astype(int))
classifier_table.head()
classifier_table['true_class'].value_counts()
classifier_table['true_class'].value_counts().sum()<jupyter_output><empty_output><jupyter_text>Save classifier table<jupyter_code>classifier_table = classifier_table.sort_index()
assert classifier_table.index.is_unique
classifier_table.head()
classifier_table.shape
classifier_table.to_csv(
    os.path.join(output_dir, 'smultixcan-mashr-classifier_data-phewas_catalog.tsv.gz'),
    sep='\t', index=False
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14776245426400225, 'fraction_numerical': 0.017027863777089782, 'mean_word_length': 7.470798569725864, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '34849', 'license': 'MIT', 'n_tokens_mistral': 2981, 'n_tokens_neox': 2793, 'n_words': 334, 'path': 'scripts/150_validation/phewas_catalog/11_roc-phewas_catalog_standard-mashr_models.ipynb', 'repo_name': 'miltondp/phenomexcan', 'size': '7106'}