欢迎关注"生信修炼手册"!
本篇文章按照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
这一步会生成以下三个文件
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.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
以上只是一个最基本的示例,更多用法和细节可以参考官方的文档。
扫描关注微信号,更多精彩内容等着你!