【康普森GS专栏】基因组选择中构建H矩阵需要设置哪些参数?

function hmatrix_julia_adjust(id_full,id_geno,Amatrix,Gmatrix,a =0.95,b=0.05,tau=1,omega=1)    print("G* = a*G + b*A22, a is ",a,"; b is ",b,"\n")    print("iG = tau*G - omega*A22, tau is ",tau,"; omega is ",omega,"\n\n")
    function diag(mat)        dd = [mat[i,i] for i in 1:size(mat,1)]        return dd    end
     id1 = id_full     idb = id_geno     A = Amatrix     G = Gmatrix     ida = setdiff(id1,idb)     A = NamedArray(A,(id1,id1))     G = NamedArray(G,(idb,idb))     reb = idb     rea = ida
    A22 = A[reb,reb]    diagA22 = diag(A22)    offdiagA22 = []    for i in 1:size(A22,1)        for j in 1:size(A22,2)            if(i != j )                push!(offdiagA22,A22[i,j])            end        end    end
    diagG = diag(G)    offdiagG = []    for i in 1:size(G,1)        for j in 1:size(G,2)            if(i != j )                push!(offdiagG,G[i,j])            end        end    end    print("\n\nStatistic of Rel Matrix G\n","\t\t","N","\t","Mean","\t","Min","\t","Max","\t","\n",      "Diagonal\t",length(diagG),"\t",round.(mean(diagG),digits=3),"\t",      round.(minimum(diagG),digits=3),"\t",      round.(maximum(diagG),digits=3),"\t",      "\nOff-diagonal\t",length(offdiagG),"\t",      round.(mean(offdiagG),digits=3),"\t",round.(minimum(offdiagG),digits=3),"\t",      round.(maximum(offdiagG),digits=3),"\t","\n\n")
    print("\n\nStatistic of Rel Matrix A22\n","\t\t","N","\t","Mean","\t","Min","\t","Max","\t","\n",    "Diagonal\t",length(diagA22),"\t",round.(mean(diagA22),digits=3),"\t",    round.(minimum(diagA22),digits=3),"\t",    round.(maximum(diagA22),digits=3),"\t",    "\nOff-diagonal\t",length(offdiagA22),"\t",    round.(mean(offdiagA22),digits=3),"\t",round.(minimum(offdiagA22),digits=3),"\t",    round.(maximum(offdiagA22),digits=3),"\t","\n\n")
    # 根据方程计算alpha和beta两个参数值    meanoffdiagG=mean(offdiagG)    meandiagG=mean(diagG)    meanoffdiagA22=mean(offdiagA22)    meandiagA22=mean(diagA22)    print("Means_off_diag\t","Means_diag:\n","G\t",meanoffdiagG,"\t",meandiagG,    "\nA22\t",meanoffdiagA22,"\t",meandiagA22,"\n")
    beta=(meandiagA22 - meanoffdiagA22)/(meandiagG - meanoffdiagG)    alpha=meandiagA22-meandiagG*beta    print("Adjust G, and the value of alpha and beta is:",alpha,beta,"\n\n\n")    G=alpha .+ beta .* G #    # 在将a和b的参数加上    G = a .* G + b .* A22
    iA22 = inv(A22)    iG =inv(G)    iA = inv(A)    x22 = tau*iG - omega*iA22
    diagx22 = diag(x22)    offdiagx22 = []    for i in 1:size(x22,1)        for j in 1:size(x22,2)            if(i != j )                push!(offdiagx22,x22[i,j])            end        end    end
    print("\n\nStatistic of Rel Matrix x22\n","\t\t","N","\t","Mean","\t","Min","\t","Max","\t","\n",      "Diagonal\t",length(diagx22),"\t",round.(mean(diagx22),digits=3),"\t",      round.(minimum(diagx22),digits=3),"\t",      round.(maximum(diagx22),digits=3),"\t",      "\nOff-diagonal\t",length(offdiagx22),"\t",      round.(mean(offdiagx22),digits=3),"\t",round.(minimum(offdiagx22),digits=3),"\t",      round.(maximum(offdiagx22),digits=3),"\t","\n\n")
   # 构建H逆矩阵   iHaa = iA[rea,rea]   iHab = iA[rea,reb]   iHba = iHab'   iHbb = iA[reb,reb] + x22   function cbind(A,B)       X = vcat(A',B')'       return(X)   end   function rbind(A,B)       X = vcat(A,B)       return(X)   end   iH = cbind(rbind(iHaa,iHba),rbind(iHab,iHbb))   return(round.(iH,digits=6))end