EMMAX是关联分析速度提升的一个代表性算法, 已广泛应用于棉花、大豆和水稻等的复杂性状关联分析。EMMAX认为，每一个SNP对复杂性状的解释率都很低，每个组件的方差在运算中都只计算一次，从而达到简化计算的目的。所以，当存在有大解释率的SNP时，违反了EMMAX的假设，这可能会产生保守的p值，导致统计能力的损失。
EMMA和EMMAX的区别
EMMA是计算每个标记自己的方差，而EMMAX假设每一个SNP的方差都很小，因此在运算过程当中只会有一个方差。
1. 下载及安装
1.1 下载地址
http://csg.sph.umich.edu//kang/emmax/download/index.html
1.2 安装
tar xvf emmax-intel-binary-20120210.tar.gz
# 调用
emmax-intel64
2. 准备工作
2.1 tped格式基因型文件
emmax识别.tped，.tfam，.nosex格式，利用plink进行转置：
./plink --vcf root.id.vcf --recode 12 transpose --out emmax_in --allow-extra-chr --chr-set 27
--vcf vcf文件名
--out 输出文件前缀
--allow-extra-chr 允许scaffold等染色体格式
--chr-set 设置染色体数目
2.2 亲缘关系 emmax_in.BN.kinf
2.2.1 emmax计算亲缘关系
./emmax-kin-intel64 emmax_in -v -d 10 -o emmax_in.BN.kinf
2.2.2 Tassel或GCTA中计算得到的亲缘关系
8. GWAS：亲缘关系——TASSEL&GCTA - 简书 (jianshu.com)
将之前计算得到的亲缘关系去掉注释、行名和列名，保存为纯数字的矩阵emmax_in.BN.kinf。
2.3 群体结构 emmax.cov.txt
Admixture中最佳分群数所对应的Q文件，手动添加行名列名，这里最佳分群为3。
随后，在R中，利用emmax_in.tfam和上一步整理好的Q文件进行整理。
> tfam <- read.table("emmax_in.tfam", header = F, stringsAsFactors = F)
> covs <- read.table("Q.txt", header = T, check.names = F, stringsAsFactors = F, skip = 1)
> covs <- covs[match(tfam$V1, covs$'<Trait>'), ]
> covs <- cbind(covs[,1], covs[,1], rep(1, nrow(covs)), covs[,2:ncol(covs)])
> write.table(covs[, -ncol(covs)], file = "emmax.cov.txt", col.names = F, row.names = F, sep = "\t", quote = F)
2.4 表型数据 emmax.trait.**.txt
emmax中一次只能输入一个性状，文件格式与.tfam类似，需要在R中进行转换。
> tfam <- read.table("emmax_in.tfam", header = F, stringsAsFactors = F)
> tr <- read.table("trait.txt", header = T, check.names = F, stringsAsFactors = F)
> head(tr)
<Trait> TL SA AD NR
1 1 46.86 10.10 0.692 22.03
2 2 32.98 8.73 0.871 15.56
3 3 37.59 9.31 0.840 18.46
4 4 24.12 6.67 0.804 20.23
5 5 32.47 7.77 0.827 19.60
6 6 27.44 8.17 0.957 23.00
> tr <- tr[match(tfam$V1, tr$'<Trait>'),]
> tr[tr == -999] <- NA
> tre <- cbind(tr[,1], tr)
> for(i in 3:ncol(tre)){
file_name <- paste0("emmax.trait.", names(tre)[i], ".txt")
write.table(tre[, c(1, 2, i)], file = file_name, col.names = F, row.names = F, sep = "\t", quote = F)
}
每个性状，各生成一个数据文件，格式如下：
3.关联分析
3.1 只利用亲缘关系Kinship进行关联分析
./emmax-intel64 -t emmax_in -o TL_emmax.k -p emmax.trait.TL.txt -k emmax_in.BN.kinf
-t：基因型文件名（除扩展名）
-o：输出文件前缀
-p：输入性状文件
-k：亲缘关系文件
3.2 K+Q进行关联分析
./emmax-intel64 -t emmax_in -o TL_emmax.qk -p emmax.trait.TL.txt -k emmax_in.BN.kinf -c emmax.cov.txt
保存.emmax.qk.ps文件，加标题，第一列ID，第二列回归系数(beta),第三列回归系数的标准差，第四列为P值。
4.在R中进行结果可视化——曼哈顿图，QQ图
将.emmax.qk.ps文件整理为pmap格式，第一列为ID，第二列为chr，第三列position，第四列为p值
> library("CMplot")
> pmap <- read.table("pmap.txt", header = T)
> head(pmap)
ID chr position p
1 Chr1__5154 Chr1 5154 0.8291300
2 Chr1__5187 Chr1 5187 0.1632200
3 Chr1__5220 Chr1 5220 0.0048982
4 Chr1__50889 Chr1 50889 0.4160400
5 Chr1__50920 Chr1 50920 0.0510200
6 Chr1__50956 Chr1 50956 0.1364600
# 阈值计算
> threshold <- 1/nrow(pmap[!is.na(pmap$position),])
# 画图
> CMplot(pmap, threshold = threshold, amplify = F, memo = "", file = "tiff", plot.type=c("m","q"))
plot.type="m"：曼哈顿图
plot.type：可以选择 "d", "c", "m", "q" or "b"
plot.type="d", SNP密度图
plot.type="c", 环形曼哈顿图
plot.type="m",曼哈顿图
plot.type="q",QQ图
plot.type="b",同时绘制环形曼哈顿图、曼哈顿图和QQ图
plot.type=c("m","q"), 绘制曼哈顿图和QQ图
引用转载请注明出处，如有错误敬请指出。
