使用Q矩阵作为协变量使用emmax进行GWAS分析
数据准备
1.上个推文中vcf转的tped文件
snp_filter.map
snp_filter.nosex
snp_filter.tfam
snp_filter.tped
2.表型数据
sample.table
3.Q矩阵
做群体结构分析的时候产生的Q矩阵
snp.3.Q
4.亲缘关系矩阵
pop.kinship
参考脚本
修改Q矩阵格式
前面加样品名，并去掉Q矩阵的最后一列
awk '{print $1}' snp_filter.tfam | \
paste - snp.3.Q | \
awk '{print $1" "$1" 1 "$2" "$3}' > pop.Qmatrix
gwas分析
/home/software/emmax/emmax-intel64 -v \
-d 10 \ #指定小数点位数
-t snp_filter \ #输入文件前缀
-p sample.table \#表型数据
-k pop.kinship \#亲缘关系矩阵
-c pop.Qmatrix \#Q矩阵
-o emmax.out \#输出文件
1> emmax.log 2>emmax.err #日志文件
结果绘图
#生成绘图文件
paste snp_filter.map emmax.out.ps | \
awk 'BEGIN{print "SNP\tCHR\tBP\tP"}{if($2==$5){print $2"\t"$1"\t"$4"\t"$NF}}' > emmax.out.ps.manht_input
#绘图
Rscript ./manhattan_cmplot.R \
emmax.out.ps.manht_input \
emmax.out.ps.manht_figure
emmax.out.ps.manht_figure_manhattan_threshold.png
emmax.out.ps.manht_figure_qqplot.png
