[lyc@200server ~]$ Rscript --no-save hist_miss.R
pdf 
  2 

接上次，我有了用R画图的可视化的缺失情况
结果在命令行的linux下看pdf产生了困难，或许我可以把它下到我的电脑上

[lyc@200server ~]$ evince histimiss.pdf
Cannot parse arguments: 无法打开显示：

https://www.jianshu.com/p/bc35bc9885f4

scp 用户名@服务器地址:文件路径 本地路径
这是从服务器下载文件。
如果从本地上传至服务器就把两个交换一下就可以了。
scp 本地文件路径 用户名@服务器地址:文件路径

天哪，这个相当于需要mac的终端去写命令然后下载服务器上的文件，而不是在服务器上把文件放去本地mac

初筛的情况

然后我先直接筛一下吧，用一条命令把那些参数都加进去

[lyc@200server ~]$ ./plink --bfile rice --geno 0.2 --maf 0.05 --mind 0.2 --make-bed --out clean1

从这里看出的我的过滤情况

7186300 variants loaded from .bim file.
141 people (0 males, 0 females, 141 ambiguous) loaded from .fam.
Ambiguous sex IDs written to clean1.nosex .
0 people removed due to missing genotype data (--mind).
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.
1467225 variants removed due to missing genotype data (--geno).
2966184 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
2752891 variants and 141 people pass filters and QC.
Note: No phenotypes present.
--make-bed to clean1.bed + clean1.bim + clean1.fam ... done.

将其转为map和ped形式，方便查看

[lyc@200server ~]$ ./plink --bfile clean1 --allow-extra-chr --recode --out clean1test

这个时候基因型频率变高了

2752891 variants loaded from .bim file.
141 people (0 males, 0 females, 141 ambiguous) loaded from .fam.
Ambiguous sex IDs written to clean1test.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.976998.
2752891 variants and 141 people pass filters and QC.
Note: No phenotypes present.
--recode ped to clean1test.ped + clean1test.map ... done.

统计剩下的个体和SNP数量

[lyc@200server ~]$ wc -l clean1test.map clean1test.ped
   2752891 clean1test.map
       141 clean1test.ped
   2753032 总用量

天哪，被 筛掉了近三分之二
等会会下载pdf了再针对性的设置一下参数
水稻一般都是自花授粉,这筛掉的太多，我都不敢hwe了

image.png

这个图本意代表，个体缺失率在多少的SNP频率是多少，大于某值的部分可以删掉，数字越大，留下的越多，可以看出设定0.2能留下大部分，如果是0.02就基本给删没了。

image.png
这个图本意代表snp缺失率在多少的个体频率是多少，大于某值的部分可以删掉，数字越大留下的越多。可以看出设定0.2，能留下大部分个体，设定0.02也是大部分

考虑其他阈值可能

但是，我这些图也太奇怪了，而且过滤的太多了，而且这句话的意思是mind没有删掉一个，geno删掉的很多，maf删掉的更多，或许可以mind来个0.15，geno来个0.4（我是指看图瞎说的，maf一般都是1%－5%，越小越宽松那我试试宽松的1％
但现在又来了一个问题，到底剩下多少才算这个质控比较合格呢，我是觉得有三分之二的位点被删掉了，貌似也太严格了，或者说，如果这就是标准，只能说，实验数据太差了，所以才会被删那么多，所以maf的值也不能放太小

image.png

所以要不mind来个0.15，geno继续保持0.2，maf搞个宽松一点的比如0.3／0.4
让我们先来看一下maf的分布

image.png
怎么得到的，和前面的方法是一样的

[lyc@200server ~]$ ls 1_QC_GWAS
1_Main_script_QC_GWAS.txt    HapMap_3_r3_1.bed  heterozygosity_outliers_list.R  inversion.txt     Relatedness.R
check_heterozygosity_rate.R  HapMap_3_r3_1.bim  hist_miss.R                     MAF_check.R
gender_check.R               HapMap_3_r3_1.fam  hwe.R                           pops_HapMap_3_r3
[lyc@200server ~]$ cp /home/lyc/1_QC_GWAS/MAF_check.R /home/lyc
[lyc@200server ~]$ ./plink --bfile rice --allow-extra-chr --freq --out MAF_check
--freq: Allele frequencies (founders only) written to MAF_check.frq .
[lyc@200server ~]$ head MAF_check.frq
 CHR                SNP   A1   A2          MAF  NCHROBS
   1      Chr1_1203_T_C    C    T       0.2837      282
   1      Chr1_1249_A_C    C    A       0.2766      282
   1      Chr1_1266_G_A    A    G       0.2664      274
   1      Chr1_1277_T_C    C    T       0.2519      270
   1      Chr1_1325_C_T    T    C       0.1653      248
   1      Chr1_1335_G_T    T    G       0.1653      248
   1      Chr1_1362_G_A    A    G       0.1536      280
   1      Chr1_1411_A_G    G    A       0.3085      282
   1      Chr1_1482_T_C    C    T       0.3191      282
[lyc@200server ~]$ Rscript --no-save MAF_check.R
null device 
          1 

从这个图看出，确实他的大部分基因频率为在0.06以下，没有分型，所以被删掉的自然就很多了，或者说改成0.04

[lyc@200server ~]$ ./plink --bfile rice --allow-extra-chr --geno 0.2 --maf 0.04 --mind 0.15 --make-bed --out clean2

就怎么说呢，还是挺严厉的

都没加别的条件，加了怕是都没有数据供我分析了，无语子
先问问师兄，需不需要继续筛了
看看这，大多数都给整没了

image.png
image.png