zhilizheng / SBayesRC

GNU General Public License v3.0
25 stars 5 forks source link

sbayesrc: "Error in if (all(dt.tune$r < 0)) { : missing value where TRUE/FALSE needed" #26

Closed bwspitzer closed 5 months ago

bwspitzer commented 6 months ago

I've installed SBayesRC and downloaded all of the required files. The [minimal example] code(https://github.com/zhilizheng/SBayesRC#minimal-example) runs fine on the example data given here. However, when I run the code on the summary statistics of one chromosome of my real data (UKBB EUR protein data), SBayesRC::sbayesrc ends in an error.

I've tried this with different chromosomes and a different set of UKBB GWAS summary statistics, with the same result. Am I missing something?

Here are excerpts from the relevant files (with names simplified for readability). I'm attaching the full files to this post.

Thank you very much, Brian

head sumstats.txt

SNP A1 A2 freq b se p N rs555672843 A T 0.00136858 -0.0941215 0.0931064 0.312063 33191 rs561879353 G C 0.000820182 -0.0869738 0.12092 0.471974 33191 rs3916645 T A 0.977441 -0.0187979 0.0226051 0.405648 33191 rs28971674 A T 0.978258 -0.016392 0.0231259 0.478441 33191 rs28805621 T A 0.98246 -0.0266427 0.0273235 0.329518 33191 rs547911002 C G 0.0895738 0.0186111 0.0118198 0.115355 33191 rs28971399 T A 0.976691 -0.0141354 0.0225909 0.531503 33191 rs28778303 C G 0.971911 -0.00177705 0.0212677 0.933409 33191 rs551729884 A T 0.0858045 0.0171966 0.0118598 0.147063 33191

script:

ma_file=...sumstats.txt ld_folder=.../HapMap3_EUR/ukbEUR_HM3 annot=.../annot_baseline2.2.txt out_prefix=.../output threads=4 export OMP_NUM_THREADS=$threads

Rscript -e "SBayesRC::tidy(mafile='$ma_file', LDdir='$ld_folder', \ output='${out_prefix}_tidy.ma', log2file=TRUE)"

Rscript -e "SBayesRC::impute(mafile='${out_prefix}_tidy.ma', LDdir='$ld_folder', \ output='${out_prefix}_imp.ma', log2file=TRUE)"

Rscript -e "SBayesRC::sbayesrc(mafile='${out_prefix}_imp.ma', LDdir='$ld_folder', \ outPrefix='${out_prefix}_sbrc', annot='$annot', log2file=TRUE)"

Warning: In SBayesRC::tidy(mafile = "... sumstats.txt", : Too many SNPs (>30%) were missing in summary data after QC. The SBayesRC results may be unreliable.

output files:

output_imp.ma output_imp.ma.log output_sbrc.annot.tmp.bin output_sbrc.log output_sbrc_tune0.95.rds output_sbrc_tune0.95.rm.snpidx output_sbrc_tune0.995.rds output_sbrc_tune0.995.rm.snpidx output_sbrc_tune0.99.rds output_sbrc_tune0.99.rm.snpidx output_sbrc_tune0.9.rds output_sbrc_tune0.9.rm.snpidx output_sbrc_tune_inter.txt output_sbrc_tune.txt output_tidy.ma output_tidy.ma.full output_tidy.ma.log

output_sbrc.log:

SBayesRC v0.2.5 ... GWAS summary statistics: .../output_imp.ma Annotation: TRUE Annotation file: .../annot_baseline2.2.txt LD input folder: .../ukbEUR_Imputed, eigen variance cutoff threshold: 0.995 Tune the best eigen cutoff: TRUE Tuning step: 0.995 0.99 0.95 0.9 Tuning iterations: 150 Tuning burn-in iterations: 100 Start hsq: 0.5 Start Pi: 0.99 0.005 0.003 0.001 0.001 Gamma: 0 0.001 0.01 0.1 1 Number of MCMC iteration: 3000 Burn-in iteration: 1000 Use 2pq to scale summary: FALSE Method to resample residual: allMixVe Threshold to resample residual: 1.1

Perform tuning...

Tune with parameter 0.995 ... Tune with parameter 0.99 ... Tune with parameter 0.95 ... Tune with parameter 0.9 ... Tuning Done Error in if (all(dt.tune$r < 0)) { : missing value where TRUE/FALSE needed Calls: Execution halted

20240508_LPA_EUR_chr21_test.sh.txt 20240508_LPA_EUR_SBayesRC_chr21_sbrc.log

20240508_EUR_LPA_sumstats_for_SBayesRC_chr21.txt

zhilizheng commented 6 months ago

Hi @bwspitzer ,

Thanks for the report. I would recommend run genome wide. Running in 1 CHR is not supported currently and is not suggested. If you need a run, you can mannually set the imputed beta to 0, and se, N to a fixed value. Or, another way: keep LD in one CHR only.

Regards, Zhili