笔记GWAS 操作流程2-3：MAF过滤

上一次我们经过去掉缺失，去掉错误的性别信息，得到的文件为：

代码语言：javascript复制

HapMap_3_r3_6.bed  HapMap_3_r3_6.fam  HapMap_3_r3_6.log
HapMap_3_r3_6.bim  HapMap_3_r3_6.hh

这里，我们根据最小等位基因频率（MAF）去筛选。

「为什么要根据MAF去筛选？」

❝最小等位基因频率怎么计算？比如一个位点有AA或者AT或者TT，那么就可以计算A的基因频率和T的基因频率，qA qT = 1，这里谁比较小，谁就是最小等位基因频率，比如qA = 0.3, qT = 0.7，那么这个位点的MAF为0.3. 之所以用这个过滤标准，是因为MAF如果非常小，比如低于0.02，那么意味着大部分位点都是相同的基因型，这些位点贡献的信息非常少，增加假阳性。更有甚者MAF为0，那就是所有位点只有一种基因型，这些位点没有贡献信息，放在计算中增加计算量，没有意义，所以要根据MAF进行过滤。 ❞

1. 去掉性染色体上的位点

「思路：」

在map文件中选择常染色体，提取snp信息
根据snp信息进行提取

「提取常染色体上的位点名称：」

因为这里是人的数据，所以染色体只需要去1~22的常染色体，提取它的家系ID和个体ID，后面用于提取。

代码语言：javascript复制

awk '{ if($1 >=1 && $1 <= 22) print $2}' HapMap_3_r3_6.bim > snp_1_22.txt

代码语言：javascript复制

wc -l snp_1_22.txt
1399306 snp_1_22.txt

常染色体上共有1399306位点。

2. 提取常染色体上的位点

这里，用到了位点提取参数--extract

代码语言：javascript复制

plink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7

「查看日志：」

代码语言：javascript复制

PLINK v1.90b6.5 64-bit (13 Sep 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to HapMap_3_r3_7.log.
Options in effect:
  --bfile HapMap_3_r3_6
  --extract snp_1_22.txt
  --make-bed
  --out HapMap_3_r3_7

515185 MB RAM detected; reserving 257592 MB for main workspace.
1431211 variants loaded from .bim file.
163 people (79 males, 84 females) loaded from .fam.
112 phenotype values loaded from .fam.
--extract: 1399306 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 51 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.998153.
1399306 variants and 163 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls.  (51 phenotypes
are missing.)
--make-bed to HapMap_3_r3_7.bed   HapMap_3_r3_7.bim   HapMap_3_r3_7.fam ...
done.

可以看到，共有163个基因型，共有1399306个SNP

3. 计算每个SNP位点的基因频率

首先，通过参数--freq，计算每个SNP的MAF频率，通过直方图查看整体分布。可视化会更加直接。

代码语言：javascript复制

plink --bfile HapMap_3_r3_7 --freq --out MAF_check

结果文件：MAF_check.frq预览：

4. 对基因频率作图

「R代码：」

代码语言：javascript复制

maf_freq <- read.table("MAF_check.frq", header =TRUE, as.is=T)
pdf("MAF_distribution.pdf")
hist(maf_freq[,5],main = "MAF distribution", xlab = "MAF")
dev.off()

「图形如下：」

可以看出，很多基因频率为0，说明没有分型，这些位点需要删掉。

4. 去掉MAF小于0.05的位点

代码语言：javascript复制

plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8

日志：

代码语言：javascript复制

PLINK v1.90b6.5 64-bit (13 Sep 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to HapMap_3_r3_8.log.
Options in effect:
  --bfile HapMap_3_r3_7
  --maf 0.05
  --make-bed
  --out HapMap_3_r3_8

515185 MB RAM detected; reserving 257592 MB for main workspace.
1399306 variants loaded from .bim file.
163 people (79 males, 84 females) loaded from .fam.
112 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 51 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.998153.
325518 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
1073788 variants and 163 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls.  (51 phenotypes
are missing.)
--make-bed to HapMap_3_r3_8.bed   HapMap_3_r3_8.bim   HapMap_3_r3_8.fam ...
done.

可以看到，325518个位点被删掉了，剩余1073788 个位点。

「结果文件：」

代码语言：javascript复制

HapMap_3_r3_8.bed  HapMap_3_r3_8.bim  HapMap_3_r3_8.fam  HapMap_3_r3_8.log

后面我们用这个文件，进行后续的质控。

相关系列：

笔记GWAS操作流程1：下载数据

笔记GWAS操作流程2-1：缺失质控

笔记GWAS操作流程2-2：性别质控

0 人点赞