Skip to content

Commit

Permalink
Accidentally added the coloc and other MR functionality to the functi…
Browse files Browse the repository at this point in the history 
…on which led to a processing speed decrease, now this can be turned on by a new function flag.
  • Loading branch information
@adriaan-vd-graaf
adriaan-vd-graaf committed Dec 3, 2024
1 parent 90e1042 commit 0799a67
Show file tree
Hide file tree
Showing 2 changed files with 42 additions and 22 deletions.
7 changes: 6 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,12 @@ options:
How many SNPs are allowed in the correlation matrix, as it can take a lot of time to
do the eigendecompositions and Rscripts otherwise.
default: 5250
--run_other_cis_mr_and_coloc
This is a flag that you can add to run other MR analyses: MR-IVW, MR-IVW LD and MR-PCA,
as well as coloc and SuSIE coloc. Requires R to be installed with the coloc and data.table
packages. Will append columns of these methods' results to the output file.
--no_normalize_sumstats
flag to _not_ normalize summary statistics
--no_exclude_hla
Expand Down
57 changes: 36 additions & 21 deletions mr_link_2_standalone.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1042,6 +1042,7 @@ def mr_link2_on_region(region: StartEndRegion,
var_explained_grid: list,
verbosity=0,
max_snp_threshold=5250,
run_other_functions = False,
) -> pd.DataFrame:
"""
Expand Down Expand Up @@ -1178,17 +1179,17 @@ def mr_link2_on_region(region: StartEndRegion,
mafs = np.zeros_like(exp_betas)
mafs[:] = 0.5


coloc_results = run_colocalization_analysis(exp_betas, out_betas, exp_pvals, out_pvals, n_exp, n_out, mafs,
correlation_mat, out_file = tmp_prepend + '_coloc_out',
in_file = tmp_prepend + '_coloc_in', ld_file = tmp_prepend + 'coloc_ld.bin')
instruments = np.zeros(exp_betas.shape[0], dtype=bool)
tmp_instruments = np.asarray(select_instruments_by_clumping(exp_pvals, correlation_mat, 5e-8, r_threshold=0.1), dtype=int)
for instrument in tmp_instruments:
instruments[instrument] = True
r_mr_results = run_external_mr_analysis(exp_betas, out_betas, exp_ses, out_ses, exp_pvals, out_pvals, n_exp, n_out, mafs,
correlation_mat,instruments= instruments, out_file = tmp_prepend + '_mr_out',
in_file = tmp_prepend + '_mr_in', ld_file = tmp_prepend + 'mr_ld.bin')
if run_other_functions:
coloc_results = run_colocalization_analysis(exp_betas, out_betas, exp_pvals, out_pvals, n_exp, n_out, mafs,
correlation_mat, out_file = tmp_prepend + '_coloc_out',
in_file = tmp_prepend + '_coloc_in', ld_file = tmp_prepend + 'coloc_ld.bin')
instruments = np.zeros(exp_betas.shape[0], dtype=bool)
tmp_instruments = np.asarray(select_instruments_by_clumping(exp_pvals, correlation_mat, 5e-8, r_threshold=0.1), dtype=int)
for instrument in tmp_instruments:
instruments[instrument] = True
r_mr_results = run_external_mr_analysis(exp_betas, out_betas, exp_ses, out_ses, exp_pvals, out_pvals, n_exp, n_out, mafs,
correlation_mat,instruments= instruments, out_file = tmp_prepend + '_mr_out',
in_file = tmp_prepend + '_mr_in', ld_file = tmp_prepend + 'mr_ld.bin')

print(' Starting MR-link2')

Expand Down Expand Up @@ -1257,8 +1258,9 @@ def mr_link2_on_region(region: StartEndRegion,

this_run['var_explained'] = var_explained

this_run = this_run.join(r_mr_results, rsuffix='')
this_run = this_run.join(coloc_results, rsuffix='')
if run_other_functions:
this_run = this_run.join(r_mr_results, rsuffix='')
this_run = this_run.join(coloc_results, rsuffix='')

results_list.append(this_run)

Expand All @@ -1269,13 +1271,15 @@ def mr_link2_on_region(region: StartEndRegion,


columns = ['region', 'var_explained', 'm_snps_overlap',
'alpha', 'se(alpha)', 'p(alpha)',
'sigma_y', 'se(sigma_y)', 'p(sigma_y)',
'sigma_x', 'function_time',
'beta_ivw', 'se_ivw', 'p_ivw',
'beta_ivw_r', 'se_ivw_r', 'p_ivw_r',
'beta_pca', 'se_pca', 'p_pca',
'PP.H1.abf', 'PP.H2.abf', 'PP.H3.abf', 'PP.H4.abf']
'alpha', 'se(alpha)', 'p(alpha)',
'sigma_y', 'se(sigma_y)', 'p(sigma_y)',
'sigma_x', 'function_time',]

if run_other_functions:
columns +=['beta_ivw', 'se_ivw', 'p_ivw',
'beta_ivw_r', 'se_ivw_r', 'p_ivw_r',
'beta_pca', 'se_pca', 'p_pca',
'PP.H1.abf', 'PP.H2.abf', 'PP.H3.abf', 'PP.H4.abf']

if 'susie_max_PP.H4.abf' in mr_results_df.columns:
columns += ['susie_max_PP.H1.abf',
Expand Down Expand Up @@ -1600,6 +1604,15 @@ def read_regions_from_file(filename, region_or_regions):
'Do the eigendecompositions and Rscripts otherwise.'
)

parser.add_argument('--run_other_cis_mr_and_coloc',
required=False,
action='store_true',
help='By setting this flag, we will run other _cis_ MR functions as well as coloc and SuSIE '
'coloc. Warning, can be slow '
)



parser.add_argument('--verbose',
default=0,
help='Set to 1 if you want to read more output, for debugging purposes ')
Expand Down Expand Up @@ -1865,7 +1878,9 @@ def read_regions_from_file(filename, region_or_regions):
tmp_prepend=tmp_dir,
verbosity=verbosity,
var_explained_grid=var_explained_grid,
max_snp_threshold=max_correlation_matrix_snps)
max_snp_threshold=max_correlation_matrix_snps,
run_other_functions=args.run_other_cis_mr_and_coloc,
)
except Exception as x:
exceptions.append((region, x))
print(f'Unable to make an MR-link2 estimate in {region} due to {x}')
Expand Down

0 comments on commit 0799a67

Please sign in to comment.