<jupyter_start><jupyter_text>Load S-MultiXcan results<jupyter_code>smultixcan_pvalues_file = os.path.join(conf.GENE_ASSOC_DIR, f'smultixcan-mashr-pvalues.pkl.xz')
display(smultixcan_pvalues_file)
smultixcan_gene_associations = pd.read_pickle(smultixcan_pvalues_file)
smultixcan_gene_associations.shape
smultixcan_gene_associations.head(5)<jupyter_output><empty_output><jupyter_text>Load fastENLOC results<jupyter_code>fastenloc_rcp_file = os.path.join(conf.GENE_ASSOC_DIR, f'fastenloc-torus-rcp.pkl.xz')
display(fastenloc_rcp_file)
fastenloc_gene_associations = pd.read_pickle(fastenloc_rcp_file)
fastenloc_gene_associations.shape
fastenloc_gene_associations.head(5)
assert fastenloc_gene_associations.min().min() >= 0
assert fastenloc_gene_associations.max().max() <= 3<jupyter_output><empty_output><jupyter_text>S-MultiXcan  Counts of significant associations<jupyter_code>smultixcan_gene_associations.shape
all_pvalues = pd.Series(smultixcan_gene_associations.values.flatten()).dropna()
display(all_pvalues.shape)
all_pvalues.describe()
(-np.log10(all_pvalues)).describe()
PVALUE_THRESHOLD = (0.05 / all_pvalues.shape[0])
display(PVALUE_THRESHOLD)
PVALUE_THRESHOLD = 5.49e-10
hits = (all_pvalues < PVALUE_THRESHOLD).sum()
display(hits)
display((hits / all_pvalues.shape[0]) * 100)
hits = (all_pvalues < 0.01).sum()
display(hits)
display((hits / all_pvalues.shape[0]) * 100)<jupyter_output><empty_output><jupyter_text>fastENLOC  Numbers of significant genes also found by fastENLOC<jupyter_code>common_genes = fastenloc_gene_associations.index.intersection(smultixcan_gene_associations.index)
display(common_genes)
fastenloc_matrix_sorted = fastenloc_gene_associations.loc[common_genes, smultixcan_gene_associations.columns]
display(fastenloc_matrix_sorted.shape)
display(fastenloc_matrix_sorted.head())
multixcan_matrix_sorted = smultixcan_gene_associations.loc[common_genes]
display(multixcan_matrix_sorted.shape)
display(multixcan_matrix_sorted.head())
assert fastenloc_matrix_sorted.shape == multixcan_matrix_sorted.shape<jupyter_output><empty_output><jupyter_text>Using stringent pvalue threshold<jupyter_code>display(PVALUE_THRESHOLD)
multixcan_signif = (multixcan_matrix_sorted < PVALUE_THRESHOLD)
display(multixcan_signif.shape)
display(multixcan_signif.head())
fastenloc_signif = (fastenloc_matrix_sorted > 0.1)
display(fastenloc_signif.shape)
display(fastenloc_signif.head())
multixcan_hits = multixcan_signif.sum().sum()
display(multixcan_hits)
display(multixcan_hits / multixcan_signif.size)
fastenloc_hits = fastenloc_signif.sum().sum()
display(fastenloc_hits)
display(fastenloc_hits / fastenloc_signif.size)
fastenloc_and_multixcan_signif = (multixcan_signif & fastenloc_signif)
assert fastenloc_and_multixcan_signif.shape == multixcan_signif.shape == fastenloc_signif.shape
fe_mu_sum = fastenloc_and_multixcan_signif.sum().sum()
display(fe_mu_sum)
_perc = fe_mu_sum / multixcan_hits
display(_perc)
display(f'{(_perc * 100):.2f}%')<jupyter_output><empty_output><jupyter_text>Conclusions 72,884 significant associations with S-MultiXcan (after taking the set of genes in common with fastENLOC, otherwise it's 72,994).From these, 22,219 associations are also colocalized (30.49%)  Export significant associations as a table  S-MultiXcan<jupyter_code>smultixcan_results = multixcan_matrix_sorted.unstack()
smultixcan_results.shape
smultixcan_results.head()
smultixcan_results.isna().sum()
smultixcan_results = smultixcan_results.dropna()
smultixcan_results.shape
assert smultixcan_results.isna().sum() == 0
(smultixcan_results < PVALUE_THRESHOLD).sum()<jupyter_output><empty_output><jupyter_text>fastENLOC<jupyter_code>fastenloc_results = fastenloc_matrix_sorted.unstack()
fastenloc_results.shape
fastenloc_results.head()
fastenloc_results.isna().sum()
fastenloc_results = fastenloc_results.dropna()
fastenloc_results.shape
final_results = pd.DataFrame(index=smultixcan_results.index.copy(), columns=['smultixcan_pval', 'fastenloc_rcp'])
final_results = final_results.assign(smultixcan_pval=smultixcan_results, fastenloc_rcp=fastenloc_results)
final_results.shape
final_results.head()
final_results.isna().sum()
final_results = final_results.dropna(how='any')
final_results.shape
PVALUE_THRESHOLD
publishable_final_results = final_results[(final_results['smultixcan_pval'] < PVALUE_THRESHOLD) & (final_results['fastenloc_rcp'] > 0.1)]
publishable_final_results.index.rename(('trait', 'gene_id'), inplace=True)
publishable_final_results.shape
publishable_final_results.head()
metadata.GENES_MAPPINGS.head()
_tmp = pd.merge(
    publishable_final_results.reset_index(),
    metadata.GENES_MAPPINGS[['gene_id', 'gene_name', 'band', 'gene_type']],
    on='gene_id',
    how='left'
).rename(columns={'band': 'gene_band'})
_tmp.head()
publishable_final_results = _tmp.set_index(['trait', 'gene_id'])[['gene_name', 'gene_band', 'gene_type', 'smultixcan_pval', 'fastenloc_rcp']]
publishable_final_results = publishable_final_results.sort_index()
assert publishable_final_results.index.is_unique
publishable_final_results.isna().sum()
publishable_final_results[publishable_final_results['gene_name'].isna()]
publishable_final_results.head()
publishable_final_results.shape
# some testing
_gene_name = 'VIP'
_gene_id = metadata.GENE_NAME_TO_ID_MAP[_gene_name]
display(_gene_id)

_result = publishable_final_results.loc[('1180-Morningevening_person_chronotype', _gene_id)]
assert _result['smultixcan_pval'] == 1.8124965657173678e-17, _result['smultixcan_pval']
assert _result['fastenloc_rcp'] == 0.25849, _result['fastenloc_rcp']
# some testing
_gene_name = 'RP11-220I1.5'
_gene_id = metadata.GENE_NAME_TO_ID_MAP[_gene_name]
display(_gene_id)

_result = publishable_final_results.loc[('1180-Morningevening_person_chronotype', _gene_id)]
assert _result['smultixcan_pval'] == 6.426619359444203e-11
assert _result['fastenloc_rcp'] == 0.1985
# for publication (tsv.gz)
output_folder = os.path.join(conf.DELIVERABLES_DIR, 'supp_tables')
os.makedirs(output_folder, exist_ok=True)
output_file = os.path.join(output_folder, 'suppl_table_S1-significant_gene_trait_associations.tsv.gz')
display(output_file)

publishable_final_results.to_csv(output_file, sep='\t', float_format='%.4e')
# some testing
_tmp = pd.read_csv(output_file, sep='\t')
assert publishable_final_results.shape == _tmp.set_index(['trait', 'gene_id']).shape
_tmp.shape
_tmp.head()
_tmp = _tmp.set_index(['trait', 'gene_id'])
_result = _tmp.loc[('1180-Morningevening_person_chronotype', metadata.GENE_NAME_TO_ID_MAP['VIP'])]
assert _result['smultixcan_pval'] == 1.8125e-17
assert _result['fastenloc_rcp'] == 0.25849

_result = _tmp.loc[('1180-Morningevening_person_chronotype', metadata.GENE_NAME_TO_ID_MAP['RP11-220I1.5'])]
assert _result['smultixcan_pval'] == 6.4266e-11
assert _result['fastenloc_rcp'] == 0.1985
assert np.allclose(
    _tmp[['smultixcan_pval', 'fastenloc_rcp']].values,
    publishable_final_results[['smultixcan_pval', 'fastenloc_rcp']].values,
    atol=1e-320,
    rtol=1e-4
)
# for publication (xlsx)
output_folder = os.path.join(conf.DELIVERABLES_DIR, 'supp_tables')
os.makedirs(output_folder, exist_ok=True)
output_file = os.path.join(output_folder, 'suppl_table_S1-significant_gene_trait_associations.xlsx')
display(output_file)

publishable_final_results.reset_index().to_excel(output_file, index=False, float_format='%.4e')
# some testing
_tmp = pd.read_excel(output_file)
assert publishable_final_results.shape == _tmp.set_index(['trait', 'gene_id']).shape
_tmp.shape
_tmp.head()
_tmp = _tmp.set_index(['trait', 'gene_id'])
_result = _tmp.loc[('1180-Morningevening_person_chronotype', metadata.GENE_NAME_TO_ID_MAP['VIP'])]
assert _result['smultixcan_pval'] == 1.8125e-17
assert _result['fastenloc_rcp'] == 0.25849

_result = _tmp.loc[('1180-Morningevening_person_chronotype', metadata.GENE_NAME_TO_ID_MAP['RP11-220I1.5'])]
assert _result['smultixcan_pval'] == 6.4266e-11
assert _result['fastenloc_rcp'] == 0.1985
assert np.allclose(
    _tmp[['smultixcan_pval', 'fastenloc_rcp']].values,
    publishable_final_results[['smultixcan_pval', 'fastenloc_rcp']].values,
    atol=1e-320,
    rtol=1e-4
)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup