笔记 GWAS 操作流程2-5：杂合率检验

一般自然群体，基因型个体的杂合度过高或者过低，都不正常，我们需要根据杂合度进行过滤。偏差可能表明样品受到污染，近亲繁殖。我们建议删除样品杂合率平均值中偏离±3 SD的个体。

我的理解：非自然群体中，比如自交系，杂交种F1，这些群体不需要过滤杂合度。

参数过滤和手动过滤
plink有个特点，所有的过滤标准，都可以生成过滤前的文件，然后可以手动过滤，也可以用参数进行过滤。

比如：--missing生成结果，也可以用--geno和--mind过滤。
比如：--hardy生成结果，可以使用--hwe过滤
比如：--freq生成结果，可以用--maf过滤
但是杂合度--het，没有过滤的函数，只能通过编程去提取ID，然后用--remove去实现。

1. 计算杂合度

plink --bfile HapMap_3_r3_9 --het --out R_check

结果文件：

.het (method-of-moments F coefficient estimates)
Produced by --het.A text file with a header line, and one line per sample with the following six fields:FID Family ID  # 家系ID
IID Within-family ID # 个体ID
O(HOM)  Observed number of homozygotes # 实际纯合个数
E(HOM)  Expected number of homozygotes # 期望纯合个数
N(NM)   Number of non-missing autosomal genotypes # 总个数
F   Method-of-moments F coefficient estimate #

所以，F的计算方法：
F=O−EN−EF = \frac{O-E}{N-E}F=N−EO−E

O: O(HOM)
E: E(HOM)
N: N(NM)

2. 杂合度可视化

R代码：

het <- read.table("R_check.het", head=TRUE)
pdf("heterozygosity.pdf")
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
hist(het$HET_RATE, xlab="Heterozygosity Rate", ylab="Frequency", main= "Heterozygosity Rate")
dev.off()

可视化：

3. 计算杂合度三倍标准差以外的个体

首先，查看哪些个体在3倍标准差以外：

het <- read.table("R_check.het", head=TRUE)
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
het_fail = subset(het, (het$HET_RATE < mean(het$HET_RATE)-3*sd(het$HET_RATE)) | (het$HET_RATE > mean(het$HET_RATE)+3*sd(het$HET_RATE)));
het_fail$HET_DST = (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE);
write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)

结果可以看出，这两个个体杂合度在3倍标准差以外：

4. 去掉这些个体

cat fail-het-qc.txt
"FID" "IID" "O.HOM." "E.HOM." "N.NM." "F" "HET_RATE" "HET_DST"
1330 "NA12342" 703770 691000 1066256 0.03402 0.33996151018142 -4.75272486494316
1459 "NA12874" 709454 695300 1072858 0.03758 0.338725162137021 -5.18288854902555

先对数据进行清洗，去掉引号，然后提取家系和个体ID

sed 's/"//g' fail-het-qc.txt |awk '{print $2}' > het_fail_ind.txt

sed 's/"//g' fail-het-qc.txt |awk '{print $1,$2}' > het_fail_ind.txt

使用remove去掉这两个个体

plink --bfile HapMap_3_r3_9 --remove het_fail_ind.txt --make-bed --out HapMap_3_r3_10

5. 结果文件

HapMap_3_r3_10.bed  HapMap_3_r3_10.bim  HapMap_3_r3_10.fam  HapMap_3_r3_10.log

笔记 GWAS 操作流程2-5：杂合率检验相关推荐

笔记 GWAS 操作流程2-4：哈温平衡检验
什么是哈温平衡? 哈迪-温伯格(Hardy-Weinberg)法则哈迪-温伯格(Hardy-Weinberg)法则是群体遗传中最重要的原理,它解释了繁殖如何影响群体的基因和基因型频率.这个法则是用H ...
笔记 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 操作流程5-2：利用GEMMA软件进行LMM+PCA+协变量
这里,我们用正常的GWAS分析,考虑所有的协变量(数值协变量+因子协变量)+ PCA协变量,然后用混合线性模型进行分析. 1. 协变量文件 c.txt文件 1 1 0 0 -0.0169445 0.0 ...
笔记 GWAS 操作流程2-3：最小等位基因频率
上一次我们经过去掉缺失,去掉错误的性别信息,得到的文件为: HapMap_3_r3_6.bed HapMap_3_r3_6.fam HapMap_3_r3_6.log HapMap_3_r3_6.bi ...
笔记 GWAS 操作流程2-6：去掉亲缘关系近的个体
这里,我们要对一些亲子关系的个体,进行一下过滤,计算类似IBS的结果. 注意: 这里讲亲子关系的个体移除,不是必须要的,比如我们分析的群体里面有亲子关系的个体,想要进行分析,不需要做这一步的筛选. 1 ...
Java计算连续自交杂合概率代系变化
在上一篇博文里我们解析了Rosalind中Independent Alleles一题.在题目背后,小编重新思考了一下题目的设定.Tom的家系如果在每代不能保证都是AaBb杂交产生的后代,而是从初代To ...
HaploMerger2: 从高杂合二倍体基因组组装中重建单倍型
本文只是按照自己的需求翻译了HaploMerger2提供的手册部分内容.HaploMerger2的帮助文档写的非常好,一定要花点时间去读啊! HaploMerger2的分析流程如下重建单倍体组装中的 ...
vcf文件（call variants得来的）怎么看变异是纯合还是杂合的
如下图片所示: 对于位置为48245131的allele来说,REF为A,ALT为C 想确定变异到底是纯合还是杂合,即两条染色体是否同时发生了变异,则看GT,GT对应的数值为0/1,说明该变异为杂合: ...
「文献」杂合基因组的策略思路
「文献」杂合基因组的策略思路文献出处: Sequencing a Juglans regia × J. microcarpa hybrid yields high-quality genome as ...