介绍如何使用SNP数据进行GWAS分析,软件选择为emmax和tassel。
软件选择
emmax
plink
qqman
CMplot
数据准备
vcf文件:all_snp.vcf
表型数据文件:sample.table
参考脚本
vcf格式转tped格式
plink --vcf ./all_snp.vcf \ # 输入文件
--maf 0.05 --geno 0.1 \ #过滤
--recode12 \#输出格式
--output-missing-genotype 0 \ #指定缺失数据用0表示
--transpose \ #指定输出的tped格式
--out snp_filter \ #输出文件前缀
--set-missing-var-ids @:# \
--allow-extra-chr
结果文件
表型数据准备
表型数据的样品id要求有两列,且id的排列顺序与tfam文件的id排列顺序一致
计算亲缘关系矩阵
/home/software/emmax/emmax-kin-intel64 -v \#打开日志模式
-d 10 \#指定精度
-o ./pop.kinship \#指定输出文件名称
snp_filter #指定输入SNP数据前缀
gwas分析
/home/software/emmax/emmax-intel64 -v -d 10 \ #
-t snp_filter \ # 输出snp数据的前缀
-p sample.table \ # 表型数据
-k pop.kinship \ # 亲缘关系矩阵
-o emmax.out \ # 输出文件前缀
1> emmax.log 2>emmax.error
结果文件为emmax.out.ps,最后一列为pvalue
gwas结果绘图
#生成绘图的输入文件
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
输入文件需要有4列,snp名称、染色体、位置、pvalue
Rscript ./manhattan_cmplot.R emmax.out.ps.manht_input emmax.out.ps.manht_figure