使用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

欢迎关注Bioinfor 生信云！