<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 S-MultiXcan results<jupyter_code>smultixcan_genes_associations_filename = os.path.join(conf.GENE_ASSOC_DIR, 'smultixcan-mashr-zscores.pkl.xz')
display(smultixcan_genes_associations_filename)

smultixcan_genes_associations = pd.read_pickle(smultixcan_genes_associations_filename)
smultixcan_genes_associations.shape
smultixcan_genes_associations.head(5)<jupyter_output><empty_output><jupyter_text>Load fastENLOC results<jupyter_code>fastenloc_genes_associations_filename = os.path.join(conf.GENE_ASSOC_DIR, 'fastenloc-torus-rcp.pkl.xz')
display(fastenloc_genes_associations_filename)

fastenloc_genes_associations = pd.read_pickle(fastenloc_genes_associations_filename)
fastenloc_genes_associations.shape
fastenloc_genes_associations.head(5)<jupyter_output><empty_output><jupyter_text>Genes in common between S-MultiXcan and fastENLOC<jupyter_code>common_genes = fastenloc_genes_associations.index.intersection(smultixcan_genes_associations.index)
display(common_genes)<jupyter_output><empty_output><jupyter_text>Load OMIM silver standard<jupyter_code>omim_silver_standard = pd.read_csv(os.path.join(conf.DATA_DIR, 'omim_silver_standard.tsv'), sep='\t')
omim_silver_standard = omim_silver_standard.dropna(subset=['ensembl_gene_id', 'trait', 'pheno_mim'])
display(omim_silver_standard.shape)
display(omim_silver_standard.head())<jupyter_output><empty_output><jupyter_text>Read gwas2gene results These results were generated with scripts in `scripts/extras/gwas2gene`.<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, _ = os.path.splitext(f_base)
    
    #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] = common_genes.intersection(set(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))<jupyter_output><empty_output><jupyter_text>Create list of UKB-OMIM traits<jupyter_code>omim_silver_standard.head()<jupyter_output><empty_output><jupyter_text>Create PrediXcan classifier table<jupyter_code>_tmp = omim_silver_standard[['trait', 'ensembl_gene_id']]
ukb_traits_common = _tmp['trait'].unique()

omim_true_classes = _tmp[['trait', 'ensembl_gene_id']].drop_duplicates()
omim_true_classes = omim_true_classes.assign(omim_value=1)
omim_true_classes = omim_true_classes.set_index(['trait', 'ensembl_gene_id'])
len(ukb_traits_common)
omim_true_classes.shape
omim_true_classes.head()
len(ukb_traits_common)
from entity import Trait
index_tuples = []

for t in ukb_traits_common:
    t_code = Trait(full_code=t).code
    if t_code not in gwas2genes_results:
        continue
    
    for g in gwas2genes_results[t_code]:
        index_tuples.append((t, g))
len(index_tuples)
index_tuples[:5]
classifier_index = pd.MultiIndex.from_tuples(
    index_tuples,
    names=['ukb_efo', 'gene']
)
len(gwas2gene_all_genes)
classifier_index.shape
predixcan_classifier_df = pd.DataFrame(index=classifier_index, columns=['score', 'predicted_class', 'true_class'])
predixcan_classifier_df = predixcan_classifier_df.sort_index()
predixcan_classifier_df.shape
predixcan_classifier_df['true_class'] = 0
predixcan_classifier_df.head()
true_classes = omim_true_classes.squeeze()
display(true_classes.shape)
display(true_classes.head())
predixcan_classifier_df.loc[predixcan_classifier_df.index.intersection(true_classes.index), 'true_class'] = 1
assert predixcan_classifier_df['true_class'].isna().sum() == 0
# some testing
predixcan_classifier_df.loc[('M41-Diagnoses_main_ICD10_M41_Scoliosis',)].head()
true_classes.loc[('M41-Diagnoses_main_ICD10_M41_Scoliosis', 'ENSG00000112234')]
'ENSG00000090263' not in true_classes.loc['M41-Diagnoses_main_ICD10_M41_Scoliosis'].index
assert predixcan_classifier_df.loc[('M41-Diagnoses_main_ICD10_M41_Scoliosis', 'ENSG00000112234'), 'true_class'] == 1.0
assert predixcan_classifier_df.loc[('M41-Diagnoses_main_ICD10_M41_Scoliosis', 'ENSG00000090263'), 'true_class'] == 0.0
len(gwas2gene_all_genes)
# score
df_score = pd.Series(index=classifier_index, dtype=float)

for trait in ukb_traits_common:
    trait_code = Trait(full_code=trait).code
    if trait_code not in gwas2genes_results:
        print(trait_code)
        continue
    trait_genes = gwas2genes_results[trait_code]
    if len(trait_genes) == 0:
        print(f'Empty: {trait}')
        continue
    
    smultixcan_zscores = smultixcan_genes_associations.loc[trait_genes, trait]
    fastenloc_rcps = fastenloc_genes_associations.loc[trait_genes, trait].fillna(0.0)
    
    smultixcan_zscores[(fastenloc_rcps < RCP_CUTOFF) & (~np.isnan(smultixcan_zscores))] = 0
    scores = smultixcan_zscores
    df_score.loc[trait] = scores.values
# for S-MultiXcan dropna
df_score = df_score.dropna().sort_index()
assert df_score.isna().sum().sum() == 0
df_score.head()
# some testing
_gene, _trait = ('ENSG00000090263', 'M41-Diagnoses_main_ICD10_M41_Scoliosis')
display(smultixcan_genes_associations.loc[_gene, _trait])
display(fastenloc_genes_associations.loc[_gene, _trait])
assert 0.0 == df_score.loc[_trait, _gene]
_gene, _trait = ('ENSG00000070061', 'O14-Diagnoses_main_ICD10_O14_Gestational_pregnancyinduced_hypertension_with_significant_proteinuria')
display(smultixcan_genes_associations.loc[_gene, _trait])
display(fastenloc_genes_associations.loc[_gene, _trait])
assert 0.0 == df_score.loc[_trait, _gene]
_gene, _trait = ('ENSG00000004534', '1200-Sleeplessness_insomnia')
display(smultixcan_genes_associations.loc[_gene, _trait])
display(fastenloc_genes_associations.loc[_gene, _trait])
assert smultixcan_genes_associations.loc[_gene, _trait] == df_score.loc[_trait, _gene]
df_score.shape
df_score.head()
df_score.describe()
predixcan_classifier_df = predixcan_classifier_df.assign(score=df_score)
# assert not predixcan_classifier_df['score'].isna().any()
from scipy import stats
_n_genes = len(gwas2gene_all_genes)
display(_n_genes)

_n_ukb_traits = len(omim_silver_standard['trait'].unique())
display(_n_ukb_traits)

display(_n_genes * _n_ukb_traits)

PVALUE_THRESHOLD = (0.05 / (_n_genes * _n_ukb_traits))
display(PVALUE_THRESHOLD)

ZSCORE_THRESHOLD = np.abs(stats.norm.ppf(PVALUE_THRESHOLD / 2))
display(ZSCORE_THRESHOLD)
def _assign_predicted_class(x):
#     trait, gene = x.name
#     smultixcan_zscore = smultixcan_genes_associations.loc[gene, trait]
    if x > ZSCORE_THRESHOLD:
        return 1
    else:
        return 0
predixcan_classifier_df = predixcan_classifier_df.assign(predicted_class=predixcan_classifier_df['score'].apply(_assign_predicted_class))
predixcan_classifier_df.shape
predixcan_classifier_df.head()
predixcan_classifier_df.loc['M41-Diagnoses_main_ICD10_M41_Scoliosis'].sort_values('true_class', ascending=False).head()<jupyter_output><empty_output><jupyter_text>Select genes per trait<jupyter_code>#selected_predixcan_classifier_df = predixcan_classifier_df.loc[predixcan_classifier_df.index.intersection(trait_genes_to_keep)]
selected_predixcan_classifier_df = predixcan_classifier_df
# some testing
selected_predixcan_classifier_df.shape
selected_predixcan_classifier_df.head()
selected_predixcan_classifier_df.sort_values('predicted_class', ascending=False).head()
_tmp = selected_predixcan_classifier_df.sort_values(['true_class', 'ukb_efo'], ascending=False)
display(_tmp.shape)
display(_tmp[_tmp['true_class'] > 0].shape)
display(_tmp[_tmp['true_class'] > 0].head())<jupyter_output><empty_output><jupyter_text>Test classes<jupyter_code>selected_predixcan_classifier_df.index.get_level_values('ukb_efo').unique().shape
selected_predixcan_classifier_df.index.get_level_values('gene').unique().shape
_pheno = 'N20-Diagnoses_main_ICD10_N20_Calculus_of_kidney_and_ureter'
_clinvar_asthma_genes = omim_silver_standard[omim_silver_standard['trait'] == _pheno]['ensembl_gene_id'].unique()
display(_clinvar_asthma_genes)
display(_clinvar_asthma_genes.shape)
_tmp = selected_predixcan_classifier_df.loc[_pheno]
_tmp.loc[_tmp.index.intersection(_clinvar_asthma_genes)]
_predixcan_asthma_genes = selected_predixcan_classifier_df.loc[_pheno]
_predixcan_asthma_genes.head()
selected_predixcan_classifier_df.shape
selected_predixcan_classifier_df['predicted_class'].value_counts()
selected_predixcan_classifier_df['true_class'].value_counts()
selected_predixcan_classifier_df.sort_values(['true_class'], ascending=[False])<jupyter_output><empty_output><jupyter_text>Save classifier table<jupyter_code># remove nans
selected_predixcan_classifier_df = selected_predixcan_classifier_df.dropna()
assert selected_predixcan_classifier_df.index.is_unique
selected_predixcan_classifier_df.head()
selected_predixcan_classifier_df.shape
# this variant is not used, do not save

# selected_predixcan_classifier_df.to_csv(
#     os.path.join(output_dir, 'combined-classifier_data.tsv.gz'),
#     sep='\t', index=False
# )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11024656426839127, 'fraction_numerical': 0.019603880355699273, 'mean_word_length': 8.94673366834171, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, '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': '52938', 'license': 'MIT', 'n_tokens_mistral': 4016, 'n_tokens_neox': 3771, 'n_words': 429, 'path': 'scripts/150_validation/omim/01_roc-omim_standard-combined-cutoff.ipynb', 'repo_name': 'miltondp/phenomexcan', 'size': '9896'}