全基因组选择中G矩阵和H矩阵构建时的计算效率问题的研究

摘要

全基因组选择，可以将表型数据、基因型数据和系谱数据整合起来，计算中需要对亲缘关系矩阵A和基因组亲缘关系矩阵G进行求逆。这篇文章利用模拟数据，研究不同算法对不同数据量的计算效率的问题。
分子数据：模拟的40k芯片数据
测序个体：30,000个个体
计算方法：
1，do循环，包括两个版本
2，利用矩阵乘法的子例程
3，广义求逆是利用LAPACK的子例程
结果显示，利用最有效的计算程序，利用并行处理，30,000个个体的40k的基因组数据求逆，需要花费3个小时的时间。

矩阵的操作效率很高，可以利用线性代数的子例程（basic Linear Algebra Subroutines, BLSA）,也有类似的开源的项目ATLAS，这些库可以指定特定的运算流程（内存和缓存的大小）

计算机存储层次机构（Computer memory hierarchy）

在现代计算机中，内存运行的趋势是通过加快缓存的交换以及降低缓存的大小。速度的限制是缓存的大小和主要内存交换的大小。但对于多线程的计算，每个线程都有自己的缓存，但是主内存是共享的。

计算方法

1，do循环是将Z矩阵直接读进内存中
2，OPTM，是Z矩阵防止外部循环中，避免间接寻址
3，OPTML，将主循环分成两个循环

计算程序

矩阵的乘法是利用BLAS的子例程DGEMM，设定的方法可以采用ATLAS或者因特尔的MKL（Intel  Math Kernel Library）

避免G矩阵奇异性的方法

G = wG + (1-w)A22，这里w是权重，主要是避免对角线为0造成的计算困难，A22为测序个体的系谱亲缘关系矩阵

测试结果

1，缓存的大小，影响运算速度，缓存有1Mb提升到6Mb时，运算速度提高3~10倍。

2，不同的数据量，三种方法运行的时间，可以看出，Optimized DGEMM，即采用MKL定制的算法，效果最好

3，采用MKL的DGEMM随着CPU数目的提升，运算速度也不断增加

结论

最佳的是采用MKL的DGEMM算法，在40K，测序个体为30,000时，计算H逆矩阵的时间大约是3小时。