OneK1K benchmark inconsistent result

[jasonbhn]

Hello!

Thank you for developing an easy to use suite of tools.

I'm trying to figure out why my OneK1K eQTL benchmark is not getting any correlations.

My process:

1. Use eQTL Catalogue QTD000620 susie results as source of causal eQTLs (cs with max PiP > 0.9)
2. Using prediction dimensions that matches cell_type == "natural killer cell" and tissue == "blood", take average across these dimensions for each prediction from decima.vep.predict_variant_effect method, with designated gene columns. (basically taken from the cell_type_match from decima-applications.
3. After plotting the predicted effect vs the observed beta I found 2 things:
   1. the predictions are mostly positive
   2. the negative effect size variants from 1k1k are almost all positive, as if they're flipped somehow.
      The result looks like this: the left are variants that are lower in pip from the cs, and the right are those with pip > 0.9.

Appreciate any insights. thanks

[avantikalal]

@MuhammedHasan please take a look

[MuhammedHasan]

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

[jasonbhn]

@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