An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression.

Code ↔ Paper

5 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 5 matches

[1] § Materials and methods › Instrument selection ↔ gwas_norm/normalise.py, lines 503–962 · score 0.60 · chi square, minor allele frequencies, minimal, GWAS
[2] § Materials and methods › Instrument selection ↔ merit/coloc/coloc_core.py, lines 1367–1445 · score 0.60 · minor allele frequencies, removing variants, squares, selection, trait, correlations
[3] § Materials and methods › Sensitivity analyses and target annotation › Colocalisation ↔ merit/coloc/coloc_core.py, lines 334–437 · score 0.58 · posterior probability, PPH3, abf, PPH4, colocalisation, traits
[4] § Materials and methods › Primary Mendelian randomisation (MR) analysis ↔ merit/mr/ivw.py, lines 21–164 · score 0.57 · inverse variance weighted, IVW, instrumented, correlated, MR, variant
[5] § Materials and methods › Sensitivity analyses and target annotation › Estimated therapeutic relevance ↔ bio_misc/drug_lookups/therapeutic_direction.py, lines 102–230 · score 0.54 · therapeutic directional, ChEMBL, activator, inhibitor, unknown, risk

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 2,566 lines · 104 KB · GPL-3.0 · 2 matches

"""
Python implementation of the coloc R package.
See the original R code here (06-04-2021 last update):
`Here <https://github.com/chr1swallace/coloc/blob/main/R/split.R>`__.
Relevant documentation:
`Here <https://chr1swallace.github.io/coloc/articles/a02_data.html>`__.
And relevant papers:
`Here <https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1004383>`__.
`Here <https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1008720>`__.
`Here <https://www.biorxiv.org/content/10.1101/2021.02.23.432421v1>`__.
Note, the `Sum of Single Effects` SuSIE version is not yet implementied
Attribution
-----------
All the credit goes towards the original authors.
All mistakes are our own.
Note:
-----
Some function were omitted, because there were not used in R coloc:
-. est.cond.nometa
-. bin2lin.nometa
-. estgeno.1.ctl
-. estgeno.1.cse
"""
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import warnings
import itertools
import scipy.stats
# from pandas.core.frame import DataFrame
from collections import OrderedDict
from merit.errors import(
_process_corr_mat,
are_columns_in_df,
are_Series_equal,
error_on_non_finite_pd_corr_matrix,
is_type,
)
from merit import constants as c
from typing import Any, List, Type, Union, Tuple
SEP = ''
"""
TODO:
Nothing ...
"""
# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Constants
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Coloc_Expected(object):
"""
Collecting some expected constants, getting around all the `SEP.join` that
would otherwise clutter the code.
Might move this to merit.constants, depending on the development of
the coloc implementation into separate scripts.
Notice:
-------
In R coloc `standard_deviation`, `sample_size` and `event_rate` are floats,
not columns. Here we have included them as columns and addapted to code to
take the average.
"""
# defining names
exp_name = 'exposure'
out_name = 'outcome'
VARIANT_ID=c.UNIVERSAL_ID.name
EXPOSURE_EFFECT_SIZE= SEP.join(
[c.EFFECT_SIZE.name, c.MERGE_SUFFIX_DELIMITER_UP, exp_name])
OUTCOME_EFFECT_SIZE= SEP.join(
[c.EFFECT_SIZE.name, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# simply se^2
EXPOSURE_VARBETA= SEP.join(
[c.COLOC_VARBETA, c.MERGE_SUFFIX_DELIMITER_UP, exp_name])
OUTCOME_VARBETA= SEP.join(
[c.COLOC_VARBETA, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# standard deviation of the gwas trait
EXPOSURE_STANDARD_DEVIATION = SEP.join(
[c.COLOC_STANDARD_DEVIATION, c.MERGE_SUFFIX_DELIMITER_UP, exp_name])
OUTCOME_STANDARD_DEVIATION = SEP.join(
[c.COLOC_STANDARD_DEVIATION, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# GWAS sample size
EXPOSURE_SAMPLE_SIZE = SEP.join(
[c.COLOC_SAMPLE_SIZE, c.MERGE_SUFFIX_DELIMITER_UP, exp_name]
)
OUTCOME_SAMPLE_SIZE = SEP.join(
[c.COLOC_SAMPLE_SIZE, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# PROPORTION OF CASES
EXPOSURE_EVENT_RATE = SEP.join(
[c.COLOC_EVENT_RATE, c.MERGE_SUFFIX_DELIMITER_UP, exp_name])
OUTCOME_EVENT_RATE = SEP.join(
[c.COLOC_EVENT_RATE, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# z-statistics
EXPOSURE_Z_STATISTICS = SEP.join(
[c.COLOC_NORMAL_DEVIATE, c.MERGE_SUFFIX_DELIMITER_UP, exp_name])
OUTCOME_Z_STATISTICS = SEP.join(
[c.COLOC_NORMAL_DEVIATE, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# Shrinkage factor: ratio of the prior variance to the total variance
EXPOSURE_SHRINKAGE_FACTOR = SEP.join(
[c.COLOC_SHRINKAGE_FACTOR, c.MERGE_SUFFIX_DELIMITER_UP, exp_name])
OUTCOME_SHRINKAGE_FACTOR = SEP.join(
[c.COLOC_SHRINKAGE_FACTOR, c.MERGE_SUFFIX_DELIMITER_UP, out_name])
# The log (approximate) Bayes Factor for colocalization
EXPOSURE_ABF = SEP.join(
[c.COLOC_LOG_APPROX_BAYES_FACTOR, c.MERGE_SUFFIX_DELIMITER_UP,
exp_name])
OUTCOME_ABF = SEP.join(
[c.COLOC_LOG_APPROX_BAYES_FACTOR, c.MERGE_SUFFIX_DELIMITER_UP,
out_name])
# minor allele frequence
# NOTE does this need to be exp/outcome stratified?
MINOR_ALLELE_FREQ = c.MINOR_ALLELE_FREQ.name
# types of gwas traits
TRAIT_TYPES = ['quant', 'bin']
# methods
METHODS=["single", "cond", "mask"]
# modes
MODES=["iterative", "allbutone"]
# colnames for dataset 1
exposure_columns = [
VARIANT_ID,
EXPOSURE_EFFECT_SIZE,
EXPOSURE_VARBETA,
EXPOSURE_STANDARD_DEVIATION,
EXPOSURE_SAMPLE_SIZE,
EXPOSURE_EVENT_RATE,
MINOR_ALLELE_FREQ
]
# names for dataset 2
outcome_columns = [
VARIANT_ID,
OUTCOME_EFFECT_SIZE,
OUTCOME_VARBETA,
OUTCOME_STANDARD_DEVIATION,
OUTCOME_SAMPLE_SIZE,
OUTCOME_EVENT_RATE,
MINOR_ALLELE_FREQ
]
# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Example datasets
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def _coloc_example_data():
"""
Returns a coloc example dataset, and corr_matrix, for five variants
"""
nvariants = 5
idx = ['snp_' + str(l) for l in range(1, 6)]
# dataset
dataset = pd.DataFrame(
{
Coloc_Expected.VARIANT_ID: idx,
Coloc_Expected.EXPOSURE_EFFECT_SIZE : [ 1.14517929e+00,
7.90320821e-01,
-1.33679243e-02,
1.78273663e-01,
-3.77606312e-01],
Coloc_Expected.EXPOSURE_VARBETA : [0.03924328, 0.03880072,
0.04455547, 0.04602953,
0.04237982],
Coloc_Expected.EXPOSURE_SAMPLE_SIZE : [1000] * nvariants,
Coloc_Expected.EXPOSURE_STANDARD_DEVIATION : [2.32924174]*nvariants,
Coloc_Expected.EXPOSURE_EVENT_RATE : [0.10] * nvariants,
Coloc_Expected.OUTCOME_EFFECT_SIZE : [ 1.07128929, 0.92979073,
-0.14149679, -0.02395321,
-0.03075124],
Coloc_Expected.OUTCOME_VARBETA : [0.04028734, 0.04478292,
0.04728015, 0.04748579,
0.04415969],
Coloc_Expected.OUTCOME_SAMPLE_SIZE : [1000] * nvariants,
Coloc_Expected.OUTCOME_STANDARD_DEVIATION : [2.33203611]*nvariants,
Coloc_Expected.OUTCOME_EVENT_RATE : [0.10, 0.11, 0.09, 0.12, 0.08],
Coloc_Expected.MINOR_ALLELE_FREQ : [0.1349454, 0.11783439,
0.13162119, 0.12019231, 0.125]
}, index = idx
)
dataset.index_name=Coloc_Expected.VARIANT_ID
# correlation matrix
corr_matrix = pd.DataFrame(
{
idx[0] : [1.000000, 0.738238, 0.634454, 0.555827, 0.478299],
idx[1] : [0.738238, 1.000000, 0.730387, 0.636565, 0.550489],
idx[2] : [0.634454, 0.730387, 1.000000, 0.735947, 0.618572],
idx[3] : [0.555827, 0.636565, 0.735947, 1.000000, 0.746561],
idx[4] : [0.478299, 0.550489, 0.618572, 0.746561, 1.000000]
}, index = idx)
# return
return dataset, corr_matrix
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# empty results dataframe
def _empty_coloc_results():
'''
Returns an empty coloc results pd.DataFrame
'''
results = pd.DataFrame(columns = [c.COLOC_HIT1, c.COLOC_HIT2, c.MR_NSNPS,
c.COLOC_PPH0, c.COLOC_PPH1, c.COLOC_PPH2,
c.COLOC_PPH3, c.COLOC_PPH4, c.COLOC_BEST1,
c.COLOC_BEST2, c.COLOC_BEST4,
c.COLOC_ZSTAT_HIT1, c.COLOC_ZSTAT_HIT2
], index=[0])
return results
# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Helper functions
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def _add_missing_columns(check_data: pd.DataFrame, complete_data: pd.DataFrame,
columns : list, index_col: Coloc_Expected.VARIANT_ID
):
'''
Add missing columns
Arguments
---------
data : pd.DataFrame
columns : list of strings
index_col : str
'''
# Shared index
check_data.set_index(check_data[index_col], drop=False, inplace=True)
complete_data.set_index(complete_data[index_col], drop=False, inplace=True)
check_data.index.name = None
complete_data.index.name = None
# mapp columns
for col in columns:
if not col in check_data.columns:
check_data[col] = complete_data[col]
# returns stuff
return check_data
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def _logsum(x):
"""
Calculate the sum of terms.
Used by the approximate Bayes factor function to calculate the posterior
probabilities. Depends on numpy (np).
Arguments
---------
x : pd.Series, numpy.array or some other kind of numeric data type
"""
max_res = np.max(x)
results = max_res + np.log(np.sum(np.exp(x - max_res)))
return results
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def _logdiff(x, y):
"""
Calculate the difference of between two vectors.
Used by the approximate Bayes factor function to calculate the posterior
probablities. Depends on numpy (np), used in combination with _logsum.
In the case of coloc its output should be a constant.
Arguments
---------
x, y : pd.Series, numpy.array or some other kind of numeric data type
Returns
-------
A pd.Series with the element wise difference of x-y after relevant
transformation. Results will be NA if x-y < 0.
"""
max_res = np.max(x + y)
results = max_res + np.log(np.exp(x - max_res) - np.exp(y - max_res))
return results
# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Functions originally in coloc/split.R
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def sd_pheno_approx(variancebeta: pd.Series, maf: pd.Series,
n: Union[pd.Series, float]):
"""
Get an estimate of the standard deviation of the phenotype distribution.
Note
----
Implements the `sdY.est` function of R coloc.
Arguments
---------
variancebeta : pd.Series
Squared standard errors of the relevant point estimates.
maf : pd.Series
The minor allele frequencies of the supplied variants.
n : pd.Series
The GWAS sample size.
Returns
-------
np.float64, reflecting an estimate of the phenotypic standard deviation
(of 'Y').
"""
# @@@ check input
is_type(variancebeta, pd.Series)
is_type(maf, pd.Series)
is_type(n,(pd.Series, float))
if not variancebeta.shape[0] == maf.shape[0]:
raise IndexError('Input does not have the same shape')
# @@@ Do the actual calculations
oneover = np.divide(1, variancebeta)
# remove accidental inf
oneover = oneover[~np.isinf(oneover)]
if isinstance(n, pd.Series):
n = n.mean()
nvx = np.multiply(2 * n, maf * np.subtract(1, maf))
# OLS to approximate the SD of the phenotype
ols_res = smf.ols(
formula="nvx ~ oneover -1",
data=pd.DataFrame({"nvx": nvx, "oneover": oneover}),
missing='drop'
)
cf = ols_res.fit().params["oneover"]
if cf < 0:
raise ValueError("The standard deviation is negative. \
Please provide your own estimate of the phenotype SD")
else:
return np.sqrt(cf)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def approx_bf(ln_abf1: pd.Series, ln_abf2: pd.Series,
p1: float=10**-4,
p2: float=10**-4,
p12: float=5*10**-6,
extra_info: bool=True,
masking_ln_abf3: pd.Series=pd.Series(dtype=float),
masking_ln_abf4: pd.Series=pd.Series(dtype=float)):
"""
Calculate the posterior probabilities of colocalized signals.
Depends on numpy (np). Its arguments include the natural logarithm of the
base factor for GWAS data set 1 and GWAS data set 2. P1-2 are constants (
or column vectors of length equal to the number of variants) representing
the probability of a random variant being associated only with GWAS
phenotype 1 or 2, but not both. P12 on the other hand is the prior
probability of colocalization.
Note
----
Implements `.f`, from `coloc.process` from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L419>`__.
Implements `combine.abf` coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/claudia.R#L115>`__.
Arguments
---------
ln_abf1, ln_abf2 : pd.Series
The log approximate base factor for phenotypes 1 or 2.
p1, p2, p12 : floats
The prior probabilities that a random variant is associated with either
trait 1 or 2, or both, respectively.
extra_info : boolean, default True
If we want addiontal detail on the best variants per hypothesis.
Similar to the information provided by R coloc.detail and the `.f`
function in R coloc.process.
masking_ln_abf3, masking_ln_abf4 : pd.Series
the log approximate Bayes factors when using `method=mask`
Returns
-------
Returns a dictionary of posterior probabilities
"""
# @@@ CHECKS
# TODO
# @@@ actual calculations
lsum = ln_abf1 + ln_abf2
# the log posterior probability
lh0_abf = 0
lh1_abf = np.add(np.log(p1), _logsum(ln_abf1))
lh2_abf = np.add(np.log(p2), _logsum(ln_abf2))
# NOTE the _logdiff here is specific to `combine.abf` and not repeated
# below when extra_info == True
lh3_abf = np.add(
np.add(np.log(p1), np.log(p2)),
_logdiff(_logsum(ln_abf1) + _logsum(ln_abf2), _logsum(lsum))
)
lh4_abf = np.add(np.log(p12), _logsum(lsum))
# returning a summary
temp_tuple = (lh0_abf, lh1_abf, lh2_abf, lh3_abf, lh4_abf)
denom = _logsum(temp_tuple)
# if extra_info is True
if extra_info == True:
best1 = ln_abf1.idxmax()
best2 = ln_abf2.idxmax()
best4 = lsum.idxmax()
# NOTE _logsum(ln_abf1) + _logsum(ln_abf2)) == lbf3
# NOTE this does NOT include an _logdiff following `coloc.process`
# see https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L423
lh3_abf = np.add(np.add(np.log(p1), np.log(p2)),
_logsum(ln_abf1) + _logsum(ln_abf2))
if masking_ln_abf3.size != 0:
lh3_abf = np.add(np.add(np.log(p1), np.log(p2)),
_logsum(masking_ln_abf3))
if masking_ln_abf4.size != 0:
lh4_abf = np.add(np.log(p12), _logsum(masking_ln_abf4))
best4 = masking_ln_abf4.idxmax()
temp_tuple = (lh0_abf, lh1_abf, lh2_abf, lh3_abf, lh4_abf)
denom = _logsum(temp_tuple)
pp_abf = np.exp(temp_tuple - denom)
# NOTE best? are the variants that match to their respective hypotheses.
results = {c.COLOC_PPH0: pp_abf[0],
c.COLOC_PPH1: pp_abf[1],
c.COLOC_PPH2: pp_abf[2],
c.COLOC_PPH3: pp_abf[3],
c.COLOC_PPH4: pp_abf[4],
c.COLOC_BEST1: best1,
c.COLOC_BEST2: best2,
c.COLOC_BEST4: best4
}
else:
# if extra_info is False
posterior = np.exp(np.subtract(temp_tuple, denom))
results = {c.COLOC_PPH0: posterior[0],
c.COLOC_PPH1: posterior[1],
c.COLOC_PPH2: posterior[2],
c.COLOC_PPH3: posterior[3],
c.COLOC_PPH4: posterior[4]
}
return results
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def vestgeno_1_ctl(f: pd.Series):
'''
Helper function to estimate the relative allele frequency of genotypes
0, 1, 2 in control subjects.
Note
----
Implements `vestgeno.1.ctl`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L543>`__.
Arguments
---------
f : pd.Series
A vector of minor allele frequencies.
Returns
-------
A pd.DF with the genotype frequencies.
'''
# @@@ check input
is_type(f, pd.Series)
if (0 > f.min()) or (f.max() > 0.50):
raise ValueError('Supplied values should range between {} and {}'.
format(0, 0.50))
# @@@ calulcating the MAF in controls
dt = pd.DataFrame(np.array([(1-f)**2,
2*f*(1-f),
f**2 ])).T
return dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def vestgeno_1_cse(G_0: pd.DataFrame, b: np.ndarray):
'''
Helper function to estimate the relative allele frequency of genotypes
0, 1, 2 in case subjects.
Implements `vestgeno.1.cse`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L559>`__.
Arguments
---------
G_0 : pd.DF
The control group genotype frequencies;
Results from `vestgeno_1_ctl`. With:
-. Column 0 recording the frequency of genotype 0,
-. Column 1 the frequency of genotype 1,
-. Column 2 the frequency of genotype 2.
b : pd.Series
A column vector of effect estimates.
Returns
-------
A pd.DF
'''
# @@@ check input
is_type(G_0, pd.DataFrame)
is_type(b, pd.Series)
b = b.to_numpy()
# @@@ Do the actual calculations
g_0 = 1
g_1 = np.exp(b - np.log(G_0.iloc[:,0]/G_0.iloc[:,1]))
g_2 = np.exp(2*b - np.log(G_0.iloc[:,0]/G_0.iloc[:,2]))
dt = pd.DataFrame({"g_0" : g_0,
"g_1" : g_1,
"g_2" : g_2})
dt1 = dt.div(g_0+g_1+g_2, axis='rows')
return dt1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def VMAF_cc(f_0: pd.Series, b: pd.Series, N_0: float, N_1: float) -> np.ndarray:
'''
Estimate the MAF for a case-control analysis, based on MAF from controls,
the log OR (`b`) and the number of controls `N_0` and cases `N_1`.
Implements `VMAF.cc`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L566>`__.
Arguments
---------
f_O : pd.Series
The minor allele frequency estimated from control subjects.
b : pd.Series
A vector of log odds ratio.
N_0 : float
The number of control subjects.
N_1 : float
The number of subjects with an event (the cases).
Returns
-------
An np.ndarray
'''
# @@@ check input
is_type(f_0, pd.Series)
is_type(b, pd.Series)
is_type(N_0, float)
is_type(N_1, float)
# @@@ the calculations
G_0 = vestgeno_1_ctl(f_0)*N_0
G_1 = vestgeno_1_cse(G_0, b)*N_1
# set same names
G_0.columns = G_1.columns
G = G_0.add(G_1)
E_2 = np.divide((G @ np.array([0,1,4]).reshape(-1,1)), (N_0+N_1))
E_1 = np.divide((G @ np.array([0,1,2]).reshape(-1,1)), (N_0+N_1))
cc = np.divide((E_2 - E_1**2) * (N_0+N_1), (N_0+N_1-1))
# to numpy
# NOTE this may break stuff - was pd.DF before
if cc.shape[1] == 1:
cc = cc.iloc[:, 0].to_numpy()
else:
raise IndexError('`cc` has more than one column')
return cc
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def map_cond(dt,
corr_matrix,
YY,
sigsnps=[],
r2_threshold = 0.80,
effect_type='quant',
variant_id=Coloc_Expected.VARIANT_ID,
effect_size=c.EFFECT_SIZE.name,
varbeta=c.COLOC_VARBETA,
event_rate=c.COLOC_EVENT_RATE,
sample_size=c.COLOC_SAMPLE_SIZE,
minor_allele_freq=Coloc_Expected.MINOR_ALLELE_FREQ,
label=""):
"""
Internal helper function for finemap.signals. Finds the next most significant
variant, after conditioning on content of `sigsnps`.
Implements `map_cond`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L256>`__.
Arguments
---------
dt : pd.DF,
A variant dataframe with columns ...
corr_matrix : pd.DF,
The correlation matrix of the varaints in `dt`.
YY : float,
The residual sum of squares.
sigsnps : list,
A list of row indices of the varaint(s) we want to condition on.
effect_size : str,
The name of the effect size column.
r2_threshold : float, default 0.80,
Pruning variants with a squared pairwise correlation above
`r2_threshold`. This is done to improve model stability, nothing else.
effect_type : str, default `quant`
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
varbeta : str,
The name of the variance of the effect size column
standard_deviation : str,
The column name of the gwas trait variance (will be estimated if
absent). Only relevant for quantitative trait GWAS.
sample_size : str,
The column name of the GWAS sample size.
minor_allele_freq : str,
The column name of the MAF.
The name of the effect size column.
label : str
...
Returns
-------
A pd.DF with the conditionally significant variants.
"""
# @@@ check input
is_type(dt, pd.DataFrame)
core_cols = [variant_id, effect_size, varbeta,
sample_size, minor_allele_freq]
are_columns_in_df(dt, core_cols)
is_type(effect_type, str)
# check effect_type
if not effect_type in Coloc_Expected.TRAIT_TYPES:
raise AttributeError("`effect_type` should be either {} or {}.".\
format(*Coloc_Expected.TRAIT_TYPES))
# Check if we have a correlation matrix
if corr_matrix is None:
raise AttributeError(" The `corr_matrix` is missing")
if corr_matrix.empty == True:
raise AttributeError(" The `corr_matrix` is missing")
is_type(corr_matrix, pd.DataFrame)
# check if corr_matrix and dt have the same entries
# and align dt and corr_matrix entries
_process_corr_mat(data=dt, corr_matrix=corr_matrix,
index_col=Coloc_Expected.VARIANT_ID
)
# @@@ The calculations
# if there are no pre-defined varaints
# take max unconditional abs(z)
if len(sigsnps) == 0:
Z = np.divide(dt[effect_size], np.sqrt(dt[varbeta])).to_numpy()
wh = np.argmax(np.abs(Z))
sig_dt = pd.DataFrame([Z[wh]], columns=[dt[variant_id].iloc[wh]])
elif len(sigsnps) > 0:
# place holder
if effect_size == 'quant':
event_rate = c.COLOC_EVENT_RATE
est = est_cond(dt=dt, corr_matrix=corr_matrix, YY=YY,
sigsnps=sigsnps, effect_type=effect_type,
r2_threshold=r2_threshold,
variant_id=variant_id,
effect_size=effect_size,
varbeta=varbeta,
event_rate=event_rate,
sample_size=sample_size,
minor_allele_freq=minor_allele_freq,
label=label
)
Z = np.divide(est[c.EFFECT_SIZE.name+"_"+label],
np.sqrt(est[c.COLOC_VARBETA+"_"+label]))
wh = np.argmax(np.abs(Z))
# NOTE the iloc is a dirty fix - should go back to this and actually
# address the index
sig_dt = pd.DataFrame([Z.iloc[wh]],
columns=[est.iloc[wh][c.UNIVERSAL_ID.name]])
# return
return sig_dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def est_cond(dt,
corr_matrix,
YY,
sigsnps,
effect_type="quant",
xtx=None,
r2_threshold = 0.80,
effect_size=c.EFFECT_SIZE.name,
variant_id=Coloc_Expected.VARIANT_ID,
varbeta=c.COLOC_VARBETA,
event_rate=c.COLOC_EVENT_RATE,
sample_size=c.COLOC_SAMPLE_SIZE,
minor_allele_freq=Coloc_Expected.MINOR_ALLELE_FREQ,
label=""):
"""
Internal helper function for est_all_cond. Estimates conditional GWAS summary
statistics.
Implements `est_cond`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L115>`__.
Arguments
---------
dt : pd.DF,
A variant dataframe with columns ...
corr_matrix : pd.DF,
The correlation matrix of the variants in `dt`.
YY : float,
The residual sum of squares.
sigsnps : list,
A list of row indices of the varaint(s) we want to condition on.
effect_type : str, default `quant`
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
xtx : pd.DF?, default None,
The "Gram matrix" matrix -- essentially the variance-covariance matrix
obtained by X.transpose(). Will be internally calculated when set to
`None`, based on the corr_matrix, MAF, beta and sample size (the last
three are supplied through `dt`).
r2_threshold : float, default 0.90,
Pruning variants with a squared pairwise correlation above
`r2_threshold`.
effect_size : str,
The name of the effect size column.
varbeta : str,
The name of the variance of the effect size column
standard_deviation : str,
The column name of the gwas trait variance (will be estimated if
absent). Only relevant for quantitative trait GWAS.
sample_size : str,
The column name of the GWAS sample size.
minor_allele_freq : str,
The column name of the MAF.
label : str
...
Returns
-------
pd.DF with conditional GWAS estimates
"""
# @@@ check input
is_type(dt, pd.DataFrame)
core_cols = [variant_id, effect_size, varbeta,
sample_size, minor_allele_freq]
are_columns_in_df(dt, core_cols)
# Check if we have a correlation matrix
if corr_matrix is None:
raise AttributeError(" The `corr_matrix` is missing")
if corr_matrix.empty == True:
raise AttributeError(" The `corr_matrix` is missing")
is_type(corr_matrix, pd.DataFrame)
# check if corr_matrix and dt have the same entries
# and align dt and corr_matrix entries
_process_corr_mat(data=dt, corr_matrix=corr_matrix,
index_col=Coloc_Expected.VARIANT_ID
)
# Guard the matrix-inverting path (np.linalg.solve, below) against
# undefined-LD NaN rows and non-positive-definite matrices. This runs
# AFTER _process_corr_mat has aligned dt/corr_matrix and BEFORE the XX /
# solve computation. It is matrix-only (no instrument matching), so it is
# safe for coloc's dt whose variant IDs live in a column with a RangeIndex.
error_on_non_finite_pd_corr_matrix(corr_matrix)
# further input
is_type(YY, float)
is_type(sigsnps, list)
if not effect_type in Coloc_Expected.TRAIT_TYPES:
raise AttributeError("`effect_type` should be either {} or {}.".\
format(*Coloc_Expected.TRAIT_TYPES))
# @@@ The actual calculations
# Extracting dt.variants we want to condition one
nuse = []
for snp in sigsnps:
nuse.append(dt[variant_id].tolist().index(snp))
# to find new beta for, conditional on nuse, excluding variants in r-squared
# > r2_threshold. With nuse because small inaccuracies can blow up, and
# we generally don't expect a second detectable signal in
# such high LD with a first signal
ld_mask = np.max(np.power(corr_matrix.iloc[nuse], 2), axis=0) > r2_threshold
nuse_ld = np.flatnonzero(np.asarray(ld_mask)).tolist()
# adding high_ld variants
nuse_comb = list(set(nuse + nuse_ld))
nuse_comb_nam = dt[variant_id].iloc[nuse_comb].to_list()
# position
use_names = [x for x in dt[variant_id].tolist() if x not in nuse_comb_nam]
# use_names = [x for x in dt[variant_id].tolist() if x not in sigsnps]
use = []
# get position
for snp in use_names:
use.append(dt[variant_id].tolist().index(snp))
# getting the variance-covariance matrix
if xtx is None:
# depending on the trait, do
if effect_type == "quant":
# if quant, assume MAF is correct estimate for sample
# expected variance of X, h_buf in GCTA, Dj/N in Yang et al
VX = 2 * dt[minor_allele_freq] * (1 - dt[minor_allele_freq])
A1 = np.array(VX).reshape(1,-1)
A2 = np.array(VX).reshape(-1,1)
# E(XtX/Ni)
P = np.sqrt(A2 @ A1)
elif effect_type == "bin":
# Correcting the MAF
# NOTE Taking the average event_rate and sample_size
VW = VMAF_cc(f_0 = dt[minor_allele_freq],
b = dt[effect_size],
N_0 = (1-dt[event_rate].mean()) * dt[sample_size].mean(),
N_1 = dt[event_rate].mean() * dt[sample_size].mean(),
)
A1 = np.array(VW).reshape(-1, 1)
A2 = np.array(VW).reshape(1,-1)
# np.sqrt(A1 @ A2)
P = np.sqrt(A1 @ A2)
else:
raise AttributeError("effect_type does not exist, either supply \
effect type as `quant` or `bin`, or supply a variance-covariance matrix through \
xtx")
else:
XX = xtx
# calculate XX
XX = dt[sample_size].mean() * corr_matrix * P
# The variance vectors
D = np.diag(XX)
D_1 = D[nuse]
D_2 = D[use]
# the covariance matrices (_1, _2, will be vectors)
XX_1 = XX.iloc[nuse, nuse]
XX_2 = XX.iloc[use, use]
XX_12 = XX.iloc[nuse, use]
XX_21 = XX.iloc[use, nuse]
# accommodate nuse/use of length 1 or mor
D_1 = np.array(D_1).reshape(1,1) if len(nuse) == 1 else np.diag(D_1)
D_2 = np.array(D_2).reshape(1,1) if sum(use) == 1 else np.diag(D_2)
# joint effects at nuse, if needed
# (X'X)^-1 @ var @ beta
b_1 = np.linalg.solve(XX_1, np.identity(len(XX_1))) @ D_1 @\
dt[effect_size].iloc[nuse].to_numpy().reshape(-1,1) # single column vector
## conditional effects 2 | 1
b_2 = np.subtract(
# term 1
dt[effect_size].iloc[use].to_numpy(),
# subtract this from term 1
np.divide(
# numerator
XX_21.to_numpy() @ np.linalg.solve(XX_1, np.identity(len(XX_1))) @ D_1 @\
np.array(dt[effect_size].iloc[nuse]),
# denominator
np.diag(D_2)
)
).reshape(1,-1)
## residual and conditional b2 variance
# row vector
# b_1 = b_1.reshape(1,-1)
# column vector
beta_mat = dt[effect_size].iloc[nuse].to_numpy().reshape(-1,1)
# Sc
# NOTE taking the mean -- this will always be a float
Sc = YY - b_1.T @ D_1 @ beta_mat
Sc = (
# numerator
(Sc - b_2 * np.diag(D_2) *
dt[effect_size].iloc[use].to_numpy().reshape(1,-1)) /
# denominator
(dt[sample_size].mean() - len(sigsnps) - 1)
)
# variance b2
vb_2 = np.divide(
# numerator
np.array(Sc) * (
np.diag(D_2) -\
np.diag(XX_21 @ np.linalg.solve(XX_1, np.identity(len(XX_1))) @\
np.array(XX_12)
)
),
# denominator
np.diag(D_2)**2)
# abs here because very occasionally can get a small negative vb2 due to
# approximation. In this case, replace by a positive vb2 of same small
# magnitude
vb_2 = np.abs(vb_2)
# results table
cond_dt = pd.concat(
[
# table 1
pd.DataFrame({
c.UNIVERSAL_ID.name: dt[variant_id].iloc[use].to_list(),
c.EFFECT_SIZE.name+"_"+label: list(itertools.chain.from_iterable(b_2)),
c.COLOC_VARBETA+"_"+label: list(itertools.chain.from_iterable(vb_2))
}) ,
# table 2
# NOTE USING nuse_comb, here to ensure multicolinear variants
# get returned as well
pd.DataFrame({
c.UNIVERSAL_ID.name: dt[variant_id].iloc[nuse_comb].to_list(),
c.EFFECT_SIZE.name+"_"+label: 0*len(nuse_comb),
c.COLOC_VARBETA+"_"+label: dt[varbeta].iloc[nuse_comb].to_list()
})
]
)
# reset index
cond_dt.reset_index(drop=True, inplace=True)
# returns stuff
return cond_dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Done
def find_best_signal(dt: pd.DataFrame,
variant_id=Coloc_Expected.VARIANT_ID,
effect_size=c.EFFECT_SIZE.name,
varbeta=c.COLOC_VARBETA,
):
"""
Internal helper function to select the most significant variant.
Implements `find.best.signal`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L285>`__.
Arguments
---------
dt : pd.DF,
A variant dataframe with columns ...
effect_size : str,
The name of the effect size column.
varbeta : str,
The name of the variance of the effect size column
Returns
-------
pd.DF with most significant variant
"""
# @@@ Cehck input
is_type(dt, pd.DataFrame)
are_columns_in_df(dt, [effect_size, varbeta])
# @@@ caclculate z-statistics
# NOTE HAVE removed the YY part, as far as I can tell it does nothing.
Z = dt[effect_size]/dt[varbeta].pow(1/2)
# largest Z
wh = np.argmax(np.abs(np.array(Z)))
# return
# NOTE perhaps this works better as a dictionary?
# NOTE 2 - what if there are ties?
Zs = Z.iloc[wh]
nam = [dt[variant_id].iloc[wh]]
return_dt = pd.DataFrame([Zs], columns=nam)
return return_dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def est_all_cond(dt: pd.DataFrame,
corr_matrix: Union[pd.DataFrame, type(None)],
fm: pd.DataFrame,
mode="iterative",
effect_type="quant",
r2_threshold = 0.80,
variant_id=Coloc_Expected.VARIANT_ID,
effect_size=c.EFFECT_SIZE.name,
varbeta=c.COLOC_VARBETA,
event_rate=c.COLOC_EVENT_RATE,
sample_size=c.COLOC_SAMPLE_SIZE,
minor_allele_freq=Coloc_Expected.MINOR_ALLELE_FREQ,
standard_deviation=c.COLOC_STANDARD_DEVIATION,
label="") -> dict:
"""
Implements `est_all_cond`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L671>`__.
Arguments
---------
dt : pd.DF,
A variant dataframe with columns ...
corr_matrix : pd.DF,
The correlation matrix of the variants in `dt`.
fm : pd.DF,
The best variants selected through `find_best_signal`.
mode : float, default `iterative`,
`iterative` or `allbutone`; to either successively condition on signals
find all putative signals and condition on all but one of them in each
analysis.
effect_type : str, default `quant`
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
r2_threshold : float, default 0.80,
Pruning variants with a squared pairwise correlation above
`r2_threshold`.
effect_size : str,
The name of the effect size column.
varbeta : str,
The name of the variance of the effect size column
event_rate : str,
The column name of the proportion of cases in the source GWAS.
sample_size : str,
The column name of the GWAS sample size.
minor_allele_freq : str,
The column name of the MAF.
standard_deviation : str,
The column name of the gwas trait variance (will be estimated if
absent). Only relevant for quantitative trait GWAS.
label : str
...
Returns
-------
A dictionary with one unconditional and x many conditional GWAS estimates
"""
# @@@ Check input
# Check if core data is present
is_type(dt, pd.DataFrame)
is_type(fm, pd.DataFrame)
is_type(corr_matrix, (pd.DataFrame, type(None)))
is_type(effect_type, str)
is_type(mode, str)
is_type(r2_threshold, float)
# Check if effect types are correct
if not effect_type in Coloc_Expected.TRAIT_TYPES:
raise AttributeError("`effect_type` should be either {} or {}.".\
format(*Coloc_Expected.TRAIT_TYPES))
# check modes
if not mode in Coloc_Expected.MODES:
raise AttributeError(
"`modes` should be a str: {}".format(Coloc_Expected.MODES)
)
# Checking columns content
cor_cols = [variant_id, effect_size, varbeta]
are_columns_in_df(dt, cor_cols)
# @@@ the actual calculations
if effect_type == 'bin':
dt, _ = bin2lin(dt,
variant_id=variant_id,
effect_size=effect_size,
varbeta=varbeta,
event_rate=event_rate,
sample_size=sample_size,
minor_allele_freq=minor_allele_freq
)
# get z-statistics
dt[c.COLOC_NORMAL_DEVIATE + label] = dt[effect_size]/dt[varbeta].pow(1/2)
# calculate YY
YY = est_YY(dt=dt, effect_type=effect_type,
standard_deviation=standard_deviation,
sample_size=sample_size,
event_rate=event_rate)
# replaced the `sigs` function with the following
snps = fm.columns.tolist()
# snps = ['snp_1', 'snp_2', 'snp_3']
if mode == "allbutone":
sigs = list(itertools.combinations(snps, len(snps)-1))
sigs.reverse()
elif mode == "iterative":
sigs = [tuple(snps[:i]) for i in range(len(snps))]
# getting the conditional dataframe
# NOTE USING ORDERED_DICT
cond_dt = OrderedDict()
for sigsnps in sigs:
# This will be True if sigsnps are not in sigs
if not sigsnps:
# making a tuple with the same structure as sigsnps
temp = dt.copy()
# making sure the index is correct
temp.set_index(temp[variant_id], drop=False, inplace=True)
temp.index.name = None
cond_dt[(None, None)] = temp
del temp
else:
# NOTE mapping sigsnps to a list in `est_cond`
temp = est_cond(dt=dt, corr_matrix=corr_matrix, YY=YY,
sigsnps=list(sigsnps),
effect_type=effect_type,
r2_threshold=r2_threshold,
variant_id=variant_id,
effect_size=effect_size,
varbeta=varbeta,
event_rate=event_rate,
sample_size=sample_size,
minor_allele_freq=minor_allele_freq,
label=label
)
# ensuring the index is correct
temp.set_index(temp[variant_id], drop=False, inplace=True)
temp.index.name = None
cond_dt[sigsnps] = temp
del temp
# return
return cond_dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def map_mask(dt: pd.DataFrame,
corr_matrix: pd.DataFrame,
r2thr=0.01,
sigsnps=[],
variant_id=Coloc_Expected.VARIANT_ID,
effect_size=c.EFFECT_SIZE.name,
varbeta=c.COLOC_VARBETA,
zstatistic=c.COLOC_NORMAL_DEVIATE
):
"""
Internal helper function for finemap.signals. Finds the next most
significant variants, masking/removing variants in `sigsnps`.
Implements `map_mask`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L72>`__.
Arguments
---------
dt : pd.DF,
A coloc dataframe with either the exposure or outcome data.
corr_matrix : pd.DF,
The correlation matrix of the variants in `dt`.
r2thr : float, default 0.01
Mask variants with a squared correlation with any of the `sigsnps`
larger than `r2thr`.
sigsnps : list,
A list of row indices of the varaint(s) we want to condition on.
effect_size : str,
The name of the effect size column.
varbeta : str,
The name of the variance of the effect size column.
zstatistic : str,
The name of the z-statistic column. Will be estimates if effect_size
and varbeta are supplied and the column is not present in `dt`.
Returns
-------
A pd.DF with the significant variants
"""
# @@@ check input
is_type(dt, pd.DataFrame)
core_cols = [variant_id, varbeta ]
are_columns_in_df(dt, core_cols)
# Check if we have a correlation matrix
if corr_matrix is None:
raise AttributeError(" The `corr_matrix` is missing")
if corr_matrix.empty == True:
raise AttributeError(" The `corr_matrix` is missing")
is_type(corr_matrix, pd.DataFrame)
# check if corr_matrix and dt have the same entries
are_Series_equal(dt[variant_id], corr_matrix.index.to_series(),
objects_names=['dt', 'corr_matrix']
)
# Files have the same entries, now make sure they are in the same order
dt.sort_values(by=[variant_id], ascending=True, inplace=True)
corr_matrix.sort_index(axis=0, ascending=True, inplace=True)
corr_matrix.sort_index(axis=1, ascending=True, inplace=True)
# @@@ actual calculations
# calculate z-statistic
if not zstatistic in dt.columns:
Z = np.divide(dt[effect_size], np.sqrt(dt[varbeta])).to_numpy()
else:
Z = dt[zstatistic].to_numpy()
uni_id = np.array(dt[variant_id])
use = np.repeat(True, dt.shape[0])
# if there are sigsnps, do
if len(sigsnps) > 0:
expectedz = np.repeat(0, dt.shape[0])
# select just the column for the sigsnps
a = np.array(corr_matrix[sigsnps]).T
friends = np.any(np.abs(a) > np.sqrt(r2thr), axis=0)
use = np.logical_xor(use, friends)
# NOTE check if this needs to be an empty df with the same rows and columns
# It seems to only fill the column name by uni_id, so presumably you
# can keep it empty.
if np.any(use) == False:
sig_dt = pd.DataFrame()
else:
# otherwise go on
idx = np.flatnonzero(
dt[variant_id].isin(sigsnps).to_numpy()
)
expectedz = np.array(corr_matrix[sigsnps]) @ Z[idx]
zdiff = np.abs(Z.T[use]) - np.abs(expectedz.T[use])
wh = np.argmax(zdiff)
sig_dt = pd.DataFrame([Z.T[use][wh]],
columns=[uni_id[use][wh]])
elif len(sigsnps) == 0:
zdiff = np.abs(Z.T)
wh = np.argmax(zdiff)
sig_dt = pd.DataFrame([Z[wh]],
columns=[dt[variant_id].iloc[wh]])
# return
return sig_dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def wmean(estimates: Union[pd.DataFrame, pd.Series, np.ndarray],
weights: Union[pd.DataFrame, pd.Series, np.ndarray, float]
) -> np.ndarray:
'''
Weighted mean -- fixed effect meta-analysis.
Arguments
---------
estimates, weights : pd.DataFrame or pd.Series
Returns
-------
A np.array
'''
# @@@ check
is_type(estimates, (pd.DataFrame, pd.Series, np.ndarray))
is_type(weights, (pd.DataFrame, pd.Series, np.ndarray, float))
# @@@ estimates
# remove indices if needed
if isinstance(estimates, ( pd.DataFrame, pd.Series )):
estimates = estimates.to_numpy()
if isinstance(weights, ( pd.DataFrame, pd.Series )):
weights = weights.to_numpy()
# make 2d
if len(weights.shape) == 1:
weights = weights.reshape(-1,1)
try:
# colsums
res = np.sum(estimates * weights, axis=0)/np.sum(weights, axis=0)
except ValueError:
# print('testing here')
res = np.sum(estimates * weights, axis=0)/np.sum(weights)
# return
return res
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def bin2lin(dt,
variant_id=Coloc_Expected.VARIANT_ID,
effect_size=c.EFFECT_SIZE.name,
varbeta=c.COLOC_VARBETA,
event_rate=c.COLOC_EVENT_RATE,
sample_size=c.COLOC_SAMPLE_SIZE,
minor_allele_freq=Coloc_Expected.MINOR_ALLELE_FREQ):
"""
Estimate beta and varbeta if a linear regression had been run on a binary
outcome, given log odds ratio and their variance + MAF in controls
Implements `bin2lin` and `bin2lin.nometa`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L633>`__.
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L574>`__.
Arguments
---------
dt : pd.DF,
A variant dataframe with columns ...
effect_size : str,
The name of the effect size column.
varbeta : str,
The name of the variance of the effect size column
event_rate : str,
The column name of the proportion of cases in the source GWAS.
sample_size : str,
The column name of the GWAS sample size.
minor_allele_freq : str,
The column name of the MAF.
Returns
Unpacks an updated `dt` and `linear_approx`. Here `dt` replaces the orignal
logOR and logOR variance with the estimated mean difference and its variance.
The `linear_approx` compares the z-statistics from both the mean difference
and the logOR, should be close to 1.
"""
# @@@ check input
is_type(dt, pd.DataFrame)
core_cols = [variant_id, effect_size, varbeta,
minor_allele_freq, event_rate, sample_size]
are_columns_in_df(dt, core_cols)
# @@@ estimate the genotype frequencies
# control frequency
g0 = vestgeno_1_ctl(dt[minor_allele_freq])
# case frequency
g1 = vestgeno_1_cse(g0, dt[effect_size])
# expected genotypes
geno_mat = np.array([[0,1,2]])
# expected proportion of case and control subjects
ex1 = geno_mat @ g1.transpose()
ex0 = geno_mat @ g0.transpose()
# squared expression
vx1 = geno_mat**2 @ g1.transpose()
vx0 = geno_mat**2 @ g0.transpose()
# combined
ex = (np.array([1 - dt[event_rate].mean()]).T @ ex0 +
np.array([dt[event_rate].mean()]).T @ ex1)
vx = (np.array([1 - dt[event_rate].mean()]).T @ vx0 +
np.array([dt[event_rate].mean()]).T @ vx1 - ex**2)
# binomial variance
vy = np.array([dt[event_rate].mean() * (1 - dt[event_rate].mean())])
# centred
vxy = vy.transpose() @ ( ex1 - ex0 )
# @@@ estimating varbeta and beta on the linear scale
# essentially getting mean differences
varbeta_mat = (vy/vx - vxy**2/vx**2) / (dt[sample_size].mean() - 1)
# row vector
# varbeta_mat = varbeta_mat.to_numpy().reshape(-1, 1)
# storing original values
beta_org = dt[effect_size].copy()
varbeta_org = dt[varbeta].copy()
if len(vxy.shape) == 1:
dt[varbeta] = varbeta_mat.to_list()
dt[effect_size] = (vxy/vx).to_list()
else:
warnings.warn("Performing meta-analysis")
# meta-analysis, axis=1 is colsum
dt[varbeta] = 1/np.sum(1/varbeta_mat, axis=1).tolist()
dt[effect_size] = wmean(vxy/vx, 1/varbeta_mat).tolist()
# @@@ diagnostic results
zstat = dt[effect_size]/np.sqrt(dt[varbeta])
zstat_org = beta_org/np.sqrt(varbeta_org)
fit = smf.ols(formula='zstat ~ zstat_org - 1',
data = pd.DataFrame(
{'zstat': zstat, 'zstat_org': zstat_org}
), missing = 'drop').fit()
linear_approx = fit.params['zstat_org']
if not (0.90 <= linear_approx <= 1.1):
warnings.warn("Poor linear approximation: {}".format(linear_approx))
# returns
return dt, linear_approx
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def est_YY(dt : pd.DataFrame,
effect_type: list,
variant_id=Coloc_Expected.VARIANT_ID,
standard_deviation=c.COLOC_STANDARD_DEVIATION,
sample_size=c.COLOC_SAMPLE_SIZE,
event_rate=c.COLOC_EVENT_RATE
) -> float:
"""
Calculates the squared residuals.
Implements code chunck from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L293>`__.
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L355>`__.
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L680>`__.
Arguments
---------
dt : pd.DF,
A coloc dataframe.
effect_type : str, default `quant`,
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
standard_deviation : str,
The column name of the gwas trait variance (will be estimated if
absent). Only relevant for quantitative trait GWAS.
sample_size : str,
The column name of the GWAS sample size.
event_rate : str,
The column name of the proportion of cases in the source GWAS. Only
applicable for a binary trait GWAS.
Returns
-------
A float
"""
# @@@ Check if core data is present
is_type(dt, pd.DataFrame)
core_cols = [variant_id, sample_size ]
are_columns_in_df(dt, core_cols)
if not effect_type in Coloc_Expected.TRAIT_TYPES:
raise AttributeError("`effect_type` should be either {} or {}.".\
format(*Coloc_Expected.TRAIT_TYPES))
# @@@ performing the actual calulcations
if effect_type == "quant":
# check the correct columns are present
are_columns_in_df(dt, [sample_size, standard_deviation])
# do the actual calculations
# NOTE taking the mean
YY = dt[sample_size].mean() * dt[standard_deviation].mean()**2
elif effect_type == "bin":
# check input
are_columns_in_df(dt, [sample_size, event_rate])
# NOTE taking the mean here, reflecting R coloc
# NOTE 2, in R `find.best.signal` divides YY by N, however
# the value is never returned - dead code, not used here.
YY = (dt[sample_size].mean() * dt[event_rate].mean() *
(1 - dt[event_rate].mean()))
# return
return YY
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def finemap_signals(dt: pd.DataFrame,
corr_matrix=None,
method="single",
r2thr=0.01,
sigsnps=[],
pthr=1e-6,
maxhits=3,
effect_type="quant",
variant_id=Coloc_Expected.VARIANT_ID,
effect_size=c.EFFECT_SIZE.name,
varbeta=c.COLOC_VARBETA,
zstatistic=c.COLOC_NORMAL_DEVIATE,
standard_deviation=c.COLOC_STANDARD_DEVIATION,
sample_size=c.COLOC_SAMPLE_SIZE,
minor_allele_freq=Coloc_Expected.MINOR_ALLELE_FREQ,
event_rate=c.COLOC_EVENT_RATE,
verbose=True
):
"""
This is an analogue to finemap.abf, adapted to find multiple signals where
they exist, via conditioning or masking - ie a stepwise procedure
Finemap multiple signals in a single dataset
Internal helper function for finemap.signals. Finds the next most
significant variants, masking/removing variants in `sigsnps`.
Implements `finemap.signals`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L328>`__.
Arguments
---------
dt : pd.DF,
A coloc dataframe.
corr_matrix : pd.DF,
The correlation matrix of the variants in `dt`. Only relevant when
method == 'cond'.
r2thr : float, default 0.01
Mask variants with a squared correlation with any of the `sigsnps`
larger than `r2thr`.
method : str, default `single`,
Use `single` to perform coloc assuming there is a single causal variant,
Use `cond` to perform a stepwise conditional analysis,
Use `mask` to perform the stepwise selection algorithm using masking
instead.
sigsnps : list,
A list of row indices of the varaint(s) we want to condition on.
NOTE:-> this does not seem to do anything in the R-version.
pthr : float, default 1 \times 10^{-6},
Stops when p-values are larger than `pthr`, where input should be
defined between 0 and 1.
maxhits : int, default 3,
Stops when the procedures finds `maxhits`.
effect_type : str, default `quant`,
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
effect_size : str,
The name of the effect size column.
varbeta : str,
The name of the variance of the effect size column
zstatistic : str,
The name of the z-statistic column. Will be estimates if effect_size
and varbeta are supplied and the column is not present in `dt`.
standard_deviation : str,
The column name of the gwas trait variance (will be estimated if
absent). Only relevant for quantitative trait GWAS.
sample_size : str,
The column name of the GWAS sample size.
minor_allele_freq : str,
The column name of the MAF.
event_rate : str,
The column name of the proportion of cases in the source GWAS. Only
applicable for a binary trait GWAS.
Returns
-------
a pd.df with the remaining significant variants after masking.
"""
# @@@ Check if core data is present
is_type(dt, pd.DataFrame)
is_type(corr_matrix, (pd.DataFrame, type(None)))
core_cols = [variant_id, sample_size ]
are_columns_in_df(dt, core_cols)
# Checking methods
if not isinstance(method, str):
raise AttributeError( "The method type should be a string")
if not method in Coloc_Expected.METHODS:
raise AttributeError(
"`method` should be a string. Pick from: {}".\
format(Coloc_Expected.METHODS)
)
# Check input for `cond` or 'mask'
if method == "cond" or method =='mask':
# Check if we have an correlation matrix
if corr_matrix is None:
raise AttributeError(" The `corr_matrix` is missing")
# check if corr_matrix and dt have the same entries
# and align dt and corr_matrix entries
_process_corr_mat(data=dt, corr_matrix=corr_matrix,
index_col=Coloc_Expected.VARIANT_ID
)
# Check input for `cond`
if method == "cond":
# Check if MAF is present
are_columns_in_df(dt, minor_allele_freq)
# Is the trait standard deviation present or do we need to calculate it
if not standard_deviation in dt.columns:
dt[standard_deviation] = sd_pheno_approx(variancebeta=dt[varbeta],
maf=dt[minor_allele_freq],
n=dt[sample_size]
)
# Check if the sigsnps and dt have the same entries
# if len(sigsnps) > 0:
# are_Series_equal(dt[variant_id], pd.Series(sigsnps),
# objects_names=['dt', 'sigsnps']
# )
# map p-values to z-values
if 0 <= pthr <= 1:
zthr = scipy.stats.norm.ppf(1- pthr/2)
else:
raise ValueError("Please supply `pthr` as a value between 0 and 1")
# Calculate mean differences
if method == "cond" and effect_type == "bin":
if verbose is True:
warnings.warn("[info] Approximating linear analysis of binary trait")
# bin2lin will also check dt columns are correct.
dt, _ = bin2lin(dt,
variant_id=variant_id,
effect_size=effect_size,
varbeta=varbeta,
event_rate=event_rate,
sample_size=sample_size,
minor_allele_freq=minor_allele_freq
)
# get YY
YY = est_YY(dt=dt, effect_type=effect_type,
standard_deviation=standard_deviation,
sample_size=sample_size,
event_rate=event_rate)
# getting variants hits
hits = pd.DataFrame()
while hits.shape[1] < maxhits:
if method == "mask":
# Masking
newhit = map_mask(dt=dt, corr_matrix=corr_matrix, r2thr=r2thr,
sigsnps=hits.columns.unique().tolist(),
variant_id=variant_id, effect_size=effect_size,
zstatistic=zstatistic,
varbeta=varbeta
)
else:
# getting conditional hits
newhit = map_cond(dt=dt, corr_matrix=corr_matrix, YY=YY,
sigsnps=hits.columns.tolist(),
effect_type=effect_type,
variant_id=variant_id,
effect_size=effect_size,
varbeta=varbeta,
event_rate=event_rate,
sample_size=sample_size,
minor_allele_freq=minor_allele_freq
)
# DO we need to exit?
if newhit is None or newhit.empty:
break
# otherwise get value
# NOTE might change newhit to a dictionary - this is just painfull
value = np.array(newhit)[0][0]
# if not sufficiently significant stop
if np.abs(value) < zthr:
break
hits = pd.concat([hits, newhit], axis=1)
# stop after one hit
if method=="single":
break
return hits
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def gethits(hits, i, mode="iterative"):
"""
Implements `gethits`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L378>`__.
Arguments
---------
hits : list of strings
A list of strings of variant ids
i : float
mode : str, default iterative,
Mode is either `iterative` or `allbutone`.
Returns
-------
A numpy.array
"""
# @@@ checking input
is_type(hits, list)
is_type(i, (np.int32,np.int16, np.int32, np.int64, int,float))
is_type(mode, str)
# check modes
if not mode in Coloc_Expected.MODES:
raise AttributeError(
"`modes` should be a str: {}".format(Coloc_Expected.MODES)
)
# @@@ algorithm
# copy, so we save the original
hits_c = hits.copy()
if mode == "allbutone":# mask everything else
# remove the ith element
del hits_c[i]
ret = hits_c
elif mode == "iterative": # mask everything earlier in list
if i == 0:
return np.empty(0, dtype='U5')
else:
# everything until i (excluding i)
ret = hits_c[:i]
# ret = hits[0:(i)]
# return
return np.setdiff1d(ret, "")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def coloc_process(coloc_res,
hits1: list=None,
hits2: list=None,
corr_exp: pd.DataFrame=None,
corr_out: pd.DataFrame=None,
r2thr=0.01,
p1=10**-4,
p2=10**-4,
p12=5*10**-6,
mode="iterative",
verbose=True):
"""Generates a coloc results table, used by coloc_signal.
Parameters
----------
coloc_res : ColocResults,
ColocResults object output by `coloc_single`.
hits1, hits2 : list of strings,
strings from `gethits` representing the hits in dataset1 and
dataset2. If the array contain more than one lead variant -> masking.
corr_(exp|out) : pd.DF,
The correlation matrix of the variants in `dt`, for the exposure and
outcome data respectively.
r2thr : float, default 0.01
Mask variants with a squared correlation with any of the `sigsnps`
larger than `r2thr`.
p1 : `float`, optional, default: `1E-4`
Prior probability for H1 (variant only associated with trait 1)
p2 : `float`, optional, default: `1E-4`
Prior probability for H2 (variant only associated with trait 2)
p12 : `float`, optional, default: `5E-6`,
Prior probability for H4 (associated with both traits).
mode : float, default: `iterative`,
`iterative` or `allbutone`; to either successively condition on signals
find all putative signals and condition on all but one of them in each
analysis.
Returns
-------
results : `pandas.DataFrame`
The coloc results.
Notes
-----
Implements `finemap.signals`, from coloc R, see the `chris wallace
<https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df8
4d7eba39d92a8c1/R/split.R#L408>`_ implementation.
The H3 probability (causal variants for distinct traits) is calculated
internally from the product of p1 and p2. See supplemental of the
`original publication` <https://journals.plos.org/plosgenetics/article
?id=10.1371/journal.pgen.1004383#s5>`_.
"""
# @@@ Check
is_type(coloc_res, ColocResults)
is_type(corr_exp, (pd.DataFrame, type(None)))
is_type(corr_out, (pd.DataFrame, type(None)))
is_type(hits1, (list, type(None)))
is_type(hits2, (list, type(None)))
is_type(r2thr, float)
is_type(p1, float)
is_type(p2, float)
is_type(p12, float)
is_type(mode, str)
# check modes
if not mode in Coloc_Expected.MODES:
raise AttributeError(
"`modes` should be a str: {}".format(Coloc_Expected.MODES)
)
# check rsquared
if not (0 <= r2thr <= 1):
raise ValueError('Please define r2thr within the range [0, 1]')
# @@@ Actual calculations
# NOTE `.f` is moved to `approx_bf`
# see R coloc `.f`: https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L455
if hits1 is None:
hits1 = ['']
#
if hits2 is None:
hits2 = ['']
# one signal per trait
if len(hits1) <= 1 and len(hits2) <= 1:
res_dict = {
c.COLOC_SUMMARY: coloc_res.results_dataframe(hit1=hits1,
hit2=hits2),
c.COLOC_RESULTS: coloc_res.variant_dataframe(),
c.COLOC_PRIORS: pd.DataFrame(
{
c.COLOC_PRIOR1 : p1,
c.COLOC_PRIOR2 : p2,
c.COLOC_PRIOR12 : p12
}, index=[0]
)
}
# break function
return res_dict
# otherwise mask
# check corr_matrix
if (len(hits1) > 1 or len(hits2) > 1) and (corr_exp is None or corr_out is None):
raise AttributeError('Please supply `corr_matrix` as a pd.DataFrame')
# - allbutone : all but one signal from each trait with >1 signal
# - iterative : successively one more signal from each trait with >1 signal
ldfriends1 = np.power(corr_exp.loc[hits1, :], 2)
ldfriends2 = np.power(corr_out.loc[hits2, :], 2)
i = [idx for idx,value in enumerate(hits1)]
j = [idx for idx,value in enumerate(hits2)]
# procudes an i*j by 2 matrix
todo=pd.DataFrame(itertools.product(i,j))
# interchanges columns, followed by a transpose
todo=todo.iloc[:, [1,0]].T
# getting pd.Series with variant names
newresult = coloc_res.variant_dataframe()[[c.UNIVERSAL_ID.name]].copy()
res = pd.DataFrame()
for r in list(range(0,todo.shape[1])):
i = todo.iloc[1, r]
j = todo.iloc[0, r]
if verbose == True:
print("[info] {0}".format(r+1))
# collecting hits
drop1, drop2 = ([] for i in range(2))
# hit 1
if len(hits1) > 1:
# int(i) to get non-numpy ints
ihits1 = gethits(hits1, int(i), mode)
if len(ihits1) >= 1:
ld_out1 = ldfriends1.loc[ihits1, ]
drop1 = ldfriends1[ld_out1 > r2thr].\
dropna(thresh=1, axis=1).columns.tolist()
# hit 2
if len(hits2) > 1:
ihits2 = gethits(hits2, int(j), mode)
if len(ihits2) >= 1:
# if ihits2 != None:
ld_out2 = ldfriends2.loc[ihits2, ]
drop2 = ldfriends2[ld_out2 > r2thr].\
dropna(thresh=1, axis=1).columns.tolist()
# dropping variants
# finding the union
dropsnps = list(set(drop1+drop2))
if verbose == True:
# COLOR R uses message -> warnings
warnings.warn("[info] dropping {0}/{1}: {2} (hits1) + {3} (hits2)"
"".format(
len(dropsnps),
coloc_res.variant_dataframe().shape[0],
len(drop1),
len(drop2)
)
)
# NOTE returns an empty dict
if len(np.setdiff1d(coloc_res.variant_dataframe()[c.UNIVERSAL_ID.name],
dropsnps)) <=1:
res_dict = {
c.COLOC_SUMMARY: pd.DataFrame(),
c.COLOC_RESULTS: pd.DataFrame(),
c.COLOC_PRIORS: pd.DataFrame(
{
c.COLOC_PRIOR1 : p1,
c.COLOC_PRIOR2 : p2,
c.COLOC_PRIOR12 : p12
}, index=[0]
)
}
return res_dict
### getting results object
# NOTE ADD some of these strings to merit.constant
df = coloc_res.variant_dataframe().copy()
# get ABF
df[Coloc_Expected.EXPOSURE_ABF] = np.where(
df[c.UNIVERSAL_ID.name].isin(drop1), -1.1,
df[Coloc_Expected.EXPOSURE_ABF])
df[Coloc_Expected.OUTCOME_ABF] = np.where(
df[c.UNIVERSAL_ID.name].isin(drop2), -1.1,
df[Coloc_Expected.OUTCOME_ABF])
# isum
df.loc[:,c.COLOC_INTERNAL_SUM] =\
np.where(df[c.UNIVERSAL_ID.name].isin(dropsnps), -1.1,
df[c.COLOC_INTERNAL_SUM].copy()).copy()
my_denom_log_abf = _logsum(df[c.COLOC_INTERNAL_SUM].values)
# new results
sum_temp1 = np.exp(
df[c.COLOC_INTERNAL_SUM].copy().values - my_denom_log_abf
).tolist()
newresult["SNP.PP.H4.row"+str(r)] = sum_temp1
newresult.loc[:,"z.df1.row"+str(r)] = np.where(
df[c.UNIVERSAL_ID.name].isin(drop1), 0,
df[Coloc_Expected.EXPOSURE_Z_STATISTICS].copy()
).copy()
newresult.loc[:, "z.df2.row"+str(r)] = np.where(
df[c.UNIVERSAL_ID.name].isin(drop2), 0,
df[Coloc_Expected.OUTCOME_ABF].copy()
).copy()
# subsetting df3
df3 = pd.DataFrame(itertools.product(
coloc_res.variant_dataframe()[c.UNIVERSAL_ID.name],
coloc_res.variant_dataframe()[c.UNIVERSAL_ID.name]),
columns=["uni_id_exposure", "uni_id_outcome"]
)
df3 = df3.merge(
df[[Coloc_Expected.EXPOSURE_ABF, c.UNIVERSAL_ID.name]],
left_on = "uni_id_exposure", right_on = c.UNIVERSAL_ID.name
)
df3 = df3.merge(
df[[Coloc_Expected.OUTCOME_ABF, c.UNIVERSAL_ID.name]],
left_on = "uni_id_outcome", right_on = c.UNIVERSAL_ID.name
)
df3 = df3[["uni_id_exposure",
"uni_id_outcome",
Coloc_Expected.EXPOSURE_ABF,
Coloc_Expected.OUTCOME_ABF
]]
# subsetting ln_abs_h3
df3["ln_abf_h3"] = df3[Coloc_Expected.EXPOSURE_ABF] +\
df3[Coloc_Expected.OUTCOME_ABF].copy()
df3["ln_abf_h3"] = np.where(
df3["uni_id_exposure"].isin(drop1), 0, df3["ln_abf_h3"])
df3["ln_abf_h3"] = np.where(
df3["uni_id_outcome"].isin(drop2), 0, df3["ln_abf_h3"])
# res object
# df3.ln_abf_h3 and df.isum are empty unless `method == mask`
ret = pd.DataFrame(approx_bf(df.ln_abf_exposure,
df.ln_abf_outcome,
masking_ln_abf3=df3.ln_abf_h3,
masking_ln_abf4=df.isum,
p1=p1,
p2=p2,
p12=p12,
extra_info=True),
index=[0])
# get best varaints
# ret[c.COLOC_BEST1] = hits1[todo.iloc[1,r]]
# ret[c.COLOC_BEST2] = hits2[todo.iloc[0,r]]
# hits
ret[c.COLOC_HIT1] = hits1[todo.iloc[1,r]]
ret[c.COLOC_HIT2] = hits2[todo.iloc[0,r]]
# number of variants
ret['nsnps'] = df.shape[0]
res = pd.concat([res, ret])
# adding hits
res_dict = {
c.COLOC_SUMMARY: res,
c.COLOC_RESULTS: newresult,
c.COLOC_PRIORS: pd.DataFrame(
{
c.COLOC_PRIOR1 : p1,
c.COLOC_PRIOR2 : p2,
c.COLOC_PRIOR12 : p12
}, index=[0]
)
}
return res_dict
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def coloc_single(inst: pd.DataFrame,
p1: float=10**-4,
p2: float=10**-4,
p12: float=5*10**-6,
effect_type: list=["quant","quant"],
variant_id: str=Coloc_Expected.VARIANT_ID,
exposure_effect_size: str=Coloc_Expected.EXPOSURE_EFFECT_SIZE,
outcome_effect_size: str=Coloc_Expected.OUTCOME_EFFECT_SIZE,
exposure_varbeta: str=Coloc_Expected.EXPOSURE_VARBETA,
outcome_varbeta: str=Coloc_Expected.OUTCOME_VARBETA,
exposure_zstatistic: str=Coloc_Expected.EXPOSURE_Z_STATISTICS,
outcome_zstatistic: str=Coloc_Expected.OUTCOME_Z_STATISTICS,
exposure_standard_deviation: str=Coloc_Expected.EXPOSURE_STANDARD_DEVIATION,
outcome_standard_deviation: str=Coloc_Expected.OUTCOME_STANDARD_DEVIATION,
exposure_sample_size: str=Coloc_Expected.EXPOSURE_SAMPLE_SIZE,
outcome_sample_size: str=Coloc_Expected.OUTCOME_SAMPLE_SIZE,
minor_allele_freq: str=Coloc_Expected.MINOR_ALLELE_FREQ,
extra_info: bool=False):
"""
Coloc main function, is similar to coloc.details in R coloc. It started
out as a port of `coloc.abf` (assuming a single causal variant).
Will calculate `sdY`, the phenotypic standard deviation if this is not
provided by `inst` and type == `quant`.
Will compute the z-statistic if not provided from the effect_size and
varbeta.
Similar to `coloc.detail`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L21>`__.
And to `coloc.abf`:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L21>`__.
Arguments
---------
inst : pd.DF
A coloc.DataFrame.
p1, p2, p12 : float,
prior probabilities for H1 (variant only associated with trait 1),
for H2 (only with trait 2), or H4 (associated with both traits).
See supplemental of:
`Here <https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1004383#s5>`__.
effect_type : list of str, default `quant`,
The effect_types or the `_exposure` and `_outcome` data.
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
effect_size : str,
The effect size column names for the `exposure` and `outcome` data.
varbeta : str,
The variance of the effect size column names for the
`exposure` and `outcome` data.
zstatistics : str,
The z-statistics column names for the `exposure` and `outcome` data.
Note these can be readily obtained from the GWAS p-values. Can be
estimated from the effect_size and varbeta columns, if absent.
standard_deviation : str,
The column name of the gwas trait variance (will be estimated if
absent).
sample_size : str,
The column name of the GWAS sample size.
minor_allele_freq : str,
The column name of the MAF.
extra_info : boolean, default False
If additional information should be returned (mimicking coloc.details),
otherwise (defaults) output will be similar to coloc.abf.
Returns
-------
A `ColocResults` object.
"""
# @@@ checks
# core data needed
core_cols = [variant_id, exposure_varbeta, outcome_varbeta,
exposure_zstatistic, outcome_zstatistic, minor_allele_freq ]
# are_columns_in_df(inst, core_cols)
# unpacking effect type
if np.all([ t in Coloc_Expected.TRAIT_TYPES for t in effect_type ]):
type_exposure, type_outcome = effect_type
# if needed estimate sdY
# check if standard deviation is provided
if (not exposure_standard_deviation in inst.columns) and\
(type_exposure == "quant"):
inst[exposure_standard_deviation] =\
sd_pheno_approx(
variancebeta=inst[exposure_varbeta],
maf=inst[minor_allele_freq],
n=inst[exposure_sample_size]
)
if (not outcome_standard_deviation in inst.columns) and\
(type_outcome == "quant"):
inst[outcome_standard_deviation] =\
sd_pheno_approx(
variancebeta=inst[outcome_varbeta],
maf=inst[minor_allele_freq],
n=inst[outcome_sample_size]
)
# set the variance prior
if type_exposure == "quant":
sd_prior_exp = 0.15 * inst[exposure_standard_deviation]
elif type_exposure == "bin":
# for case control (binary data)
sd_prior_exp = 0.20
else:
raise ValueError("Please specify `effect_type` two element list with {0} or {1}".\
format(*Coloc_Expected.TRAIT_TYPES))
if type_outcome == "quant":
sd_prior_out = 0.15 * inst[outcome_standard_deviation]
elif type_outcome == "bin":
# for case-control (binary data)
sd_prior_out = 0.20
else:
raise ValueError("Please specify `effect_type` two element list with {0} or {1}".\
format(*Coloc_Expected.TRAIT_TYPES))
# if needed calculate z normal deviate
if not exposure_zstatistic in inst.columns:
inst[exposure_zstatistic] =\
np.divide( inst[exposure_effect_size], np.sqrt(inst[exposure_varbeta])
)
if not outcome_zstatistic in inst.columns:
inst[outcome_zstatistic] =\
np.divide( inst[outcome_effect_size], np.sqrt(inst[outcome_varbeta])
)
# extract core columns
dt = inst[core_cols].copy()
# Shrinkage factor: ratio of the prior variance to the total variance
dt[Coloc_Expected.EXPOSURE_SHRINKAGE_FACTOR] = \
np.divide(
np.power(sd_prior_exp, 2),
(np.power(sd_prior_exp, 2) + dt[exposure_varbeta])
)
dt[Coloc_Expected.OUTCOME_SHRINKAGE_FACTOR] = \
np.divide(
np.power(sd_prior_out, 2),
(np.power(sd_prior_out, 2) + dt[outcome_varbeta])
)
# appending the log (approximate) Bayes Factor for colocalization
dt[Coloc_Expected.EXPOSURE_ABF] =\
(0.5 *
(
np.log(1 - dt[Coloc_Expected.EXPOSURE_SHRINKAGE_FACTOR]) +
dt[Coloc_Expected.EXPOSURE_SHRINKAGE_FACTOR] *
np.power( dt[Coloc_Expected.EXPOSURE_Z_STATISTICS], 2)
)
)
dt[Coloc_Expected.OUTCOME_ABF] =\
(0.5 *
(
np.log(1 - dt[Coloc_Expected.OUTCOME_SHRINKAGE_FACTOR]) +
dt[Coloc_Expected.OUTCOME_SHRINKAGE_FACTOR] *
np.power( dt[Coloc_Expected.OUTCOME_Z_STATISTICS], 2)
)
)
# INTERNAL SUM Just wrapping in tuple to get better code-alignment
dt[c.COLOC_INTERNAL_SUM] = (dt[Coloc_Expected.EXPOSURE_ABF] +
dt[Coloc_Expected.OUTCOME_ABF]
)
# posterior prob that each SNP is THE causal variant for a shared signal
dt[c.COLOC_PPH4] = np.exp(dt[c.COLOC_INTERNAL_SUM] -
_logsum(dt[c.COLOC_INTERNAL_SUM]))
# coloc summary following classic coloc
if extra_info is False:
summary = approx_bf(ln_abf1=dt[Coloc_Expected.EXPOSURE_ABF],
ln_abf2=dt[Coloc_Expected.OUTCOME_ABF],
p1=p1, p2=p2, p12=p12, extra_info=extra_info
)
# results
res = {c.COLOC_SUMMARY: summary, c.COLOC_VARIANT_RESULTS: dt}
elif extra_info is True:
summary = approx_bf(ln_abf1=dt[Coloc_Expected.EXPOSURE_ABF],
ln_abf2=dt[Coloc_Expected.OUTCOME_ABF],
p1=p1, p2=p2, p12=p12, extra_info=extra_info
)
summary[c.COLOC_BEST1] = dt.loc[ summary[c.COLOC_BEST1], variant_id]
summary[c.COLOC_BEST2] = dt.loc[ summary[c.COLOC_BEST2], variant_id]
summary[c.COLOC_BEST4] = dt.loc[ summary[c.COLOC_BEST4], variant_id]
# results
res = {c.COLOC_SUMMARY: summary, c.COLOC_VARIANT_RESULTS: dt}
else:
raise ValueError("`extra_info` is a boolean")
# returning
return ColocResults(**res)
# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
class ColocResults(object):
"""Return Results object."""
ALL_ARGS = [c.COLOC_SUMMARY, c.COLOC_VARIANT_RESULTS]
# SET_ARGS = ['uni_id', 'b', 'se_fixed', 'pvalue_fixed', 'instrument']
IGNORE_ARGS = [c.COLOC_VARIANT_RESULTS]
# OUT_ARGS = [el for el in ALL_ARGS if el not in IGNORE_ARGS]
OUT_ARGS = [c.COLOC_SUMMARY]
def __init__(self, **kwargs):
"""
Initialise
"""
# Check that we have not supplied any arguments that can't be set
for k in kwargs.keys():
if k not in self.__class__.ALL_ARGS:
raise AttributeError("unrecognised argument '{0}'".format(k))
for s in self.__class__.ALL_ARGS:
try:
setattr(self, s, kwargs[s])
except KeyError:
warnings.warn("argument '{0}' is set to 'None'".format(s))
setattr(self, s, None)
def __repr__(self):
values = ",".join(["{0}={1}".format(s, getattr(self, s))
for s in self.__class__.ALL_ARGS if s not in
self.__class__.IGNORE_ARGS])
return "{0}({1})".format(self.__class__.__name__, values)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@property
def nsnps(self):
"""Return the number of variants."""
return len(self.variant_results[c.UNIVERSAL_ID.name])
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def results_dataframe(self, **kwargs):
"""Return dataframe with summary results."""
pp_dict = [getattr(self, s) for s in self.__class__.OUT_ARGS][0]
# check if is a dict or a pd.DF
# if a pd.DF - assume it is already formatted.
if isinstance(pp_dict, dict):
dt = pd.DataFrame([[float(self.nsnps)] + list(pp_dict.values())],
columns=[c.MR_NSNPS] + list(pp_dict.keys()))
for k, v in kwargs.items():
dt.loc[:, k] = v
elif isinstance(pp_dict, pd.DataFrame):
dt = pp_dict
# return stuff
return dt
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def variant_dataframe(self, **kwargs):
"""Return dataframe with individual variant results."""
return self.variant_results
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# def results_H3_dataframe(self, **kwargs):
# """Return dataframe with individual variant results."""
# return self.results_H3
# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# MAIN ENTRY POINT
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DONE
def coloc_signals(input_data,
corr_matrix,
method=["single", "single"],
mode="iterative",
p1=1e-4,
p2=1e-4,
p12=5*10**-6,
maxhits=3,
r2thr=0.01,
pthr = 1e-06,
verbose=True,
effect_type=["quant","quant"],
variant_id=Coloc_Expected.VARIANT_ID,
exposure_effect_size=Coloc_Expected.EXPOSURE_EFFECT_SIZE,
outcome_effect_size=Coloc_Expected.OUTCOME_EFFECT_SIZE,
exposure_varbeta=Coloc_Expected.EXPOSURE_VARBETA,
outcome_varbeta=Coloc_Expected.OUTCOME_VARBETA,
exposure_zstatistic: str=Coloc_Expected.EXPOSURE_Z_STATISTICS,
outcome_zstatistic: str=Coloc_Expected.OUTCOME_Z_STATISTICS,
exposure_standard_deviation=Coloc_Expected.EXPOSURE_STANDARD_DEVIATION,
outcome_standard_deviation=Coloc_Expected.OUTCOME_STANDARD_DEVIATION,
exposure_sample_size=Coloc_Expected.EXPOSURE_SAMPLE_SIZE,
outcome_sample_size=Coloc_Expected.OUTCOME_SAMPLE_SIZE,
exposure_event_rate=Coloc_Expected.EXPOSURE_EVENT_RATE,
outcome_event_rate=Coloc_Expected.OUTCOME_EVENT_RATE,
minor_allele_freq=Coloc_Expected.MINOR_ALLELE_FREQ
):
"""
The main entry point for coloc to consider multiple causal signals.
Implements `coloc.signals`, from coloc R:
`Here <https://github.com/chr1swallace/coloc/blob/eca81cdd0e5cc9fc0c5dd4df84d7eba39d92a8c1/R/split.R#L760>`__.
Arguments
---------
input_data : pd.DF,
A coloc pd.dataframe, with columns ... . The original coloc `data1` and
`data2` are combined in this single dataframe, with column suffixes
`_exposure` and `_outcome` respectivly.
corr_matrix : pd.DF,
The correlation matrix of the variants in `input_data`. If the expousre
and outcome GWAS are from distinct populations please supply a list of
two matrices [exposure_cor, outcome_cor]. Note is only required when
`method` is `cond` or 'mask'.
method : list of two strings, default `single`,
The method to apply to the `_exposure` and `_outcome` data.
Use `single` to perform coloc assuming there is a single causal
variant - equivalent to the original coloc implementation,
Use `cond` to perform a stepwise conditional analysis,
Use `mask` to perform the stepwise selection algorithm using masking
instead.
mode : float, default `iterative`,
`iterative` or `allbutone`; to either successively condition on signals
find all putative signals and condition on all but one of them in each
analysis.
p1, p2, p12 : float,
prior probabilities for H1 (any random variant only associated with
trait 1), for H2 (only with trait 2), or H4 (associated with both traits).
H3 (causal variants for distinct traits) is simply the product of
p1 and p2. See supplemental of:
`Here <https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1004383#s5>`__.
Note that when the data is sufficiently large, p1 and p2 can often be
estimated from the emperical data under consideration.
pthr : float, default 1 \times 10^{-6},
Stops when p-values are larger than `pthr`, where input should be
defined between 0 and 1.
maxhits : int, default 3,
Stops when the procedures finds `maxhits`.
r2thr : float, default 0.01
Mask variants with a squared correlation with any of the `sigsnps`
larger than `r2thr`.
effect_type : list of str, default `quant`,
The effect_types or the `_exposure` and `_outcome` data.
If the GWAS trait was `quant`: continuous/quantative, or `bin`:
case-control/binary.
Returns
-------
Resutns a pd.DF of coloc results.
"""
# @@@ CHECKING INPUT
# Check if core data is present
is_type(input_data, pd.DataFrame)
is_type(corr_matrix, (pd.DataFrame, type(None), list))
is_type(effect_type, list)
is_type(method, list)
is_type(mode, str)
# expected columns
Exposure_col = [variant_id, exposure_effect_size, exposure_varbeta,
exposure_zstatistic, exposure_standard_deviation,
exposure_sample_size, exposure_event_rate,
minor_allele_freq]
Outcome_col = [variant_id, outcome_effect_size, outcome_varbeta,
outcome_zstatistic, outcome_standard_deviation,
outcome_sample_size, outcome_event_rate,
minor_allele_freq]
# Check if effect types are correct
if np.all([ t in Coloc_Expected.TRAIT_TYPES for t in effect_type ]):
effect_type_exp = effect_type[0]
effect_type_out = effect_type[1]
else:
raise ValueError(
"`effect_type` should be a lis with two entries. Pick from: {}".\
format(Coloc_Expected.TRAIT_TYPES)
)
# Checking methods
if len(method) != 2:
raise AttributeError( "`method` should be a list with two entries")
if not np.all([ t in Coloc_Expected.METHODS for t in method ]):
raise ValueError(
"`method` should be a list with two entries. Pick from: {}".\
format(Coloc_Expected.METHODS)
)
else:
method_exp=method[0]
method_out=method[1]
# check modes
if not mode in Coloc_Expected.MODES:
raise ValueError(
"`modes` should be a str: {}".format(Coloc_Expected.MODES)
)
# Checking columns content
# are_columns_in_df(input_data, list(set(Coloc_Expected.exposure_columns +
# Coloc_Expected.outcome_columns)))
# Check if all SNPs are in the corr_matrix
if np.any(method_exp in ['cond', 'mask']) or np.any(method_out in ['cond', 'mask']):
if isinstance(corr_matrix, list):
if len(corr_matrix) == 2:
corr_exp = corr_matrix[0]
corr_out = corr_matrix[1]
else:
raise IndexError("`corr_matrix`, please supply a single pd.DataFrame, or a list of two pd.DataFrame's")
else:
corr_exp = corr_matrix.copy()
corr_out = corr_matrix.copy()
print('test1')
# align with dt
_process_corr_mat(data=input_data, corr_matrix=corr_exp,
index_col=Coloc_Expected.VARIANT_ID
)
_process_corr_mat(data=input_data, corr_matrix=corr_out,
index_col=Coloc_Expected.VARIANT_ID
)
else:
# if single just assing whatever is corr_matrix
corr_exp = corr_matrix
corr_out = corr_matrix
# add standard deviation of Y
if not Coloc_Expected.EXPOSURE_STANDARD_DEVIATION in input_data.columns:
input_data[exposure_standard_deviation] =\
sd_pheno_approx(
variancebeta=input_data[exposure_varbeta],
maf=input_data[minor_allele_freq],
n=input_data[exposure_sample_size]
)
if not Coloc_Expected.OUTCOME_STANDARD_DEVIATION in input_data.columns:
input_data[outcome_standard_deviation] =\
sd_pheno_approx(
variancebeta=input_data[outcome_varbeta],
maf=input_data[minor_allele_freq],
n=input_data[outcome_sample_size]
)
# map p-value to z-statistics
if 0 <= pthr <= 1:
zthr = scipy.stats.norm.ppf(1- pthr/2)
else:
raise ValueError("Please supply `pthr` as a value between 0 and 1")
# @@@ DOING ACTUAL ANALYSES
# Finemaping
if verbose is True:
print("[info] Fine-mapping exposure")
# NOTE the .copy() is essential here because it overwrites some columns
# when set to effect_type = 'bin'
fm1 = finemap_signals(
input_data.loc[:,input_data.columns.isin(Exposure_col)].copy(),
corr_matrix=corr_exp, method=method_exp, maxhits=maxhits,
r2thr=r2thr, pthr=pthr, effect_type=effect_type_exp,
variant_id=variant_id,
effect_size=exposure_effect_size,
varbeta=exposure_varbeta,
zstatistic=exposure_zstatistic,
standard_deviation=exposure_standard_deviation,
sample_size=exposure_sample_size,
event_rate=exposure_event_rate,
minor_allele_freq=minor_allele_freq
)
if verbose is True:
print("[info] Variants independently associated in the exposure dataset")
print(fm1)
print('')
print("[info] Fine-mapping outcome")
fm2 = finemap_signals(
input_data.loc[:,input_data.columns.isin(Outcome_col)].copy(),
corr_matrix=corr_out, method=method_out, maxhits=maxhits,
r2thr=r2thr, pthr=pthr, effect_type=effect_type_out,
variant_id=variant_id,
effect_size=outcome_effect_size,
varbeta=outcome_varbeta,
zstatistic=outcome_zstatistic,
standard_deviation=outcome_standard_deviation,
sample_size=outcome_sample_size,
event_rate=outcome_event_rate,
minor_allele_freq=minor_allele_freq
)
if verbose is True:
print("[info] Variants independently associated in the outcome dataset")
print(fm2)
if fm1.empty:
fm1 = find_best_signal(input_data,
variant_id=variant_id,
effect_size=exposure_effect_size,
varbeta=exposure_varbeta,
)
if fm2.empty:
fm2 = find_best_signal(input_data,
variant_id=variant_id,
effect_size=outcome_effect_size,
varbeta=outcome_varbeta,
)
# conditionals if needed
if (not fm1.empty) and (method_exp == "cond"):
# NOTE the copy() is intended, the function overwrites columns
cond1 = est_all_cond(
dt=input_data.loc[:,input_data.columns.isin(Exposure_col)].copy(),
corr_matrix=corr_exp,
fm=fm1, mode=mode, effect_type=effect_type_exp,
variant_id=variant_id,
effect_size=exposure_effect_size,
varbeta=exposure_varbeta,
event_rate=exposure_event_rate,
sample_size=exposure_sample_size,
standard_deviation=exposure_standard_deviation,
minor_allele_freq=minor_allele_freq,
label=Coloc_Expected.exp_name
)
if (not fm2.empty) and (method_out == "cond"):
cond2 = est_all_cond(
dt=input_data.loc[:,input_data.columns.isin(Outcome_col)].copy(),
corr_matrix=corr_out,
fm=fm2, mode=mode, effect_type=effect_type_out,
variant_id=variant_id,
effect_size=outcome_effect_size,
varbeta=outcome_varbeta,
event_rate=outcome_event_rate,
sample_size=outcome_sample_size,
standard_deviation=outcome_standard_deviation,
minor_allele_freq=minor_allele_freq,
label=Coloc_Expected.out_name
)
# doublemask
doublemask_method = ['single', 'mask']
if not np.all([ t in Coloc_Expected.METHODS for t in doublemask_method ]):
raise AttributeError(
"`doublemask_method` should be a list with two entries. Pick from: {}".\
format(Coloc_Expected.METHODS)
)
if (method_exp in doublemask_method) and (method_out in doublemask_method):
# coloc_single == coloc.detail
coloc_res = coloc_single(inst =input_data.copy(),
p1=p1, p2=p2, p12=p12,
effect_type=effect_type,
variant_id=variant_id,
exposure_effect_size=exposure_effect_size,
outcome_effect_size=outcome_effect_size,
exposure_varbeta=exposure_varbeta,
outcome_varbeta=outcome_varbeta,
exposure_zstatistic=exposure_zstatistic,
outcome_zstatistic=outcome_zstatistic,
exposure_standard_deviation=exposure_standard_deviation,
outcome_standard_deviation=outcome_standard_deviation,
exposure_sample_size=exposure_sample_size,
outcome_sample_size=outcome_sample_size,
minor_allele_freq=minor_allele_freq,
extra_info=True)
# get results
proc_res_1 = coloc_process(coloc_res,
hits1=fm1.columns.tolist(),
hits2=fm2.columns.tolist(),
corr_exp = corr_exp,
corr_out = corr_out,
r2thr=r2thr, mode=mode,
p1=p1, p2=p2, p12=p12,
verbose=verbose
)
# get summary
proc_res = proc_res_1[c.COLOC_SUMMARY]
# double cond
if (method_exp in ["cond"]) and (method_out in ["cond"]):
i = [idx for idx, value in enumerate(cond1)]
j = [idx for idx, value in enumerate(cond2)]
todo = pd.DataFrame(itertools.product(i,j), columns=["i", "j"])
# resuts pd.DF
proc_res = pd.DataFrame()
proc_res_variants = pd.DataFrame(
{c.UNIVERSAL_ID.name: input_data[variant_id]})
# looping over rows
for k in list(range(todo.shape[0])):
cond_dt = pd.merge( cond1[list(cond1.keys())[todo.loc[k,"i"]]],
cond2[list(cond2.keys())[todo.loc[k,"j"]]])
### Check if all the columns are there
column_check = [exposure_sample_size, outcome_sample_size,
exposure_standard_deviation,
outcome_standard_deviation,
minor_allele_freq]
cond_dt = _add_missing_columns(cond_dt, input_data.copy(),
columns=column_check,
index_col=variant_id)
# Getting coloc results
coloc_res = coloc_single(inst=cond_dt,
p1=p1, p2=p2, p12=p12,
effect_type=effect_type,
variant_id=variant_id,
exposure_effect_size=exposure_effect_size,
outcome_effect_size=outcome_effect_size,
exposure_varbeta=exposure_varbeta,
outcome_varbeta=outcome_varbeta,
exposure_zstatistic=exposure_zstatistic,
outcome_zstatistic=outcome_zstatistic,
exposure_standard_deviation=exposure_standard_deviation,
outcome_standard_deviation=outcome_standard_deviation,
exposure_sample_size=exposure_sample_size,
outcome_sample_size=outcome_sample_size,
minor_allele_freq=minor_allele_freq,
extra_info=True)
# Processing
proc_res_1 = coloc_process(coloc_res,
hits1=[fm1.columns.tolist()[todo.loc[k,"i"]]],
hits2=[fm2.columns.tolist()[todo.loc[k,"j"]]],
corr_exp = corr_exp,
corr_out = corr_out,
r2thr=r2thr, mode=mode, p1=p1, p2=p2,
p12=p12,
verbose=verbose
)
# Getting Summary reults
proc_res= pd.concat([proc_res, proc_res_1[c.COLOC_SUMMARY]])
# merging — use iteration-specific suffixes to avoid
# duplicate column names across signal-pair iterations
proc_res_variants = pd.merge(
proc_res_variants, coloc_res.variant_dataframe(),
left_on=c.UNIVERSAL_ID.name, right_on=c.UNIVERSAL_ID.name,
suffixes=("", "_{0}".format(k))
)
# cond mask/-
if method_exp in ["cond"] and method_out in ["mask", "single"]:
proc_res = pd.DataFrame()
# Looping over the keys
for k in list(range(len(cond1))):
# this merges cond1 = exposure with outcome data
cond_dt = pd.merge(
cond1[list(cond1.keys())[k]],
input_data.loc[:,input_data.columns.isin(Outcome_col)].copy(),
)
### adding back potentially missing data
column_check = [exposure_sample_size, outcome_sample_size,
exposure_standard_deviation,
outcome_standard_deviation,
minor_allele_freq]
cond_dt = _add_missing_columns(cond_dt, input_data,
columns=column_check,
index_col=variant_id)
# run coloc
coloc_res = coloc_single(inst=cond_dt,
p1=p1, p2=p2, p12=p12,
effect_type=effect_type,
variant_id=variant_id,
exposure_effect_size=exposure_effect_size,
outcome_effect_size=outcome_effect_size,
exposure_varbeta=exposure_varbeta,
outcome_varbeta=outcome_varbeta,
exposure_zstatistic=exposure_zstatistic,
outcome_zstatistic=outcome_zstatistic,
exposure_standard_deviation=exposure_standard_deviation,
outcome_standard_deviation=outcome_standard_deviation,
exposure_sample_size=exposure_sample_size,
outcome_sample_size=outcome_sample_size,
minor_allele_freq=minor_allele_freq,
extra_info=True)
# get summary
proc_res_1 = coloc_process(coloc_res,
hits1=[fm1.columns.tolist()[k]],
hits2=fm2.columns.tolist(),
corr_exp = corr_exp,
corr_out = corr_out,
r2thr=r2thr,
mode=mode, p1=p1, p2=p2, p12=p12,
verbose=verbose
)
proc_res = pd.concat([proc_res, proc_res_1[c.COLOC_SUMMARY]])
## mask/- cond
if method_exp in ["mask", "single"] and method_out in ["cond"]:
proc_res = pd.DataFrame()
# NOTE change with `enumerate`?
for k in list(range(0,len(cond2))):
# this merges cond2 = outcomes with exposure data
cond_dt = pd.merge(
cond2[list(cond2.keys())[k]],
input_data.loc[:,input_data.columns.isin(Exposure_col)].copy()
)
### adding back potentially missing data
column_check = [exposure_sample_size, outcome_sample_size,
exposure_standard_deviation,
outcome_standard_deviation,
minor_allele_freq]
cond_dt = _add_missing_columns(cond_dt, input_data,
columns=column_check,
index_col=variant_id)
# COLOC results
coloc_res = coloc_single(inst=cond_dt,
p1=p1, p2=p2, p12=p12,
effect_type=effect_type,
variant_id=variant_id,
exposure_effect_size=exposure_effect_size,
outcome_effect_size=outcome_effect_size,
exposure_varbeta=exposure_varbeta,
outcome_varbeta=outcome_varbeta,
exposure_zstatistic=exposure_zstatistic,
outcome_zstatistic=outcome_zstatistic,
exposure_standard_deviation=exposure_standard_deviation,
outcome_standard_deviation=outcome_standard_deviation,
exposure_sample_size=exposure_sample_size,
outcome_sample_size=outcome_sample_size,
minor_allele_freq=minor_allele_freq,
extra_info=True)
# get summary
proc_res_1 = coloc_process(coloc_res,
hits1=fm1.columns.tolist(),
hits2=[fm2.columns.tolist()[k]],
corr_exp = corr_exp,
corr_out = corr_out,
r2thr=r2thr,
mode=mode, p12=p12, p1=p1, p2=p2,
verbose=verbose
)
proc_res = pd.concat([proc_res,
proc_res_1[c.COLOC_SUMMARY]])
# returning results objects
proc_res.reset_index(inplace=True)
proc_res[c.MR_NSNPS] = proc_res[c.MR_NSNPS].astype(str)
results = proc_res[
[c.COLOC_HIT1, c.COLOC_HIT2,
c.MR_NSNPS, c.COLOC_PPH0, c.COLOC_PPH1,
c.COLOC_PPH2, c.COLOC_PPH3, c.COLOC_PPH4,
c.COLOC_BEST1, c.COLOC_BEST2, c.COLOC_BEST4]
]
# add zstats
try:
results[c.COLOC_ZSTAT_HIT1] = fm1[results[c.COLOC_HIT1]].apply(pd.Series)
results[c.COLOC_ZSTAT_HIT2] = fm2[results[c.COLOC_HIT2]].apply(pd.Series)
except:
results[c.COLOC_ZSTAT_HIT1] = fm1[results[c.COLOC_HIT1]].iloc[0].to_list()
results[c.COLOC_ZSTAT_HIT2] = fm2[results[c.COLOC_HIT2]].iloc[0].to_list()
# returning ColocResults objects
res_obj = {c.COLOC_SUMMARY: results,
c.COLOC_VARIANT_RESULTS: proc_res_1[c.COLOC_RESULTS]}
# return stuff
# return results
return ColocResults(**res_obj)

coloc_core.py at commit 0e05762, under GPL-3.0 · at the source

Overview

Authors: Abigail R. ter Kuile¹, Chris Finan^2,3, Sandesh Chopade², Marion van Vugt^2,4,5, Nikita Hukerikar^2,6, Serena Barral⁷, Argyris Stringaris^1,8, Amand F. Schmidt^2,3,5, Karoline Kuchenbaecker⁸, Jean-Baptiste Pingault^1,9

ORCID iDs: Abigail R. ter Kuile, Chris Finan, Sandesh Chopade, Marion van Vugt, Nikita Hukerikar, Serena Barral, Argyris Stringaris, Amand F. Schmidt, Karoline Kuchenbaecker, Jean-Baptiste Pingault

Department of Clinical, Educational and Health Psychology, University College London,London, UK
Institute of Cardiovascular Science, University College London,London, UK
The National Institute for Health Research University College London Hospitals Biomedical Research Centre, University College London,London, UK
Department of Cardiology, Division of Heart and Lungs, University Medical Centre Utrecht, Utrecht University,Utrecht, the Netherlands
Department of Cardiology, Amsterdam Cardiovascular Sciences, Amsterdam University Medical Centre, University of Amsterdam,Amsterdam, the Netherlands
Institute of Health Informatics, University College London,London, UK
Developmental Neurosciences, Zayed Centre for Research into Rare Disease in Children, UCL GOS Institute of Child Health, University College London,London, UK
Division of Psychiatry, University College London,London, UK
Social Genetic and Developmental Psychiatry, King’s College London,London, UK

Institutions: University College London (United Kingdom); Utrecht University (Netherlands); University of Amsterdam (Netherlands); King's College London (United Kingdom)

Journal: Translational psychiatry, volume 16, issue 1, article 392

Dates: received 26 September 2025; accepted 22 May 2026; published online 2 June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41398-026-04137-9 · PMID 42230543 · PMCID PMC13443190 · OpenAlex W7163184988

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), depression (population), clinical / translational (subfield)

Methods: Statistics

Keywords: Genomics, Depression, Clinical pharmacology

MeSH: Antidepressive Agents*, Drug Repositioning*, Major Depressive Disorder*, Mendelian Randomization Analysis*, Genome-Wide Association Study, Humans, Quantitative Trait Loci (* major topic)

Topic: Genetic Associations and Epidemiology (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: EC | EU Framework Programme for Research and Innovation H2020 | H2020 Priority Excellent Science | H2020 European Research Council (H2020 Excellent Science - European Research Council) (I-IRISK grant agreement No. 863981); Great Ormond Street Hospital Children’s Charity; British Heart Foundation (PG/18/5033837, PG/22/10989); National Institute for Health Research University College London Hospitals Biomedical Research Centre

Citations: not cited yet (Europe PMC); 101 references in the paper

Abstract

Major depression (MD) treatments have limited efficacy and target few mechanisms, highlighting the need for innovative drug discovery. Drugs targeting genetically supported proteins are 2.6 times more likely to succeed in drug development. Here, we use genetic methods to identify and prioritise MD drug targets, leveraging genome-wide association study (GWAS) summary statistics from >525,000 MD cases. We derived exposure data from 10 datasets measuring protein quantitative trait loci (pQTLs) and gene expression levels (eQTLs) in blood, cerebrospinal fluid, and brain tissues. We performed cis-Mendelian randomisation (MR) on 3469 druggable targets (genes encoding proteins targeted by existing compounds or experimentally predicted to be druggable). To strengthen causal inference, we implemented robust MR estimators, colocalisation, external replication, and assessed directional consistency across tissues. We integrated cis-MR effect directions with drug mechanisms and clinical annotations to infer potential therapeutic effects. Validation analyses showed that 82% of drugs approved for depression/anxiety had ≥1 significant MR target, compared to 51% for compounds in clinical trials. For repurposing, we prioritised 54 targets of compounds developed for other conditions with estimated beneficial effects on MD (e.g., an inhibitor for a risk-increasing target). Ten high-priority targets of brain-penetrating compounds included ACE and NISCH (cardiovascular drugs), NDUFA2, NDUFB6, and NDUFS1 (metformin), CDK4, NTRK3, and MET (oncology inhibitors), and GLS and NOS2 (enzyme inhibitors). We found genetic evidence for established and novel MD targets across the drug development pipeline. Novel targets point to mechanisms beyond monoaminergic systems, most with approved drugs for other conditions, offering immediate repurposing opportunities.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.

cfinan/merit

License: GPL-3.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 0e05762f86400756a0152cb5aa645edb81d8d831, 4 September 2026

Languages: Python (398), R (28), Jupyter (26), Shell (22)

Size: 1,787 files, 474 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, environment (pyproject.toml, requirements.txt, resources/conda/conda_build_config.yaml), tests, continuous integration, documentation, 27 notebooks

Not found: CITATION.cff

Tools: NumPy (194 files), pandas (164 files), SciPy (28 files), pysam (24 files), statsmodels (9 files), Matplotlib (4 files), seaborn (4 files), Numba (3 files), rpy2 (3 files), BCFtools (1 file), data.table (1 file), psych (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

476 files

docs/source/conf.py, Python, 98 lines
merit/__init__.py, Python, 16 lines
merit/coloc/__init__.py, Python, 17 lines
merit/coloc/coloc_core.py, Python, 2,566 lines, 2 matches
merit/coloc/coloc_sensitivity.py, Python, 359 lines
merit/common.py, Python, 24 lines
merit/constants.py, Python, 794 lines
merit/errors.py, Python, 1,702 lines
merit/example_data/__init__.py, Python, 1 line
merit/example_data/downloads.py, Python, 79 lines
merit/example_data/dummy_ipd_data.py, Python, 312 lines
merit/example_data/dummy_summary_data.py, Python, 1,126 lines
merit/example_data/example_datasets/test_file_pools/cohort_c/chr_hint/create_bgen.sh, Shell, 1 line
merit/example_data/examples.py, Python, 1,213 lines
merit/gwas/__init__.py, Python, 9 lines
merit/gwas/dataset.py, Python, 2,897 lines
merit/gwas/utils.py, Python, 3,684 lines
merit/ipd/__init__.py, Python, 80 lines
merit/ipd/admission.py, Python, 753 lines
merit/ipd/common.py, Python, 1,766 lines
merit/ipd/corr_factory.py, Python, 2,127 lines
merit/ipd/corr_matrix.py, Python, 4,026 lines
merit/ipd/dosage_cache.py, Python, 1,313 lines
merit/ipd/files.py, Python, 2,856 lines
merit/ipd/manifest.py, Python, 346 lines
merit/ipd/regions.py, Python, 3,606 lines
merit/ld/__init__.py, Python, 31 lines
merit/ld/clumping.py, Python, 940 lines
merit/ld/factory_clumping.py, Python, 1,022 lines
merit/ld/handlers.py, Python, 781 lines
merit/ld/variants.py, Python, 81 lines
merit/misc/__init__.py, Python, 9 lines
merit/misc/sojo.py, Python, 550 lines
merit/misc/steiger_filtering.py, Python, 1,247 lines
merit/mr/__init__.py, Python, 47 lines
merit/mr/base.py, Python, 282 lines
merit/mr/conmix.py, Python, 229 lines
merit/mr/egger.py, Python, 226 lines
merit/mr/ivw.py, Python, 207 lines, 1 match
merit/mr/median.py, Python, 321 lines
merit/mr/mendelian_randomization.py, Python, 47 lines
merit/mr/mode.py, Python, 549 lines
merit/mr/mvmr.py, Python, 432 lines
merit/mr/selection.py, Python, 127 lines
merit/mr/utils.py, Python, 71 lines
merit/mr/wald.py, Python, 118 lines
merit/pipelines/__init__.py, Python, 306 lines
merit/pipelines/_outcome_loader.py, Python, 482 lines
merit/pipelines/assembly.py, Python, 1,066 lines
merit/pipelines/cli/__init__.py, Python, 10 lines
merit/pipelines/cli/_args.py, Python, 3,107 lines
merit/pipelines/cli/_log.py, Python, 353 lines
merit/pipelines/cli/run_coloc.py, Python, 331 lines
merit/pipelines/cli/run_extract.py, Python, 562 lines
merit/pipelines/cli/run_mr.py, Python, 741 lines
merit/pipelines/cli/run_mr_hpc.py, Python, 754 lines
merit/pipelines/cli/run_mr_parallel.py, Python, 779 lines
merit/pipelines/coloc.py, Python, 1,383 lines
merit/pipelines/columns.py, Python, 61 lines
merit/pipelines/constants.py, Python, 370 lines
merit/pipelines/data/__init__.py, Python, 1 line
merit/pipelines/extract.py, Python, 3,557 lines
merit/pipelines/features.py, Python, 1,043 lines
merit/pipelines/filtering.py, Python, 1,690 lines
merit/pipelines/formatting.py, Python, 320 lines
merit/pipelines/hpc/__init__.py, Python, 27 lines
merit/pipelines/hpc/jaf.py, Python, 401 lines
merit/pipelines/hpc/planner.py, Python, 1,171 lines
merit/pipelines/hpc/task_script.py, Python, 512 lines
merit/pipelines/io.py, Python, 766 lines
merit/pipelines/ld_reference.py, Python, 587 lines
merit/pipelines/maf.py, Python, 1,586 lines
merit/pipelines/methods/__init__.py, Python, 25 lines
merit/pipelines/methods/_config.py, Python, 267 lines
merit/pipelines/methods/_protocol.py, Python, 190 lines
merit/pipelines/methods/_registry.py, Python, 272 lines
merit/pipelines/methods/_support.py, Python, 170 lines
merit/pipelines/methods/egger.py, Python, 198 lines
merit/pipelines/methods/ivw.py, Python, 180 lines
merit/pipelines/methods/median.py, Python, 189 lines
merit/pipelines/methods/mode.py, Python, 223 lines
merit/pipelines/methods/mvmr.py, Python, 148 lines
merit/pipelines/methods/wald.py, Python, 111 lines
merit/pipelines/mr.py, Python, 4,608 lines
merit/pipelines/parallel/__init__.py, Python, 40 lines
merit/pipelines/parallel/executor.py, Python, 581 lines
merit/pipelines/parallel/options.py, Python, 164 lines
merit/pipelines/parallel/publish.py, Python, 58 lines
merit/pipelines/parallel/runtime.py, Python, 566 lines
merit/pipelines/parallel/worker.py, Python, 331 lines
merit/pipelines/parameters.py, Python, 613 lines
merit/pipelines/region_plan.py, Python, 167 lines
merit/pipelines/result.py, Python, 57 lines
merit/pipelines/sample_size.py, Python, 81 lines
merit/pipelines/validation.py, Python, 283 lines
merit/stats.py, Python, 49 lines
resources/analysis/fn_12_missingness_failure_prediction/__init__.py, Python, 14 lines
resources/analysis/fn_12_missingness_failure_prediction/_paths.py, Python, 52 lines
resources/analysis/fn_12_missingness_failure_prediction/definedness.py, Python, 154 lines
resources/analysis/fn_12_missingness_failure_prediction/features.py, Python, 243 lines
resources/analysis/fn_12_missingness_failure_prediction/index_build.py, Python, 369 lines
resources/analysis/fn_12_missingness_failure_prediction/index_parse.py, Python, 102 lines
resources/analysis/fn_12_missingness_failure_prediction/io_util.py, Python, 126 lines
resources/analysis/fn_12_missingness_failure_prediction/metrics.py, Python, 71 lines
resources/analysis/fn_12_missingness_failure_prediction/model.py, Python, 388 lines
resources/analysis/fn_12_missingness_failure_prediction/pipeline.py, Python, 433 lines
resources/analysis/fn_12_missingness_failure_prediction/plot_roc.py, Python, 294 lines
resources/analysis/fn_12_missingness_failure_prediction/predictors.py, Python, 201 lines
resources/analysis/fn_12_missingness_failure_prediction/project.py, Python, 212 lines
resources/analysis/fn_12_missingness_failure_prediction/report.py, Python, 532 lines
resources/analysis/fn_12_missingness_failure_prediction/roc.py, Python, 171 lines
resources/analysis/fn_12_missingness_failure_prediction/run.py, Python, 714 lines
resources/analysis/fn_12_missingness_failure_prediction/sampling.py, Python, 55 lines
resources/analysis/fn_12_missingness_failure_prediction/simulate_production.py, Python, 639 lines
resources/analysis/fn_12_missingness_failure_prediction/subsets.py, Python, 617 lines
resources/analysis/fn_12_missingness_failure_prediction/windows_io.py, Python, 126 lines
resources/analysis/fn_51_missingness_vs_d/__init__.py, Python, 14 lines
resources/analysis/fn_51_missingness_vs_d/_paths.py, Python, 53 lines
resources/analysis/fn_51_missingness_vs_d/calibrate.py, Python, 339 lines
resources/analysis/fn_51_missingness_vs_d/io_util.py, Python, 137 lines
resources/analysis/fn_51_missingness_vs_d/metrics.py, Python, 233 lines
resources/analysis/fn_51_missingness_vs_d/pairwise.py, Python, 121 lines
resources/analysis/fn_51_missingness_vs_d/production_parity.py, Python, 149 lines
resources/analysis/fn_51_missingness_vs_d/sampling.py, Python, 55 lines
resources/analysis/fn_51_missingness_vs_d/select_windows.py, Python, 292 lines
resources/analysis/fn_51_missingness_vs_d/summarize.py, Python, 367 lines
resources/analysis/fn_51_missingness_vs_d/thresholds.py, Python, 108 lines
resources/analysis/fn_51_missingness_vs_d/vcf_io.py, Python, 389 lines
resources/analysis/fn_51_missingness_vs_d/verify_pairwise.py, Python, 383 lines
resources/analysis/fn_51_missingness_vs_d/windows.py, Python, 453 lines
resources/bin/demo_corr_factory.py, Python, 652 lines
resources/bin/demo_dosage_cache.py, Python, 462 lines
resources/bin/download_1kgp.sh, Shell, 470 lines
resources/bin/gen_corr_golden.py, Python, 114 lines
resources/bin/probe_query_edge_truncation.py, Python, 215 lines
resources/ci_cd/before_script.sh, Shell, 30 lines
resources/ci_cd/build.sh, Shell, 11 lines
resources/ci_cd/pages.sh, Shell, 27 lines
resources/ci_cd/run_docker.sh, Shell, 68 lines
resources/ci_cd/run_pages.sh, Shell, 32 lines
resources/conda/conda-build.sh, Shell, 4 lines
resources/examples/1kgp_01_setup.ipynb, Jupyter, 178 lines
resources/examples/1kgp_02_variant_queries.ipynb, Jupyter, 220 lines
resources/examples/1kgp_03_dosage_caches.ipynb, Jupyter, 185 lines
resources/examples/1kgp_04_ld_analysis.ipynb, Jupyter, 351 lines
resources/examples/coloc.ipynb, Jupyter, 238 lines
resources/examples/drug_target_mr.ipynb, Jupyter, 316 lines
resources/examples/ld_clumping.ipynb, Jupyter, 346 lines
resources/examples/mendelian_randomization.ipynb, Jupyter, 237 lines
resources/misc/2026-sort/merit/ipd/cohort.py, Python, 159 lines
resources/misc/2026-sort/merit/ipd/common.py, Python, 1,889 lines
resources/misc/2026-sort/merit/ipd/engine.py, Python, 184 lines
resources/misc/2026-sort/merit/ipd/file_container.py, Python, 386 lines
resources/misc/2026-sort/merit/ipd/file_pool.py, Python, 276 lines
resources/misc/2026-sort/merit/ipd/readers.py, Python, 1,613 lines
resources/misc/2026-sort/merit/ipd/workers.py, Python, 70 lines
resources/misc/code_other_packages/MendelianRandomization/mr_conmix-methods.R, R, 181 lines
resources/misc/code_other_packages/MendelianRandomization/mr_mode.R, R, 57 lines
resources/misc/code_other_packages/TwoSampleMR/mr/mr_wald_ratio.R, R, 38 lines
resources/misc/code_other_packages/coloc/coloc_baseline.R, R, 391 lines
resources/misc/code_other_packages/coloc/coloc_sensitivity.R, R, 114 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/phenocov_mvmr.R, R, 66 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/pl2_mvmr_one.R, R, 52 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/pl2_mvmr_two.R, R, 39 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/pl_mvmr.R, R, 73 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/qhet_mvmr.R, R, 218 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/qmin.R, R, 112 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/snpcov_mvmr.R, R, 91 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/strength_mvmr.R, R, 213 lines
resources/misc/code_other_packages/mvmr_wes_spiller/code/strhet_mvmr.R, R, 220 lines
resources/misc/code_other_packages/mvmr_wes_spiller/docs/MVMR.rmd, R, 174 lines
resources/misc/code_other_packages/sojo/load_data_and_arguments.py, Python, 68 lines
resources/misc/code_other_packages/sojo/sojo_newest.R, R, 311 lines
resources/misc/code_other_packages/sojo/sojo_pytest.r, R, 288 lines
resources/misc/dead_code/analysis_tools_dead_code.py, Python, 180 lines
resources/misc/dead_code/cohorts.py, Python, 2,042 lines
resources/misc/dead_code/config.py, Python, 884 lines
resources/misc/dead_code/config2.py, Python, 154 lines
resources/misc/dead_code/dataset_dead_code.py, Python, 774 lines
resources/misc/dead_code/emerald.py, Python, 208 lines
resources/misc/dead_code/ensembl/__init__.py, Python, 1 line
resources/misc/dead_code/ensembl/constants.py, Python, 500 lines
resources/misc/dead_code/ensembl/coords.py, Python, 167 lines
resources/misc/dead_code/ensembl/ensembl_errors.py, Python, 129 lines
resources/misc/dead_code/ensembl/rest.py, Python, 1,183 lines
resources/misc/dead_code/ensembl/variants.py, Python, 2,577 lines
resources/misc/dead_code/examples/allele_norm.ipynb, Jupyter, 101 lines
resources/misc/dead_code/examples/chrpos_parsing.ipynb, Jupyter, 150 lines
resources/misc/dead_code/examples/cohort_init_from_genomic_config.ipynb, Jupyter, 60 lines
resources/misc/dead_code/examples/drug_target_mr_example.ipynb, Jupyter, 374 lines
resources/misc/dead_code/examples/liftover.ipynb, Jupyter, 125 lines
resources/misc/dead_code/examples/run_simple_prs.ipynb, Jupyter, 454 lines
resources/misc/dead_code/examples/run_simple_prs_example.ipynb, Jupyter, 303 lines
resources/misc/dead_code/examples/sojo.ipynb, Jupyter, 278 lines
resources/misc/dead_code/interface.py, Python, 317 lines
resources/misc/dead_code/ipd/__init__.py, Python, 1 line
resources/misc/dead_code/ipd/agnostic_handler.py, Python, 1 line
resources/misc/dead_code/ipd/bed_handler.py, Python, 1 line
resources/misc/dead_code/ipd/bgen_handler.py, Python, 2,624 lines
resources/misc/dead_code/ipd/ipd_tools.py, Python, 982 lines
resources/misc/dead_code/ipd/vcf_handler.py, Python, 1 line
resources/misc/dead_code/mrbase/__init__.py, Python, 1 line
resources/misc/dead_code/mrbase/mrbase_api.py, Python, 1,377 lines
resources/misc/dead_code/pools/__init__.py, Python, 1 line
resources/misc/dead_code/pools/bgen_pool.py, Python, 43 lines
resources/misc/dead_code/pools/cohort_pool.py, Python, 1 line
resources/misc/dead_code/pools/file_pool.py, Python, 415 lines
resources/misc/dead_code/pools/vcf_pool.py, Python, 88 lines
resources/misc/dead_code/prs/__init__.py, Python, 1 line
resources/misc/dead_code/prs/scores.py, Python, 1,257 lines
resources/misc/dead_code/prs/simple_prs.py, Python, 1,263 lines
resources/misc/dead_code/prs/validation.py, Python, 465 lines
resources/misc/dead_code/prs/validation_calibration.py, Python, 509 lines
resources/misc/dead_code/qctool_dead_code.py, Python, 439 lines
resources/misc/dead_code/refpops.py, Python, 1,531 lines
resources/misc/dead_code/utils/liftover.py, Python, 787 lines
resources/misc/dead_code/utils/plink.py, Python, 140 lines
resources/misc/dead_code/utils/vcf.py, Python, 72 lines
resources/misc/dead_code/visual/visualisation.py, Python, 178 lines
resources/misc/dead_code/wrappers/__init__.py, Python, 1 line
resources/misc/dead_code/wrappers/bgenix.py, Python, 534 lines
resources/misc/dead_code/wrappers/plink.py, Python, 523 lines
resources/misc/dead_code/wrappers/qctool.py, Python, 745 lines
resources/misc/dead_code/wrappers/wrapper_errors.py, Python, 30 lines
resources/misc/examples/COLOC_example.ipynb, Jupyter, 382 lines
resources/misc/examples/LD_tests.ipynb, Jupyter, 347 lines
resources/misc/examples/MR_example.ipynb, Jupyter, 227 lines
resources/misc/examples/allele_conversion_testing.ipynb, Jupyter, 202 lines
resources/misc/examples/coloc.ipynb, Jupyter, 153 lines
resources/misc/examples/example_path.py, Python, 7 lines
resources/misc/examples/mrbase_query.ipynb, Jupyter, 121 lines
resources/misc/examples/mrbase_query.py, Python, 113 lines
resources/misc/helper_scripts/LD_clumping_tests.ipynb, Jupyter, 211 lines
resources/misc/helper_scripts/build_druggable_list.sh, Shell, 71 lines
resources/misc/helper_scripts/dup_test.ipynb, Jupyter, 18 lines
resources/misc/helper_scripts/emerald.ipynb, Jupyter, 187 lines
resources/misc/helper_scripts/ensembl.py, Python, 817 lines
resources/misc/helper_scripts/ensembl_regulatory_activity.sh, Shell, 52 lines
resources/misc/helper_scripts/ensembl_regulatory_peaks.sh, Shell, 56 lines
resources/misc/helper_scripts/extract_eqtl.sh, Shell, 108 lines
resources/misc/helper_scripts/extract_somalogic.sh, Shell, 83 lines
resources/misc/helper_scripts/fetch_gene_coords.sh, Shell, 37 lines
resources/misc/helper_scripts/fetch_on_asid.sh, Shell, 35 lines
resources/misc/helper_scripts/fetch_targets.sh, Shell, 22 lines
resources/misc/helper_scripts/generate_plink_testset.sh, Shell, 69 lines
resources/misc/helper_scripts/get_chembl_indication.sh, Shell, 57 lines
resources/misc/helper_scripts/make_physical_regions.sh, Shell, 10 lines
resources/misc/helper_scripts/run_simple_prs.sh, Shell, 42 lines
resources/misc/helper_scripts/test_bgenix.py, Python, 18 lines
resources/misc/helper_scripts/test_corr_plink.py, Python, 66 lines
resources/misc/helper_scripts/test_emerald.py, Python, 28 lines
resources/misc/helper_scripts/test_filter.py, Python, 99 lines
resources/misc/helper_scripts/test_ld_matrix.py, Python, 34 lines
resources/misc/helper_scripts/test_park_amd.py, Python, 66 lines
resources/misc/helper_scripts/test_plink.py, Python, 49 lines
resources/misc/helper_scripts/test_simple_prs.py, Python, 96 lines
resources/misc/helper_scripts/test_vcf_write.py, Python, 18 lines
resources/misc/issue_debug/51_issue_debug.py, Python, 49 lines
resources/misc/issue_debug/53_issue_debug.py, Python, 90 lines
resources/misc/issue_debug/57_issue_debug.py, Python, 62 lines
resources/misc/issue_debug/63_issue_debug.ipynb, Jupyter, 120 lines
resources/misc/issue_debug/63_issue_debug.py, Python, 45 lines
resources/misc/meta_analysis/__init__.py, Python, 1 line
resources/misc/meta_analysis/meta_analysis.py, Python, 159 lines
resources/misc/mv_mr_examples/multivariatemr_test.R, R, 35 lines
resources/misc/old_tests/test_bgen_handler.py, Python, 1,131 lines
resources/misc/old_tests/test_ipd_tools.py, Python, 234 lines
resources/misc/old_tests/tests/example_data/test_examples.py, Python, 129 lines
resources/misc/old_tests/tests/mrbase/test_mrbase_api.py, Python, 720 lines
resources/misc/old_tests/tests/prs/conftest.py, Python, 580 lines
resources/misc/old_tests/tests/prs/test_prs.py, Python, 303 lines
resources/misc/old_tests/tests/prs/test_scores.py, Python, 429 lines
resources/misc/old_tests/tests/prs/test_simple_prs.py, Python, 335 lines
resources/misc/old_tests/tests/prs/test_validation.py, Python, 118 lines
resources/misc/old_tests/tests/prs/test_validation_calibration.py, Python, 229 lines
resources/misc/old_tests/tests/test_errors.py, Python, 554 lines
resources/misc/old_tests/tests/test_refpops.py, Python, 1,051 lines
resources/misc/old_tests/tests/utils/test_pyplink.py, Python, 96 lines
resources/misc/old_tests/tests/utils/test_steiger_filtering.py, Python, 375 lines
resources/misc/old_tests/tests/utils/test_vcf.py, Python, 19 lines
resources/misc/old_tests/tests/wrappers/test_bgenix.py, Python, 101 lines
resources/misc/old_tests/tests/wrappers/test_plink.py, Python, 150 lines
resources/misc/old_tests/tests/wrappers/test_qctool.py, Python, 485 lines
resources/misc/tests_to_check/genotypes/conftest.py, Python, 32 lines
resources/misc/tests_to_check/genotypes/test_bgen_pool.py, Python, 32 lines
resources/misc/tests_to_check/genotypes/test_file_pool.py, Python, 151 lines
resources/misc/tests_to_check/genotypes/test_vcf_pool.py, Python, 22 lines
resources/misc/tests_to_check/test_checks.py, Python, 1,032 lines
resources/misc/tests_to_check/test_config.py, Python, 180 lines
resources/misc/tests_to_check/test_data/test_mr_data.py, Python, 99 lines
resources/misc/tests_to_check/test_dummy_data.py, Python, 298 lines
resources/misc/tests_to_check/test_ld_calcs.py, Python, 161 lines
resources/misc/tests_to_check/test_ld_clumping.py, Python, 111 lines
resources/misc/tests_to_check/test_mrbase.py, Python, 55 lines
resources/misc/tests_to_check/test_plink_utils.py, Python, 169 lines
resources/misc/tests_to_check/test_read_gwas.py, Python, 1,535 lines
resources/misc/tests_to_check/test_refpops.py, Python, 997 lines
resources/misc/tests_to_check/test_vcf_utils.py, Python, 175 lines
resources/misc/tests_to_check/utils/test_liftover.py, Python, 124 lines
resources/misc/work_in_progress/bma_mr/summary_mvMR_BF.R, R, 451 lines
resources/misc/work_in_progress/bma_mr/summary_mvMR_SSS.R, R, 503 lines
resources/misc/work_in_progress/mr_sim/mr_twosample_sim.py, Python, 134 lines
resources/misc/work_in_progress/mvmr_lasso/mvMR_glmnet.R, R, 1,586 lines
resources/misc/work_in_progress/pleio_lasso/mr_lasso-methods.R, R, 97 lines
resources/misc/work_in_progress/pleio_lasso/mr_mvlasso-methods.R, R, 106 lines
resources/misc/work_in_progress/selection/selection_pointestimate.py, Python, 111 lines
resources/misc/work_in_progress/selection/selection_precision.py, Python, 132 lines
resources/qa-tests/fn-51-missingness-vs-d/run.sh, Shell, 148 lines
resources/qa-tests/fn-57-external-ld-corpus/run.sh, Shell, 297 lines
resources/test_data/ld_clumping/build_clumping_corpus.py, Python, 5,003 lines
resources/test_data/ld_clumping/fixtures/naive_clump_r_input.R, R, 78 lines
resources/test_data/ld_clumping/naive_clump.R, R, 517 lines
resources/test_data/ld_clumping/naive_clump.py, Python, 440 lines
resources/test_data/ld_clumping/preflight.py, Python, 1,850 lines
resources/test_data/ld_clumping/roundtrip_numeric.R, R, 29 lines
resources/test_data/ld_clumping/run_f1_parity.py, Python, 149 lines
resources/test_data/ld_clumping/test_agreement_gate.py, Python, 533 lines
resources/test_data/ld_clumping/test_artifact_emission.py, Python, 178 lines
resources/test_data/ld_clumping/test_auxiliary_fixtures.py, Python, 85 lines
resources/test_data/ld_clumping/test_build_clumping_corpus.py, Python, 146 lines
resources/test_data/ld_clumping/test_naive_clump.R, R, 264 lines
resources/test_data/ld_clumping/test_naive_clump.py, Python, 143 lines
resources/test_data/ld_clumping/test_parameter_sweeps.py, Python, 514 lines
resources/test_data/ld_clumping/test_sumstats_generation.py, Python, 170 lines
resources/test_data/ld_clumping/validate_clumping_corpus.py, Python, 1,448 lines
resources/test_data/ld_reference/build_reference_corpus.py, Python, 3,927 lines
resources/test_data/ld_reference/reference_corr.R, R, 307 lines
resources/test_data/ld_reference/validate_reference_corpus.py, Python, 1,319 lines
tests/benchmarks/r_code/test_mendelian_randomization.py, Python, 570 lines
tests/benchmarks/r_code/test_merit_vs_r.py, Python, 322 lines
tests/benchmarks/r_code/test_steiger_filtering.py, Python, 275 lines
tests/coloc/test_coloc_core.py, Python, 1,461 lines
tests/coloc/test_coloc_sensitivity.py, Python, 102 lines
tests/conftest.py, Python, 60 lines
tests/example_data/test_dummy_ipd_data.py, Python, 127 lines
tests/example_data/test_dummy_summary_data.py, Python, 553 lines
tests/example_data/test_plink_write.py, Python, 115 lines
tests/gwas/test_dataset.py, Python, 3,377 lines
tests/gwas/test_utils.py, Python, 2,534 lines
tests/ipd/admission_fixtures.py, Python, 188 lines
tests/ipd/conftest.py, Python, 788 lines
tests/ipd/corr_test_oracle.py, Python, 671 lines
tests/ipd/helpers.py, Python, 216 lines
tests/ipd/property_oracle.py, Python, 340 lines
tests/ipd/test_admission.py, Python, 1,093 lines
tests/ipd/test_benchmark_corrcoef.py, Python, 205 lines
tests/ipd/test_cohort_init.py, Python, 259 lines
tests/ipd/test_cohort_reader.py, Python, 124 lines
tests/ipd/test_cohort_sample_subsetting.py, Python, 511 lines
tests/ipd/test_common.py, Python, 475 lines
tests/ipd/test_common_sidecar.py, Python, 668 lines
tests/ipd/test_corr_factory.py, Python, 3,036 lines
tests/ipd/test_corr_matrix.py, Python, 5,486 lines
tests/ipd/test_corr_matrix_cleanup_primitives.py, Python, 1,301 lines
tests/ipd/test_corr_matrix_default_goldens.py, Python, 149 lines
tests/ipd/test_corr_matrix_e2e.py, Python, 255 lines
tests/ipd/test_corr_test_oracle.py, Python, 478 lines
tests/ipd/test_corrcoef.py, Python, 674 lines
tests/ipd/test_dosage_cache.py, Python, 1,223 lines
tests/ipd/test_external_ld_reference_e2e.py, Python, 932 lines
tests/ipd/test_external_ld_reference_generation.py, Python, 1,456 lines
tests/ipd/test_files.py, Python, 2,525 lines
tests/ipd/test_manifest.py, Python, 676 lines
tests/ipd/test_metamorphic_admission.py, Python, 852 lines
tests/ipd/test_metamorphic_common.py, Python, 462 lines
tests/ipd/test_missingness_calibration.py, Python, 417 lines
tests/ipd/test_missingness_failure_prediction.py, Python, 1,190 lines
tests/ipd/test_package_exports.py, Python, 103 lines
tests/ipd/test_property_admission.py, Python, 1,064 lines
tests/ipd/test_property_common.py, Python, 794 lines
tests/ipd/test_property_manifest.py, Python, 497 lines
tests/ipd/test_property_regions.py, Python, 542 lines
tests/ipd/test_random_samples.py, Python, 180 lines
tests/ipd/test_regions.py, Python, 3,843 lines
tests/ipd/test_tools.py, Python, 266 lines
tests/ld/__init__.py, Python, 1 line
tests/ld/conftest.py, Python, 552 lines
tests/ld/external_clumping_support.py, Python, 295 lines
tests/ld/test_clumping.py, Python, 832 lines
tests/ld/test_cross_validation.py, Python, 414 lines
tests/ld/test_external_clumping_agreement_e2e.py, Python, 262 lines
tests/ld/test_external_clumping_auxiliary_e2e.py, Python, 314 lines
tests/ld/test_external_clumping_e2e.py, Python, 1,137 lines
tests/ld/test_external_clumping_factories.py, Python, 417 lines
tests/ld/test_external_clumping_file_routes.py, Python, 373 lines
tests/ld/test_external_clumping_identity.py, Python, 180 lines
tests/ld/test_external_clumping_selection_e2e.py, Python, 194 lines
tests/ld/test_external_clumping_semantics.py, Python, 1,152 lines
tests/ld/test_external_clumping_sumstats_e2e.py, Python, 302 lines
tests/ld/test_external_clumping_thresholds.py, Python, 1,264 lines
tests/ld/test_external_clumping_unresolved.py, Python, 623 lines
tests/ld/test_factory_clumping.py, Python, 408 lines
tests/ld/test_file_clumping.py, Python, 323 lines
tests/ld/test_handlers.py, Python, 216 lines
tests/ld/test_helpers.py, Python, 85 lines
tests/ld/test_multiprocess.py, Python, 337 lines
tests/ld/test_variants.py, Python, 84 lines
tests/misc/test_sojo.py, Python, 83 lines
tests/misc/test_steiger_filtering.py, Python, 70 lines
tests/mr/__init__.py, Python, 1 line
tests/mr/test_api_exports.py, Python, 24 lines
tests/mr/test_egger.py, Python, 126 lines
tests/mr/test_ivw.py, Python, 166 lines
tests/mr/test_mendelian_randomization.py, Python, 1,066 lines
tests/mr/test_mvmr.py, Python, 143 lines
tests/pipelines/__init__.py, Python, 1 line
tests/pipelines/cli/__init__.py, Python, 1 line
tests/pipelines/cli/conftest.py, Python, 49 lines
tests/pipelines/cli/test_args.py, Python, 921 lines
tests/pipelines/cli/test_command_skeletons.py, Python, 106 lines
tests/pipelines/cli/test_log.py, Python, 178 lines
tests/pipelines/cli/test_output_hash_reproducibility.py, Python, 102 lines
tests/pipelines/cli/test_run_coloc.py, Python, 675 lines
tests/pipelines/cli/test_run_extract.py, Python, 921 lines
tests/pipelines/cli/test_run_mr.py, Python, 2,011 lines
tests/pipelines/cli/test_run_mr_hpc.py, Python, 680 lines
tests/pipelines/cli/test_run_mr_parallel.py, Python, 809 lines
tests/pipelines/cli/test_shared_cli.py, Python, 1,781 lines
tests/pipelines/cli/test_subprocess_contracts.py, Python, 1,594 lines
tests/pipelines/cli/test_typed_ld_reference_e2e.py, Python, 676 lines
tests/pipelines/conftest.py, Python, 639 lines
tests/pipelines/hpc/__init__.py, Python, 1 line
tests/pipelines/hpc/test_jaf.py, Python, 412 lines
tests/pipelines/hpc/test_planner.py, Python, 589 lines
tests/pipelines/hpc/test_task_script.py, Python, 426 lines
tests/pipelines/methods/__init__.py, Python, 1 line
tests/pipelines/methods/test_config.py, Python, 343 lines
tests/pipelines/methods/test_egger.py, Python, 370 lines
tests/pipelines/methods/test_ivw.py, Python, 331 lines
tests/pipelines/methods/test_median.py, Python, 391 lines
tests/pipelines/methods/test_min_snps.py, Python, 61 lines
tests/pipelines/methods/test_mode.py, Python, 541 lines
tests/pipelines/methods/test_registry.py, Python, 613 lines
tests/pipelines/methods/test_support.py, Python, 280 lines
tests/pipelines/methods/test_wald.py, Python, 206 lines
tests/pipelines/mr_results_report_contract.py, Python, 177 lines
tests/pipelines/parallel/__init__.py, Python, 1 line
tests/pipelines/parallel/test_executor.py, Python, 964 lines
tests/pipelines/parallel/test_publish.py, Python, 132 lines
tests/pipelines/parallel/test_runtime.py, Python, 477 lines
tests/pipelines/parallel/test_worker.py, Python, 551 lines
tests/pipelines/test_apply_maf.py, Python, 816 lines
tests/pipelines/test_assembly.py, Python, 569 lines
tests/pipelines/test_coloc.py, Python, 1,579 lines
tests/pipelines/test_coloc_e2e.py, Python, 1,281 lines
tests/pipelines/test_columns.py, Python, 58 lines
tests/pipelines/test_constants.py, Python, 197 lines
tests/pipelines/test_extract.py, Python, 1,785 lines
tests/pipelines/test_features.py, Python, 649 lines
tests/pipelines/test_filtering.py, Python, 2,344 lines
tests/pipelines/test_fn17_e2e.py, Python, 228 lines
tests/pipelines/test_formatting.py, Python, 278 lines
tests/pipelines/test_human_chromosomes.py, Python, 68 lines
tests/pipelines/test_io.py, Python, 685 lines
tests/pipelines/test_io_integration.py, Python, 149 lines
tests/pipelines/test_ld_reference.py, Python, 1,238 lines
tests/pipelines/test_maf_crosscheck.py, Python, 289 lines
tests/pipelines/test_maf_source.py, Python, 525 lines
tests/pipelines/test_mr.py, Python, 5,018 lines
tests/pipelines/test_mr_e2e.py, Python, 2,111 lines
tests/pipelines/test_mr_results_contract.py, Python, 89 lines
tests/pipelines/test_mr_results_report.py, Python, 1,106 lines
tests/pipelines/test_outcome_loader.py, Python, 543 lines
tests/pipelines/test_outcome_loader_isolation.py, Python, 263 lines
tests/pipelines/test_package_exports.py, Python, 226 lines
tests/pipelines/test_parameters.py, Python, 599 lines
tests/pipelines/test_result.py, Python, 62 lines
tests/test_ci_image_config.py, Python, 89 lines
tests/test_constants.py, Python, 52 lines
tests/test_errors_formatting.py, Python, 80 lines
tests/test_imports.py, Python, 62 lines
tests/test_pyproject_metadata.py, Python, 127 lines
tests/test_stats.py, Python, 36 lines
tests/test_version_config.py, Python, 35 lines
LICENCE.txt, License, 674 lines
README.md, Text, 54 lines

At the source: gitlab.com/cfinan/merit

cfinan/bio-misc

License: GPL-3.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 4060630a617911862a403c065d2cd0f5bebe0cc3, 15 February 2024

Languages: Python (87), Shell (14), Jupyter (1)

Size: 268 files, 102 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, environment (requirements.txt, setup.py, resources/conda/envs/py310/conda_create.yml, resources/conda/envs/py310/conda_update.yml, resources/conda/envs/py38/conda_create.yml, resources/conda/envs/py38/conda_update.yml, resources/conda/envs/py39/conda_create.yml, resources/conda/envs/py39/conda_update.yml), tests, continuous integration, documentation, 1 notebook

Not found: CITATION.cff

Tools: NumPy (19 files), pandas (8 files), pysam (3 files), Biopython (2 files), Numba (2 files), BCFtools (1 file), SciPy (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

104 files

_version.py, Python, 1 line
bio_misc/__init__.py, Python, 4 lines
bio_misc/annotation/__init__.py, Python, 1 line
bio_misc/annotation/boolean_matrix.py, Python, 428 lines
bio_misc/annotation/genecards.py, Python, 540 lines
bio_misc/annotation/vcf_id_mapper.py, Python, 239 lines
bio_misc/common/__init__.py, Python, 1 line
bio_misc/common/arguments.py, Python, 394 lines
bio_misc/common/constants.py, Python, 35 lines
bio_misc/common/log.py, Python, 248 lines
bio_misc/common/utils.py, Python, 717 lines
bio_misc/deprecated.py, Python, 8 lines
bio_misc/drug_lookups/__init__.py, Python, 1 line
bio_misc/drug_lookups/bnf_queries.py, Python, 67 lines
bio_misc/drug_lookups/chembl_queries.py, Python, 403 lines
bio_misc/drug_lookups/common.py, Python, 2 lines
bio_misc/drug_lookups/search_effects.py, Python, 595 lines
bio_misc/drug_lookups/target_effects.py, Python, 489 lines
bio_misc/drug_lookups/target_summary.py, Python, 766 lines
bio_misc/drug_lookups/therapeutic_direction.py, Python, 441 lines, 1 match
bio_misc/example_data/__init__.py, Python, 1 line
bio_misc/example_data/example_datasets/__init__.py, Python, 1 line
bio_misc/example_data/example_datasets/clump_data/__init__.py, Python, 1 line
bio_misc/example_data/example_datasets/ipd_data/__init__.py, Python, 1 line
bio_misc/example_data/example_datasets/overlaps/__init__.py, Python, 1 line
bio_misc/example_data/example_datasets/zarr_test/__init__.py, Python, 1 line
bio_misc/example_data/examples.py, Python, 1,296 lines
bio_misc/format/__init__.py, Python, 1 line
bio_misc/format/chr_pos.py, Python, 231 lines
bio_misc/format/common.py, Python, 268 lines
bio_misc/format/uni_id.py, Python, 246 lines
bio_misc/id_map/__init__.py, Python, 1 line
bio_misc/id_map/ensembl_gene.py, Python, 287 lines
bio_misc/id_map/ensembl_human_gene_uniprot.py, Python, 261 lines
bio_misc/id_map/ensembl_human_uniprot_gene.py, Python, 237 lines
bio_misc/id_map/ensembl_uniprot.py, Python, 334 lines
bio_misc/id_map/somalogic_seq_ids.py, Python, 232 lines
bio_misc/ipd/__init__.py, Python, 1 line
bio_misc/ipd/common.py, Python, 1,171 lines
bio_misc/ipd/corr_matrix.py, Python, 600 lines
bio_misc/ipd/corr_matrix_single.py, Python, 345 lines
bio_misc/ipd/corr_matrix_zarr.py, Python, 294 lines
bio_misc/ipd/dosage_zarr.py, Python, 341 lines
bio_misc/ipd/extract_allele_counts.py, Python, 311 lines
bio_misc/ipd/index_bgen.py, Python, 108 lines
bio_misc/ipd/ld_clump.py, Python, 1,478 lines
bio_misc/ipd/ld_factory.py, Python, 159 lines
bio_misc/ipd/ld_handlers.py, Python, 1,078 lines
bio_misc/ipd/readers.py, Python, 1,613 lines
bio_misc/liftover/__init__.py, Python, 1 line
bio_misc/liftover/crossmap.py, Python, 270 lines
bio_misc/liftover/quick_lift.py, Python, 357 lines
bio_misc/overlaps/__init__.py, Python, 1 line
bio_misc/overlaps/genes_within.py, Python, 293 lines
bio_misc/overlaps/nearest_genes.py, Python, 240 lines
bio_misc/overlaps/site_overlap.py, Python, 693 lines
bio_misc/overlaps/trans_qtl_link.py, Python, 432 lines
bio_misc/overlaps/trans_within.py, Python, 307 lines
bio_misc/pubmed/__init__.py, Python, 1 line
bio_misc/pubmed/pubmed_util.py, Python, 1,258 lines
bio_misc/pubmed/query_genes.py, Python, 694 lines
bio_misc/reactome/__init__.py, Python, 1 line
bio_misc/reactome/common.py, Python, 247 lines
bio_misc/reactome/drug_path_dist.py, Python, 226 lines
bio_misc/reactome/pairwise_path_dist.py, Python, 132 lines
bio_misc/reactome/set_path_dist.py, Python, 174 lines
bio_misc/transforms/__init__.py, Python, 1 line
bio_misc/transforms/effect_flipper.py, Python, 544 lines
bio_misc/transforms/log_pvalue.py, Python, 272 lines
bio_misc/umls/__init__.py, Python, 1 line
bio_misc/umls/mesh_mapping.py, Python, 318 lines
docs/source/conf.py, Python, 103 lines
resources/bin/cat-bgen-wrap.sh, Shell, 22 lines
resources/bin/liftover-bgen.sh, Shell, 280 lines
resources/ci_cd/before_script.sh, Shell, 26 lines
resources/ci_cd/pages.sh, Shell, 32 lines
resources/ci_cd/run_docker.sh, Shell, 57 lines
resources/ci_cd/run_pages.sh, Shell, 29 lines
resources/conda/build/py310/build.sh, Shell, 3 lines
resources/conda/build/py310/run_test.sh, Shell, 5 lines
resources/conda/build/py38/build.sh, Shell, 3 lines
resources/conda/build/py38/run_test.sh, Shell, 5 lines
resources/conda/build/py39/build.sh, Shell, 3 lines
resources/conda/build/py39/run_test.sh, Shell, 5 lines
resources/examples/bgen_reader_test.ipynb, Jupyter, 196 lines
resources/tasks/cat-bgen.task.sh, Shell, 43 lines
resources/tasks/liftover-bgen.task.sh, Shell, 76 lines
resources/test_scripts/make_clump_test.py, Python, 202 lines
resources/test_scripts/test_bgen_read.py, Python, 110 lines
resources/test_scripts/test_corr_dp.py, Python, 504 lines
setup.py, Python, 220 lines
tests/test_imports.py, Python, 65 lines
tests/test_ipd/conftest.py, Python, 246 lines
tests/test_ipd/test_common.py, Python, 82 lines
tests/test_ipd/test_corr_matrix.py, Python, 374 lines
tests/test_ipd/test_corrcoef.py, Python, 113 lines
tests/test_ipd/test_ld_clump.py, Python, 235 lines
tests/test_ipd/test_np_commands.py, Python, 210 lines
tests/test_overlaps/conftest.py, Python, 65 lines
tests/test_overlaps/test_genes_within.py, Python, 67 lines
tests/test_overlaps/test_site_overlaps.py, Python, 90 lines
tests/test_overlaps/test_trans_within.py, Python, 68 lines
LICENSE.txt, License, 674 lines
README.md, Text, 79 lines

At the source: gitlab.com/cfinan/bio-misc

cfinan/gwas-norm

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 04b05b9eeb43a67266e4350a0e27f75ec0f2e0a3, 26 July 2026

Languages: Python (110), Shell (20), Jupyter (2)

Size: 659 files, 132 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, environment (pixi.lock, pixi.toml, pyproject.toml, resources/conda/conda_build_config.yaml), tests, continuous integration, documentation, 2 notebooks

Not found: CITATION.cff

Tools: Biopython (10 files), NumPy (10 files), pysam (9 files), pandas (4 files), SciPy (4 files), BCFtools (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

134 files

docs/source/conf.py, Python, 86 lines
gwas_norm/__init__.py, Python, 1 line
gwas_norm/_version.py, Python, 1 line
gwas_norm/columns.py, Python, 1,083 lines
gwas_norm/common.py, Python, 1,167 lines
gwas_norm/constants.py, Python, 318 lines
gwas_norm/controllers.py, Python, 2,368 lines
gwas_norm/crossmap.py, Python, 647 lines
gwas_norm/errors.py, Python, 28 lines
gwas_norm/example_data/__init__.py, Python, 1 line
gwas_norm/example_data/examples.py, Python, 817 lines
gwas_norm/files.py, Python, 1,314 lines
gwas_norm/genome_assemblies.py, Python, 196 lines
gwas_norm/gwas_norm_run.py, Python, 1,093 lines
gwas_norm/handlers.py, Python, 1,619 lines
gwas_norm/hpc/__init__.py, Python, 1 line
gwas_norm/hpc/hpc_gwas_norm.py, Python, 712 lines
gwas_norm/metadata/__init__.py, Python, 1 line
gwas_norm/metadata/analysis.py, Python, 1,086 lines
gwas_norm/metadata/base.py, Python, 550 lines
gwas_norm/metadata/cohort.py, Python, 1,397 lines
gwas_norm/metadata/column.py, Python, 251 lines
gwas_norm/metadata/convert.py, Python, 964 lines
gwas_norm/metadata/file.py, Python, 1,633 lines
gwas_norm/metadata/gwas_data.py, Python, 503 lines
gwas_norm/metadata/info.py, Python, 247 lines
gwas_norm/metadata/phenotype.py, Python, 955 lines
gwas_norm/metadata/study.py, Python, 1,505 lines
gwas_norm/metadata/test.py, Python, 623 lines
gwas_norm/normalise.py, Python, 3,150 lines, 1 match
gwas_norm/parsers.py, Python, 890 lines
gwas_norm/stats.py, Python, 182 lines
gwas_norm/tests.py, Python, 4,684 lines
resources/bin/list-norm-files.sh, Shell, 99 lines
resources/ci_cd/before_script.sh, Shell, 23 lines
resources/ci_cd/pages.sh, Shell, 32 lines
resources/ci_cd/pypi-upload.sh, Shell, 8 lines
resources/ci_cd/run_docker.sh, Shell, 57 lines
resources/ci_cd/run_pages.sh, Shell, 32 lines
resources/examples/metadata_api.ipynb, Jupyter, 291 lines
resources/misc/cohort.py, Python, 108 lines
resources/misc/convert_xmls.sh, Shell, 115 lines
resources/misc/extract-xml-info.py, Python, 372 lines
resources/misc/make-test-data.py, Python, 729 lines
resources/misc/old/old_code/1kg.py, Python, 485 lines
resources/misc/old/old_code/alfa.py, Python, 274 lines
resources/misc/old/old_code/arrays.py, Python, 847 lines
resources/misc/old/old_code/batch_proc_old.py, Python, 342 lines
resources/misc/old/old_code/define_match.py, Python, 119 lines
resources/misc/old/old_code/liftover_test.ipynb, Jupyter, 69 lines
resources/misc/old/old_code/mapper_file_code.py, Python, 283 lines
resources/misc/old/old_code/mappers.py, Python, 3,020 lines
resources/misc/old/old_code/mapping_files.py, Python, 1,248 lines
resources/misc/old/old_code/mapping_utils.py, Python, 388 lines
resources/misc/old/old_code/merge_1kg_freq.sh, Shell, 10 lines
resources/misc/old/old_code/process_alfa.sh, Shell, 29 lines
resources/misc/old/old_code/process_ensembl_dbsnp.sh, Shell, 33 lines
resources/misc/old/old_code/process_ukbb.sh, Shell, 17 lines
resources/misc/old/old_code/test_config_obj.py, Python, 388 lines
resources/misc/old/old_code/test_config_parser.py, Python, 250 lines
resources/misc/old/old_code/vcf.py, Python, 12 lines
resources/misc/old/tests/utils/conftest.py, Python, 49 lines
resources/misc/old/tests/utils/test_upgrade_norm.py, Python, 221 lines
resources/misc/old/v0.2-code/config.py, Python, 604 lines
resources/misc/old/v0.2-code/database/__init__.py, Python, 1 line
resources/misc/old/v0.2-code/database/add_to_db.py, Python, 594 lines
resources/misc/old/v0.2-code/external.py, Python, 56 lines
resources/misc/old/v0.2-code/gwas_norm_premap.py, Python, 935 lines
resources/misc/old/v0.2-code/log.py, Python, 248 lines
resources/misc/old/v0.2-code/map_file/__init__.py, Python, 1 line
resources/misc/old/v0.2-code/map_file/dbsnp/__init__.py, Python, 1 line
resources/misc/old/v0.2-code/map_file/dbsnp/clinvar.py, Python, 391 lines
resources/misc/old/v0.2-code/map_file/dbsnp/download.py, Python, 813 lines
resources/misc/old/v0.2-code/map_file/dbsnp/parse.py, Python, 2,450 lines
resources/misc/old/v0.2-code/old/config_obj.py, Python, 1,316 lines
resources/misc/old/v0.2-code/old/config_parser.py, Python, 2,105 lines
resources/misc/old/v0.2-code/old/crossmap.py, Python, 975 lines
resources/misc/old/v0.2-code/old/crossmap_api.py, Python, 273 lines
resources/misc/old/v0.2-code/old/gtex_xml_writer.py, Python, 281 lines
resources/misc/old/v0.2-code/old/readers.py, Python, 746 lines
resources/misc/old/v0.2-code/old/writers.py, Python, 1 line
resources/misc/old/v0.2-code/processors.py, Python, 1,436 lines
resources/misc/old/v0.2-code/sort.py, Python, 676 lines
resources/misc/old/v0.2-code/utils/__init__.py, Python, 1 line
resources/misc/old/v0.2-code/utils/extract_pvalue.py, Python, 334 lines
resources/misc/old/v0.2-code/utils/genome_chunks.py, Python, 748 lines
resources/misc/old/v0.2-code/utils/make_extract_pvalue_job_array.py, Python, 133 lines
resources/misc/old/v0.2-code/utils/quick_lift.py, Python, 343 lines
resources/misc/old/v0.2-code/utils/upgrade_norm.py, Python, 1,658 lines
resources/misc/old/v0.2-code/utils/upgrade_top_hits.py, Python, 1,745 lines
resources/misc/test-norm-structure.py, Python, 32 lines
resources/misc/test_decompress.py, Python, 16 lines
resources/misc/test_variants/create_test_mapper.py, Python, 44 lines
resources/misc/test_variants/ensembl-cmd-line.sh, Shell, 24 lines
resources/misc/test_variants/other-impute.sh, Shell, 41 lines
resources/misc/test_variants/scan-cmd-line.sh, Shell, 40 lines
resources/misc/test_variants/scan-suhre.sh, Shell, 41 lines
resources/misc/test_variants/tabix-cmd-line.sh, Shell, 42 lines
resources/misc/test_variants/test_mapper.py, Python, 33 lines
resources/misc/update-xml-cohorts.py, Python, 528 lines
resources/misc/update-xmls.py, Python, 789 lines
resources/tasks/gwas-norm.task.sh, Shell, 196 lines
tests/__init__.py, Python, 1 line
tests/conftest.py, Python, 491 lines
tests/data/dummy_master/rename.sh, Shell, 5 lines
tests/data/ref_genome_tests/old/ref_vcf_setup.sh, Shell, 102 lines
tests/data/ref_genome_tests/ref_vcf_setup.sh, Shell, 177 lines
tests/metadata/conftest.py, Python, 193 lines
tests/metadata/test_analysis.py, Python, 417 lines
tests/metadata/test_base.py, Python, 53 lines
tests/metadata/test_cohort.py, Python, 539 lines
tests/metadata/test_column.py, Python, 112 lines
tests/metadata/test_convert.py, Python, 111 lines
tests/metadata/test_file.py, Python, 453 lines
tests/metadata/test_gwas_data.py, Python, 97 lines
tests/metadata/test_metadata_info.py, Python, 86 lines
tests/metadata/test_phenotypes.py, Python, 159 lines
tests/metadata/test_study.py, Python, 826 lines
tests/metadata/test_test.py, Python, 562 lines
tests/test_columns.py, Python, 342 lines
tests/test_controllers.py, Python, 600 lines
tests/test_crossmap.py, Python, 140 lines
tests/test_gwas_norm.py, Python, 504 lines
tests/test_info_columns_without_info_element.py, Python, 185 lines
tests/test_info_handlers.py, Python, 220 lines
tests/test_merge_log_index_error.py, Python, 176 lines
tests/test_normalise.py, Python, 1,190 lines
tests/test_parsers.py, Python, 655 lines
tests/test_safe_writerow.py, Python, 160 lines
tests/test_stats.py, Python, 122 lines
tests/test_target_genome_assembly_normalise.py, Python, 110 lines
tests/test_test_data.py, Python, 431 lines
LICENCE.txt, License, 21 lines
README.md, Text, 35 lines

At the source: gitlab.com/cfinan/gwas-norm

chembl/tractability_pipeline_v2

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 84f53f49664455f89a802f0e2caf7927f7ff06c9, 9 February 2025

Languages: Python (14)

Size: 67 files, 14 scripts

Software Heritage: archived

Found in: “Code availability”

Holds: README, license file, environment (requirements.txt, setup.py, ot_tractability_pipeline_v2/SpaCy_NER_PROTAC_model_packaged/en_NER_PROTAC-0.2.5/setup.py)

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (5 files), pandas (5 files)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

16 files

ot_tractability_pipeline_v2/SpaCy_NER_PROTAC_model_packaged/en_NER_PROTAC-0.2.5/en_NER_PROTAC/__init__.py, Python, 12 lines
ot_tractability_pipeline_v2/SpaCy_NER_PROTAC_model_packaged/en_NER_PROTAC-0.2.5/setup.py, Python, 64 lines
ot_tractability_pipeline_v2/__init__.py, Python, 1 line
ot_tractability_pipeline_v2/bin/__init__.py, Python, 1 line
ot_tractability_pipeline_v2/bin/run_pipeline.py, Python, 825 lines
ot_tractability_pipeline_v2/buckets_ab.py, Python, 759 lines
ot_tractability_pipeline_v2/buckets_othercl.py, Python, 346 lines
ot_tractability_pipeline_v2/buckets_protac.py, Python, 1,318 lines
ot_tractability_pipeline_v2/buckets_sm.py, Python, 753 lines
ot_tractability_pipeline_v2/queries_ab.py, Python, 49 lines
ot_tractability_pipeline_v2/queries_othercl.py, Python, 52 lines
ot_tractability_pipeline_v2/queries_protac.py, Python, 106 lines
ot_tractability_pipeline_v2/queries_sm.py, Python, 373 lines
setup.py, Python, 26 lines
LICENSE, License, 21 lines
README.md, Text, 256 lines

At the source: github.com/chembl/tractability_pipeline_v2

Code availability

Cis-Mendelian randomisation and colocalisation analyses were performed using purpose-built Python (v3.9) packages (Mendelian Randomisation for Identifying Targets [MeRIT], v0.3.3a0; https://gitlab.com/cfinan/merit), which incorporates the coloc R package for colocalisation (https://github.com/chr1swallace/coloc). Estimated therapeutic relevance and known depression-related drug effects were derived using bio-misc (v0.2.0a0; https://gitlab.com/cfinan/bio-misc). GWAS summary statistics were processed and standardised using gwas-norm (https://gitlab.com/cfinan/gwas-norm). Druggable genome targets were defined using Open Targets tractability assessments (https://github.com/chembl/tractability_pipeline_v2; https://platform-docs.opentargets.org/target/tractability). Protein-protein interaction network analysis was conducted using STRING (v12.0; https://string-db.org).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
722 scripts, each with its path and the digest of its content;
5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

ukbiobank.ac.uk/enable-your-research/apply-for-access, at UK Biobank; found in “Data availability”
zenodo:19135140, at Zenodo; found in DataCite
zenodo:19135141, at Zenodo; found in DataCite

Data availability

GWAS summary statistics for MD are available from the Psychiatric Genomics Consortium at https://pgc.unc.edu/for-researchers/download-results/. GWAS summary statistics including 23andMe data require an approved application through 23andMe available to qualified researchers under an agreement with 23andMe that protects the privacy of the 23andMe participants (visit https://research.23andme.com/dataset-access/). UK Biobank data are available through application at https://www.ukbiobank.ac.uk/enable-your-research/apply-for-access/. Open Targets tractability data (version dated 2024-05-23) can be downloaded from http://ftp.ebi.ac.uk/pub/databases/opentargets/platform/latest/input/target/tractability/. QTL datasets are available from: Blood plasma pQTL data: deCODE (https://www.decode.com/summarydata/), UKB-PPP (https://www.synapse.org/#!Synapse:syn51364943/), INTERVAL (http://www.phpc.cam.ac.uk/ceu/proteins/), and Gudjonsson (https://www.ebi.ac.uk/gwas/publications/35078996). Brain pQTL data: ROSMAP and Banner (https://www.synapse.org/#!Synapse:syn24172458). CSF pQTL data: Yang (https://dss.niagads.org/datasets/ng00102/). Blood eQTL data: eQTLGen (https://www.eqtlgen.org/). Brain eQTL data: MetaBrain (https://www.metabrain.nl/). Compound mechanism of action and drug effect data were obtained from ChEMBL (v33; https://www.ebi.ac.uk/chembl/). Cis-MR results and genetic instruments for all targets tested, including non-significant results, are available on Zenodo (10.5281/zenodo.19135141).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 3 keywords, 7 MeSH terms, 4 funders, 101 references.

Cite

This paper

ter Kuile, A. R., Finan, C., Chopade, S., van Vugt, M., Hukerikar, N., Barral, S., Stringaris, A., Schmidt, A. F., Kuchenbaecker, K., & Pingault, J.-B. (2026). An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression. Translational psychiatry, 16(1), 392. https://doi.org/10.1038/s41398-026-04137-9

BibTeX

@article{terkuile2026integrative, author = {ter Kuile, Abigail R. and Finan, Chris and Chopade, Sandesh and van Vugt, Marion and Hukerikar, Nikita and Barral, Serena and Stringaris, Argyris and Schmidt, Amand F. and Kuchenbaecker, Karoline and Pingault, Jean-Baptiste}, title = {{An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression}}, journal = {Translational psychiatry}, year = {2026}, month = jun, volume = {16}, number = {1}, pages = {392}, publisher = {Nature Publishing Group}, issn = {2158-3188}, doi = {10.1038/s41398-026-04137-9}, url = {https://doi.org/10.1038/s41398-026-04137-9}, pmid = {42230543}, pmcid = {PMC13443190} }

RIS

TY - JOUR AU - ter Kuile, Abigail R. AU - Finan, Chris AU - Chopade, Sandesh AU - van Vugt, Marion AU - Hukerikar, Nikita AU - Barral, Serena AU - Stringaris, Argyris AU - Schmidt, Amand F. AU - Kuchenbaecker, Karoline AU - Pingault, Jean-Baptiste TI - An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression T2 - Translational psychiatry J2 - Transl Psychiatry PY - 2026 DA - 2026/06/02 VL - 16 IS - 1 SP - 392 SN - 2158-3188 PB - Nature Publishing Group DO - 10.1038/s41398-026-04137-9 UR - https://doi.org/10.1038/s41398-026-04137-9 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41398-026-04137-9", "type": "article-journal", "title": "An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression", "container-title": "Translational psychiatry", "author": [ { "family": "ter Kuile", "given": "Abigail R." }, { "family": "Finan", "given": "Chris" }, { "family": "Chopade", "given": "Sandesh" }, { "family": "van Vugt", "given": "Marion" }, { "family": "Hukerikar", "given": "Nikita" }, { "family": "Barral", "given": "Serena" }, { "family": "Stringaris", "given": "Argyris" }, { "family": "Schmidt", "given": "Amand F." }, { "family": "Kuchenbaecker", "given": "Karoline" }, { "family": "Pingault", "given": "Jean-Baptiste" } ], "container-title-short": "Transl Psychiatry", "volume": "16", "issue": "1", "page": "392", "DOI": "10.1038/s41398-026-04137-9", "PMID": "42230543", "PMCID": "PMC13443190", "ISSN": "2158-3188", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41398-026-04137-9", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 2 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41562-026-02476-7 [code]: Genome-wide meta-analysis of quantitatively measured generalized anxiety symptoms in individuals of European ancestry.
Journal: Nature human behaviour
In common: BCFtools, psych, data.table, clinical / translational, genetics / omics, 5 references, author Abigail R. ter Kuile
[2] doi:10.1038/s41467-026-76676-0 [code]: Determinants of functional burden pleiotropy and gene dosage responses across human traits.
Journal: Nature communications
In common: Numba, data.table, statsmodels, 5 other tools, ukbiobank.ac.uk/enable-your-research/apply-for-access, genetics / omics, 1 reference
[3] doi:10.1038/s41592-026-03211-w [code]: Spatial isoform sequencing at single-cell resolution reveals cell-type-specific spatial isoform variability in multiple brain cell types.
Journal: Nature methods
In common: pysam, Biopython, Numba, 7 other tools, genetics / omics
[4] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: BCFtools, Biopython, Numba, 7 other tools
[5] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pysam, rpy2, Numba, 6 other tools, genetics / omics
[6] doi:10.1038/s41592-026-03057-2 [code]: CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: pysam, Biopython, Numba, 6 other tools, genetics / omics
[7] doi:10.1016/j.celrep.2026.117110 [code]: Single-nucleus multiome analysis in the human prefrontal cortex identifies gene expression and cis-regulatory elements associated with aging.
Journal: Cell reports
In common: BCFtools, pysam, statsmodels, 5 other tools, genetics / omics, 1 reference
[8] doi:10.1038/s43587-026-01106-1 [code]: Repurposing drugs for the prevention of vascular dementia using evidence from drug target Mendelian randomization.
Journal: Nature aging
In common: genetics / omics, 7 references
[9] doi:10.1038/s41562-026-02486-5 [code]: Genome-wide association studies of infant and toddler temperament in European and multi-ancestry populations.
Journal: Nature human behaviour
In common: Biopython, psych, data.table, 5 other tools, genetics / omics, 1 reference
[10] doi:10.7554/elife.108109 [code]: Multimodal MRI marker of cognition explains the association between cognition and mental health in the UK Biobank.
Journal: eLife
In common: psych, data.table, seaborn, 4 other tools, ukbiobank.ac.uk/enable-your-research/apply-for-access

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

cfinan/merit

cfinan/bio-misc

cfinan/gwas-norm

chembl/tractability_pipeline_v2

Code availability

Tracing map

Data

Datasets cited

Data availability

Versions

Version 1, 27 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity