项目文件列表
如何运行
cmke . source /home/asc/intel/compilers_and_libraries/linux/bin/compilervars.sh make其中第二句是声明intel icpc 和 intel icc 编译器，以便make的时候使用，否则会出错。
然而这样也会出现一个问题，诸如malloc is undefined 等符号未定义问题，原因在开会过程中一直找，从gcc 版本，cmake版本，intel icpc icc编译器等方面入手，但是在开会过程中并没有找到，最后由另一位成员发现，更改gcc版本为较新的7.5，在编译的时候使用新的gcc可以解决这个问题.
以下是成功编译的截图
从工程的角度来讲，一个软件不应该对其它的软件库太过依赖，另外去除软件对特殊类库的依赖也有利于软件的推广何使用。对于我们实现的gwas软件，依赖于intel mkl比较多，而intel mkl 的使用需要安装Intel Parallel Studio，而且这个软件的安装还需要认证激活，如果是学生，还可以进行学生认证使用，但是其它科研工作者使用就并不是那么简单了，而且安装使用激活过程对于非专业人员也是有一定难度的，并且还有一定的经费支出。
综上所述，gwas软件去掉intel数学库无论在工程上还是软件适配性角度来说，都是有必要的。
下面的cmake编译配置说明了，当前的gwas使用的是intel的编译器
SET(CMAKE_CXX_COMPILER icpc) SET(CMAKE_C_COMPILER icc) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -g -mkl -static-intel -no-intel-extensions -static-libstdc++ -march=native -qopenmp")这些行声明了c++,c编译器分别是icpc、icc并且编译参数中使用了-mkl -static-intel来添加相关依赖
1.分析依赖库有那些
通过通读代码，发现依赖库主要是在计算部分，如下图：
#include <mkl.h> //X2 = Xt * X cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, M, M, N, 1, X, N, X, N, 0, X2, M);简单分析可知依赖库是intel mkl 和 BLAS
//solve(Xt * X) int inv_success = LAPACKE_dgesv(LAPACK_COL_MAJOR, M, M, X2, M, ipiv, B, M);依赖库是LAPACKE
info = LAPACKE_dgesdd(LAPACK_COL_MAJOR, 'A', M, N, X, M, d, U, M, Vt, N); //transpose mkl_dimatcopy ('C', 'T', N, N, 1.0, Vt, N, N);依赖库是LAPACKE intel MKL
2.如何编译安装依赖库
经过1中的分析，和对代码的理解，初步的gcc通用编译步骤就是使用LAPACKE、BLAS和理解并重写mkl_dimatcopy矩阵转置函数。
LAPACKE (https://github.com/Reference-LAPACK/lapack/tree/master/LAPACKE)
BLAS采用OpenBLAS (http://www.openblas.net/)
以上类库的安装和编译过程在这里就不多说了，都是采用基本的makefile或cmake进行编译安装
但是依赖库的动态链接需要在cmake中进行配置
如下图为配置号的cmake文件
设置编译器为g++,gcc
去掉CMAKE_CXX_FLAGS中的类库依赖
添加liblapack.a静态链接库 linopenblas_haswellp-r0.3.10.so 动态链接库
解决编译过程中的gfortran相关函数未定义的错误，添加libgfortran.a静态链接库
3.部分类库函数的实现
最后一个类库依赖是mkl_dimatcopy，首先需要清楚这个函数的输入输出参数，然后才能根据调用方式进行代码函数重构：
在intel mkl document 可以找到这个函数
在理解和读懂文档和函数的作用之后，重构函数，对项目进行重新编译
矩阵就地转置 重新编译 运行测试
1.为什么要进行分割
我们拿到的数据为UKB的50w人的21号染色体上snp的数据（以下简称bgen数据），单从大小上来讲的话就有37GB，snp的个数也是达到了1261158个，而且我们使用的机器是4核CPU，内存也只有15GB，另外我们还需要对数据据进行解析存储到c++的各种容器中，需要对数据中的值进行处理，在计算部分也需要进行矩阵运算，即使有多线程的加持，这也将是一大笔的资源消耗，因此数据切分变得尤为重要，这对数据读取和广义线性模型计算部分都将有很大的优化必要，而且利用不同大小的数据可以方便进行代码调试和性能统计分析。
2.如何进行分割
分割工具安装使用bgenix 对bgen数据进行切分，很方便的获得bgen数据的子集
https://enkre.net/cgi-bin/code/bgen/dir?ci=trunk
直接在文件目录中执行以下命令即可编译出bgenix可执行文件
./waf configure ./wafbuild/apps/目录下可以看到可执行文件bgenix
分割脚本书写
要想完整的运行整个数据集，对数据的有效分割是必不可少的
大体的思路是：按照块的大小整个数据划分成 总数/块大小 个文件（也需要处理不足块大小的数据），这样方便后续的文件读入和计算，这也是并行化处理和计算的需要，因为要想使得整个程序的处理过程更快，可以分布式进行bgen数据处理和计算，因此将数据按照块大小划分可以更好的在集群上运行，集群上的每一台机器采取作业调度的管理方式。
完整的bgen数据按照块大小block_size进行划分的脚本
#!/bin/bash #user configurations bgen_file="./ukb_imp_chr21_v3.bgen" out_dir="./split" block_size=10000 if [ "${bgen_file##*.}" != "bgen" ] then echo "${bgen_file} is not a bgen file" exit 0 fi out_prefix="`basename -- ${bgen_file%.*}`" bgi_file="${bgen_file}.bgi" echo ${out_prefix} #checking bgen file if [ ! -e "${bgen_file}" ] then echo "${bgen_file}" echo "bgen file not exist!" echo "Please provide bgen file!" exit 0 fi #rebuild index? if [ ! -e "${bgi_file}" ] then echo "${bgi_file}" echo "index not exist, rebuilding" ./bgenix -g ${bgen_file} -index fi total_count=`sqlite3 -header -csv ${bgi_file} "SELECT COUNT(*) FROM Variant" | tail -n 1` echo "total SNPs: $total_count" n_blocks=`expr $total_count / $block_size` block_tail=`expr $total_count % $block_size` echo "n_blocks: $n_blocks" echo "block_tail: $block_tail" limit=$block_size offset=0 #dealing with body echo "dealing with body" for((i=0;i<n_blocks;i++)) do echo "limit: $limit" echo "offset: $offset" file_name="${out_prefix}_${i}.bgen" echo ${out_dir}/${file_name} ./bgenix -g ${bgen_file} -limit $limit -offset $offset > ${out_dir}/${file_name} offset=`expr $offset + $limit` done #dealing with tails echo "dealing with tail" echo "limit: $limit" echo "offset: $offset" file_name="${out_prefix}_${i}.bgen" echo ${out_dir}/${file_name} ./bgenix -g ${bgen_file} -limit $limit -offset $offset > ${out_dir}/${file_name}脚本逻辑
首先必须提供ukb_imp_chr21_v3.bgen 输入的数据文件
其次用户可选的提供 ukb_imp_chr21_v3.bgen.bgi 文件，这个文件是索引文件，用于快速找到需要切割的数据的位置，加快切割速度。如果用户不进行提供的话，使用bgenix进行索引的重构
#rebuild index file ./bgenix -g ${bgen_file} -index 首先统计所有的数据的条数 sqlite3 -header -csv ${bgi_file} "SELECT COUNT(*) FROM Variant" | tail -n 1其中sqlite3是数据库可执行文件(下载安装参考 https://www.sqlite.org/download.html)，
-header 和 -csv 表示输出的格式是带有列名的csv文件
然后再对这个文件从Variant表统计总的个数，tail -n 1 即是吧统计结果的最后一行拿到，这个值就是总的snp的个数。
其次计算出块的个数和剩余的数据个数 n_blocks=`expr $total_count / $block_size` block_tail=`expr $total_count % $block_size`循环对数据进行切分，计算出偏移位置offset，然后取之后的limit条数据，定向输出到名为file_name的文件中
for((i=0;i<n_blocks;i++)) do echo "limit: $limit" echo "offset: $offset" file_name="${out_prefix}_${i}.bgen" echo ${out_dir}/${file_name} ./bgenix -g ${bgen_file} -limit $limit -offset $offset > ${out_dir}/${file_name} offset=`expr $offset + $limit` done最后把剩下的尾巴做同样的输出处理
./bgenix -g ${bgen_file} -limit $limit -offset $offset > ${out_dir}/${file_name}2020-06-30 snptest gwas 分析测试对比 为了测试所设计的gwas软件的性能如何，我安装运行snptest来进行对比 1.snptest 安装 https://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.html#download
2.snptest 测试对比 切分500条数据并放到机器上进行测试，文件列表是这样的：
使用以下脚本进行运行，指定输出文件夹和使用的sample_file， -pheno lung表示要分析snp与lung上的表型数据的关系 -con_names 指的是协变量 -chunk 指的是划分整个数据为及部分来进行分析处理 最后日志输出到run.log中 另外脚本读入第一个命令行参数为*.bgen文件，即切分的500条数据
#!/bin/bash output_dir="result" sample_file="ukb_imp_chr_v3.sample" time ./snptest -chunk 10 -data $1 $sample_file - o./$output_dir/$1.out -frequentist 1 -method score -pheno lung -cov_names pc1 pc2 pc3 pc4 pc5 pc6 pc7 pc8 pc9 pc10 &> run.log运行结果如下图，由于结果比较长，这里只贴出首尾的部分： 使用head,tail 命令打印首尾20行输出，可以看到开始时间和结束时间： Analysis: “SNPTEST 2.5.4-beta3 analysis, started 2020-06-30 15:10:31” Completed successfully at 2020-06-30 17:43:37 总的时间是约2小时33分钟6秒，即就是9186秒
对于我们实现的gwas软件呢，队友给我的单线程运行时间是1000个snp，运行时间是138.9s，相比于snptest 的500个snp的测试时间9186秒有了极大的速度提升，预估1000个snp来进行snptest时间将大于18000s,也就是说我们实现的gwas提升了将近130倍，下图展示了运行时间上的差异：
可以看到再四线程的环境下，竟然是280倍之多，其中提升速度的主要部分是计算部分glm（广义线性模型的使用）。
