biona001 / ghostknockoffgwas

Knockoff-based analysis of GWAS summary statistics data

GhostKnockoffGWAS

Documentation	Build Status	Code Coverage

This is a package for analyzing summary statistics data from genome-wide association studies (GWAS) under the statistical knockoff framework. Compared to marginal association testing which controls the FWER, the knockoff framework conducts conditional independence testing while controlling the FDR. As a consequence, GhostKnockoffGWAS can be both more precise and powerful than current state-of-the-art GWAS+fine-mapping methods. Its detailed evaluations can be found in our companion paper.

New users

To get started, please refer to the documentation.

In GhostKnockoffGWAS, the main working assumption is that we do not have access to individual level genotype or phenotype data. Rather, for each SNP, we have its Z-scores with respect to some phenotype from a GWAS, and access to LD (linkage disequilibrium) data. The user is expected supply the Z-scores, while we supply pre-processed LD files freely downloadable from the cloud.

Advantages/disadvantages of GhostKnockoffGWAS

Compared to existing knockoff methods for GWAS, the main advantages of GhostKnockoffGWAS is (1) its ease of use and (2) its computational efficiency. The only user-provided input is marginal Z-scores. Computationally, running a knockoff-based GWAS pipeline took approximately 15 minutes on 650,000 SNPs. The main limitation of GhostKnockoffGWAS is that it relies on the availability of pre-processed LD files suitable for the user's target samples.

Bug fixes and user support

If you encounter a bug or need user support, please open a new issue on Github. Please provide as much detail as possible for bug reports, ideally a sequence of reproducible code that lead to the error.

PRs and feature requests are welcomed!

Citation

If you use GhostKnockoffGWAS in your research, please cite the following references:

He Z, Chu BB, Yang J, Gu J, Chen Z, Liu L, Morrison T, Bellow M, Qi X, Hejazi N, Mathur M, Le Guen Y, Tang H, Hastie T, Ionita-laza, I, Sabatti C, Candes C. "In silico identification of putative causal genetic variants", bioRxiv, 2024.02.28.582621; doi: https://doi.org/10.1101/2024.02.28.582621.

ghostknockoffgwas's Issues

Improve documentation

• Readme should describe GhostKnockoffGWAS in relation to mainstream GWAS methods, as most users won't be familiar with knockoff based inference.
• The developer documentation is a little obsolete and should be revamped. In particular, the syntax for calling BlockGroupGhostMatrix (which depended on R package ghostbasil) should be updated to the syntax used by Ghostbasil.jl.
• Gallery: We should include a few manhattan plots comparing marginal vs knockoff analysis

Error in CMplot

As reported by @GuJQ5, an error appears when drawing the manhattan plot:

Error in CMplot(x1t, plot.type = "m", LOG10 = FALSE, col = c("grey30",  :

  unused argument (bin.range = c(0, 500))

The reason is that the CMplot package has been update on Jan 19, 2024, where the argument bin.range is removed in the latest version.

FileIO error

As reported by @GuJQ5, there is some problem with reading .h5 formatted LD files.

$ app_linux_x86/bin/GhostKnockoffGWAS --zfile UCB_GWAS.txt.gz --LD-files EUR --N 44161 --genome-build 19 --out UCB_GWAS_results

┌ Warning: Error requiring `UnPack` from `JLD2`
│ exception = 
│ IOError: stat("/home/groups/sabatti/.julia/packages/JLD2/ryhNR/ext/UnPackExt.jl"): permission denied (EACCES) 
│ Stacktrace: 
│ [1] uv_error 
│ @ ./libuv.jl:100 [inlined] 
│ [2] stat(path::String) 
│ @ Base.Filesystem ./stat.jl:152 
│ [3] isfile(path::String) 
│ @ Base.Filesystem ./stat.jl:461 
│ [4] macro expansion 
│ @ /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/Requires.jl:37 [inlined] 
│ [5] top-level scope 
│ @ /home/groups/sabatti/.julia/packages/JLD2/ryhNR/src/JLD2.jl:607 
│ [6] eval 
│ @ ./boot.jl:370 [inlined] 
│ [7] eval 
│ @ /home/groups/sabatti/.julia/packages/JLD2/ryhNR/src/JLD2.jl:1 [inlined] 
 [8] (::JLD2.var"#131#134")() 
│ @ JLD2 /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:101 
│ [9] macro expansion 
│ @ ./timing.jl:393 [inlined] 
│ [10] err(f::Any, listener::Module, modname::String, file::String, line::Any) 
│ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:47 
│ [11] (::JLD2.var"#130#133")() 
│ @ JLD2 /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:100 
│ [12] withpath(f::Any, path::String) 
│ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:37 
│ [13] (::JLD2.var"#129#132")() 
│ @ JLD2 /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:99 
│ [14] listenpkg(f::Any, pkg::Base.PkgId) 
│ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:20 
│ [15] macro expansion 
│ @ /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:98 [inlined] 
│ [16] __init__() 
│ @ JLD2 /home/groups/sabatti/.julia/packages/JLD2/ryhNR/src/JLD2.jl:606 
└ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:51 

┌ Warning: Error requiring `FileIO` from `HDF5` 
│ exception = 
│ IOError: stat("/home/groups/sabatti/.julia/packages/HDF5/HtnQZ/src/fileio.jl"): permission denied (EACCES) 
│ Stacktrace: 
│ [1] uv_error 
│ @ ./libuv.jl:100 [inlined] 
│ [2] stat(path::String) 
│ @ Base.Filesystem ./stat.jl:152 
│ [3] isfile(path::String) 
│ @ Base.Filesystem ./stat.jl:461 
│ [4] top-level scope 
│ @ /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/Requires.jl:37 
│ [5] eval 
│ @ ./boot.jl:370 [inlined] 
│ [6] eval 
│ @ /home/groups/sabatti/.julia/packages/HDF5/HtnQZ/src/HDF5.jl:1 [inlined] 
│ [7] (::HDF5.var"#115#121")() 
│ @ HDF5 /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:101 
│ [8] macro expansion 
│ @ ./timing.jl:393 [inlined] 
│ [9] err(f::Any, listener::Module, modname::String, file::String, line::Any) 
│ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:47 
│ [10] (::HDF5.var"#114#120")() 
│ @ HDF5 /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:100 
│ [11] withpath(f::Any, path::String) 
│ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:37 
│ [12] (::HDF5.var"#113#119")() 
│ @ HDF5 /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:99 
│ [13] listenpkg 
│ @ /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:20 [inlined] 
│ [14] macro expansion 
│ @ /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:98 [inlined] 
│ [15] __init__() 
│ @ HDF5 /home/groups/sabatti/.julia/packages/HDF5/HtnQZ/src/HDF5.jl:119 
└ @ Requires /home/groups/sabatti/.julia/packages/Requires/Z8rfN/src/require.jl:51 

Welcome to GhostKnockoffGWAS analysis! 
You have specified the following options: 
zfile = /oak/stanford/groups/zihuai/UCB data/Mar_13_James/UCB_GWAS.txt.gz 
LD_files = /oak/stanford/groups/zihuai/UCB data/Mar_13_James/EUR 
N (sample size) = 44161 
hg_build = 19 
outdir = /oak/stanford/groups/zihuai/UCB data/Mar_13_James/ 
outfile = /oak/stanford/groups/zihuai/UCB data/Mar_13_James/UCB_GWAS_results seed = 2023 
verbose = true 
random_shuffle = true 
skip_shrinkage_check = false 

count_matchable_snps processed chr 1, cumulative SNPs = 46145 
count_matchable_snps processed chr 2, cumulative SNPs = 94450 
count_matchable_snps processed chr 3, cumulative SNPs = 135010 
count_matchable_snps processed chr 4, cumulative SNPs = 173294 
count_matchable_snps processed chr 5, cumulative SNPs = 209902 
count_matchable_snps processed chr 6, cumulative SNPs = 253171 
count_matchable_snps processed chr 7, cumulative SNPs = 285925 
count_matchable_snps processed chr 8, cumulative SNPs = 316982 
count_matchable_snps processed chr 9, cumulative SNPs = 343328 
count_matchable_snps processed chr 10, cumulative SNPs = 373313 
count_matchable_snps processed chr 11, cumulative SNPs = 403207 
count_matchable_snps processed chr 12, cumulative SNPs = 431788 
count_matchable_snps processed chr 13, cumulative SNPs = 452614 
count_matchable_snps processed chr 14, cumulative SNPs = 471730 
count_matchable_snps processed chr 15, cumulative SNPs = 489924 
count_matchable_snps processed chr 16, cumulative SNPs = 510232 
count_matchable_snps processed chr 17, cumulative SNPs = 528617 
count_matchable_snps processed chr 18, cumulative SNPs = 546334 
count_matchable_snps processed chr 19, cumulative SNPs = 561544 
count_matchable_snps processed chr 20, cumulative SNPs = 576703 
count_matchable_snps processed chr 21, cumulative SNPs = 585317 
count_matchable_snps processed chr 22, cumulative SNPs = 594385 
Error encountered while load FileIO.File{FileIO.DataFormat{:HDF5}, String}("EUR/chr1/LD_start100826405_end102041015.h5"). 

Fatal error: 
ERROR: HDF5 load error: neither load nor fileio_load is defined 
due to FileIO.SpecError(HDF5, :load) 
Will try next loader. 

Stacktrace: 
[1] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::FileIO.Formatted; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:209 
[2] action 
@ /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:196 [inlined] 
[3] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::Symbol, ::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185 
[4] action 
@ /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185 [inlined] 
[5] load(::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:113 
[6] load 
@ /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:109 [inlined] 
[7] macro expansion 
@ /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/ghostbasil_parallel.jl:131 [inlined] 
[8] macro expansion 
@ ./timing.jl:393 [inlined] 
[9] ghostknockoffgwas(LD_files::String, z::Vector{Float64}, chr::Vector{Int64}, pos::Vector{Int64}, effect_allele::Vector{String}, non_effect_allele::Vector{String}, N::Int64, hg_build::Int64, outdir::String; outname::String, seed::Int64, target_chrs::Vector{Int64}, A_scaling_factor::Float64, kappa_lasso::Float64, LD_shrinkage::Bool, target_fdrs::Vector{Float64}, verbose::Bool, skip_shrinkage_check::Bool, random_shuffle::Bool) 
@ GhostKnockoffGWAS /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/ghostbasil_parallel.jl:130 
[10] macro expansion 
@ /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/app.jl:43 [inlined] 
[11] macro expansion 
@ ./timing.jl:393 [inlined] 
[12] julia_main() 
@ GhostKnockoffGWAS /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/app.jl:42 
[13] top-level scope 
@ none:1

Stacktrace: 
[1] handle_error(e::FileIO.LoaderError, q::Base.PkgId, bt::Vector{Union{Ptr{Nothing}, Base.InterpreterIP}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/error_handling.jl:61 
[2] handle_exceptions(exceptions::Vector{Tuple{Any, Union{Base.PkgId, Module}, Vector}}, action::String) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/error_handling.jl:56
[3] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::FileIO.Formatted; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:228 
[4] action 
@ /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:196 [inlined] 
[5] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::Symbol, ::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185 
[6] action 
@ /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185 [inlined] 
[7] load(::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) 
@ FileIO /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:113 
[8] load 
@ /home/groups/sabatti/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:109 [inlined] 
[9] macro expansion 
@ /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/ghostbasil_parallel.jl:131 [inlined] 
[10] macro expansion 
@ ./timing.jl:393 [inlined] 
[11] ghostknockoffgwas(LD_files::String, z::Vector{Float64}, chr::Vector{Int64}, pos::Vector{Int64}, effect_allele::Vector{String}, non_effect_allele::Vector{String}, N::Int64, hg_build::Int64, outdir::String; outname::String, seed::Int64, target_chrs::Vector{Int64}, A_scaling_factor::Float64, kappa_lasso::Float64, LD_shrinkage::Bool, target_fdrs::Vector{Float64}, verbose::Bool, skip_shrinkage_check::Bool, random_shuffle::Bool) 
@ GhostKnockoffGWAS /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/ghostbasil_parallel.jl:130 
[12] macro expansion 
@ /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/app.jl:43 [inlined] 
[13] macro expansion 
@ ./timing.jl:393 [inlined] 
[14] julia_main() 
@ GhostKnockoffGWAS /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/app.jl:42 
[15] top-level scope 
@ none:1

manhattan plots not utilizing within-group knockoff filter

very slow manhattan plots

It seems like when there are a lot of discoveries (e.g. > 100), it takes a really long time (>20h) to generate the resulting manhattan plots

Explicitly check `CHROM` field is Int

According to the docs

--chr input: This MUST be an integer and it must match the CHROM field in your VCF file. Thus, if your VCF file has CHROM field like chr1, CHR1, or CHROM1 etc, each record must be renamed into 1

However currently there's no check of this right now.

Using your method with our data

Hello Benjamin!

We are interested in using your method in our study, but we have found that only 20% of the SNPs overlap between our data and the reference information you provided. We have the raw genotype data available.

1. Is it possible for us to generate the corresponding LD files based on our own data?
2. Can we specify the groups ourselves for our analysis, or do we need to use the groups you have provided?

We want to ensure that we can properly apply your method to our dataset. Any guidance you can provide would be greatly appreciated.

Thank you for your time!

Improve coverage

We need more unit tests. We can think about doing the following:

1. Upload a small and freely distributable VCF file into this repo
2. Use solve_blocks() function to create .h5 formatted LD files and test its properties
3. Download GhostKnockoffGWAS from zenodo
4. Call GhostKnockoffGWAS on .h5 files from step 2 and test its properties

VCF: run `solve_blocks()` on typed SNPs only

Some VCF files might have imputed SNPs.

In ~/knockoffs/400phenotypes-part1.ipynb on Sherlock, I have some code to compute group representatives only on the typed variants. Putting this here in case I forget.

Tutorial for `solveblock`

Currently, there exist documentation for the solveblock function as well as some high-level overview for what it is supposed to do. However, we should probably write a step-by-step tutorial.

Custom Reference Panel in PLINK Format

Hello Benjamin,

I am attempting to use the GhostKnockoffGWAS tool with my own reference panel for chr17 in PLINK format (bed, bim, fam files). I have converted my PLINK files to VCF format and attempted to use the solve_block function (https://biona001.github.io/GhostKnockoffGWAS/dev/man/solveblocks/), but I am unable to locate the solve_block executable in the bin directory and receive an error.

Could you provide guidance on:
How to properly integrate a custom reference panel in PLINK format?
The correct way to locate or generate the solve_block executable?

Thank you for your time and for developing this valuable tool.

wget `example_zfile.txt` broken

2 things to fix in the quick start of documentation

1. In step (1), the download link to the software is fixed at v0.1.0. It should be the latest version
2. The 3rd step actually doesn't work.

Roadmap until 1.0 release

Features to add:

1. Add knockoff q-values to output (see eq19 of this paper)
2. Ability to read different columns in zfile, based on user input
3. Add code to make manhattan plot and show this in documentation
4. Within-group knockoff filter (see this paper)

Docs for `solve_blocks()`

Parse LD files with either hg19 or hg38, but not both

The EUR LD files contain both pos_hg19 and pos_hg38.

But in the vcf branch we will soon support custom LD files constructed from individual level VCF data, which will only have either pos_hg19 or pos_hg38.

Set `random-shuffle=false` by default

Consider turning this option off due to the following:

• In our simulations, we never used this option, and FDR was always under control.
• In our flagship paper, the method section did not mention randomly permutate Z-scores and its knockoffs.

Non-linux support?

Rather than compiling binaries for GhostKnockoffGWAS via PackageCompiler.jl, consider providing a docker image, see this example. Dockers provide an alternative way to run GhostKnockoffGWAS, even on non-linux machines without Julia installed

`manhattan.R` does not work on output of `solveblock`

I get an error

Error in x1_sug_t[1, "BP"] - ra_w : 
  non-numeric argument to binary operator

`solveblock`: matrix contains Infs or NaNs

When there are only 1 SNP in the user specified region, cor(X) is NaN. Must throw error here. Along the same lines, probably we should add a check regarding minimum sample size.

Roadmap for v0.2.0 release

The following should work before releasing v0.2.0

• Implement solve_block() function and compile it into solveblock binary executable (#13)
• Add documentation & tutorial for solveblock (#15)
• Update ghostknockoffgwas so it runs on custom built LD files (#11)
• Test (on EUR LD files) we get the same result as v0.1.2
• Test (on custom built LD files) the result is reasonable

`solveblock` throws crazy error when no SNPs is between `start_bp` and `end_bp`

MWE (the first SNP is chr1 POS 58814)

using GhostKnockoffGWAS
vcffile = "/oak/stanford/groups/zihuai/paisa/VCF/chr1.vcf.gz"
chr = 1
start_bp = 1 
end_bp = 10583 
outdir = "/oak/stanford/groups/zihuai/paisa/LD_files"
hg_build = 19
julia> solve_blocks(vcffile, chr, start_bp, end_bp, outdir, hg_build)

 ** On entry to DGETRI parameter number  6 had an illegal value
ArgumentError: invalid argument #6 to LAPACK call

Stacktrace:
  [1] chklapackerror
    @ /share/software/user/open/julia/1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/lapack.jl:38 [inlined]
  [2] getri!(A::Matrix{Float64}, ipiv::Vector{Int64})
    @ LinearAlgebra.LAPACK /share/software/user/open/julia/1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/lapack.jl:1045
  [3] inv!
    @ /share/software/user/open/julia/1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:508 [inlined]
  [4] _inv(A::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ LinearAlgebra /share/software/user/open/julia/1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/symmetric.jl:632
  [5] inv
    @ /share/software/user/open/julia/1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/symmetric.jl:648 [inlined]
  [6] choose_group_reps(Σ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}}, groups::Vector{Int64}; threshold::Float64, prioritize_idx::Nothing, Σinv::Nothing)
    @ Knockoffs /home/groups/sabatti/.julia/packages/Knockoffs/AqvY6/src/group.jl:1974
  [7] choose_group_reps
    @ /home/groups/sabatti/.julia/packages/Knockoffs/AqvY6/src/group.jl:1929 [inlined]
  [8] macro expansion
    @ /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/make_hdf5.jl:269 [inlined]
  [9] macro expansion
    @ ./timing.jl:393 [inlined]
 [10] solve_blocks(vcffile::String, chr::Int64, start_bp::Int64, end_bp::Int64, outdir::String, hg_build::Int64; snps_to_keep::Nothing, tol::Float64, min_maf::Float64, min_hwe::Float64, force_block_diag::Bool, method::String, linkage::String, force_contiguous::Bool, group_cor_cutoff::Float64, group_rep_cutoff::Float64, verbose::Bool)
    @ GhostKnockoffGWAS /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/make_hdf5.jl:266
 [11] solve_blocks(vcffile::String, chr::Int64, start_bp::Int64, end_bp::Int64, outdir::String, hg_build::Int64)
    @ GhostKnockoffGWAS /home/groups/sabatti/.julia/dev/GhostKnockoffGWAS/src/make_hdf5.jl:228
 [12] top-level scope
    @ In[94]:20

ERROR: id is missing

On the data here: http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz

I get an error with solveblock: id is missing. The issue is caused by missing SNP name for one of the records. A simple fix is to return . for such SNPs.

Optimize VCF reading time

In the documentation, it reads

For optimal performance, it is best to filter the VCF file down to records between start_bp and end_bp (e.g. with bcftools) before running solveblock.

We can automate this step internally within solve_block() (if the input VCF file is indexed) by calling bcftools using the bcftools_jll.jl package.