This function is a wrapper from the rrBLUP package to be used when a mixed model including markers to perform GWAS is specified and once the variance components have been estimated the fixed effects are obtained as B= (X'V-X)-X'V-y and the score calculation is obtained with the F statistic as F = Beta^2 / Var(Beta) where Var(Beta) = SSe/(n-p) * [XH-X']-, and quantile value for the beta distribution is calculated as q = (n-p) / (n-p + 1 * F) which once obtained, the -log10 for such value is the score value.