LDSC fails after harmonization

[ampregnall]

Hello! Thank you again for creating and maintaining this package. I am running into a potential bug. If I use the harmonize() function and then try to perform LDSC to account for test statistic inflation I get an error. My code is as follows:

sumstats.harmonize(
    ref_seq=args.fasta,
    ref_rsid_vcf=args.dbsnp, 
    ref_infer=args.popvcf, # Ancestry specific. Logic handled by Nextflow
    ref_alt_freq="AF", 
    threads=args.threads, # pass threads from nextflow process, 
    sweep_mode=True 
    )

# Perform LDSC correction
sumstats_hapmap3 = sumstats.filter_hapmap3(inplace=False)
sumstats_hapmap3.estimate_h2_by_ldsc(ref_ld = args.ldsc,  w_ld = args.ldsc)


if np.float64(sumstats_hapmap3.ldsc_h2['Intercept'][0]) > 1:
     # Perform correction

The error traceback is:

Traceback (most recent call last): 
File "ampregnall/tools/nf-meta-gwas/bin/munge_sumstats.py", line 43, 
in <module> sumstats_hapmap3.estimate_h2_by_ldsc(ref_ld = args.ldsc, w_ld = args.ldsc) 
File "/opt/conda/lib/python3.12/site-packages/gwaslab/g_Sumstats.py", line 1544, 
in estimate_h2_by_ldsc self.ldsc_h2, self.ldsc_h2_results = _estimate_h2_by_ldsc(insumstats=insumstats,
    File "/opt/conda/lib/python3.12/site-packages/gwaslab/qc/qc_decorator.py", line 241, in wrapper result = func(*args, **kwargs)
    File "/opt/conda/lib/python3.12/site-packages/gwaslab/util/util_ex_ldsc.py", line 404, 
    in _estimate_h2_by_ldsc summary = estimate_h2(sumstats, args = default_kwargs, log = log) 
    File "/opt/conda/lib/python3.12/site-packages/gwaslab/extension/ldsc/ldsc_sumstats.py", line 331, 
    in estimate_h2 M_annot, w_ld_cname, ref_ld_cnames, sumstats, novar_cols = _read_ld_sumstats( 
    File "/opt/conda/lib/python3.12/site-packages/gwaslab/extension/ldsc/ldsc_sumstats.py", line 254, 
    in _read_ld_sumstats sumstats = _merge_and_log(ref_ld, sumstats, 'reference panel LD', log) 
    File "/opt/conda/lib/python3.12/site-packages/gwaslab/extension/ldsc/ldsc_sumstats.py", line 239, 
    in _merge_and_log raise ValueError(msg.format(N=len(sumstats), F=noun)) 
ValueError: -After merging with reference panel LD, 0 SNPs remain.

However, if I save the summary statistics from the harmonize step and call a separate script that loads that file and performs LDSC everything works. From my pipeline design perspective I think separating the processes is actually better, but something weird seems to be happening. I have attached my log file. Thanks!

gwaslab-logs.txt

[Cloufield]

Hi,
Thanks a lot for the feedback!
The merge failed because dropna() was run on the whole table before ldsc.
There might be some extra columns being all NA, so every row was dropped before the LD merge, which then showed 0 SNPs.
When you called a separate script, maybe you also dropped those empty columns, and the dropna did not drop every row.
I just fixed it in v4.1.5.

[ampregnall]

Thank you! I appreciate your responsiveness in identifying the issue. Do you have a set of parameters that you would recommend for the harmonize() function or for filtering post-harmonization to setup a within and/or across-population IVW meta-analysis. I have cobbled together some manual filtering but am curious what you think the best way to use the tool is. Happy to close the issue with the update.

sumstats.basic_check(
    remove=True,
    remove_dup=True,
    normalize=True,
    remove_dup_kwargs={"mode": "md", "keep": "first", "keep_col": "P"},
)

# Harmonize summary statistics
sumstats.harmonize(
    basic_check=False,
    ref_seq=args.fasta,
    ref_rsid_vcf=args.dbsnp,
    ref_infer=args.popvcf,  # Ancestry specific. Logic handled by Nextflow
    ref_alt_freq="AF",
    threads=args.threads,  # pass threads from nextflow process,
    sweep_mode=True,
)

status_str = sumstats.data["STATUS"].astype(str)
mask = ~(status_str.str[6].isin(["7", "8"]) | status_str.str[5].isin(["6"]))
sumstats.data = sumstats.data[mask].copy()

if args.type == "binary":
    sumstats.data["MAC"] = sumstats.data["MAF"] * sumstats.data["N_CASE"] * 2
else:
    sumstats.data["MAC"] = sumstats.data["MAF"] * sumstats.data["N"] * 2

[Cloufield]

Hi,
Thanks for the feedback.
I think QC creteria vary a lot, and depend on what you plan to do after the harmonization.
My suggestion is to run sumstats.lookup_status() to check the status and determine what variants you want to keep and exclude based on your datasets. (you already did that by filtering out variants with certain status).
Regarding cross-population meta-analysis, it usually requires more strict filters than within population meta-analysis. For example, sometimes you may need additional filtering like MAF>0.01.