「主要介绍」
❝
--geno筛选个体；--mind筛选SNP ❞
GWAS分析时，拿到基因型数据，拿到表型数据，要首先做以下几点：
- 1，查看自己的表型数据，是否有问题
- 2，查看自己的基因型数据，是否有问题
然后再进行建模，得到显著性SNP以及可视化结果。
清洗数据的时间占80%的时间，有句话这样讲：“Garbage in, Garbage out（垃圾进，垃圾出）”，所以清洗数据非常重要，今天学习一下基因组数据如何清洗。
「为何对缺失数据进行筛选？」
❝无论是测序还是芯片，得到的基因型数据要进行质控，而对缺失数据进行筛选，可以去掉低质量的数据。如果一个个体，共有50万SNP数据，发现20%的SNP数据（10万）都缺失，那这个个体我们认为质量不合格，如果加入分析中可能会对结果产生负面的影响，所以我们可以把它删除。同样的道理，如果某个SNP，在500个样本中，缺失率为20%（即该SNP在100个个体中都没有分型结果），我们也可以认为该SNP质量较差，将去删除。当然，这里的20%是过滤标准，可以改变质控标准。下文中的质控标准是2%。 ❞
1. plink数据格式转化
数据使用上一篇的数据，因为数据是plink的bfile格式，二进制不方便查看，我们将其转化为文本map和ped的格式。
代码语言：javascript复制plink --bfile HapMap_3_r3_1 --recode --out test
结果生成：test.map test.ped
2. 查看基因型个体和SNP数量
代码语言：javascript复制wc -l test.map test.ped
可以看出，共有165个基因型个体，共有1447897个SNP数据。
「预览一下ped文件：」
「预览一下map文件：」
3. 查看一下个体缺失的位点数，每个SNP缺失的个体数
看一下描述：
代码语言：javascript复制--missing: Sample missing data report written to plink.imiss, and variant-based
missing data report written to plink.lmiss.
结果生成两个文件，分别是一个个体ID上SNP缺失的信息，另一个是每个SNP在个体ID中缺失的信息。
- 个体缺失位点的统计在
plink.imiss中 - 单个SNP缺失的个体数在
plink.lmiss.中
「个体缺失位点统计预览：」第一列为家系ID，第二列为个体ID，第三列是否表型缺失，第四列缺失的SNP个数，第五列总SNP个数，第六列缺失率。
「SNP缺失的个体数文件预览：」第一列为染色体，第二列为SNP名称，第三列为缺失个数，第四列为总个数，第五列为缺失率
「R语言做直方图」
代码的意思是读取这两个文件，然后用频率的那一列作图，将图保存为pdf输出。
代码语言：javascript复制indmiss<-read.table(file="plink.imiss", header=TRUE)
snpmiss<-read.table(file="plink.lmiss", header=TRUE)
# read data into R
pdf("histimiss.pdf") #indicates pdf format and gives title to file
hist(indmiss[,6],main="Histogram individual missingness") #selects column 6, names header of file
pdf("histlmiss.pdf")
hist(snpmiss[,5],main="Histogram SNP missingness")
dev.off() # shuts down the current device
因为，当前目录已经有hist_miss.R脚本，所以可以直接在目录下运行下面代码，会生成两个pdf文件。
Rscript hist_miss.R
「单个位点缺失率统计：」
「基因型个体缺失率统计：」
4. 对个体及SNP缺失率进行筛选
- 1， 如果一个SNP在个体中2%都是缺失的，那么就删掉该SNP，参数为：
--mind 0.02 - 2，如果一个个体，有2%的SNP都是缺失的，那么就删掉该个体，参数为：
--geno 0.02
4. 1 对个体缺失率进行筛选
「先过滤个体缺失率高于2%的SNP」
代码语言：javascript复制 plink --bfile HapMap_3_r3_1 --geno 0.02 --make-bed --out HapMap_3_r3_2
「转化为map和ped的形式：」
代码语言：javascript复制plink --bfile HapMap_3_r3_2 --recode --out test
查看一下过滤后的行数，
「之前的为：」
代码语言：javascript复制 1457897 test.map
165 test.ped
「现在的为：」
代码语言：javascript复制 1430443 test.map
165 test.ped
可以看出，过滤了2万多个位点。
从当时的log日志里也可以看出这一点：
代码语言：javascript复制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
Logging to HapMap_3_r3_2.log.
Options in effect:
--bfile HapMap_3_r3_1
--geno 0.02
--make-bed
--out HapMap_3_r3_2
515185 MB RAM detected; reserving 257592 MB for main workspace.
1457897 variants loaded from .bim file.
165 people (80 males, 85 females) loaded from .fam.
112 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 53 nonfounders present.
Calculating allele frequencies... done.
Warning: 225 het. haploid genotypes present (see HapMap_3_r3_2.hh ); many
commands treat these as missing.
Total genotyping rate is 0.997378.
27454 variants removed due to missing genotype data (--geno).
1430443 variants and 165 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls. (53 phenotypes
are missing.)
--make-bed to HapMap_3_r3_2.bed HapMap_3_r3_2.bim HapMap_3_r3_2.fam ...
done.
可以看到--geno，过滤了27454个位点。
4. 2 对SNP缺失率进行筛选
「过滤SNP缺失率高于2%的个体」
代码语言：javascript复制plink --bfile HapMap_3_r3_2 --mind 0.02 --make-bed --out HapMap_3_r3_3
「查看日志：」
代码语言：javascript复制Options in effect:
--bfile HapMap_3_r3_2
--make-bed
--mind 0.02
--out HapMap_3_r3_3
515185 MB RAM detected; reserving 257592 MB for main workspace.
1430443 variants loaded from .bim file.
165 people (80 males, 85 females) loaded from .fam.
112 phenotype values loaded from .fam.
0 people removed due to missing genotype data (--mind).
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 53 nonfounders present.
Calculating allele frequencies... done.
Warning: 179 het. haploid genotypes present (see HapMap_3_r3_3.hh ); many
commands treat these as missing.
Total genotyping rate is 0.997899.
1430443 variants and 165 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls. (53 phenotypes
are missing.)
--make-bed to HapMap_3_r3_3.bed HapMap_3_r3_3.bim HapMap_3_r3_3.fam ...
done.
「没有过滤掉个体，剩余：」个体：165 SNP：1430443
5 同时对个体和SNP的缺失率进行筛选
「两步合在一起，即过滤位点，又过滤个体：」
代码语言：javascript复制plink --bfile HapMap_3_r3_1 --geno 0.02 --mind 0.02 --make-bed --out HapMap_3_r3_5
plink日志：
代码语言：javascript复制Options in effect:
--bfile HapMap_3_r3_1
--geno 0.02
--make-bed
--mind 0.02
--out HapMap_3_r3_5
515185 MB RAM detected; reserving 257592 MB for main workspace.
1457897 variants loaded from .bim file.
165 people (80 males, 85 females) loaded from .fam.
112 phenotype values loaded from .fam.
1 person removed due to missing genotype data (--mind).
ID written to HapMap_3_r3_5.irem .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 52 nonfounders present.
Calculating allele frequencies... done.
Warning: 225 het. haploid genotypes present (see HapMap_3_r3_5.hh ); many
commands treat these as missing.
Total genotyping rate in remaining samples is 0.997486.
26686 variants removed due to missing genotype data (--geno).
1431211 variants and 164 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls. (52 phenotypes
are missing.)
--make-bed to HapMap_3_r3_5.bed HapMap_3_r3_5.bim HapMap_3_r3_5.fam ...
done.
可以看出，两者最终结果是一样的。
