笔记 GWAS 操作流程2-3：最小等位基因频率

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

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，后面用于提取。

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

wc -l snp_1_22.txt
1399306 snp_1_22.txt

常染色体上共有1399306位点。

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

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

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

查看日志：

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_7515185 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频率，通过直方图查看整体分布。可视化会更加直接。

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

结果文件：MAF_check.frq
预览：

4. 对基因频率作图

R代码：

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的位点

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

日志：

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_8515185 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 个位点。

结果文件：

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

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

笔记 GWAS 操作流程2-3：最小等位基因频率相关推荐

笔记 GWAS 操作流程5-1：根红苗正的GWAS分析软件：GEMMA
笔记 GWAS 操作流程5-1:根红苗正的GWAS分析软件:GEMMA 1. GEMMA软件介绍这个肯定厉害了,是大家闺秀,是名门望族,是根红苗正的GWAS分析软件. GEMMA名称来源: G: G ...
笔记 GWAS 操作流程6-2：手动计算GWAS分析中的GLM和Logistic模型
1. 名词解释 GWAS 全基因组关联分析手动计算使用R语言编程GLM模型和Logistic模型,提取Effect和Pvalue GLM 一般线性模型 Logistic 主要分析广义线性模型,Y变 ...
笔记 GWAS 操作流程2-4：哈温平衡检验
什么是哈温平衡? 哈迪-温伯格(Hardy-Weinberg)法则哈迪-温伯格(Hardy-Weinberg)法则是群体遗传中最重要的原理,它解释了繁殖如何影响群体的基因和基因型频率.这个法则是用H ...
笔记 GWAS 操作流程5-2：利用GEMMA软件进行LMM+PCA+协变量
这里,我们用正常的GWAS分析,考虑所有的协变量(数值协变量+因子协变量)+ PCA协变量,然后用混合线性模型进行分析. 1. 协变量文件 c.txt文件 1 1 0 0 -0.0169445 0.0 ...
笔记 GWAS 操作流程2-6：去掉亲缘关系近的个体
这里,我们要对一些亲子关系的个体,进行一下过滤,计算类似IBS的结果. 注意: 这里讲亲子关系的个体移除,不是必须要的,比如我们分析的群体里面有亲子关系的个体,想要进行分析,不需要做这一步的筛选. 1 ...
plink描述性统计--等位基因频率、缺失值
0.3 描述性统计 0.3.1 等位基因频率 --freq,产生的文件后缀为.frq,该文件包含基因型的等位基因和最小等位基因频率(MAF)和每个SNP的等位基因编码的信息. plink --bfil ...
plink, vcftool计算等位基因频率（allele frequency，vcf）
计算等位基因频率有两种方式,第一种用vcftool计算: /path/to/vcftools --vcf file.vcf --freq --chr 1 --out filefreq 很简单的一个命令 ...
Genome Aggregation Database (gnomAD) 简介 | 参考人群等位基因频率数据库
Genome Aggregation Database (gnomAD) 这是一个关于什么的数据库?broad institute开发的,整合了目前几乎所有的公共的WES和WGS测序数据,并对数据做了 ...
关于 minor allele frequency（次等位基因频率）的理解
引用自NCBI的概念(https://www.ncbi.nlm.nih.gov/projects/SNP/docs/rs_attributes.html#gmaf) Global minor alle ...