全基因组关联分析（GWAS）软件：gemma
版本：Version 0.94.1
输入文件
GEMMA 需要四个主要的输入文件，包括基因型、表型、相关矩阵和协变量（可选）。
基因型和表型文件可以是两种格式，都是 PLINK binary ped 格式或者都 BIMBAM 格式。不可混用。
输入文件
接受比较多的输入类型，这是用的 plink 格式：
GEMMA 识别基因型和表型的 PLINK binary ped 文件格式.
该格式需要三个文件：*.bed, *.bim, *.bam，所有文件都有相同的前缀。
用户可以用 PLINK 将标准 ped 文件转换成 binary ped 文件：
命令：
plink --file [file_prefix] --make-bed --out [bedfile_prefix]
对于 *.fam 文件，GEMMA 仅读取第二列（individual id）和第六列（phenotype）.
用户可以定义其它列作为 phenotype 列，使用 ‘-n [num]’ 参数指定要分析的表型，’-n 1’ 使用原始的第六列作为 phenotypes，’-n 2’ 使用第七列，以此类推。
GEMMA 用 0/1 编码等位基因。
*.bim 文件的第 5 列是 minor allele ，编码为 1；第 6 列为 major allele ，编码为 0.
因此，第 5 列的 minor allele 是效应基因。
GEMMA 将按照提供的方式读入表型，’-9‘ 和 ‘NA’ 作为缺失表型。
pink --file ms_0.4_NAU --maf 0.05 --geno 0.2 --recode --make-bed --out snp --noweb
- snp.bim
- snp.bed
- snp.fam：gemma 仅读取第 2 行（材料 id）和第 6 行（性状）
计算相关矩阵
./gemma -bfile [prefix] -gk [num] -o [prefix]
- -gk：数字表示不同的类型，一般用 1
- -bfile：输入 bed 文件
注意：用 genotype 计算相关矩阵。相关矩阵的计算不依赖于具体的 phenotype，但是如果表型值有缺失，会在计算的时候剔除相关材料，所以可以用任意不含缺失值的数据计算相关矩阵，对于同一套 genotype ，可以用一个相关矩阵
Association Tests with Multivariate Linear Mixed Models
./gemma -bfile [prefix] -k [filename] -lmm [num] -n [num1] [num2] [num3] -o [prefix]
- -bfile：bed 格式 genotype
- -k：相关矩阵
- -lmm：检测方法。默认
1，Wald test;2: likelihood ratio test;3：score test - -n：表型。可以指定多个表型，进行多表型效应分析。也可以仅指定一个表型，此时与单变量混合线性模型一样。
为了用 -n 参数方便循环，在此使用多变量混合线性模型进行检测，但是仅指定一个变量（实际同单变量混合线性模型相同）
注意: 似乎 -n 的最大为 46（即一个 .fam 文件中最多放 46 个性状），多了就无法识别了。
运行 gemma
#!/bin/bash
# gemma 软件路径
gemma=gemma
# 输入文件，此处为 bed 格式
bed_file="snp"
# 相关矩阵
rel_ma="snp.cXX.txt"
# 性状数量
num="8"
# missingness
ms=0.2
#maf
maf=0.05
for i in `seq 1 ${num}`
do
nohup ${gemma} -bfile ${bed_file} -k ${rel_ma} -lmm 1 -miss ${ms} -maf ${maf} -n ${i} -o ${i} &
done
输出文件
共有两个输出文件，都在当前目录下的文件夹中。
‘preix.log.txt’ 文件有一些关于运行参数和计算时间的细节信息。除此之外，还包含 PBE 检测和它的标准误。
’prefix.assoc.txt’ 文件包含结果。例如：
chr rs ps n_mis n_obs allele1 allele0 af beta se p_wald
1 rs3683945 3197400 0 1410 A G 0.443 -1.586575e-01 3.854542e-02 4.076703e-05
1 rs3707673 3407393 0 1410 G A 0.443 -1.563903e-01 3.855200e-02 5.252187e-05
1 rs6269442 3492195 0 1410 A G 0.365 -2.349908e-01 3.905200e-02 2.256622e-09
1 rs6336442 3580634 0 1410 A G 0.443 -1.566721e-01 3.857380e-02 5.141944e-05
1 rs13475700 4098402 0 1410 A C 0.127 2.209575e-01 5.644804e-02 9.497902e-05
这 11 列是：染色体、SNP id、碱基对在染色体上的位置、相应 SNP 的缺失个体数量、相应 SNP 的无缺个体数量、最小等位基困、最大等位基因、等位基因频率、beta 检测、beta 标准误、Wald 检验的 P 值。