基因组选择中如何清洗基因组数据

原创

育种数据分析之放飞自我 2022-02-16 16:24:30 ©著作权

文章标签 数据缺失率 r语言 文章分类 代码人生

©著作权归作者所有：来自51CTO博客作者育种数据分析之放飞自我的原创作品，请联系作者获取转载授权，否则将追究法律责任

打油诗

题目：plink666

不用不知道

一用忘不掉

朝思暮也想

真的好奇妙

1. 背景

一个朋友之前问过这个问题，问题可以分为：

基因组选择中如何将ATCG的数据转化为012的形式
如何进行基因组数据的过滤筛选
有没有示例代码可以演示

今天我将用一份模拟的芯片下机数据，演示一下如何进行基因组数据的筛选。

2. 数据筛选的几个标准

1，去除缺失率大于10%的SNP位点

这句话的意思是，比如有50K个SNP位点，有1000个体，如果某一个SNP位点在1000个体中，缺失率达到10%以上，即有超过100的个体在这个SNP位点上没有分型结果，就把这个位点删除。因为缺失率大的SNP，后期填充准确性不高。

样本检出率大于90%

这里的样本检出率，比如共有50K的芯片，需要个体的芯片分型数要大于45K，否则去掉该样本

3，最小等位基因频率（MAF）大于1%

我们在进行基因组选择（SSGBLUP）时，需要构建G矩阵和H矩阵，构建G矩阵是利用SNP的分型信息。构建好的G矩阵描述的是个体间的关系，如果一个位点在999个个体上都为A，在1个个体上为G，即MAF为0.1%，这种SNP位点在计算中意义不大，去掉可以节约计算资源。

去除哈迪温伯格平衡检验(HWE)P值小于10的-6次方的SNP

位点进行哈温平衡检验时，去除P值小于10的-6次方的SNP

3. PLINK格式

3.1 map文件

格式说明链接: http://zzz.bwh.harvard.edu/plink/data.shtml#map

map格式的文件, 图谱信息文件, 主要包括染色体名称, SNP名称，所在染色体的遗传坐标和物理坐标.

map文件没有行头

map文件包括四列: 染色体, SNP名称, SNP位置, 碱基对坐标

染色体编号为数字, 未知为0
SNP名称为字符或数字, 如果不重要, 可以从1编号, 注意要和bed文件SNP列对应
染色体的摩尔位置(如果未知，可以用0表示)
SNP物理坐标

如果只有SNP名称, 可以手动构建map文件, 第二列为SNP名称, 其它几列如果未知，可以写为0.

3.2 ped文件

格式说明链接:http://zzz.bwh.harvard.edu/plink/data.shtml#ped

bed格式的文件, 主要包括SNP的信息, 包括个体ID, 系谱信息, 表型和SNP的分型信息.

数据没有行头, 空格或者tab隔开的文件

前六列必须有, 包括系谱信息, 表型信息

第一列: Family ID # 如果没有, 可以用个体ID代替
第二列: Individual ID # 个体ID编号
第三列: Paternal ID # 父本编号
第四列: Maternal ID # 母本编号
第五列: Sex (1=male; 2=female; other=unknown) # 性别, 如果未知, 用0表示
第六列: Phenotype # 表型数据, 如果未知, 用0表示
第七列以后: 为SNP分型数据, 可以是AT CG或11 12, 或者A T C G或1 1 2 2

前六列, 必须要有, 如果没有相关数据, 用0表示.

4. PLINK过滤命令

4.1 缺失大于10%

--geno 0.1

4.2 样本call rate大于90%

--mind 0.1

4.3 最小等位基因频率大于0.01

--maf 0.01

4.4 哈温平衡P值大于1e-6

--hwe 1e-6

4.5 输出过滤后的012形式的文件

--recodeA

4.6 将过滤后的文件命名为c

--out c

5. 测试数据

测试数据为QMSim生成的模拟数据，经过处理之后生成b.map和b.ped文件，处理流程如下

5.1 预览文件

(base) [dengfei@localhost test]$ ls
b.map  b.ped

5.2 使用plink进行处理

plink --file b --chr-set 30 --geno 0.1 --mind 0.1 --maf 0.01 --hwe 1e-6 --recodeA --out c

5.3 处理结果日志

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
Note: --recodeA flag deprecated.  Use 'recode A ...'.
Logging to c.log.
Options in effect:
  --chr-set 30
  --file b
  --geno 0.1
  --hwe 1e-6
  --maf 0.01
  --mind 0.1
  --out c
  --recode A

63985 MB RAM detected; reserving 31992 MB for main workspace.
.ped scan complete (for binary autoconversion).
Performing single-pass .bed write (50000 variants, 1200 samples).
--file: c-temporary.bed + c-temporary.bim + c-temporary.fam written.
50000 variants loaded from .bim file.
1200 samples (0 males, 0 females, 1200 ambiguous) loaded from .fam.
Ambiguous sex IDs written to c.nosex .
0 samples removed due to missing genotype data (--mind).
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1200 founders and 0 nonfounders present.
Calculating allele frequencies... done.
0 variants removed due to missing genotype data (--geno).
--hwe: 2 variants removed due to Hardy-Weinberg exact test.
0 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
49998 variants and 1200 samples pass filters and QC.
Note: No phenotypes present.
--recode A to c.raw ... done.