Hello,
I know that your package report standard errors of variance components and heritability in log or console of R studio.
So, I wonder how can I get or save standard errors of variance components and heritability as an output for example in a data frame or a table?

Regards,
Siavash

Running the following I'm getting based off of the newest github commit.

I also sometimes get
recursive gc invocation but can't quite get the code down to reproduce it yet...

Warning: stack imbalance in '=', 2 then -74

library("hibayes")

set.seed(42)
n = 10000
m = 1000
k = 10
h2 = 0.5

n = 1000
m = 1000
Z = matrix(sample(c(0, 1, 2), n*m, prob=c(0.25, 0.5, 0.25), replace=TRUE), n)
Z = apply(Z, 2, function(x){x-mean(x)})
k = 10
qtl = sort(sample(1:m, k))
betaTrue = array(0,m)
betaTrue[qtl] = rnorm(k)
h2 = 0.5

g = Z%*%betaTrue
vg = var(g)
ve = (1-h2)/h2 * vg
y = g + rnorm(n,0,sqrt(ve))

fit = bayes(y=y, Z=Z, model="R", niter=1000, nburn=500)

Hi there,

Thanks for making this available!
I am wondering if I can use LD matrix downloaded from (https://zenodo.org/record/3350914#.X94SHOn0lBx ) instead of computing by myself. For example, can I specify the path to this LD matrices and run SBayesR using hibayes package?

ldm1 = read("xx") # path to downloaded LD matrix
fit = sbayes(sumstat=sumstat, ldm=ldm1, model="SBayesR", niter=20000, nburn=10000, outfreq=10, verbose=TRUE)

Thanks!

Ying

Hi, it seems you are making a very useful and interesting R package here.

You may know my packages {bigstatsr} and {bigsnpr}.
I think your package could make use of those.

It would be also beneficial if we could agree on some common framework so that people don't see your package as just another tool, but rather as something that fits well in the existing ecosystem and comes as a complementary tool.

What do you think?

老师，您好
请问在GWAS分析时，结果可以像gemma一样输出位点的P值吗？

Hei, quick question, as I couldn't find immediate information in the manual etc: does hibayes enable parallelization? Or does it always run on one core only?
Thanks!

Hi,

I am trying to use a summary-level Bayes model (i.e. creating a PRS model using GWAS sumstat and 1000Genome reference panel with HAPMAP3 snplist). So far, I created an LD variance-covariance matrix with this command: ldm4 = ldmat(geno, map, ldchr=FALSE, chisq=5, threads=4); however, the final model seems to be unstable.

I tried to create an LD variance-covariance matrix with this command instead: ldm3 = ldmat(geno, map, ldchr=FALSE, threads=4). The problem is this step requires a very large memory even when I only limited the snplist to HAPMAP3 variants (jobs always got killed because of not enough memory). Would you please advise on how to solve this problem?

Thank you!

hibayes可以添加基因型填补功能么，vcf文件基因型缺失无法运行程序

Hello,
I applied the BSLMM algorithm to predict genomic breeding values, estimate variance components, and... in hibayes.
This is the script that I used:
fit <- bayes(y = pheno[, 8], M = geno, model ="BSLMM" , niter = 100000, nburn = 20000,threads=10)
After running this script, I faced with an error which is:
Error in make_grm(M, lambda = lambda, inverse = TRUE, verbose = verbose) :
system is computationally singular: reciprocal condition number = 4.452607e-26,
try to specify parameter 'lambda' with a small value, eg: 0.001 or bigger
I could solve it by adding lambda=0.01 as your software pointed out.
I wonder why when I used other algorithms such as BayesA, BayesCpi and ... I did not face this error?
Is there any difrences between GRM (genomic relationship matrices) construction and their inversion in different methods?
Also, I wonder how can I solve this problem without giving a lambda value in BSLMM?
Regards,
Siavash

I have experienced the error below when I run the example code using bayes results with hibayes package.

data <- cbind(map[, 1:3], (1 - fit[["pip"]]))
chr5 <- data[data[, 2] == 5, ]

visualize the results

CMplot(chr5, plot.type = "m", width = 9, height = 5, threshold = 0.01,

• ylab = expression(-log[10](1 - italic(PIP))), LOG10 = TRUE,
• amplify = FALSE)

Rectangular-Manhattan Plotting (1 - fit[["pip"]]).
Error in jpeg(paste("Rectangular-Manhattan.", taxa[i], ".jpg", sep = ""), :
unable to start jpeg() device
Ek olarak: Warning messages:
1: In jpeg(paste("Rectangular-Manhattan.", taxa[i], ".jpg", sep = ""), :
unable to open file 'Rectangular-Manhattan.(1 - fit[["pip"]]).jpg' for writing
2: In jpeg(paste("Rectangular-Manhattan.", taxa[i], ".jpg", sep = ""), :
opening device failed

Is it possible to output the MCMC estimates of gebvs every, e.g. 100th sample, of the MCMC chain? Sometimes it is useful to see the posterior distribution of each individual's gebv. I realise this would result in large files, and perhaps not everybody needs them, but it would be nice to have the option. Perhaps I missed it somewhere. This looks like a really nice package anyway - thank you!

您好，请问我在服务器安装hibayes包时，总是无法编译成功，报错如下：ERROR: compilation failed for package ‘hibayes’。请问应该如何解决呢？R语言版本3.6.1.

Hello!
The packages "hibayes" and "G2P" were performed to GS of weight trait, and the predictive accuracies were significantly different. The accuracies were 0.82 and 0.15 by BayesB in hibayes and G2P, respectively. I don't know the reasons.

Hi Li

I am getting this error: Error in as.vector(data) : no method for coercing this S4 class to a vector; when I try to run

pve <- apply(geno,2,var) * (fit$g^2) / var(pheno[,1]).

Can you please advise.

在进行GS分析时，能够加入固定效应吗？说明书没有具体的说明

Hello,
Thanks for your excellent package! I wondered what the default thinning interval parameter in the software is. And how can we change this parameter? Does the package include all the samples to infer the posterior distribution features?

Beware that if you're linking to the code from {Rcpp}, which you do, you have to use a GPL license.
Right, @eddelbuettel?

Hello,

When using variable selection method on markers, how do you calculate variance explained both genetic and phenotypic of all snps? I tried to follow the example in vignette but I am getting a lower sum PEV explained of all markers than snp heritability

Hello, When I have multiple observations of some genotypes (for example over several locations which I model as fixed effect terms), how do I format my genetic data? Do I need multiple rows of genetic data, one row per observation of a genotype, or not?