使用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.phehapmap1.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这一步会生成以下三个文件
- hapmap1.bed
- hapmap1.bim
- 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.003939miss_stat.lmiss ： 每个突变位点基因型缺失的频率，部分内容如下
CHR SNP N_MISS N_GENO F_MISS
1 rs6681049 0 89 0
1 rs4074137 0 89 0
1 rs7540009 0 89 0miss_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_stat6. 关联分析
进行疾病和突变位点基因型之间的关联分析，命令如下
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以上只是一个最基本的示例，更多用法和细节可以参考官方的文档。
腾讯云开发者
