1. 基因组选择中H矩阵的构建
这里的1为非测序个体, 2为测序个体. A11, A12, A21, A22可以由系谱构建的A矩阵提取. G为基因组构建的矩阵. H矩阵构建的相关代码见: 【GS专栏】全基因组选择中如何构建H矩阵.
2. 直接构建H^{-1}矩阵
一步法混合线性方程组中, 直接利用的是H逆矩阵, 因此直接构建H逆矩阵更加方便.
H逆矩阵构建方法:
3. 构建H逆矩阵需要考虑的参数
3.1 G矩阵矫正到A22矩阵的尺度
原因:
G矩阵构建是由测序个体的基因频率构建的,它的频率可能与A矩阵的基因频率(可以追溯到基础群体base population)不一样, 这导致G矩阵和A22矩阵存在尺度上的差异, 因此需要矫正.
矫正方法是建立两个方程组:
原因: 对于某些性状, G矩阵可能无法解释所有的变异, 这时可以设定A矩阵解释的比例.
对于某些性状, 适当设置tau 和 omega, 可以矫正估算的无偏性.
4. 完整的H逆矩阵构建方法
5. 构建H逆矩阵Julia代码展示
Julia是新一代的计算科学语言, 速度非常快, 比R语言、Python要快很多, 内存占用也比较少. 比C语言、Fortran要简单很多,是一门未来的明星语言. 是时候尝试一下新东西了.
代码思路:
矩阵
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