Extract genome-wide data for multivariate analysis

Extract genome-wide data for multivariate analysis#

Extracts strong, random and null sets of effects from genome-wide summary statistics and writes them as a single MASH-format input.

Overview#

The workflow produces three z-score matrices used to fit a multivariate prior:

  • Strong effects (Z_s) — extracted from genome-wide cis fine-mapping results. The top locus per condition is taken (credible-set coverage threshold 0.7) and the z-scores are merged into one data frame.

  • Null effects (Z_n) — up to M candidate SNPs per region satisfying |z| ≤ 2 are overlapped with the independent-SNP list to keep only independent variants; the union is taken.

  • Random effects (Z_r) — variants randomly sampled from the supplied independent-variant list.

These matrices feed Ultimate Deconvolution / MASH to learn patterns of effect sharing across conditions (e.g. cell types or tissues).

Requires genome-wide SuSiE fine-mapping outputs. See the Important Note below for toy-data limitations.

When to run it. After per-region fine-mapping or tensorQTL association scans have produced results across conditions, and before mash_fit.

Input#

susie_to_mash - builds the MASH input from fine-mapped SuSiE results:

  • --fine_mapping_meta input/finemapping/protocol_example.fine_mapping_meta.tsv (index of the fine-mapped RDS files, one row per region and context. Columns chrom, start, end, region_id, gene_id, context, susie_path, sumstats_path, conditions.)

chrom	start	end	region_id	gene_id	context	susie_path	sumstats_path	...
chr22	15528191	15529138	chr22:15528191-15529138	ENSG00000130538	De_Jager	protocol_example.susie_fit.rds	protocol_example.sumstats_db.rds	...

  • --finemapping_column susie_path (which column of the meta table holds the fine-mapping RDS path.)

  • --sig_p_cutoff 0.1 (p-value threshold for calling a variant strong. The toy signal is weak, so loosen this or every strong variant is filtered out.)

random_null_tensorqtl - extracts random and null effects directly from tensorQTL summary statistics:

  • --region_file input/finemapping/protocol_example.region (genome region partition defining the analysis units, columns chrom, start, end, gene_id.)

chrom	start	end	gene_id
chr11	187771	3551668	ENSG00000142089

  • --sum_files input/protocol_example.sumstats_list.txt (list of the per-trait tensorQTL summary-statistics files. This file does not ship with the repository, so no preview is shown.)

  • --traits bulk_rnaseq (trait or context labels matching the summary-statistics files.)

  • --independent_variant_list (optional; LD-pruned independent SNPs, used to keep only independent null and random variants.)

Note on susie_to_mash inputs#

susie_to_mash consumes one FineMappingResult RDS per region (the S4 output of fine_mapping.R), listed in --fine_mapping_meta under --finemapping_column. For each region it takes the lead variant (max PIP) of every credible set in each condition as the strong set, and samples the random / null background from the region’s full variant universe (pecotmr::mashInput).

Toy MWE. The legacy per-region SuSiE fixtures can be converted once to a QtlFineMappingResult — see tests/fixtures/mash/build_qtl_fine_mapping_result.R, which builds protocol_example.QtlFineMappingResult.rds from the committed toy fine-mapping tables. Point --fine_mapping_meta at that RDS to run the toy demo. The toy signal is weak, so loosen --sig_p_cutoff (e.g. 0.1) or every strong variant is filtered out.

For a genome-wide run, provide a --fine_mapping_meta listing every region’s FineMappingResult RDS.

  • --cwd output/mash_preprocessing (working directory all outputs are written under. Defaults to output.)

  • --name protocol_example (tag used in the output filename and in the per-gene cache directory)

Output#

  • <cwd>/<name>.mash_input.rds (susie_to_mash, or random_null_tensorqtl_2) (a list of 10 elements containing the strong.b, random.b and null.b data frames of z-scores across all conditions, ready for Ultimate Deconvolution or MASH. Both entry points write this same file.)

List of 10
 $ strong.b:'data.frame':	62 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:62] 0.0402 0.0177 0.1852 -0.2188 0.1352 ...
  ..$ Ast_De_Jager_eQTL  : num [1:62] -0.74 -0.215 0.165 -0.263 0.113 ...
  ..$ Oli_De_Jager_eQTL  : num [1:62] -0.0564 0.195 0.1379 -0.1245 -0.061 ...
  ..$ OPC_De_Jager_eQTL  : num [1:62] 0.1748 0.2064 0.1182 -0.0751 0.1023 ...
  ..$ Exc_De_Jager_eQTL  : num [1:62] -0.56004 0.04682 0.2279 -0.13765 0.00667 ...
  ..$ Inh_De_Jager_eQTL  : num [1:62] -0.6306 0.1502 0.1777 -0.2118 -0.0562 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:62] -0.0379 0.1315 0.0361 -0.315 -0.1005 ...
  ..$ PCC_De_Jager_eQTL  : num [1:62] -0.0403 0.0209 0.0614 -0.1755 0.0422 ...
  ..$ AC_De_Jager_eQTL   : num [1:62] -0.13649 -0.00365 0.08409 -0.19717 0.08202 ...
  ..$ Mic_Kellis_eQTL    : num [1:62] -0.396 0 0 0 0 ...
  ..$ Ast_Kellis_eQTL    : num [1:62] -0.9948 0.1155 0.0732 -0.2074 -0.044 ...
  ..$ Oli_Kellis_eQTL    : num [1:62] -0.559 0.387 0.1339 -0.3099 -0.0334 ...
  ..$ OPC_Kellis_eQTL    : num [1:62] -0.11556 0.17907 0.06332 0.06868 0.00266 ...
  ..$ Exc_Kellis_eQTL    : num [1:62] -0.5978 0.0797 0.1463 -0.1659 0.049 ...
  ..$ Inh_Kellis_eQTL    : num [1:62] -0.40719 0.04581 0.12052 -0.2968 0.00677 ...
  ..$ Ast.10_Kellis_eQTL : num [1:62] 0 0 0 0 0 0 0 0 0 0 ...
 $ random.b:'data.frame':	200 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:200] -0.014 -0.0691 -0.0426 -0.1391 0.1639 ...
  ..$ Ast_De_Jager_eQTL  : num [1:200] -0.006636 -0.051634 0.064082 0.000642 -0.040912 ...
  ..$ Oli_De_Jager_eQTL  : num [1:200] -0.0702 -0.0912 0.0985 -0.0492 0.0375 ...
  ..$ OPC_De_Jager_eQTL  : num [1:200] 0.0517 0.0577 -0.0443 0.0677 0.019 ...
  ..$ Exc_De_Jager_eQTL  : num [1:200] 0.0404 -0.02711 0.02452 -0.01076 0.00808 ...
  ..$ Inh_De_Jager_eQTL  : num [1:200] -0.09783 -0.0174 0.01247 0.00424 0.0561 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:200] -0.005966 0.006705 0.000253 -0.004673 0.001346 ...
  ..$ PCC_De_Jager_eQTL  : num [1:200] -0.0113 -0.0248 -0.0239 0.0113 0.0307 ...
  ..$ AC_De_Jager_eQTL   : num [1:200] -0.01516 0.00196 -0.05586 0.01852 0.0188 ...
  ..$ Mic_Kellis_eQTL    : num [1:200] 0.0641 -0.1708 -0.0346 -0.0661 0.0293 ...
  ..$ Ast_Kellis_eQTL    : num [1:200] 0.05854 -0.03466 -0.03777 -0.00221 -0.01986 ...
  ..$ Oli_Kellis_eQTL    : num [1:200] -0.01408 0.02196 0.02674 0.00619 -0.03527 ...
  ..$ OPC_Kellis_eQTL    : num [1:200] -0.018068 -0.026661 -0.051184 0.000499 -0.018648 ...
  ..$ Exc_Kellis_eQTL    : num [1:200] 0.00046 -0.08625 0.02572 0.00127 -0.07891 ...
  ..$ Inh_Kellis_eQTL    : num [1:200] -0.07466 -0.03319 -0.00879 0.06033 -0.02189 ...
  ..$ Ast.10_Kellis_eQTL : num [1:200] 0 0 0 0 0 0 0 0 0 0 ...
 $ null.b  :'data.frame':	200 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:200] -0.0252 -0.0296 -0.1086 -0.0394 -0.2567 ...
  ..$ Ast_De_Jager_eQTL  : num [1:200] -0.0597 0.0402 -0.03 0.0394 -0.1582 ...
  ..$ Oli_De_Jager_eQTL  : num [1:200] -0.0376 0.0668 0.1072 0.0819 -0.3407 ...
  ..$ OPC_De_Jager_eQTL  : num [1:200] 0.058 0.0907 -0.0941 0.0193 -0.0585 ...
  ..$ Exc_De_Jager_eQTL  : num [1:200] 0.01518 -0.00338 0.01137 -0.03883 -0.02836 ...
  ..$ Inh_De_Jager_eQTL  : num [1:200] 0.00196 -0.05773 0.06773 -0.08281 0.29775 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:200] 0.000521 -0.004139 -0.010968 -0.019917 0.081309 ...
  ..$ PCC_De_Jager_eQTL  : num [1:200] -0.0193 0.0273 -0.0237 0.0212 -0.0884 ...
  ..$ AC_De_Jager_eQTL   : num [1:200] 0.001409 0.016122 -0.049309 -0.000438 -0.045112 ...
  ..$ Mic_Kellis_eQTL    : num [1:200] 0.0506 0.0105 0.0444 0.0287 0.4847 ...
  ..$ Ast_Kellis_eQTL    : num [1:200] 1.18e-02 1.34e-06 1.17e-01 1.28e-01 -5.49e-01 ...
  ..$ Oli_Kellis_eQTL    : num [1:200] 0.0742 0.0675 -0.025 -0.0363 0.2675 ...
  ..$ OPC_Kellis_eQTL    : num [1:200] -0.0333 -0.0226 -0.0231 0.0717 -0.0915 ...
  ..$ Exc_Kellis_eQTL    : num [1:200] 0.03593 0.00316 0.01524 -0.0089 0.08042 ...
  ..$ Inh_Kellis_eQTL    : num [1:200] -0.043 -0.0471 0.0418 0.0791 -0.0831 ...
  ..$ Ast.10_Kellis_eQTL : num [1:200] 0 0 0 0 0 0 0 0 0 0 ...
 $ strong.s:'data.frame':	62 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:62] 0.3036 0.0925 0.0569 0.095 0.0626 ...
  ..$ Ast_De_Jager_eQTL  : num [1:62] 0.2298 0.0777 0.0476 0.08 0.053 ...
  ..$ Oli_De_Jager_eQTL  : num [1:62] 0.2397 0.0912 0.0556 0.0944 0.0623 ...
  ..$ OPC_De_Jager_eQTL  : num [1:62] 0.3184 0.1084 0.0672 0.1122 0.0741 ...
  ..$ Exc_De_Jager_eQTL  : num [1:62] 0.1051 0.0577 0.0336 0.0586 0.039 ...
  ..$ Inh_De_Jager_eQTL  : num [1:62] 0.198 0.0861 0.0513 0.0872 0.0581 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:62] 0.0511 0.0366 0.0225 0.0349 0.0249 ...
  ..$ PCC_De_Jager_eQTL  : num [1:62] 0.0844 0.0551 0.0326 0.0543 0.0367 ...
  ..$ AC_De_Jager_eQTL   : num [1:62] 0.0729 0.0483 0.0302 0.0484 0.0343 ...
  ..$ Mic_Kellis_eQTL    : num [1:62] 0.319 1000 1000 1000 1000 ...
  ..$ Ast_Kellis_eQTL    : num [1:62] 0.3174 0.0941 0.054 0.0926 0.0604 ...
  ..$ Oli_Kellis_eQTL    : num [1:62] 0.2342 0.0776 0.0471 0.0783 0.0527 ...
  ..$ OPC_Kellis_eQTL    : num [1:62] 0.3137 0.0972 0.0568 0.0953 0.0631 ...
  ..$ Exc_Kellis_eQTL    : num [1:62] 0.1788 0.0754 0.0442 0.0756 0.0495 ...
  ..$ Inh_Kellis_eQTL    : num [1:62] 0.2481 0.0783 0.0454 0.0762 0.0506 ...
  ..$ Ast.10_Kellis_eQTL : num [1:62] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
 $ null.s  :'data.frame':	200 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:200] 0.0633 0.0667 0.0866 0.0795 0.3139 ...
  ..$ Ast_De_Jager_eQTL  : num [1:200] 0.0486 0.0502 0.0665 0.0603 0.2429 ...
  ..$ Oli_De_Jager_eQTL  : num [1:200] 0.0498 0.0522 0.0671 0.0621 0.2443 ...
  ..$ OPC_De_Jager_eQTL  : num [1:200] 0.0663 0.0688 0.0905 0.0842 0.3302 ...
  ..$ Exc_De_Jager_eQTL  : num [1:200] 0.0225 0.0234 0.0308 0.0282 0.1132 ...
  ..$ Inh_De_Jager_eQTL  : num [1:200] 0.0421 0.044 0.0577 0.0529 0.2107 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:200] 0.0127 0.0133 0.0181 0.0163 0.0726 ...
  ..$ PCC_De_Jager_eQTL  : num [1:200] 0.0201 0.0213 0.0289 0.0259 0.104 ...
  ..$ AC_De_Jager_eQTL   : num [1:200] 0.0183 0.0191 0.0269 0.0237 0.1017 ...
  ..$ Mic_Kellis_eQTL    : num [1:200] 0.0602 0.0619 0.0864 0.078 0.305 ...
  ..$ Ast_Kellis_eQTL    : num [1:200] 0.0543 0.0565 0.0765 0.0707 0.2917 ...
  ..$ Oli_Kellis_eQTL    : num [1:200] 0.0421 0.0444 0.0597 0.0545 0.2213 ...
  ..$ OPC_Kellis_eQTL    : num [1:200] 0.0539 0.0566 0.077 0.0708 0.2842 ...
  ..$ Exc_Kellis_eQTL    : num [1:200] 0.0329 0.0344 0.0457 0.0427 0.1859 ...
  ..$ Inh_Kellis_eQTL    : num [1:200] 0.0461 0.047 0.0634 0.0589 0.2358 ...
  ..$ Ast.10_Kellis_eQTL : num [1:200] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
 $ random.s:'data.frame':	200 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:200] 0.0603 0.096 0.0694 0.0657 0.0769 ...
  ..$ Ast_De_Jager_eQTL  : num [1:200] 0.0466 0.074 0.0521 0.0507 0.0578 ...
  ..$ Oli_De_Jager_eQTL  : num [1:200] 0.048 0.0752 0.0526 0.0517 0.0585 ...
  ..$ OPC_De_Jager_eQTL  : num [1:200] 0.0635 0.1011 0.0711 0.0693 0.0789 ...
  ..$ Exc_De_Jager_eQTL  : num [1:200] 0.0219 0.0345 0.0243 0.0237 0.0266 ...
  ..$ Inh_De_Jager_eQTL  : num [1:200] 0.0405 0.0648 0.0457 0.0439 0.0501 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:200] 0.0125 0.0196 0.0145 0.0139 0.016 ...
  ..$ PCC_De_Jager_eQTL  : num [1:200] 0.0198 0.0305 0.0234 0.021 0.0251 ...
  ..$ AC_De_Jager_eQTL   : num [1:200] 0.0176 0.0284 0.0209 0.0196 0.0229 ...
  ..$ Mic_Kellis_eQTL    : num [1:200] 0.0598 0.0894 0.065 0.0641 0.0709 ...
  ..$ Ast_Kellis_eQTL    : num [1:200] 0.0523 0.0838 0.0587 0.0559 0.0625 ...
  ..$ Oli_Kellis_eQTL    : num [1:200] 0.041 0.0651 0.0453 0.0444 0.0509 ...
  ..$ OPC_Kellis_eQTL    : num [1:200] 0.053 0.0829 0.0582 0.0571 0.0627 ...
  ..$ Exc_Kellis_eQTL    : num [1:200] 0.033 0.0506 0.0348 0.0343 0.0387 ...
  ..$ Inh_Kellis_eQTL    : num [1:200] 0.045 0.0694 0.0492 0.0481 0.0545 ...
  ..$ Ast.10_Kellis_eQTL : num [1:200] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
 $ strong.z:'data.frame':	62 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:62] 0.132 0.191 3.253 -2.304 2.161 ...
  ..$ Ast_De_Jager_eQTL  : num [1:62] -3.22 -2.77 3.46 -3.28 2.13 ...
  ..$ Oli_De_Jager_eQTL  : num [1:62] -0.235 2.139 2.48 -1.319 -0.98 ...
  ..$ OPC_De_Jager_eQTL  : num [1:62] 0.549 1.905 1.759 -0.67 1.38 ...
  ..$ Exc_De_Jager_eQTL  : num [1:62] -5.33 0.811 6.79 -2.35 0.171 ...
  ..$ Inh_De_Jager_eQTL  : num [1:62] -3.185 1.744 3.466 -2.428 -0.967 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:62] -0.741 3.598 1.601 -9.025 -4.042 ...
  ..$ PCC_De_Jager_eQTL  : num [1:62] -0.477 0.379 1.885 -3.231 1.151 ...
  ..$ AC_De_Jager_eQTL   : num [1:62] -1.8726 -0.0757 2.7865 -4.0712 2.3891 ...
  ..$ Mic_Kellis_eQTL    : num [1:62] -1.24 0 0 0 0 ...
  ..$ Ast_Kellis_eQTL    : num [1:62] -3.134 1.227 1.355 -2.239 -0.729 ...
  ..$ Oli_Kellis_eQTL    : num [1:62] -2.386 4.987 2.843 -3.956 -0.634 ...
  ..$ OPC_Kellis_eQTL    : num [1:62] -0.3683 1.8424 1.1147 0.7209 0.0422 ...
  ..$ Exc_Kellis_eQTL    : num [1:62] -3.343 1.056 3.308 -2.194 0.991 ...
  ..$ Inh_Kellis_eQTL    : num [1:62] -1.641 0.585 2.653 -3.893 0.134 ...
  ..$ Ast.10_Kellis_eQTL : num [1:62] 0 0 0 0 0 0 0 0 0 0 ...
 $ random.z:'data.frame':	200 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:200] -0.232 -0.72 -0.614 -2.116 2.132 ...
  ..$ Ast_De_Jager_eQTL  : num [1:200] -0.1425 -0.6982 1.2298 0.0127 -0.7073 ...
  ..$ Oli_De_Jager_eQTL  : num [1:200] -1.463 -1.214 1.873 -0.953 0.641 ...
  ..$ OPC_De_Jager_eQTL  : num [1:200] 0.813 0.571 -0.624 0.976 0.24 ...
  ..$ Exc_De_Jager_eQTL  : num [1:200] 1.844 -0.786 1.009 -0.455 0.304 ...
  ..$ Inh_De_Jager_eQTL  : num [1:200] -2.4158 -0.2683 0.273 0.0964 1.1196 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:200] -0.4757 0.3416 0.0175 -0.3357 0.084 ...
  ..$ PCC_De_Jager_eQTL  : num [1:200] -0.572 -0.813 -1.023 0.536 1.221 ...
  ..$ AC_De_Jager_eQTL   : num [1:200] -0.8634 0.0689 -2.6709 0.9463 0.8218 ...
  ..$ Mic_Kellis_eQTL    : num [1:200] 1.073 -1.91 -0.533 -1.03 0.413 ...
  ..$ Ast_Kellis_eQTL    : num [1:200] 1.1182 -0.4138 -0.6432 -0.0395 -0.3177 ...
  ..$ Oli_Kellis_eQTL    : num [1:200] -0.343 0.337 0.59 0.139 -0.693 ...
  ..$ OPC_Kellis_eQTL    : num [1:200] -0.34097 -0.32177 -0.8792 0.00874 -0.29735 ...
  ..$ Exc_Kellis_eQTL    : num [1:200] 0.0139 -1.705 0.7393 0.0371 -2.0365 ...
  ..$ Inh_Kellis_eQTL    : num [1:200] -1.658 -0.478 -0.179 1.255 -0.402 ...
  ..$ Ast.10_Kellis_eQTL : num [1:200] 0 0 0 0 0 0 0 0 0 0 ...
 $ null.z  :'data.frame':	200 obs. of  16 variables:
  ..$ Mic_De_Jager_eQTL  : num [1:200] -0.398 -0.444 -1.254 -0.495 -0.818 ...
  ..$ Ast_De_Jager_eQTL  : num [1:200] -1.229 0.8 -0.452 0.654 -0.651 ...
  ..$ Oli_De_Jager_eQTL  : num [1:200] -0.755 1.28 1.597 1.319 -1.394 ...
  ..$ OPC_De_Jager_eQTL  : num [1:200] 0.876 1.318 -1.039 0.229 -0.177 ...
  ..$ Exc_De_Jager_eQTL  : num [1:200] 0.675 -0.145 0.37 -1.379 -0.251 ...
  ..$ Inh_De_Jager_eQTL  : num [1:200] 0.0466 -1.3118 1.1748 -1.5654 1.4129 ...
  ..$ DLPFC_De_Jager_eQTL: num [1:200] 0.041 -0.31 -0.605 -1.221 1.121 ...
  ..$ PCC_De_Jager_eQTL  : num [1:200] -0.962 1.283 -0.819 0.82 -0.85 ...
  ..$ AC_De_Jager_eQTL   : num [1:200] 0.077 0.8421 -1.835 -0.0184 -0.4436 ...
  ..$ Mic_Kellis_eQTL    : num [1:200] 0.841 0.17 0.514 0.368 1.589 ...
  ..$ Ast_Kellis_eQTL    : num [1:200] 2.17e-01 2.37e-05 1.53 1.81 -1.88 ...
  ..$ Oli_Kellis_eQTL    : num [1:200] 1.762 1.519 -0.418 -0.665 1.209 ...
  ..$ OPC_Kellis_eQTL    : num [1:200] -0.619 -0.4 -0.3 1.013 -0.322 ...
  ..$ Exc_Kellis_eQTL    : num [1:200] 1.091 0.092 0.333 -0.208 0.433 ...
  ..$ Inh_Kellis_eQTL    : num [1:200] -0.932 -1.002 0.659 1.344 -0.353 ...
  ..$ Ast.10_Kellis_eQTL : num [1:200] 0 0 0 0 0 0 0 0 0 0 ...
 $ XtX     : num [1:16, 1:16] 2.182 1.09 1.481 0.759 2.405 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:16] "Mic_De_Jager_eQTL" "Ast_De_Jager_eQTL" "Oli_De_Jager_eQTL" "OPC_De_Jager_eQTL" ...
  .. ..$ : chr [1:16] "Mic_De_Jager_eQTL" "Ast_De_Jager_eQTL" "Oli_De_Jager_eQTL" "OPC_De_Jager_eQTL" ...

  • <cwd>/<name>_cache/<name>.<gene_id>.qss.rds (random_null_tensorqtl_1) (per-gene intermediate written by the tensorQTL entry point before the effects are pooled; not produced by susie_to_mash)

Minimal Working Example#

This pipeline has two entry-point workflows. Run whichever matches your input; both write the strong / null / random z-score matrices for downstream multivariate analysis. Both write {name}.mash_input.rds, so run one or the other rather than both in sequence.

Build the MASH input from fine-mapped SuSiE results#

The susie_to_mash workflow requires --name, --fine_mapping_meta (the index of fine-mapped RDS files) and --independent_variant_list; credible-set coverage and per-chunk size are tunable.

Timing: ~2 min (on toy dataset)

sos run pipeline/mash_preprocessing.ipynb susie_to_mash \
    --name protocol_example_mash \
    --fine_mapping_meta input/finemapping/protocol_example.fine_mapping_meta.tsv \
    --finemapping_column susie_path \
    --sig_p_cutoff 0.1 \
    --cwd output/mash_preprocessing

Extract random and null effects from tensorQTL summary statistics#

The random_null_tensorqtl workflow requires --name, --region_file, --independent_variant_list, and the list of summary-statistics files and traits.

Timing: TBD (on toy dataset)

sos run pipeline/mash_preprocessing.ipynb random_null_tensorqtl \
    --name protocol_example_mash \
    --region_file input/finemapping/protocol_example.region \
    --sum_files input/protocol_example.sumstats_list.txt \
    --traits bulk_rnaseq \
    --cwd output/mash_preprocessing

Command Interface#

sos run pipeline/mash_preprocessing.ipynb -h

usage: sos run pipeline/mash_preprocessing.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  susie_to_mash
  random_null_tensorqtl

Global Workflow Options:
  --name VAL (as str, required)
  --cwd output (as path)
                        Path to work directory where output locates
  --modular-script-dir code/script (as path)
  --seed 999 (as int)
  --n-random 10 (as int)
                        Random / null background rows sampled per region
  --n-null 10 (as int)
  --exclude-condition  (as list)
                        Conditions (columns) to drop from every partition
  --coverage 0.95 (as float)
                        Credible-set coverage for FineMappingResult strong
                        selection (susie_to_mash)
  --sig-p-cutoff 1e-06 (as float)
                        Strong-partition significance cutoff. Production default
                        1E-6; loosen (e.g. 0.1) for weak-signal toy data or
                        every strong variant is filtered out.
  --[no-]z-only (default to False)
                        Emit only the z-score matrices (drop the effect-size .b
                        / .s matrices)
  --container ''
                        Containers that contains the necessary packages
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 1h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 1 (as int)
                        Number of threads
  --independent-variant-list . (as path)
                        LD-pruned independent SNP list (optional). When
                        provided, restricts the random/null background to
                        matching variants (allele-aware via matchVariants);
                        strong is unaffected. Default path(".") means not
                        provided.

Sections
  susie_to_mash:        Build the MASH input from per-region FineMappingResult
                        RDS (e.g. fine_mapping.R output). Strong = the lead
                        variant (max PIP) of each credible set per condition;
                        random / null are sampled from each region's full
                        variant universe. mashInput merges across regions,
                        cleans via filterInvalidSummaryStat, and appends the
                        strong XtX cross-product.
    Workflow Options:
      --fine-mapping-meta VAL (as path, required)
                        Index of per-region fine-mapping results (one row per
                        region)
      --finemapping-column 'susie_path'
                        Column in --fine_mapping_meta holding each region's
                        FineMappingResult RDS path
      --region-column 'region_id'
                        Column holding the region id (labels rows across
                        regions)
  random_null_tensorqtl_1:
    Workflow Options:
      --sum-files  paths

      --region-file VAL (as path, required)
      --traits  paths

  random_null_tensorqtl_2:

Workflow implementation#

The pipeline defines two workflows. susie_to_mash (steps susie_to_mash_1 and susie_to_mash_2) builds the MASH input from SuSiE fine-mapping outputs. random_null_tensorqtl extracts the random and null effects directly from tensorQTL summary statistics. The shared [global] block below declares the common parameters.

[global]
parameter: name = str
# Path to work directory where output locates
parameter: cwd = path("./output")
parameter: modular_script_dir = path('code/script')
parameter: seed = 999
# Random / null background rows sampled per region
parameter: n_random = 10
parameter: n_null = 10
# Conditions (columns) to drop from every partition
parameter: exclude_condition = []
# Credible-set coverage for FineMappingResult strong selection (susie_to_mash)
parameter: coverage = 0.95
# Strong-partition significance cutoff. Production default 1E-6; loosen (e.g.
# 0.1) for weak-signal toy data or every strong variant is filtered out.
parameter: sig_p_cutoff = 1E-6
# Emit only the z-score matrices (drop the effect-size .b / .s matrices)
parameter: z_only = False
# Containers that contains the necessary packages
parameter: container = ""
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "1h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 1
# LD-pruned independent SNP list (optional). When provided, restricts the
# random/null background to matching variants (allele-aware via matchVariants);
# strong is unaffected. Default path(".") means not provided.
parameter: independent_variant_list = path(".")

# Build the MASH input from per-region FineMappingResult RDS (e.g. fine_mapping.R
# output). Strong = the lead variant (max PIP) of each credible set per
# condition; random / null are sampled from each region's full variant universe.
# mashInput merges across regions, cleans via filterInvalidSummaryStat, and
# appends the strong XtX cross-product.
[susie_to_mash]
# Index of per-region fine-mapping results (one row per region)
parameter: fine_mapping_meta = path
# Column in --fine_mapping_meta holding each region's FineMappingResult RDS path
parameter: finemapping_column = 'susie_path'
# Column holding the region id (labels rows across regions)
parameter: region_column = 'region_id'
import csv
with open(fine_mapping_meta) as _f:
    _rows = list(csv.DictReader(_f, delimiter = '\t'))
_fmr_paths = [r[finemapping_column] for r in _rows]
_region_ids = ",".join(str(r[region_column]) for r in _rows)
input: _fmr_paths
output: f"{cwd}/{name}.mash_input.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_preprocessing.R \
        --objects ${_input} \
        --region-ids ${_region_ids} \
        --coverage ${coverage} \
        --n-random ${n_random} \
        --n-null ${n_null} \
        --exclude-condition "${','.join(exclude_condition)}" \
        --sig-p-cutoff ${sig_p_cutoff} \
        ${'--z-only' if z_only else ''} \
        ${('--independent-variant-list ' + str(independent_variant_list)) if independent_variant_list.is_file() else ''} \
        --seed ${seed} \
        --output ${_output}

Get the random and null effects per analysis unit#

[random_null_tensorqtl_1]
parameter: sum_files = paths
parameter: region_file = path
parameter: traits = paths
import re
meta = []
with open(region_file) as regions:
    regions.readline()  # header
    for line in regions:
        parts = line.strip().split("\t")
        chr_id, start, end, gene_id = parts[0], parts[1], parts[2], parts[3]
        pattern = re.compile(r"\.%s\." % chr_id[3:])
        matched = []
        for sf in sum_files:
            with open(str(sf)) as af:
                for aline in af:
                    if pattern.search(aline):
                        matched.append(aline.strip())
        meta.append(dict(gene_id = gene_id, paths = ",".join(matched),
                         region = f"{chr_id}:{start}-{end}"))
input: for_each = 'meta'
output: f'{cwd:a}/{name}_cache/{name}.{_meta["gene_id"]}.qss.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_sumstats_construct.R \
        --tensorqtl-paths ${_meta['paths'].replace(',', ' ')} \
        --conditions ${','.join([str(t) for t in traits])} \
        --region ${_meta['region']} \
        --study ${name} \
        --output ${_output}

[random_null_tensorqtl_2]
input: group_by = 'all'
output: f"{cwd}/{name}.mash_input.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_preprocessing.R \
        --objects ${_input} \
        --n-random ${n_random} \
        --n-null ${n_null} \
        --exclude-condition "${','.join(exclude_condition)}" \
        --sig-p-cutoff ${sig_p_cutoff} \
        ${'--z-only' if z_only else ''} \
        ${('--independent-variant-list ' + str(independent_variant_list)) if independent_variant_list.is_file() else ''} \
        --seed ${seed} \
        --output ${_output}