LDpred2-auto: automatic determination of use_MLE and shrink_corr

[sritchie73]

Hi,

I've just seen the updated tutorial on LDpred2-auto (last time I ran this, the use_MLE, allow_jump_sign, and shrink_corr arguments didn't exist), and I'm wondering if its possible to provide further guidance on (1) how and when to reduce shrink_corr from 0.95 to 0.4, and (2) what constitutes a "low power GWAS" that would cause convergence issues, ultimately with a goal of programmatically determining both parameters.

For use_MLE, is there a typical (effective) sample size below which you find you have to set use_MLE = FALSE?

For shrink_corr, I'm wondering if it makes sense to set this based on the proportion of GWAS summary statistics passing the QC filters when doing the SNP matching. Is there a cutoff below which this proportion passing QC indicates sufficient LD divergence that the shrink_corr should be set to 0.4 instead of 0.95? Or should the shrink_corr scale between 0.4 and 0.95 somehow proportionally to the % of GWAS summary statistics passing SNP QC?

[privefl]

That's very important questions.
I wanted to investigate automatic values and decisions for this based on some tests, but I never had time to do this (and will be probably never since I'm leaving academia in one month).

What I thought about setting use_MLE to TRUE and then rerunning with FALSE when

• too many chains diverged
• had predicted r2 at low value; predicted r2 is good proxy for GWAS power, but is not directly computed within LDpred2 (some code afterwards); alternatively, an upper bound for r2 is uniroot(function(r2) r2 - h2 / (1 + (1 - r2) * M_c / (N * h2)), interval = c(0, h2))$root, where M_c is the number of causal variants.

I think setting shrink_corr based on the proportions removed by QC can be a good idea.
What we did in some pipeline was to rerun with some lower value (0.05 or 0.1 lower each time down to 0.4) as long as not enough chains were kept.

[sritchie73]

Thanks @privefl ,

So the GWAS power relates to the power of the GWAS to predict the target trait in the data in which LDpred2 was run, rather than the power of the original GWAS?

When you say too many chains diverged / not enough chains kept, do you have a rule of thumb here on what proportion of chains should be kept?

[privefl]

For me, everything is basically the same, a monotonic transformation of N*h2/p.

I would aim to keep at half of them at least.

[sritchie73]

Thanks for the advice,

My plan is to iterate over shrink_corr = c(0.95, 0.9, 0.85, ..., 0.40) and use_MLE = c(TRUE, FALSE) at each shrink_corr, terminating in the first combination of shrink_corr and use_MLE that yield >50% chain convergence.

[privefl]

I'm testing a few things right now.
You are interested mainly in predictive performance, right?

[privefl]

I couldn't find something working for everything.
Basically shrink_corr of 0.8 or less seems good for external LD (maybe less when it is a bit multi-ancestry like with the CAD 2015 data), and 0.95 for internal LD.
use_MLE should be disabled for low power, and usually you get very wide range for the alpha estimates, or complete divergence.

[sritchie73]

Hi Florian,

Thanks for investigating. Yes - I'm primarily interested in predictive performance. We're in the progress of revising this preprint (https://www.medrxiv.org/content/10.1101/2024.08.22.24312440v3) to retrain the "metaPRS" (elastic-net of multiple LDpred2-derived PRSs across traits) in multiple distinct ancestries.

We've just finished a test run using the iteration approach I described in an earlier comment and have found some cases where the LDpred2-auto with use_MLE = FALSE appears to be getting stuck in some kind of cycle when estimating the proportion of causal variants:

This is with shrink_corr = 0.95, but the plots looks very similar for all shrink_corr tested down to 0.4.

With use_MLE = TRUE LD-pred2 very quickly aborts without converging:

Again, plots all look very similar to this regardless of the shrink_corr used.

This is for a EUR-ancestry GWAS being applied to AFR-ancestry samples.

Notably, the chains which do converge look identical to the chains which don't converge - does it make sense in these instances to relax the metric used to assess chain convergence to keep these chains? Or should we exclude LDpred2-auto in these instances?

Ultimately the PRS generated by LDpred2-auto is going to be compared to PRS generated by LDpred2-infinitesimal, LDpred2-lassosum2, and LDpred2-grid to find the most predictive PRS for the GWAS in some independent samples, then plugged into elasticnet with the optimal LDpred2 PRS for a suite of other traits; so non-predictive/nonsense PRS will ultimately be discarded.

[privefl]

This is quite odd.
Do you happen to have some code to download and process the sumstats?
Then, do you then use them with the LD I provide?

[sritchie73]

No, we compute the LD directly from the training dataset.

---

The GWAS summary statistics are on the GWAS Catalog with accession GCST90475667; http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST90475001-GCST90476000/GCST90475667/harmonised/GCST90475667.h.tsv.gz
has the file in the GWAS-catalog harmonized format on build GRCh38.

The following should hopefully get the summary statistics into a format usable by snp_match():

library(data.table)  # For fast read of table
library(R.utils) # enable read directly from gzipped file

# Download summary statistics
if (!file.exists("GCST90475667.h.tsv.gz")) {
  system("wget ftp://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST90475001-GCST90476000/GCST90475667/harmonised/GCST90475667.h.tsv.gz", wait=TRUE)
}

# Load summary statistics:
gwas <- fread("GCST90475667.h.tsv.gz")

# SE not reported, so need to compute from 95% CI for odds ratio
Zconst <- qnorm(p=0.05/2, lower.tail=FALSE)
gwas[, beta_se := (log(ci_upper) - log(odds_ratio))/Zconst] # Rearrangement of 95% CI = beta +/- 1.96 * se

# Arrange into format recognised by snp_match():
gwas <- gwas[,.(chr=chromosome, pos_b38=base_pair_location, EA=effect_allele, OA=other_allele, beta_se,
                             neg_log10_p=-log10(p_value), EAF=effect_allele_frequency)]

These have been filtered to the set of HapMap3+ variants with >1% MAF in the training dataset (a random 10K AFR participants from All of Us in this case).

LD has been computed from the training dataset directly; with the 10K sample size determined so that SNPs with MAF down to 1% are included after using the QC threshold you recommend in the LDpred2 tutorial when computing the correlation matrix (1/sqrt(sample size)).

After filtering SNPs and matching to the training dataset using snp_match, we also apply the recommended per SNP QC filters:

# N-eff is computed for the overall study as per-SNP case/control numbers (or sample sizes) are not provided with the 
# GWAS
n_eff <- 4 / (1/154194 + 1/278454) # 4/(1/cases + 1/controls)

# Check whether standard deviations of genotypes in the GWAS summary statistics 
# are consistent with the ones in the training dataset, filtering out SNPs which are too different
# Here, a1freq_train is computed for each SNP from the training data genotype dosage matrix
gwas[, sd_train := sqrt(2 * a1freq_train * (1 - a1freq_train))] 

# Infer standard deviation of the genotypes in the gwas
gwas[, sd_gwas := 2 / sqrt(n_eff * beta_se^2 + beta^2)] 

# Set up column with flag information about QC exclusion reason:
gwas[, snp_qc := ""]

# Remove SNPs where the estimated genotype dosage SD in the GWAS diverges too far from the SD of the 
# allele frequencies in the training dataset
gwas[snp_qc == "" & sd_gwas < 0.05, snp_qc := "Low MAF in GWAS (Estimated genotype dosage SD < 0.05)"]
gwas[snp_qc == "" & sd_train < 0.05, snp_qc := "Low MAF in training data (Genotype dosage SD < 0.05)"]
gwas[snp_qc == "" & (sd_gwas < (0.5 * sd_train) | sd_gwas > (sd_train + 0.1)), 
        snp_qc := "Estimated genotype dosage SD in GWAS too different from genotype dosage in training data"]

# Where the GWAS summary statistics include allele frequencies, also remove SNPs that are too divergent based
# directly on the observed allele frequencies
gwas[, sd_gwas_af := sqrt(2 * a1freq_gwas * (1 - a1freq_gwas))]
gwas[, a1freq_diff := a1freq_train - a1freq_gwas]
gwas[snp_qc == "" & sd_gwas_af < 0.05, snp_qc := "Low MAF in GWAS (Reported genotype dosage SD < 0.05)"]
gwas[snp_qc == "" & (sd_gwas < (0.7 * sd_gwas_af) | sd_gwas > (sd_gwas_af + 0.1)), 
        snp_qc := "Estimated genotype dosage SD in GWAS diverged from SD calculated from reported allele frequencies"]
gwas[snp_qc == "" & abs(a1freq_diff) > 0.07, 
        snp_qc := "GWAS allele frequency too divergent from training dataset allele frequency (abs diff > 0.07)"]

# Generate plots as desired, then exclude SNPs failing qc:
gwas <- gwas[snp_qc == ""]

For reference, here are the SNP-QC plots we generate on our end for this GWAS:

[sritchie73]

Also noting that INFO scores haven't been incorporated into the SNP QC as the training data are WGS not imputed genotypes.

[privefl]

The SD plot is awful.

One thing I don't understand; you have a GWAS from Europeans but are using LD from Africans, is that right?

[sritchie73]

Yes -

The intuition here (which may very well be wrong) is the even though the SD plot looks awful, the SNPs that are kept will be the ones that are most transferable between ancestries. The (simplified) idea then is to combine EUR GWAS -> AFR LD, AFR GWAS -> AFR LD, EUR GWAS -> AFR LD, EUR GWAS -> EUR LD (and other combinations of ancestry-specific training data and GWAS along with multi-ancestry GWAS), fit the four resulting PRS in ancestry-specific elasticnets in an independent set of EUR and AFR training data, then combine into a single PRS; i.e. to essentially up-weight the most transferrable SNPs when combined into a single PRS.

[privefl]

I see, but I don't think this would work as LD should be quite similar to the input GWAS data.

[sritchie73]

Thanks - we'll test what happens, skipping LDpred2 auto when <50% of the chains converge.