2021-03-26 QC1.0基于缺失率的过滤
https://blog.csdn.net/qq_40605470/article/details/113183901
虽然我的起始文件是vcf不是bfile，但是转换了格式就都一样
目前我的目录里是这些
[lyc@200server ~]$ ls
eee.log plink prettify rice.nosex
eee.map plink2 rice.bed Rice.recode.vcf
eee.nosex plink2_linux_x86_64_20210302.zip rice.bim test
eee.ped plink2.log rice.fam
LICENSE plink_linux_x86_64_20201019.zip rice.log
文本plink数据包含两个文件：
包含有关个体及其基因型的信息(.ped)；
包含有关遗传标记的信息(.map)
二进制plink数据由三个文件组成：
包含单个标识符（ID）和基因型(.bed)
包含有关个体(.fam)
遗传标记(.bim)
map
[lyc@200server ~]$ head -n 5 eee.map
1 Chr1_1203_T_C 0 1203
1 Chr1_1249_A_C 0 1249
1 Chr1_1266_G_A 0 1266
1 Chr1_1277_T_C 0 1277
1 Chr1_1325_C_T 0 1325
ped
[lyc@200server ~]$ head -n 5 eee.ped
PB-362 PB-362 0 0 0 -9 C C C C A A C C 0 0 0 0 G G G G C C A A T T T T T T A A A A T T T T A A A A T T A A G G G G G G T T C C A A T T C C G G T T T T T T
[lyc@200server ~]$ wc -l eee.map eee.ped
7186300 eee.map
141 eee.ped
7186441 总用量
可以看出共有141个基因型个体，共有7186300个SNP数据。
如果每次都不知道plink的命令
bash: plink: 未找到命令...
相似命令是： 'link'
plink前面加上./就行，好像加地址也可以
我想要查看个体缺失的位点数，每个SNP缺失的个体数
发现bfile的用法,也就是只用前缀名称的意思
(--bfile expects a filename *prefix*;
'.bed', '.bim', and '.fam' are automatically appended.)
[lyc@200server ~]$ ./plink --bfile rice --allow-extra-chr --missing
PLINK v1.90b6.21 64-bit (19 Oct 2020) www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang GNU General Public License v3
Logging to plink.log.
Options in effect:
--allow-extra-chr
--bfile rice
--missing
773821 MB RAM detected; reserving 386910 MB for main workspace.
7186300 variants loaded from .bim file.
141 people (0 males, 0 females, 141 ambiguous) loaded from .fam.
Ambiguous sex IDs written to plink.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 141 founders and 0 confounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.878065.
--missing: Sample missing data report written to plink.imiss, and variant-based
missing data report written to plink.lmiss.
大概意思就是说141个体，无性别（本来就是水稻）plink.nosex 写的是不清楚的性别ID，总的基因型频率是0.878065，样本缺失数据报告被写入plink.imiss，基于变异的数据缺失报告写入plink.lmiss
所以我现在文件是
eee.log plink2_linux_x86_64_20210302.zip plink-temporary.bed rice.log
eee.map plink2.log plink-temporary.bim rice.nosex
eee.nosex plink.imiss plink-temporary.fam Rice.recode.vcf
eee.ped plink_linux_x86_64_20201019.zip prettify test
LICENSE plink.lmiss rice.bed
plink plink.log rice.bim
plink2 plink.nosex rice.fam
输出plink.imiss和 plink.lmiss两个文件。
plink.imiss：个体缺失位点的统计
plink.lmiss：单个SNP缺失的个体数
个体缺失位点统计
第一列为家系ID，第二列为个体ID，第三列是否表型缺失，第四列是缺失的SNP个数，第五列总SNP个数，第六列缺失率
[lyc@200server ~]$ head -n 10 plink.imiss
FID IID MISS_PHENO N_MISS N_GENO F_MISS
PB-362 PB-362 Y 759143 7186300 0.1056
PB-363 PB-363 Y 889516 7186300 0.1238
PB-364 PB-364 Y 840674 7186300 0.117
PB-365 PB-365 Y 859383 7186300 0.1196
PB-366 PB-366 Y 745249 7186300 0.1037
PB-367 PB-367 Y 968356 7186300 0.1348
PB-368 PB-368 Y 936545 7186300 0.1303
PB-369 PB-369 Y 869221 7186300 0.121
PB-370 PB-370 Y 806860 7186300 0.1123
SNP缺失的个体数
第一列为染色体，第二列为SNP名称，第三列为缺失个数，第四列为总个数，第五列为缺失率。
[lyc@200server ~]$ head -n 10 plink.lmiss
CHR SNP N_MISS N_GENO F_MISS
1 Chr1_1203_T_C 0 141 0
1 Chr1_1249_A_C 0 141 0
1 Chr1_1266_G_A 4 141 0.02837
1 Chr1_1277_T_C 6 141 0.04255
1 Chr1_1325_C_T 17 141 0.1206
1 Chr1_1335_G_T 17 141 0.1206
1 Chr1_1362_G_A 1 141 0.007092
1 Chr1_1411_A_G 0 141 0
1 Chr1_1482_T_C 0 141 0
好吧，缺失结果可视化我没搞出来
(题外话有空可以看看王健康老师写的基因定位与育种设计）
在“每个被试”基础上实施QC，后在“每个SNP”基础上进行QC，以最大限度地提高研究中剩余的SNP数。这种方法可防止由于小部分基因分型差的个体而错误地去除某个SNP，但是可能会由于小部分基因分型差的SNP而错误地去除一些个体。
https://www.jianshu.com/p/382a1967cc35
所以说要先搞个体再搞SNP
http://blog.sina.com.cn/s/blog_bcc268080102xdmy.html
1、单个样本的质量控制
（1）样本检出率（sample call rate）:是指对于某个样本而言，通过测序并成功判型的SNPs与所有检出的SNPs的比值，通常的标准应当在80%或90%以上。
（2）杂合性程度（heterozygosity）:这个参数过高即被排除，因为过高的杂合说明样本可能被污染，从而导致杂合基因型数目不相称，通常标准应该控制在23%-30%之间。
2、单核苷酸多态性的质量控制
（1）SNP检出率（SNP call rate）:指对于某一个SNP位点，被成功检测到的样本与所有样本的比值，一般要求在90%以上。
（2）较小等位基因频率（minor allele frequency, MAF）：对于那些MAF较小的SNPs，能得到的信息量较少，而且目前GWAS对这些SNP的检出效能也不高，通常要求MAF在3%以上。
（3）哈代-温伯格平衡（Hardy-Weinberg equilibrium, HWE）检出，HWE有助于确定那些有明显基因分型错误的SNPs，因此一般要求位点SNP的等位基因频率符合哈代-温伯格平衡。
我在文献中看到说SNP的缺失的过滤，应该先于个体的缺失的过滤，相反的程序可能会导致不必要个体的缺失
❝ 无论是测序还是芯片，得到的基因型数据要进行质控，而对缺失数据进行筛选，可以去掉低质量的数据。如果一个个体，共有50万SNP数据，发现20%的SNP数据（10万）都缺失，那这个个体我们认为质量不合格，如果加入分析中可能会对结果产生负面的影响，所以我们可以把它删除。同样的道理，如果某个SNP，在500个样本中，缺失率为20%（即该SNP在100个个体中都没有分型结果），我们也可以认为该SNP质量较差，将去删除。当然，这里的20%是过滤标准，可以改变质控标准。下文中的质控标准是2%。
❞
https://blog.csdn.net/weixin_35698867/article/details/112482182这句话就很清楚明了
弄的这一大套东西或许都是一篇文献里的内容，但是英语太差就只能瞎搜搜
A tutorial on conducting genome‐ wide association studies:Quality control and statistical analysis
如果我需要在服务器上连网
su root
输入密码
然后vpn-connect
exit退回自己的账户就可以用了
https://www.jianshu.com/p/8340055d3646?utm_campaign=haruki思路来源
我觉得还是得可视化一下缺失结果，这样能更好的确定阈值选取多少，然后我就打算去摸索这个hist_miss.R的脚本
打不开https://github.com/MareesAT/GWA_tutorial
于是直接上代码
git clone https://github.com/MareesAT/GWA_tutorial.git
然后拥有了
GWA_tutorial
[lyc@200server ~]$ cp /home/lyc/GWA_tutorial/1_QC_GWAS.zip /home/lyc
用这个把a文件复制到b目录下
然后解压缩
[lyc@200server ~]$ unzip 1_QC_GWAS.zip
然后发下还是不能用那个脚本
image.png
image.png
cp /home/lyc/1_QC_GWAS/hist_miss.R /home/lyc
这回应该可以了
