OneK1K benchmark inconsistent result #52

Activity

1. avantikalal

   assigned

   MuhammedHasan

   on Mar 17, 2026
2. avantikalal commented on Mar 17, 2026

   

   avantikalal

   on Mar 17, 2026

   Collaborator

   More actions

   @MuhammedHasan please take a look
3. MuhammedHasan commented on Mar 17, 2026

   

   MuhammedHasan

   on Mar 17, 2026

   Collaborator

   More actions

   @jasonbhn Can you provide code or jupyter notebook to reproduce the plots above?

   👍React with 👍1Reacted by jasonbhn
4. jasonbhn commented on Mar 18, 2026

   

   jasonbhn

   on Mar 18, 2026

   Author

   More actions

   @MuhammedHasan thanks! I used the GRCh38.primary_assembly.genome.fa fasta from gencode I believe.

   !wget https://ftp.ebi.ac.uk/pub/databases/spot/eQTL/susie/QTS000038/QTD000620/QTD000620.credible_sets.tsv.gz
   !wget https://ftp.ebi.ac.uk/pub/databases/spot/eQTL/sumstats/QTS000038/QTD000620/QTD000620.cc.tsv.gz
   

   python blocks to get result

   import seaborn as sns
   from matplotlib import pyplot as plt
   import pandas as pd
   import torch,gc
   import h5py
   import scanpy
   import polars as pl
   import glob
   from mygene import MyGeneInfo
   from decima.vep import predict_variant_effect
   from decima import DecimaResult
   
   device = "cuda" if torch.cuda.is_available() else "cpu"
   
   %matplotlib inline
   NK_Susie_eQTLCatalogue = pl.read_csv('QTD000620.credible_sets.tsv.gz',separator='\t')
   
   NK_Susie_eQTLCatalogue_high_conf_cs = NK_Susie_eQTLCatalogue.with_columns(
       pl.col('pip').max().over('gene_id','cs_id').alias('max_pip')
   ).filter(
       (pl.col('max_pip')>0.9)
   ).sort('gene_id','cs_id','pip')
   
   genes = NK_Susie_eQTLCatalogue_high_conf_cs['gene_id'].unique().to_list()
   mg = MyGeneInfo()
   res = mg.querymany(
       genes,
       scopes="ensembl.gene",
       fields="symbol,name,entrezgene,ensembl.gene",
       species="human",
       as_dataframe=True,
   )['symbol'].reset_index().rename(columns={'query':'gene_id'})
   
   NK_Susie_eQTLCatalogue_high_conf_cs_decima_in = NK_Susie_eQTLCatalogue_high_conf_cs.with_columns([
       (pl.col('variant').str.split('_').list[0]).alias('chrom'), 
       ((pl.col('variant').str.split('_').list[1].cast(pl.Int64))-1).alias('pos'), 
       (pl.col('variant').str.split('_').list[2]).alias('ref'), 
       (pl.col('variant').str.split('_').list[3]).alias('alt'), 
   ]).select('chrom','pos','ref','alt',pl.exclude(['chrom','pos','ref','alt'])).join(pl.from_pandas(res),on='gene_id',how='left')
   
   decima_old_res = DecimaResult.load()              # or DecimaResult.load("rep0") / model name
   allowed = set(decima_old_res.genes)               # the model’s gene universe (~18,457)
   NK_Susie_eQTLCatalogue_high_conf_cs_decima_in = NK_Susie_eQTLCatalogue_high_conf_cs_decima_in.filter(pl.col('symbol').is_not_null())
   missing = sorted(set(NK_Susie_eQTLCatalogue_high_conf_cs_decima_in["symbol"]) - allowed)
   NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed = NK_Susie_eQTLCatalogue_high_conf_cs_decima_in.filter(~pl.col('symbol').is_in(missing))
   
   
   NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv = NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed.filter((pl.col('ref').str.len_chars()==1)&(pl.col('alt').str.len_chars()==1))
   
   Onek1k_eQTLCatalogue_nk_susie_result = predict_variant_effect(NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv.to_pandas(),
                                   gene_col="symbol",
                                   genome='../borzoi-pytorch/notebooks/GRCh38.primary_assembly.genome.fa')
   
   precomp_result = DecimaResult.load()
   meta = precomp_result.anndata  # AnnData with cell-type labels
   tracks = meta.obs[
   (meta.obs["cell_type"]=='natural killer cell')&(meta.obs["tissue"]=='blood')].index
   
   Onek1k_eQTLCatalogue_nk_susie_result_processed = pd.concat([Onek1k_eQTLCatalogue_nk_susie_result.iloc[:,:14],
                            Onek1k_eQTLCatalogue_nk_susie_result[[i for i in Onek1k_eQTLCatalogue_nk_susie_result.columns if i in tracks]].mean(axis=1)],axis=1).\
   rename(columns={0:'Decima_LFC'})
   
   qtl_labels = NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv.rename({'symbol':'gene'}).\
   select('chrom','pos','ref','alt','gene','cs_id','pip','max_pip','beta').to_pandas()
   
   NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv_compare = Onek1k_eQTLCatalogue_nk_susie_result_processed.merge(qtl_labels,how='left',on=['chrom','pos','ref','alt','gene'])
   NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv_compare['inferred_causal'] = [True if i>0.9 else False for i in NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv_compare['pip']]
   
   # plot
   p = sns.lmplot(
       data=NK_Susie_eQTLCatalogue_high_conf_cs_decima_in_allowed_snv_compare,
       x='beta',
       y='Decima_LFC',col='inferred_causal')
   #plt.xlim(-2, 2) # Limit x-axis from -5 to 5
   #plt.ylim(-1, 1) # Limit y-axis from -1 to 1
   p.refline(x=0, color="black", linestyle="--", linewidth=1)
   p.refline(y=0, color="black", linestyle="--", linewidth=1)
   
   for ax in p.axes.flat:
       ax.set_xlabel('OneK1K beta')
       ax.set_ylabel('Decima NK Cells tracks')
   plt.tight_layout()
   

   cheers
5. MuhammedHasan commented on Jun 4, 2026

   

   MuhammedHasan

   on Jun 4, 2026

   Collaborator

   More actions

   Hi @jasonbhn,

   Following merge request fixes this issue: #57 and it should reproduce the results.
6. MuhammedHasan

   closed this as completedon Oct 1, 2026