Change Block Size

[Genutty]

When running the julia implementation, we are getting large mean shrinkage parameters (>0.5). We would like to try reducing the block sizes. Would it be possible to alter this by changing "delta_strong_size = 500" to "delta_strong_size = 250" in the ghostbasil_parallel.jl? Or, is there another way to effectively do this?

[biona001]

@Genutty Thanks for asking the question. You cannot change delta_strong_size, which is a parameter used in Lasso regression. It does not have anything to do with block sizes and we did not test other values besides 500.

If you are getting large shrinkage values (>0.5), you basically only have 2 options:

1. Download a different LD file more suitable for your study
2. Create such LD files yourself, using the solveblock executable. This would require individual level data.

Before you do either of these steps, it might be helpful to first figure out why you are getting large shrinkage (>0.5 is quite a lot). Could you let me know the following?

• Which LD panel did you use?
• What is the ancestry of the original GWAS?
• How did you compute Z-scores?
• How many samples/SNPs are present in your study?
• Is your HG build 19 or 38?

[Genutty]

Hi @biona001. Thanks for your response. Our GWAS analysis is including only Hispanic/Latino individuals, for which there is no available appropriate LD reference panel. We did receive this error when using the EUR reference panel. Subsequently, we have tried to use solveblock to create our own LD files, which still produce large shrinkage parameters. However, this was run using the start and end positions from EUR ldetect cut-offs. We have additionally been trying to calculate appropriate LD minima start/end for our HL partcipants, but have encountered additional errors using ld_split.R (see below).

Which LD panel did you use? See above
What is the ancestry of the original GWAS? HL/AMR
How did you compute Z-scores? Beta/SE
How many samples/SNPs are present in your study? The GWAS summary statistics are based on 16M SNPs and up to 233,413 (cases + controls)
Is your HG build 19 or 38? build 38

Error messages from running ld_split.R to get start-bp and end-bp of AMR before running solveblock vary by chromosome, but include: 1) split is NULL, 2) NaNs in corr matrix, 3) daig of corr matrix is not 0. Below are the options we have tried.

Example 1:

Rscript --vanilla /path to script/ld_split.R ${chr} \\
/PathToBed/File_chr${chr}.bed \\
/PathToFBM/chr${chr}_hl_temp \\
/PathToOut/chr${chr}_hl_LDblocks.txt \\
0.01 0.3

We tested these additional thresholds for snp_ldsplit():

# thr_r2 0.05 to 0.2 by 0.05
# min_size = 10
# max_size = start with 100 and go up to 20,000
# max_K = keep 500
# max_r2 = 0.5

[biona001]

Hi @Genutty, thanks for following up on this! Firstly, note that the start/end positions of the EUR ldetect regions are based on HG19 not HG38, so you probably need to find other start/end positions.

To answer your questions on generating your own LD files,

1. If the split is NULL, it means that the specified parameters thr_r2, min_size, ...etc, cannot successfully split the data into independent regions.
2. I'm not super sure about NaNs in corr matrix, but my guess is you have SNPs that are the same for every sample (all 0s or all 2s).
3. I'm also not sure about the diag of corr matrix not being 1 -- which correlation matrix are you talking about?

For (1), Since you have 16M SNPs, the max_size parameter being 20000 probably is not large enough, it might need to go up to 100000 (which is obviously too large, so don't do it). 16M SNPs implies that you are using either WGS/WES or imputed SNPs. In our papers, we have not yet explored how our model performs in WGS/WES data. But if your data was imputed from array data, e.g. using a reference panel, then we recommend filtering out all imputed SNPs and run the pipeline only on the original SNPs from the array data. The reason is because the knockoff methodology tests conditional independence hypothesis, i.e. a SNPs is selected only if it offers unique information on the phenotype that no other SNPs provide. By definition, all imputed SNPs are functions of the typed SNPs, and therefore they are null by definition (they should never be selected by the knockoff method).

Thus, my recommendation would be (1) filter the 16M variants to only array SNPs (should give you <1M SNPs) and do basic quality controls (e.g. filter out singletons), (2) separate data by chromosomes, (3) run ld_split.R on each chromosome separately, (4) run solveblock on each chromosome separately, (4) run standard GWAS on all SNPs, and finally (5) run GhostKnockoffGWAS

[Genutty]

Hi @biona001, thanks for your response. We used the ldetect cut-offs in build 38 provided here: https://github.com/jmacdon/LDblocks_GRCh38?tab=readme-ov-file.

Regarding 1) Yes, we are using SNPs imputed to the TOPMED R2 reference panel. Additionally, our GWAS summary statistics are from a meta-analysis. So, restricting to only genotyped SNPs would be difficult given that different studies used different genotyping arrays. Even within our own study that we are trying to estimate LD, there different arrays used over time, and have not looked at using sequencing data yet. However, we could try restricting to HapMap3 variants (~3M) or a union set across the arrays that were used across studies.

Regarding 2 and 3, we did filter our variants based on MAF, but we can look closer to see if low freq variants are causing issue by increasing the thresholds.

[biona001]

Oh okay. I wasn't aware of these HG38 blocks from ldetect, they should be fine to use provided there isn't high shrinkage.

I think it is okay to leave the meta-analysis summary statistics file as it is, so maybe it is better to focus on generating the LD files that only includes typed SNPs. If you have been running solveblock, you must have individual-level data right? Does this individual level data also have 16M variants? If so, I'd recommend filtering out imputed SNPs from these, and then run the ld_split.R and solveblock on the filtered data. Once this is done, you'll get LD files that includes only typed SNPs. Of course , if you use these with the meta-analysis summary statistics, many SNPs will be missed in the analysis, but I think it's a reasonable compromise.

[Genutty]

Yes, we have the individual-level imputed data, so we can subset to those variants that were present across arrays. We will test this and see how the results change. Thank you for your help!