这里使用rMVP包进行GWAS分析,分析之前需要做如下准备,
需要有基因型文件和表型文件
,
从表型文件中获取样本ID
cut -f1 phe > ALL-id.txt
sed -i '1d' ALL-id.txt # 删除文件第一列的表头
####提取VCF子集
bcftools view -S ALL-id.txt SE.vcf > ALL.vcf
# ped、map
plink --vcf ALL.vcf --recode 12 --out ALL --allow-extra-chr
cat ALL.map | awk '{printf("%s\t%s_%s\t%s\t%s\n",$1,$1,$4,$3,$4)}' > test.map
cp ALL.map bk.ALL.map
mv test.map ALL.map
# bed
plink --noweb --file ALL --make-bed --out ALL --allow-extra-chr
# 生成.raw 文件,即为 012 矩阵,缺失的点标记为 NA:
plink --bfile ALL --recodeA --out ALL --allow-extra-chr
#将 2-6 列删除,并格式化 RS ID:
cat ALL.raw | cut -d" " -f1,7- | sed 's/_[A-Z]//g' >ALL.genotype.txt
# 矩阵转置
#awk '{i=1;while(i <= NF){col[i]=col[i] $i " ";i=i+1}} END {i=1;while(i<=NF){print col[i];i=i+1}}' genotype.txt | sed 's/[ \t]*$//g' > genotype_transpose.txt
# 获取SNP位置信息
echo "snpid chr pos" >ALL.snpsloc.txt
awk '{print $2,"chr"$1,$4}' ALL.map >>ALL.snpsloc.txt
曼哈顿图
结果路径
网友评论