对从ADNI数据库下载的SNP数据及进行质控（QC）

• 简介
• 一、先查看数据中的个体和SNP缺失情况
• • 1.查看
  • 2.生成绘图以可视化缺失结果。
• 二、QC第一步：删除缺失度大于某个数值的SNP和个体
• • 1.删除SNP
  • 2.删除个体
  • 3.重复1操作
  • 4.重复2操作
• 三、QC第二步：检查性别（非必要步骤）
• • 1.检查性别
  • 2.可视化结果
• 3.以下两个脚本可用于处理性别差异的个体。（可选，但是我没做）
• • 1.删除存在性别差异的个体
  • • 1）找出并生成相应文本
    • 2）根据文本删除个体数据
  • 2.修改 input-sex
• 四、QC第三步：最小等位基因频率（MAF）
• • • 定义：
    • 意义：
    • 功能：
    • 解释：
  • 1.仅包含常染色体的SNP（可选，但是我做的全染色体）
  • • 1.生成1-22染色体文档
    • 2.根据文档提取数据
    • 3.生成MAF数据
    • 4.生成MAF图
  • 2.全染色体的SNP的MAF图
  • 3.筛除MAF频率较低的SNP
  • • 1.筛除命令
    • 2.查看结果
• 五、QC第四步：Hardy-Weinberg平衡（HWE）
• • 1.Hardy-Weinberg平衡（HWE）
  • • 名词解释：
    • 如何检测：
    • 意义：
    • 运用：
  • 1.查看所有SNP的HWE p值分布
  • 2.选择HWEp值低于0.00001的SNP，并生成图
  • • 1.筛选p值低于0.00001的SNP，并显示结果
    • 2.生成图
  • 3.使用–hwe进行质量控制
  • • 1.先用于对照组（control），筛除p值低于1e-6的SNP
    • 2.然后对于全体（case），筛除p值低于1e-10的SNP
• 六、QC第五步：控制杂合率偏差
• • 1.杂合率
  • • 名词解释：
    • 意义：
    • 计算：
  • 2.检查SNP的杂合率
  • • 1.筛除和修剪
    • 2.查看prune.in和prune.out文件
  • 3.对符合要求的SNP进行–het操作
  • • 1.–het操作
    • 2.查看R-check文件
    • 3.生成杂合率图
  • 4.筛除超过设定杂合率阈值的SNP
  • • 1.使用教程提供的heterozygosity_outliers_list.R 来生成列表
  • 2.根据列表筛除数据
• 七、QC第六步：清除可能存在的亲缘样本
• • 1.relatedness
  • • 名词解释：
    • 作用：
    • 运用：
  • 2.分析的数据集之间的隐含关联性
  • • 1.分析大于阈值0.2的个体
  • 2.查看文件
  • 3.可视化亲缘关系
  • 3.删除所有关联性
  • • 1.过滤非初代个体
  • 4.手动删除响应率高的相关性个体

简介

时间过于仓促，所以本博客过于简陋，后续会分小节补充完整。

本流程参考github上的教程。

https://github.com/MareesAT/GWA_tutorial/

参考论文
https://www.ncbi.nlm.nih.gov/pubmed/29484742

我们下载了ADNI Omni2.5M microarray SNP data
共三个二进制文件，为了方便操作，将其同一命名为a
a.bed,a.bim,a.fam

bed 包括了个体标识符和基因型
bim包含遗传标记
fam包含个体有关信息

一、先查看数据中的个体和SNP缺失情况

1.查看

plink –bfile a –missing

输出了四个文件
plink.hh plink.imiss plink.limss plink.log
由于未指定文件名，所以输出默认plink。
其中imiss文件记录了个体缺失的基因的信息
如下图：
FID父系
IID个体编号
MISS_PHENO表型是否缺失，Y表示缺失
N_MISS 缺失基因个数
N_GENO 基因总数
F_MISS 缺失率

lmiss记录各SNP缺失信息。
CHR 染色体编号
SNP 单核酸多态性编号

2.生成绘图以可视化缺失结果。

Rscript –no-save hist_miss.R
我们来看下R脚本里面的代码

indmiss<-read.table(file="plink.imiss", header=TRUE)
snpmiss<-read.table(file="plink.lmiss", header=TRUE)
#read data into R，传输数据给R包，
#read.table()函数从具有多列表格形式的文件中读取数据。
#使用好它可以简单的从文本文件或CSV这种文件中读取数据。 
#header 逻辑值。数据文件中是否有标题行。
#如果header设置为TRUE，则要求第一行要比数据列的数量少一列。
pdf("histimiss.pdf")#indicates pdf format and gives title to file生成的文件格式pdf，并赋予标题
hist(indmiss[,6],main="Histogram individual missingness") 
#hist表示画出直方图
#selects column 6, names header of file选择第六列数据为x值，然后赋予文件头pdf("histlmiss.pdf") 
hist(snpmiss[,5],main="Histogram SNP missingness")  
dev.off() 

二、QC第一步：删除缺失度大于某个数值的SNP和个体

1.删除SNP

plink –bfile a –geno 0.2 –make-bed –out a_2

这回生成四个文件

过程：
2379855 variants loaded from .bim file.
1094 variants removed due to missing genotype data (–geno).
2378761 variants and 812 people pass filters and QC.
删除了1094个SNP

2.删除个体

plink –bfile a_2 –mind 0.2 –make-bed –out a_3

同样生成四个文件
切记，两次顺序不能调换，必须先进行删除SNP的操作。

2378761 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
0 people removed due to missing genotype data (–mind).
说明没有个体数据被删除。

3.重复1操作

plink –bfile a_3 –geno 0.02 –make-bed –out a_4

2378761 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
61766 variants removed due to missing genotype data (–geno).
2316995 variants and 812 people pass filters and QC.
删除了61766个SNP
生成a_4

4.重复2操作

plink –bfile a_4 –mind 0.02 –make-bed –out a_5

2316995 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
0 people removed due to missing genotype data (–mind).
没有个体数据被删除
生成a_5

三、QC第二步：检查性别（非必要步骤）

1.检查性别

根据先验知识确定女性受试者的F值必须小于0.2，确定男性受试者的F值必须大于
0.8。这个F值是基于X染色体近亲交配(纯合性)估计的。
#不符合这些要求的受试者会被PLINK标记为“问题”
plink –bfile a_5 –check-sex

不指定输出名字，默认为plink前缀名字。
结果：
–check-sex: 46816 Xchr and 0 Ychr variant(s) scanned, no problems detected.
Report written to plink.sexcheck .
没有性别出错个体
并生成plink.sexcheck

2.可视化结果

虽然没有个体出错，但是基于入门学习，还是要做一下，了解代码。
Rscript –no-save gender_check.R

在了解代码前，先看一下plink.checksex的文件内部
FID IID PEDSEX SNPSEX STATUS F
1 067_S_0056 2 2 OK -0.1409
2 123_S_0108 1 1 OK 0.995
3 067_S_0059 2 2 OK -0.1285
4 024_S_0985 1 1 OK 0.9408
5 109_S_2237 1 1 OK 0.9382
6 131_S_0123 1 1 OK 0.9368
7 098_S_0160 1 1 OK 0.938

查看R脚本内部代码

gender <- read.table("plink.sexcheck", header=T,as.is=T)
#header=T，表示将第一行当作变量名。
#as.is是设置特定列的数据类型的，比如numeric,character.
#如果不设置，就是默认的把原来的字符转成factor型的。会对数据处理产生影响的。
pdf("Gender_check.pdf")#文件格式和名称
hist(gender[,6],main="Gender", xlab="F")
#产生数值向量x为第六列的直方图，xlab="F"表示x轴的命名
dev.off()pdf("Men_check.pdf")
male=subset(gender, gender$PEDSEX==1)
hist(male[,6],main="Men",xlab="F")
dev.off()pdf("Women_check.pdf")
female=subset(gender, gender$PEDSEX==2)
hist(female[,6],main="Women",xlab="F")
dev.off()

生成的结果

3.以下两个脚本可用于处理性别差异的个体。（可选，但是我没做）

1.删除存在性别差异的个体

1）找出并生成相应文本

grep “PROBLEM” plink.sexcheck| awk ‘{print$1,$2}’> sex_discrepancy.txt

此命令将生成具有“PROBLEM"状态的个人列表
管道命令 | 执行完前面一个命令后自动执行下个命令，即将前一个命令的输出作为后一
个命令的输入
grep
用于查找文件里符合条件的字符串。
awk
Awk是一个强大的文本分析工具，它每次读入一条记录，并把每条记录切分成字
段后进行分析。Awk官方文档是非常好的学习材料，通过man awk查看。

2）根据文本删除个体数据

plink –a_5 –remove sex_discrepancy.txt –make-bed –out a_6

此命令将删除具有"PROBLEM"状态的个人列表

2.修改 input-sex

plink –bfile a_5 –impute-sex –make-bed –out a_6

这会将基于基因型信息的性别输入到您的数据集中
就是根据基因型的判断信息修改错误的个体性别信息

四、QC第三步：最小等位基因频率（MAF）

定义：

最小等位基因频率通常是指在给定人群中的不常见的等位基因发生频率

意义：

MAF广泛应用于复杂疾病的全基因组关联研究。在关联研究中，较小的MAF将会使统计效能降低，从而造成假阴性的结果。为了研究人群中罕见突变与疾病的关联，通常通过加大样本量的方法来弥补MAF降低所带来的统计效能的损失。在研究中，这一突变位点往往属于同一基因或同一通路上一组罕见突变中的一个。通常认为一个位点的MAF＞0.05才有研究的意义，如果满足不了0.05的话，说明这个位点的突变频率很低

功能：

让数据仅包括高于设定MAF阈值的SNP

解释：

MAF较低的SNP很少见，所以这个操作并不能有效检测SNP表型关联能力。MAF阈值取决于样本量，较大的样本量可以使用较低的MAF阈值。例如
大样本N=100,000 MAF阈值设定为0.01
中样本N=10000 MAF阈值设定为0.05

1.仅包含常染色体的SNP（可选，但是我做的全染色体）

1.生成1-22染色体文档

awk ‘{if($1 >=1&& $1 <=22) print $2 }’ a_5.bim > snp_1_22.txt

2.根据文档提取数据

plink –bfile a_5 –extract snp_1_22.txt –make-bed -out a_6

extract 根据文本内容提取相应的数据生成新的数据。
结果
–extract: 2255862 variants remaining.
–make-bed to a_6.bed + a_6.bim + a_6.fam … done.

3.生成MAF数据

plink –bfile a_6 –freq –out MAF_check
#–freq表示计算等位基因频率
#freq和out之间加上–chr 1 表示在1号染色体里面计算等位基因的频率

结果
–freq: Allele frequencies (founders only) written to MAF_check.frq .

4.生成MAF图

Rscript –no-save MAF_check.R

查看MAF_check.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()

结果

2.全染色体的SNP的MAF图

plink –bfile a_5 –freq –out MAF_check
Rscript –no-save MAF_check.R

详细解释在常染色体过程中有
根据上面两个代码就生成了MAF图

结果：

3.筛除MAF频率较低的SNP

1.筛除命令

plink –bfile a_5 –maf 0.05 –make-bed –out a_6

–maf 0.05
设定0.05为MAF阈值，低于0.05的SNP都将被删除

结果：
2316995 variants loaded from .bim file.
1073372 variants removed due to minor allele threshold(s)

2.查看结果

plink –bfile a_6 recode –out a_6
#我们将b（二进制）文件转换为plink文件
wc -l a_6.ped a_6.map
#查看文件内数据的数目

结果：
812 a_6.ped
1243623 a_6.map
1244435 total
我们可以看到a_5的SNP值确实被筛除了，a_6的SNP数目加上被筛除的数目等于a_5的数目。

小知识：
常规GWAS的传统MAF阈值在0.01%或0.05%之间，这取决于样本大小。

五、QC第四步：Hardy-Weinberg平衡（HWE）

1.Hardy-Weinberg平衡（HWE）

名词解释：

“哈迪­温伯格定律”是指在理想状态下，各等位基因的频率在遗传中是稳定不变的，即保持着基因平衡。该定律运用在生物学、生态学、遗传学。条件：①种群足够大；②种群个体间随机交配；③没有突变；④没有选择；⑤没有迁移；⑥没有遗传漂变。

如何检测：

「卡方检验！」，一个群体是否符合这种状况，即达到了遗传平衡，也就是一对等位基因的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），符合遗传平衡定律。

意义：

卡方检验就是统计样本的实际观测值与理论推断值之间的偏离程度，实际观测值与理论推断值之间的偏离程度就决定卡方值的大小，如果卡方值越大，二者偏差程度越大；反之，二者偏差越小；若两个值完全相等时，卡方值就为0，表明理论值完全符合。

运用：

对于二元性状： 建议排除hwe的p值<1e-10的数据，对照标准<1e-6，不太严格的阈值有助于避免丢弃可能与疾病相关的SNP。
对于定量性状： 建议阈值为<1e-6。

1.查看所有SNP的HWE p值分布

plink –bfile a_6 hardy
head -10 plink.hwe
#展示前10个数据

结果：
–hardy: Writing Hardy-Weinberg report (founders only) to plink.hwe … done.
我们得到了plink.hwe

plink.hwe的前10个数据
CHR SNP TEST A1 A2 GENO O(HET) E(HET) P
0 rs35144699 ALL(NP) T A 88/355/369 0.4372 0.4401 0.8733
0 rs2571902 ALL(NP) A T 0/435/364 0.5444 0.3962 1.076e-36
0 rs2751692 ALL(NP) A T 0/434/365 0.5432 0.3957 2.462e-36
0 rs10015934 ALL(NP) T C 84/353/374 0.4353 0.4361 1
0 rs10155688 ALL(NP) A G 23/179/610 0.2204 0.2387 0.03807
0 rs1040254 ALL(NP) T C 10/114/688 0.1404 0.1514 0.05757
0 rs1042882 ALL(NP) G A 88/333/383 0.4142 0.4327 0.2222
0 rs1045074 ALL(NP) A G 114/385/310 0.4759 0.4707 0.8226
0 rs10502267 ALL(NP) A G 18/179/614 0.2207 0.23 0.2829

其中：
STATUS 是否正常
TEST 类型
A1 minor 位点
A2 major 位点
GENO 基因型分布：A1A1, A1A2, A2A2
O(HET) 观测杂合度频率
E(HET) 期望杂合度频率
P 哈温平衡的卡方检验P-value值

2.选择HWEp值低于0.00001的SNP，并生成图

其中：允许放大严重偏离的SNP

1.筛选p值低于0.00001的SNP，并显示结果

awk ‘{ if($9 <0.00001) print $0 }’ plink.hwe>plinkzoomhwe.hwe
#第九列适配选择选项，然后输出去掉第一行
head -10 plinkzoomhwe.hwe
#显示前10行数据

0 rs2571902 ALL(NP) A T 0/435/364 0.5444 0.3962 1.076e-36
0 rs2751692 ALL(NP) A T 0/434/365 0.5432 0.3957 2.462e-36
0 rs12557859 ALL(NP) G A 71/81/660 0.09975 0.2369 6.022e-44
0 rs2301187 ALL(NP) G T 30/46/735 0.05672 0.1222 3.312e-26
0 rs2306837 ALL(NP) A G 0/431/368 0.5394 0.3939 6.105e-36
0 rs2435805 ALL(NP) G A 0/432/365 0.542 0.3951 2.713e-36
0 rs2437438 ALL(NP) C A 0/430/366 0.5402 0.3943 5.934e-36
0 rs2496945 ALL(NP) G A 0/435/364 0.5444 0.3962 1.076e-36
0 rs2520748 ALL(NP) A G 29/139/630 0.1742 0.2164 6.722e-07
0 rs2522574 ALL(NP) G A 0/432/364 0.5427 0.3954 2.447e-36

2.生成图

Rscript –no-save hwe.R
wc -l plink.hwe
#查看plink.hwe的数据的数目
wc -l plinkzoomhwe.hwe
#同理

hwe.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值低于0.00001的SNP。这样可以感受严重偏离的SNP数据。

结果：

整体SNP的p值图

严重偏离的SNPp值图

3.使用–hwe进行质量控制

1.先用于对照组（control），筛除p值低于1e-6的SNP

代码：
plink –bfile a_6 –hwe 1e-6 –make-bed –out a_hwe_filter_step_1

结果：
1243623 variants loaded from .bim file.
–hwe: 1657 variants removed due to Hardy-Weinberg exact test.
1241966 variants and 812 people pass filters and QC.

2.然后对于全体（case），筛除p值低于1e-10的SNP

用较低的阈值，可有助于避免丢弃可能与疾病相关的SNP。
plink –bfile a_6 –hwe 1e-10 –hwe-all –make-bed –out a_7

结果：
1241966 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
–hwe: 0 variants removed due to Hardy-Weinberg exact test.
1241966 variants and 812 people pass filters and QC.

六、QC第五步：控制杂合率偏差

1.杂合率

名词解释：

杂合性是指某一个位点上含有一对及其以上的不同的等位基因。包括同系合性和同种合性。群体遗传多态性的均匀度的度量常采用杂合度作为参数。杂合性是在同源染色体上的一个或多个位点上有不同等位基因存在的状态，是种群的基本属性之一。

意义：

杂合性是研究自然群体遗传变异的主要研究对象。个体的杂合率是杂合基因的比例。个体内杂合度高可能表明样品质量低，杂合度低可能是因为近亲繁殖所致。

计算：

杂合性（H）的计算公式为：杂合性=每个基因座为杂合子的频率总和/基因座总数

2.检查SNP的杂合率

1.筛除和修剪

代码：
plink –bfile a_7 –exclude inversion.txt –range –indep-pairwise 50 5 0.2 –out indepSNP

注释：
这里用到了inversion.txt文件，这是作者提供的有关LD高区域的SNP编号。具体使用可根据自己为基因数据做出修改。由于我刚入门不太懂LD于SNP之间的相关性，所以我没有改，按部就班的挪用了作者的inversion.txt。
exclude 筛除
参数参数“ 50 5 0.2”分别代表：窗口大小，每步移动窗口经过的SNP数量（步长），以及一个SNP在所有其他SNP上同时回归的多重相关系数（multiple correlation coefficient for a SNP being regressed on all other SNPs simultaneously）(即r^2 threshold)。
命令–indep-pairwise修剪SNP，修剪时不用考虑SNP的p值。

结果：
1241966 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
–exclude range: 16019 variants excluded.
–exclude range: 1225947 variants remaining.
1225947 variants and 812 people pass filters and QC.

Note: No phenotypes present.
–indep-pairwise: Ignoring 2755 chromosome 0 variants.
Pruned 81637 variants from chromosome 1, leaving 11364.
Pruned 88216 variants from chromosome 2, leaving 11158.
Pruned 74905 variants from chromosome 3, leaving 9840.
Pruned 70204 variants from chromosome 4, leaving 9045.
Pruned 66212 variants from chromosome 5, leaving 8687.
Pruned 61754 variants from chromosome 6, leaving 8291.
Pruned 59383 variants from chromosome 7, leaving 7949.
Pruned 55617 variants from chromosome 8, leaving 7021.
Pruned 48742 variants from chromosome 9, leaving 6783.
Pruned 54677 variants from chromosome 10, leaving 7488.
Pruned 52396 variants from chromosome 11, leaving 7109.
Pruned 51783 variants from chromosome 12, leaving 7285.
Pruned 39073 variants from chromosome 13, leaving 5392.
Pruned 35271 variants from chromosome 14, leaving 4947.
Pruned 33086 variants from chromosome 15, leaving 4833.
Pruned 35394 variants from chromosome 16, leaving 5378.
Pruned 28475 variants from chromosome 17, leaving 4878.
Pruned 32410 variants from chromosome 18, leaving 4961.
Pruned 20820 variants from chromosome 19, leaving 4078.
Pruned 25938 variants from chromosome 20, leaving 4160.
Pruned 15402 variants from chromosome 21, leaving 2370.
Pruned 16372 variants from chromosome 22, leaving 2592.
Pruned 25668 variants from chromosome 23, leaving 3595.
Pruned 235 variants from chromosome 24, leaving 19.
Pruned 114 variants from chromosome 25, leaving 156.
Pruned 18 variants from chromosome 26, leaving 11.
Pruning complete. 1073802 of 1223192 variants removed.
Marker lists written to indepSNP.prune.in and indepSNP.prune.out .

这里我们可以看出每个染色体修剪的情况。
最后生成两个文件prune.in（符合LD阈值的不相关SNP）和prune.out（不符合）。

2.查看prune.in和prune.out文件

head -10 indepSNP.prune.in
head -10 indepSNP.prune.out

结果：

rs11240777
kgp8975187
kgp2139601
kgp11435978
kgp8336649
kgp5005988
kgp7172368
kgp11531255
kgp4294569
kgp4467256

rs3094315
rs3131972
kgp3709449
kgp5225889
kgp3392992
kgp7224902
kgp12144557
kgp3674608
kgp5074587
kgp10636897

我们可以得知这两个文件存储的时SNP的位点信息。in文件中包含的就是通过筛选条件的、独立的SNP位点，out文件记录具有LD的SNP位点。

3.对符合要求的SNP进行–het操作

1.–het操作

代码：
plink –bfile a_7 –extract indepSNP.prune.in –het –out R_check

注释：
extract 提取，即提取a_7符合in文件的SNP位点信息的SNP
–het 计算杂合度

结果：
1241966 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
–extract: 149390 variants remaining.

149390 variants and 812 people pass filters and QC.
Note: No phenotypes present.
–het: 145765 variants scanned, report written to R_check.het .

2.查看R-check文件

head -10 R_check.het

结果：
FID IID O(HOM) E(HOM) N(NM) F
1 067_S_0056 102882 9.852e+04 145578 0.09272
2 123_S_0108 98819 9.85e+04 145547 0.006776
3 067_S_0059 103230 9.849e+04 145538 0.1007
4 024_S_0985 99025 9.858e+04 145676 0.00938
5 109_S_2237 100256 9.851e+04 145555 0.0372
6 131_S_0123 98497 9.858e+04 145660 -0.001666
7 098_S_0160 99138 9.834e+04 145314 0.01699
8 027_S_0256 98962 9.846e+04 145493 0.01062
9 116_S_1243 98842 9.857e+04 145655 0.005731

注释：
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系数估计
杂合度计算 杂合度=（N-O）/N

3.生成杂合率图

代码：
Rscript –no-save check_heterozygosity_rate.R

check_heterozygosity_rate.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")
#变量het$HET_RATE,x轴名称，y轴名称
dev.off()

结果：

4.筛除超过设定杂合率阈值的SNP

1.使用教程提供的heterozygosity_outliers_list.R 来生成列表

代码：
Rscript –no-save heterozygosity_outliers_list.R

heterozygosity_outliers_list.R 代码：

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)));
#筛选大于平均值3倍标准差和小于平均值3倍标准差的SNP
#subset 是对数据的某些字段进行操作
#subset(dataset , subset , select )
het_fail$HET_DST = (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE);
#z-score标准化，用来算离数据中心的偏差的
write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)

查看生成的文件fail-het-qc.txt
cat fail-het-qc.txt
wc -l fail-het-qc.txt

结果：
“FID” “IID” “O.HOM.” “E.HOM.” “N.NM.” “F” “HET_RATE” “HET_DST”
1 “067_S_0056” 102882 98520 145578 0.09272 0.293286073445163 -4.31018561529304
3 “067_S_0059” 103230 98490 145538 0.1007 0.290700710467369 -4.73649630913152
44 “072_S_2070” 103145 98580 145673 0.09686 0.291941540299163 -4.53189098636989
49 “128_S_0200” 105319 98590 145680 0.1429 0.277052443712246 -6.98701282075642
119 “041_S_4427” 103067 98370 145365 0.09989 0.290977883259382 -4.69079219474187
187 “029_S_4279” 104605 98550 145627 0.1285 0.281692268604036 -6.22193380249188
251 “031_S_4476” 105947 98570 145665 0.1566 0.272666735317338 -7.71018956803152
304 “052_S_1250” 101889 98550 145621 0.07095 0.300313828362667 -3.15135140101021
367 “109_S_4531” 105157 98570 145658 0.1398 0.278055444946381 -6.82162399306853
384 “013_S_4395” 104268 98380 145362 0.1254 0.282701118586701 -6.05558055160113
496 “073_S_4403” 106851 98370 145366 0.1804 0.264951914477938 -8.98231679083646
611 “109_S_4260” 102149 98550 145624 0.07646 0.298542822611657 -3.4433795212354
614 “014_S_4577” 102819 98360 145335 0.09502 0.292537929610899 -4.43355000157678
621 “109_S_2200” 102412 98480 145526 0.08357 0.296263210697745 -3.8192737167699
677 “035_S_0292” 105183 98320 145279 0.1461 0.275993089159479 -7.16169396931848
681 “127_S_0684” 104892 98550 145625 0.1347 0.279711587982833 -6.548536038677
685 “029_S_0878” 104972 98440 145460 0.139 0.278344562078922 -6.77395032928023
690 “023_S_0061” 102715 98510 145551 0.08948 0.294302340760283 -4.14260929121193
754 “031_S_4021” 102176 98270 145210 0.08329 0.296357000206597 -3.80380839490325
796 “037_S_0539” 104444 98560 145632 0.1251 0.282822456602944 -6.03557264774465

21 fail-het-qc.txt
共21个个体数据不符合要求
注释：
HET_RATE 杂合度
HET_DST z-score分数

2.根据列表筛除数据

sed ‘s/"// g’ fail-het-qc.txt | awk ‘{print$1, $2}’> het_fail_ind.txt
plink –bfile a_7 –remove het_fail_ind.txt –make-bed –out a_8
cat het_fail_ind.txt

注释：
通过从文件中删除所有引号并仅选择前两列，使该文件与PLINK兼容。
–remove 删除a_7中符合txt文件中的个体数据

结果：
1241966 variants loaded from .bim file.
812 people (448 males, 364 females) loaded from .fam.
–remove: 792 people remaining.
1241966 variants and 792 people pass filters and QC.

het_fail_ind.txt文件内容
FID IID
1 067_S_0056
3 067_S_0059
44 072_S_2070
49 128_S_0200
119 041_S_4427
187 029_S_4279
251 031_S_4476
304 052_S_1250
367 109_S_4531
384 013_S_4395
496 073_S_4403
611 109_S_4260
614 014_S_4577
621 109_S_2200
677 035_S_0292
681 127_S_0684
685 029_S_0878
690 023_S_0061
754 031_S_4021
796 037_S_0539
去掉了双引号，以及只保留了第一列和第二列。

最后还剩792个个体数据。

七、QC第六步：清除可能存在的亲缘样本

1.relatedness

名词解释：

表明一对样本在遗传基因上的关联性。常规的WGAS受试者都是无关的，即不存在亲缘关系。隐含的相关性可能会干扰关联分析。如果，有基于家庭的样本（例，父母，后代），则不需要删除相关对，但是在统计分析应该考虑家庭之间的相关性。
IBD全称Identity By Descent, 又叫做血缘同源，指的是两个个体中共有的等位基因来源于共同祖先
IBS全称Identity By State, 又叫做状态同源，指的是两个个体中共有的等位基因序列相同。
对于某个等位基因，IBS state只要求allel的个数相同即可，而IBD state则进一步要求相同的allel来自于共同祖先。
来自https://www.cnblogs.com/esctrionsit/p/13415050.html

作用：

进行样本间的亲缘关系估计，并将清除未标记明显亲子关系的个体亲缘关系过近的个体。

运用：

通过descent（IBD）计算所有样品对的同一性。设置阈值pihat>0.2并创建关联性高于所选阈值的个人列表。

2.分析的数据集之间的隐含关联性

1.分析大于阈值0.2的个体

代码：
plink ‐‐bfile HapMap_3_r3_10 ‐‐extract indepSNP.prune.in ‐‐genome ‐‐min 0.2 ‐‐out pihat_min0.2

注释：
–genome –min 0.2 排查阈值大于0.2的个体。

结果：
1241966 variants loaded from .bim file.
792 people (438 males, 354 females) loaded from .fam.
–extract: 149390 variants remaining.

149390 variants and 792 people pass filters and QC.
Note: No phenotypes present.
Excluding 3625 variants on non-autosomes from IBD calculation.
IBD calculations complete.
Finished writing pihat_min0.2.genome .
这里没有计算非常染色体的SNP。

生成了三个文件pihat_min0.2.genome，pihat_min0.2.hh，pihat_min0.2.log

2.查看文件

cat pihat_min0.2.genome

结果：
FID1 IID1 FID2 IID2 RT EZ Z0 Z1 Z2 PI_HAT PHE DST PPC RATIO
66 024_S_2239 203 024_S_4084 UN NA 0.2942 0.5499 0.1559 0.4308 -1 0.835167 1.0000 8.2500

196 137_S_4466 214 021_S_0159 UN NA 0.1903 0.5921 0.2175 0.5136 -1 0.855272 1.0000 14.7976

215 023_S_0058 719 023_S_4035 UN NA 0.2736 0.4905 0.2359 0.4811 -1 0.850079 1.0000 9.5729

216 031_S_4032 704 031_S_4203 UN NA 0.2059 0.5033 0.2908 0.5424 -1 0.865546 1.0000 13.0035

注释：
来源https://www.cnblogs.com/esctrionsit/p/13415050.html
RT Relationship type inferred from .fam/.ped file，从fam/ped文件推断的样本关系
EZ IBD sharing expected value, based on just .fam/.ped relationship
Z0 P(IBD=0)
Z1 P(IBD=1)

Z2 P(IBD=2)
PI_HAT Proportion IBD, i.e. P(IBD=2) + 0.5P(IBD=1)，IBD比例
PHE Pairwise phenotypic code (1, 0, -1 = case-case, case-ctrl, and ctrl-ctrl pairs, respectively)
DST IBS distance, i.e. (IBS2 + 0.5IBS1) / (IBS0 + IBS1 + IBS2)

PPC IBS binomial test
RATIO HETHET : IBS0 SNP ratio (expected value 2)

RT取值：
列名 含义
FS full siblings，全同胞，指有着相同父母的个体之间的关系
HS half siblings，半同胞，指同父异母或同母异父的个体之间的关系
PO parent-offspring，亲子代关系

OT other，其他
UN 不相关的个体

3.可视化亲缘关系

awk ‘{ if ($8 >0.9) print $0 }’ pihat_min0.2.genome>zoom_pihat.genome
#如果第8列即Z1>0.9，输出该行
cat zoom_pihat.genome
Rscript –no-save Relatedness.R

查看Relatedness.R

pdf("relatedness.pdf")
relatedness = read.table("pihat_min0.2.genome", header=T)
par(pch=16, cex=1)
#par函数用于设定或询问绘图参数。
#参数设定可通过par(参数名 = 取值)或par(赋值参数列表)的形式进行
#pch：指定绘制点所使用的符号，取值范围[0, 24]，其中4是“差号”，20是“点”
#cex：指定符号的大小。cex是一个数值，表示pch的倍数，默认是1.5倍with(relatedness,plot(Z0,Z1,main="pihat_min0.2.genome", xlim=c(0,1), ylim=c(0,1), type="n"))
#(文件，画图plot(x,y，x轴范围，y轴范围，变量类型))这里，我加了标题。mainwith(subset(relatedness,RT=="PO") , points(Z0,Z1,col=3 #col=4))
with(subset(relatedness,RT=="UN") , points(Z0,Z1,col=4 #col=3))
#这里脚本应该是出错了，两个颜色反了，pihat_min0.2.genome文件里只有4个UN,所以我修改了
#points,点的位置
legend(1,1, xjust=1, yjust=1, legend=levels(relatedness$RT), pch=16, col=c(4,3))
#绘制图例legend()，(legend:字符或表达式向量)
#(定位图例,xjust,yjust)
#图例位置由x，y决定。legend 默认坐标为左上角的坐标（因为设定xjust=0，yjust=1），
#因此可以通过xjust，yjust来调节。xjust=1，yjust=0，即为右下角坐标。
#xjust=1，yjust=1，即为右上角坐标。pdf("zoom_relatedness.pdf")
relatedness_zoom = read.table("zoom_pihat.genome", header=T)
par(pch=16, cex=1)
with(relatedness_zoom,plot(Z0,Z1, main="zoom_pihat.genome"，xlim=c(0,0.02), ylim=c(0.98,1), type="n"))
with(subset(relatedness_zoom,RT=="PO") , points(Z0,Z1,col=3))
with(subset(relatedness_zoom,RT=="UN") , points(Z0,Z1,col=4))
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() 

结果：
FID1 IID1 FID2 IID2 RT EZ Z0 Z1 Z2 PI_HAT PHE DST PPC RATIO

无有亲缘关系样本对
原因，我的数据集是ADNI的SNP数据。

生成了三个pdf文件

3.删除所有关联性

1.过滤非初代个体

为了证明大多数相关性是由于父母后代所致，我们仅包括初代（数据集中没有父母的个人）。
代码：
plink –bfile a_8 –filter-founders –make-bed –out a_9

结果：
1241966 variants loaded from .bim file.
792 people (438 males, 354 females) loaded from .fam.
0 people removed due to founder status (–filter-founders)

1241966 variants and 792 people pass filters and QC.
Note: No phenotypes present.
–make-bed to a_9.bed + a_9.bim + a_9.fam … done.

2.再次寻找pihat>0.2的个体
plink –bfile a_9 –extract indepSNP.prune.in –genome –min 0.2 –out pihat_min0.2_in_founders
cat pihat_min0.2_in_founders.genome

结果：
1241966 variants loaded from .bim file.
792 people (438 males, 354 females) loaded from .fam.
–extract: 149390 variants remaining.

149390 variants and 792 people pass filters and QC.
Note: No phenotypes present.
Excluding 3625 variants on non-autosomes from IBD calculation.
IBD calculations complete.
Finished writing pihat_min0.2_in_founders.genome .

pihat_min0.2_in_founders.genome
FID1 IID1 FID2 IID2 RT EZ Z0 Z1 Z2 PI_HAT PHE DST PPC RATIO
66 024_S_2239 203 024_S_4084 UN NA 0.2942 0.5499 0.1559 0.4308 -1 0.835167 1.0000 8.2500

196 137_S_4466 214 021_S_0159 UN NA 0.1903 0.5921 0.2175 0.5136 -1 0.855272 1.0000 14.7976

215 023_S_0058 719 023_S_4035 UN NA 0.2736 0.4905 0.2359 0.4811 -1 0.850079 1.0000 9.5729

216 031_S_4032 704 031_S_4203 UN NA 0.2059 0.5033 0.2908 0.5424 -1 0.865546 1.0000 13.0035

共四个个体

3.仅含初代算缺失率
plink –bfile a_9 –missing
cat plink.imiss
head -10 plink.imiss
wc -l plink.imiss

结果：
1241966 variants loaded from .bim file.
792 people (438 males, 354 females) loaded from .fam.

–missing: Sample missing data report written to plink.imiss, and variant-based
missing data report written to plink.lmiss.

cat:
814 073_S_4259 Y 514 1241712 0.0004139
815 037_S_4302 Y 411 1241966 0.0003309
816 153_S_2109 Y 1109 1241966 0.0008929
817 082_S_4244 Y 471 1241966 0.0003792
818 135_S_4356 Y 381 1241966 0.0003068

head:
FID IID MISS_PHENO N_MISS N_GENO F_MISS
2 123_S_0108 Y 2170 1241966 0.001747
4 024_S_0985 Y 822 1241966 0.0006619
5 109_S_2237 Y 1705 1241966 0.001373
6 131_S_0123 Y 993 1241966 0.0007995
7 098_S_0160 Y 3966 1241966 0.003193
8 027_S_0256 Y 2108 1241966 0.001697
9 116_S_1243 Y 935 1241966 0.0007528
10 073_S_2264 Y 1914 1241966 0.001541
11 094_S_2216 Y 495 1241966 0.0003986

wc:
793 plink.imiss
共792个个体数据

4.手动删除响应率高的相关性个体

由于这一部分，我和教程使用数据集不一样，所以，我就不拿这个数据集做此次步骤，下面是我用教程里的数据集做的步骤。

使用UNIX文本编辑器（例如vi（m））检查“相关对”中哪个响应率最高。 ＃生成Pihat大于0.2的个人的FID和IID列
表，检查响应率较低的个体。
＃在我们的数据集中，单个13291 NA07045的响应率较低。

vi 0.2_low_call_rate_pihat.txt

该命令是使用vi 创建并书写txt文件。
vi创建txt文档后，我们进入输入I进入input状态，

然后输入自己的内容。
输入完毕按键盘上的esc文件退出input状态
在输入：x退出vi
i
13291 NA07045
Press esc on keyboard!
按键盘上的esc键
然后输入
：x
如果有多个“相关”对，可以使用与我们单个的“相关”对相同的方法扩展上面生成的列表。
＃删除pihat> 0.2的“相关”对中响应率最低的个人
plink ‐‐bfile HapMap_3_r3_11 ‐‐remove 0.2_low_call_rate_pihat.txt ‐‐make‐bed ‐‐out HapMap_3_r3_12

到此第一个教程就结束了。