[GWAS]基于plink的亲缘关系质控

本篇我们要对一下具有亲子关系的个体进行过滤，然后计算类似于IBS的结果
血缘同源（Identity By Descent，IBD），子代中共有的等位基因来源于同一祖先。
状态同源（identical by state ，IBS），两个个体拥有相同的等位基因（不一定来源以同一祖先）。
不太了解的同学👉IBD/IBS

这里讲亲子关系的移除其实不是必须要的，比如我们分析的群体里面有亲子关系的个体，想要进行分析，不需要做这一步的筛选

首先我们计算pihat>0.2的组合

说明文档：

–genome invokes an IBS/IBD computation, and then writes a report with the following fields to plink.genome:
FID1 👉 Family ID for first sample
IID1 👉 Individual ID for first sample
FID2 👉 Family ID for second sample
IID2 👉 Individual ID for second sample
RT 👉 Relationship type inferred from .fam/.ped file
EZ 👉 IBD sharing expected value, based on just .fam/.ped relationship
Z0 👉 P(IBD=0)
Z1 👉 P(IBD=1)
Z2 👉 P(IBD=2)

接着我们提取Z1大于0.9的个体

awk '{if($8>0.9) print $0}' pihat_min0.2.genome > zoom_pihat.genome
wc -l zoom_pihat.genome

我们用R语言可视化这些个体

# 先要设置工作路径

relatedness = read.table("pihat_min0.2.genome", header=T)
with(relatedness,plot(Z0,Z1, xlim=c(0,1), ylim=c(0,1), type="n"))
with(subset(relatedness,RT=="PO") , points(Z0,Z1,col=4))
with(subset(relatedness,RT=="UN") , points(Z0,Z1,col=3))
# legend(1,1, xjust=1, yjust=1, legend=levels(relatedness$RT), pch=16, col=c(4,3))

relatedness_zoom = read.table("zoom_pihat.genome", header=T)
with(relatedness_zoom,plot(Z0,Z1, xlim=c(0,0.02), ylim=c(0.98,1), type="n"))
with(subset(relatedness_zoom,RT=="PO") , points(Z0,Z1,col=4))
with(subset(relatedness_zoom,RT=="UN") , points(Z0,Z1,col=3))
# legend(0.02,1, xjust=1, yjust=1, legend=levels(relatedness$RT), pch=16, col=c(4,3))

pdf("hist_relatedness.pdf")
relatedness = read.table("pihat_min0.2.genome", header=T)
hist(relatedness[,10],main="Histogram relatedness", xlab= "Pihat")
dev.off()

PO是指：亲子关系
UN是指：非亲缘关系
然后我们删除亲子关系的个体：

plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11

我们可以看到其中52个个体被移除了

然后我们来查看数据
这里的文件是bed二进制文件，不方便查看，我们将其转化为ped文件和map文件
注意：这里我们使用的是ped和map格式，如果ped文件中有表型数据（第六列），如果想指定表型数据，就使用–pheno参数，包括三列：家系，个体，表型值

plink --bfile HapMap_3_r3_11 --recode --out test

从统计数据中得出：

• 基因型个体：110个（ped文件中）
• SNP个数：1073181（map文件中）

希望能帮到你~