How to get the residuals?

[moniquevdor]

Hi there, I am having some issues trying to get the residuals after running my models.
I am using this tutorial as an example, and I am running BSLMM. The file has steps for using GEMMA to fit a BSLMM and then get predicted values of the phenotype from BSLMM estimates. Then my plan is to subtract predicted values of phenotype from observed values of phenotype, which would get me the residuals.

I made my kinship matrix with a neutral set of SNPs:

./gemma -bfile "$GENO_DIR/caespitosa_gemma_neutr" \
           -gk 1 \
           -o "caespitosa_kinship"

Then ran BSLMM on each separate trait

PHENO_FILE="input_data/gemma_pheno/gemma_pheno_caespitosa.txt" 
BFILE="../gen_data/caespitosa_gemma"

for TRAIT in 1 2 3 4 5 ; do
    echo "Running trait $TRAIT with LMM"
    nvalue=$((TRAIT + 2)) 
    ./gemma -bfile ${BFILE} \
            -k output/caespitosa_kinship.cXX.txt \
            -p ${PHENO_FILE} \
            -n $nvalue \
            -bslmm \
            -o caespitosa_pheno_T${TRAIT}
done

And then am trying to get the predicted values the same way:

PHENO_FILE="input_data/gemma_pheno/gemma_pheno_caespitosa.txt" 
BFILE="../gen_data/caespitosa_gemma"

for TRAIT in 1 2 3 4 5; do
    echo "Running trait $TRAIT with predict"
    nvalue=$((TRAIT + 2))

    # Run GEMMA prediction
    ./gemma -bfile ${BFILE} \
        -p ${PHENO_FILE}.txt \
        -n ${nvalue} \
        -epm output/caespitosa_pheno_T${TRAIT}.param.txt \
        -emu output/caespitosa_pheno_T${TRAIT}.log.txt \
        -predict \
        -o caespitosa_pheno_T${TRAIT}_prdt
done

But I am getting the same error every time

Running trait 1 with predict
GEMMA 0.98.5 (2021-08-25) by Xiang Zhou, Pjotr Prins and team (C) 2012-2021
Reading Files ...
ERROR: Enforce failed for Problem reading FAM file (phenotypes do not match geno file) in src/gemma_io.cpp at line 600 in ReadFile_fam

I don't understand why it does not work, as it works fine in my bslmm. This is what my input files look like

(outlier_analysis) [mvan0045@m3-login1 gemma]$ head -n 10 input_data/gemma_pheno/gemma_pheno_caespitosa.txt
GCMa05  GCMa05  NA      18.9859945586842        1.00176543053093        0.626840007778175       0
CHMa01  CHMa01  NA      18.6030023090797        1.15752251830044        0.715276050768868       0
CPMa10  CPMa10  NA      19.2448561913047        1.05030718781501        0.682523328403877       0
GLMa02  GLMa02  NA      18.6357921316338        1.18760930990084        0.801529994231572       0
PIFa10  PIFa10  9.21044036697652        NA      0.671512768604101       0.693483505944199       1
CPMa01  CPMa01  NA      19.6791417364033        0.921480575339401       0.383187016894056       0
GLMa01  GLMa01  NA      18.2463273687124        1.15752251830044        0.844676718586457       0
FPMa06  FPMa06  NA      20.1854115425048        1.0885856600443 0.715276050768868       0
GLFa07  GLFa07  8.5173931714189 NA      0.351086219275306       0       1
FPMa02  FPMa02  NA      20.1144598066585        0.833861106417453       0.569030896263964       0
(outlier_analysis) [mvan0045@m3-login1 gemma]$ head ../gen_data/caespitosa_gemma.fam
GCMa05 GCMa05 0 0 0 -9
CHMa01 CHMa01 0 0 0 -9
CPMa10 CPMa10 0 0 0 -9
GLMa02 GLMa02 0 0 0 -9
PIFa10 PIFa10 0 0 0 -9
CPMa01 CPMa01 0 0 0 -9
GLMa01 GLMa01 0 0 0 -9
FPMa06 FPMa06 0 0 0 -9
GLFa07 GLFa07 0 0 0 -9
FPMa02 FPMa02 0 0 0 -9
(outlier_analysis) [mvan0045@m3-login1 gemma]$ head ../gen_data/caespitosa_gemma.bim
chromosome_10   chromosome_10_265148    0       265148  G       C
chromosome_10   chromosome_10_299552    0       299552  A       T
chromosome_10   chromosome_10_427034    0       427034  C       A
chromosome_10   chromosome_10_427051    0       427051  C       T
chromosome_10   chromosome_10_597889    0       597889  T       A
chromosome_10   chromosome_10_598222    0       598222  T       C
chromosome_10   chromosome_10_610851    0       610851  G       A
chromosome_10   chromosome_10_611080    0       611080  A       G
chromosome_10   chromosome_10_611081    0       611081  T       C
chromosome_10   chromosome_10_611101    0       611101  A       T
(outlier_analysis) [mvan0045@m3-login1 gemma]$ head output/caespitosa_pheno_T1.param.txt
chr     rs      ps      n_miss  alpha   beta    gamma
chromosome_10   chromosome_10_265148    265148  0       1.269917e-04    5.626064e-02    2.340000e-03
chromosome_10   chromosome_10_299552    299552  0       -2.212754e-04   -1.234974e-01   1.490000e-03
chromosome_10   chromosome_10_427034    427034  0       -5.300191e-05   -1.037424e-01   3.100000e-04
chromosome_10   chromosome_10_427051    427051  0       6.667576e-05    3.782375e-02    1.600000e-03
chromosome_10   chromosome_10_597889    597889  0       -1.533992e-04   -3.244584e-01   2.500000e-04
chromosome_10   chromosome_10_598222    598222  0       -2.044246e-04   -1.249578e-01   1.080000e-03
chromosome_10   chromosome_10_610851    610851  0       -1.236232e-04   -5.875985e-02   2.440000e-03
chromosome_10   chromosome_10_611080    611080  0       -1.276230e-04   -4.073762e-02   2.530000e-03
chromosome_10   chromosome_10_611081    611081  0       3.461184e-04    1.834051e-01    1.420000e-03
(outlier_analysis) [mvan0045@m3-login1 gemma]$ head -n 50 output/caespitosa_pheno_T1.log.txt
##
## GEMMA Version    = 0.98.5 (2021-08-25)
## Build profile    = /gnu/store/8rvid272yb53bgascf5c468z0jhsyflj-profile
## GCC version      = 7.5.0
## GSL Version      = 2.6
## OpenBlas         = OpenBLAS 0.3.9  - OpenBLAS 0.3.9 DYNAMIC_ARCH NO_AFFINITY Prescott MAX_THREADS=128
##   arch           = Prescott
##   threads        = 1
##   parallel type  = threaded
##
## Command Line Input = ./gemma -bfile ../gen_data/caespitosa_gemma -k output/caespitosa_kinship.cXX.txt -p input_data/gemma_pheno/gemma_pheno_caespitosa.txt -n 3 -bslmm -o caespitosa_pheno_T1
##
## Date = Wed Oct  8 08:49:43 2025
##
## Summary Statistics:
## number of total individuals = 152
## number of analyzed individuals = 66
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs/var = 27037
## number of analyzed SNPs/var = 27031
## REMLE log-likelihood in the null model = -87.2542
## MLE log-likelihood in the null model = -85.7597
## pve estimate in the null model = 1.76384e-06
## se(pve) in the null model = 0.701472
## vg estimate in the null model = 0
## ve estimate in the null model = 0
## beta estimate in the null model =
## se(beta) =
## estimated mean = 10.5162
##
## MCMC related:
## initial value of h = 1.76384e-06
## initial value of rho = 1
## initial value of pi = 0.000369946
## initial value of |gamma| = 10
## random seed = -1
## acceptance ratio = 0.300477
##
## Computation Time:
## total computation time = 1.75702 min
## computation time break down:
##      time on calculating relatedness matrix = 0 min
##      time on eigen-decomposition = 1.38e-05 min
##      time on calculating UtX = 0.00058885 min
##      time on proposal = 0.0137673 min
##      time on mcmc = 1.75345 min
##      time on Omega = 0.88554 min
##

If someone could help me see why it is not working I would very much appreciate it. If there is a different way to get the residuals of the model I would also love to hear about it.
Thank you!