GWAS（2）——gwas相关的文件处理命令总结

筛选需要的ID

当我们发现有些样本需要保留或剔除时，可以用vcftools进行处理。参考链接：https://www.omicsclass.com/article/647

• 例子
  假设我们根据需要保留一些样本，先将需要的样本ID存为txt文件（如yes.txt），再用命令：

vcftools --vcf test_vcf.vcf --recode --recode-INFO-all --stdout --keep yes.txt >
out.vcf

得到了需要的样本vcf文件。

计算缺失率、杂合率和最小等位基因频率（maf、het、missing）

maf——plink --noweb --vcf out --freq --out base.original.freq
het——plink --vcf out.vcf --freqx --out Het.snp.vcf --allow-extra-chr
missing——vcftools --vcf out.vcf --missing-site –out 

得到文件可以进行三种频率的计算。

基因型过滤

1. 首先进行 maf0.05 的过滤，将 maf 小于等于 0.05 的基因型过滤掉。

vcftools --vcf out.vcf --maf 0.05 --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out 1out.vcf

计算过滤后新文件的三种频率，作图。
计算过滤后新文件的每条染色体的标记数。

2. 接着进行缺失率 0.15 的过滤，将 missing 大于 0.15 的标记过滤掉。

plink --bfile out --geno 0.15 --make-bed --out clean
#或者
plink --bfile out --geno 0.15 --out clean.vcf

计算过滤后新文件的三种频率，作图。
计算过滤后新文件的每条染色体的标记数。

最后进行杂合率0.15的过滤，将het大于0.15的标记过滤掉。
需要在 R 中操作，命令大致如下：
①首先 read.table 读取步骤 2 中得到的杂合率结果文件；
②接着用公式，

het$HET_RATE = het$C.HET./(het$C.HOM.A1.+het$C.HET.+het$C.HOM.A2.)

计算位点杂合率（实际上上一步应该已经得到了）
③筛选，用 data <- [which(het$ HET_RATE > 0.15), ],得到 het 大于 0.15 的标记汇总。
④读写成 txt 文件后，在服务器上运行以下命令：

vcftools --vcf miss0.15.recode.vcf --recode --recode-INFO-all --stdout --exclude no.txt > 1.vcf

就得到了过滤好的最终文件。

计算PCA

计算前10个PCA结果

plink --bfile wgas1 --pca 10 --out wgas1_pca 

PCA作图R代码

rm(list = ls())
library(ggplot2)
library(RColorBrewer)
data <- read.table("wgas1_pca.eigenvec",sep = "\t",header = T)
ggplot(data,aes(x=pc1,y=pc2,color=species))+ geom_point()
#######
tiff(filename = "pca.tif", width = 30, height = 16, units = "cm", compression = "lzw", res = 600)
ggplot(data,aes(x=pc1,y=pc2,color=species))+ 
  geom_point()+ 
  theme_bw() +
  theme(panel.border=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.line= element_line(colour = "black"))
dev.off()