使用plink进行case/control关联分析

欢迎关注"生信修炼手册"！

本篇文章按照plink官方提供的教程，进行一个实际操作。可以看做是官方教程的一个翻译版本。官方教程的链接如下

http://zzz.bwh.harvard.edu/plink/tutorial.shtml

1.  下载测试数据

wget http://zzz.bwh.harvard.edu/plink/hapmap1.zip
unzip hapmap1.zip

文件列表如下

├── hapmap1.map
├── hapmap1.ped
├── pop.phe
└── qt.phe

hapmap1.map和hapmap1.ped 是plink需要的两个基础文件，所有样本分成了case和control两组，表型信息是一种疾病的状态。pop.phe和qt.phe 分别表示样本的人群分布和数量性状。

2. 查看输入文件的基本信息

plink运行时，会联网检查软件是否是最新版，如果不想进行这一操作，可以添加--noweb选项。plink 需要两个输入文件，分别为.ped和.map格式。当这两个文件名称相同时，直接用--file参数指定文件名就可以了。如果不同，也可以用--ped和--map分别指定这两个文件。在不加其他参数的情况下，可以查看输入文件的基本信息，包括样本个数和突变位点个数等信息。

命令如下

plink --file hapmap1 --noweb

需要注意的是，plink默认情况下，会对输入数据进行过滤，主要是过滤突变位点和和样本。主要包括以下几个参数
--mind : 对样本进行过滤，去除缺失基因型频率大于给定阈值的样本
--maf： 对SNP位点进行过滤，去除MAF小于给定阈值的SNP位点
--geno :   对SNP位点进行过滤， 去除缺失基因型频率大于给定阈值的SNP位点
--hwe :  对SNP位点进行过滤， 去除不符合哈温伯格平衡的SNP位点。
--me :   只对家系数据，根据Mendel error rate进行过滤。

3. 创建二进制的PED 文件

ped文件的大小随着样本个数和突变位点个数的增多，会变得非常大。此时可以转换成二进制的文件，加快处理速度。命令如下

plink --file hapmap1 --make-bed --out hapmap1 --noweb

这一步会生成以下三个文件

1. hapmap1.bed

2. hapmap1.bim

3. hapmap1.fam

这里的bed指的就是binary  ped, bim和fam是纯文本文件。替换成二进制之后，原始的ped和map中的信息，用bed, bim, fam三个文件进行存储。

4. 加载二进制文件

加载二进制的ped文件，对应的参数需要替换成--bfile,命令如下

plink --bfile hapmap1 --noweb

和第一步相比，只是换成了二进制文件，处理速度更快，除此之外，没有其他不同的地方。

5. 统计基本信息

分别统计样本和突变位点基因型缺失的比例，命令如下

plink --bfile hapmap1 --missing --out miss_stat --noweb

这一步会生成以下3个文件

miss_stat.imiss ： 每个样本基因型缺失的频率，部分内容如下

FID  IID MISS_PHENO   N_MISS   N_GENO   F_MISS
  HCB181    1          N      671    83534 0.008033
  HCB182    1          N     1156    83534  0.01384
  HCB183    1          N      498    83534 0.005962
  HCB184    1          N      412    83534 0.004932
  HCB185    1          N      329    83534 0.003939

miss_stat.lmiss ： 每个突变位点基因型缺失的频率，部分内容如下

CHR         SNP   N_MISS   N_GENO   F_MISS
   1   rs6681049        0       89        0
   1   rs4074137        0       89        0
   1   rs7540009        0       89        0

miss_stat.log ： log 文件

当突变位点较多时，可以按照染色体分别统计，命令如下

plink --bfile hapmap1 --chr 1 --out res1 --missing

统计突变位点的MAF, 命令如下

plink --bfile hapmap1 --freq --out freq_stat --noweb

生成的文件内容如下

CHR         SNP   A1   A2          MAF  NCHROBS
   1   rs6681049    1    2       0.2135      178
   1   rs4074137    1    2      0.07865      178
   1   rs7540009    0    2            0      178
   1   rs1891905    1    2       0.4045      178
   1   rs9729550    1    2       0.1292      178
   1   rs3813196    1    2      0.02809      178

如果样本还要其他的metadata, 比如人群分布信息，也可以按照不同人群分别进行统计，命令如下

plink --bfile hapmap1 --freq --within pop.phe --out freq_stat --noweb

输出文件内容如下

CHR         SNP     CLST   A1   A2      MAF    MAC  NCHROBS
   1   rs6681049        1    1    2   0.2333     21       90
   1   rs6681049        2    1    2   0.1932     17       88
   1   rs4074137        1    1    2      0.1      9       90
   1   rs4074137        2    1    2  0.05682      5       88
   1   rs7540009        1    0    2        0      0       90
   1   rs7540009        2    0    2        0      0       88
   1   rs1891905        1    1    2   0.4111     37       90
   1   rs1891905        2    1    2   0.3977     35       88

由于所有样本来自两个不同人群，所以每个突变位点出现了两次，分别代表不同人群中的MAF。

当只关注某个SNP位点时，用法如下

plink --bfile hapmap1 --snp rs1891905 --freq --within pop.phe --out snp1_frq_stat

6. 关联分析

进行疾病和突变位点基因型之间的关联分析，命令如下

plink --bfile hapmap1 --assoc --out as1  --noweb

输出结果如下

CHR SNP BP A1 F_A F_U A2 CHISQ P OR 
1 rs6681049 1 1 0.1591 0.2667 2 3.067 0.07991 0.5203 
1 rs4074137 2 1  0.07955  0.07778 2 0.001919 0.9651 1.025 
1 rs7540009 3 0 0 0 2 NA NA NA 
1 rs1891905 4 1 0.4091 0.4 2 0.01527 0.9017 1.038

理论上来说，p值显著的突变位点就是我们想要寻找的和性状相关的突变位点。为了降低假阳性率，可以对多重假设检验的p值进行校正，命令如下

plink --bfile hapmap1 --assoc --adjust --out as2

以上只是一个最基本的示例，更多用法和细节可以参考官方的文档。

扫描关注微信号，更多精彩内容等着你！