专栏首页育种数据分析之放飞自我笔记 | GWAS 操作流程2-4:哈温平衡检验

笔记 | GWAS 操作流程2-4:哈温平衡检验

「什么是哈温平衡?」

❝哈迪-温伯格(Hardy-Weinberg)法则 哈迪-温伯格(Hardy-Weinberg)法则是群体遗传中最重要的原理,它解释了繁殖如何影响群体的基因和基因型频率。这个法则是用Hardy,G.H (英国数学家) 和Weinberg,W.(德国医生)两位学者的姓来命名的,他们于同一年(1908年)各自发现了这一法则。他们提出在一个不发生突变、迁移和选择的无限大的随机交配的群体中,基因频率和基因型频率将逐代保持不变。---百度百科 ❞

「怎么做哈温平衡检验?」

「卡方适合性检验!」 ,一个群体是否符合这种状况,即达到了遗传平衡,也就是一对等位基因的3种基因型的比例分布符合公式:p2+2pq+q2=1,p+q=1,(p+q)2=1.基因型MM的频率为p2,NN的频率为q2,MN的频率为2pq。MN:MN:NN=P2:2pq:q2。MN这对基因在群体中达此状态,就是达到了遗传平衡。如果没有达到这个状态,就是一个遗传不平衡的群体。但随着群体中的随机交配,将会保持这个基因频率和基因型分布比例,而较易达到遗传平衡状态。应用Hardy-Weinberg遗传平衡吻合度检验方法,把计算得到的基因频率代入,计算基因型平衡频率,再乘以总人数,求得预期值(e)。把观察数(O)与预期值(e)作比较,进行χ2检验。病例组和对照组的基因型分布的观察值和预期值差异无显著性(P>0.05),符合遗传平衡定律. ❞

「哈温平衡过滤和MAF过滤的区别?」

❝之前,我对这两个概念有点混淆,后来明白过来了。这两个概念一个是对基因频率进行的筛选,一个是对基因型频率进行的筛选。对于一个位点“AA AT TT”,其中A的频率为基因频率,AA为基因型频率。MAF直接是对基因频率进行筛选,而哈温平衡检验,则是根据基因型推断出理想的(AA,AT,TT)的分布,然后和实际观察的进行适合性检验,然后得到P值,根据P值进行筛选。即P值越小,说明该位点越不符合哈温平衡。 ❞

「两个目的:」

  • 计算所有位点的哈温检测结果
  • 删除SNP中不符合哈温平衡的位点

1. 计算所有位点的HWE的P值

plink --bfile HapMap_3_r3_8 --hardy

plink.hwe的数据格式:

  • CHR 染色体
  • SNP SNP的ID
  • TEST 类型
  • A1 minor 位点
  • A2 major 位点
  • GENO 基因型分布:A1A1, A1A2, A2A2
  • O(HET) 观测杂合度频率
  • E(HET) 期望杂合度频率
  • P 哈温平衡的卡方检验P-value值

结果预览:

2. 提取哈温p值小于0.0001的位点

这里我们使用awk:

awk '{if($9 < 0.0001) print $0}' plink.hwe >plinkzoomhwe.hwe

共有123个位点,其中UNAFF为45个位点。

3. 设定过滤标准1e-4

plink --bfile HapMap_3_r3_8 --hwe 1e-4 --make-bed --out HapMap_3_r3_9

日志:

Options in effect:
  --bfile HapMap_3_r3_8
  --hwe 1e-4
  --make-bed
  --out HapMap_3_r3_9

515185 MB RAM detected; reserving 257592 MB for main workspace.
1073788 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.998136.
--hwe: 45 variants removed due to Hardy-Weinberg exact test.
1073743 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_9.bed + HapMap_3_r3_9.bim + HapMap_3_r3_9.fam ...
done.

可以看到,共有45个SNP根据哈温的P值过滤掉了,和上面手动计算的一样。

4. 可视化

R代码:

hwe<-read.table (file="plink.hwe", header=TRUE)
pdf("histhwe.pdf")
hist(hwe[,9],main="Histogram HWE")
dev.off()

hwe_zoom<-read.table (file="plinkzoomhwe.hwe", header=TRUE)
pdf("histhwe_below_theshold.pdf")
hist(hwe_zoom[,9],main="Histogram HWE: strongly deviating SNPs only")
dev.off()

哈温的P值直方图:

过滤掉SNP位点的P值:

过滤后的结果文件

HapMap_3_r3_9.bed  HapMap_3_r3_9.bim  HapMap_3_r3_9.fam  HapMap_3_r3_9.log

本文分享自微信公众号 - 育种数据分析之放飞自我(R-breeding),作者:邓飞2013

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-04-17

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 笔记 | GWAS 操作流程2-1:缺失质控

    清洗数据的时间占80%的时间,有句话这样讲:“Garbage in, Garbage out(垃圾进,垃圾出)”,所以清洗数据非常重要,今天学习一下基因组数据如...

    邓飞
  • 笔记GWAS 操作流程2-3:MAF过滤

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

    邓飞
  • 混合线性模型学习笔记2

    这一个章节主要是介绍混线性模型的应用,其实我们很多本科时候学的统计学知识(大都是一般线性模型,回归分析,方差分析等等)都可以放在混合线性模型的框架下进行分析,就...

    邓飞
  • 怎样将IP地址转换为数字

    首先我们选择一个IP地址,这里我找了一个,PING一下,是通的。好就用这个IP地址做测试吧。

    似水的流年
  • 《程序员的自我修养》第二章学习笔记

    第二章 编译和链接 2.1被隐藏了的过程 我们知道,一个程序由源代码到可执行文件往往由这几步构成: 预处理(Prepressing)-> 编译(Compilat...

    xcywt
  • 我的小工具-远程读卡器web客户端(PHP+LUA)

    本工具是在浏览器中以Lua脚本的形式对CPU卡,M1卡就行读、写等各种操作,配和使用改造过后的E711读卡器。 远端把读卡器接到电脑上,并运行读写卡服务...

    特立独行的猫a
  • 为什么有些程序员不愿意缩进代码?

    作为已经写了十几年代码的老程序员,虽然在编写代码的时候大部分情况还是遵循编码规范,但在这基础上会展示自己一些特性,有些程序员不喜欢缩进代码也是源于此,如同一个人...

    程序员互动联盟
  • GPS Tracking Devices for Security Purposes

    Today, the utilization of GPS has become a piece of our regular day to day exist...

    用户7437050
  • Satpy基础系列教程(3)-Satpy总览

    Satpy is designed to provide easy access to common operations for processing met...

    zhangqibot
  • JavaWeb高级编程(下篇)

    JSP标签语法中包含一些简写可以帮助轻松编写JSP。这些简写中第一个就是taglib指令。

    范中豪

扫码关注云+社区

领取腾讯云代金券