这是我个人使用的plink软件备忘录,将我所有会用到但是又没有完全记住的选项和功能记录在这里,方便查看。
windows中使用plink
-
进入plink官网,下载 plink
-
将 plink.exe 加入到环境变量 (我的做法是建立快捷方式,命名为 plink ,放到 D:\shortcuts 中)
输入输出
输入
compound格式
使用 –compound-genotypes 选项,读取 compound格式的plink文件(一个SNP的两个碱基之间没有空格)
1 | plink --file compound --compound-genotypes --recode --out new |
二进制文件(bed)
查看plink2.0文章的算法部分,发现plink2.0 采用了按位运算,大大提高了计算效率。但是按位运算的前提是需要先把基因型文件转化为二进制格式。因此这里,需要介绍一下plink 生成的bed文件格式:
官方说明:http://zzz.bwh.harvard.edu/plink/binary.shtml
这个说明很详细,看这个就懂了。
首先,plink生成二进制文件时,同时会生成bim和fam文件。因此,plink 读取bed文件到的时候,同时也会读取 bim 和 fam 文件,获取snp和样本数目信息。
一个SNP的碱基信息会转化为2 bits。规则如下:
1 | 00 Homozygote "1"/"1" (first allele in .bim file) |
bed 文件将8个bits 作为一个 block (即一个字节)。一开始的三个block是由plink指定的,从第四个开始才开始记录基因型数据。
前两个 block 称为 magic number ,作用是用于plink软件来区分这个bed是不是真的bed文件。第三个 block 称为 mode ,表明这个bed文件是SNP-major(00000001) 还是 individual-major(0000000),这个意思就是说数据按什么顺序排序,默认都是按照SNP的顺序,即先把bim文件中第一个SNP的所有基因型排列完,然后是第二个SNP,第三个SNP……(一直到最后,如果样本不能整除4,那么最后一个block多余的bits设为0,确保每一个新行都从一个新的block开始)。
1 | |-magic number--| |-mode-| |--genotype data---------| |
这里其相应的 bim 文件内容如下(A1 为等位基因频率较小的碱基 )
1 | 1 snp1 0 1 G A |
fam 文件如下
1 | 1 1 0 0 1 0 |
这里还有一个小问题是,在读取每个字节数据的时候是从后向前读的。如果把一个字节的8个位置从后到前标记伪 A 到 H,举例如下
1 | 01101100 |
第一个SNP的前4个基因型则为 AB, CD,EF,GH 四个位置对应的分型,如下:
1 | 01101100 |
上面例子中的 3 个位点的解析如下,第一个SNP的 A1 和 A2 碱基为 G 和 A ,因此 00 对应 GG,01 对应 GA ,11 对应 GG ,10 对应缺失。由于这里总共是 6 个样本,因此一个 SNP 需要占用 2 个字节,第二个字节多余的位置使用 00 填充。
因此第一个SNP的第一个样本是GG,之后是AA,00, AA, AA, AA。
1 |
|
VCF 文件
联合使用 –vcf filename 和 –const-fid 来读取 vcf 文件,这里 –const-fid 是将所有个体的家系号均设为 0 。
注意 VCF 文件作为 plink 的输入文件,其位点必须严格按照物理位置顺序排序,不然就会遇到下面的报错。
1 | Error: .bim file has a split chromosome. Use --make-bed by itself to |
你需要先用 –make-bed 生成标准的二进制文件(排序好的),再进行其他处理。
使用 –bcf 选项可以读取 BCF2 文件。
输出
二进制格式
使用 --make-bed
输出二进制格式,使用--make-just-bim
和 --make-just-fam
只输出 bim 或 fam 文件。
compound格式
使用 –recode compound-genotypes 选项,生成 compound格式的plink文件(一个SNP的两个碱基之间没有空格)
1 | plink --allow-extra-chr --chr-set 95 --file plink59 --recode compound-genotypes --out plink59_compound |
012格式
使用 --recodeA
选项可以输出 012 格式的基因型文件,默认是按照 A1 (最小等位基因) 进行计数。
可以通过添加 --recode-allele
选项指定每一个 SNP 是计数哪一个碱基,其后面接着输入文件名称,这个输入文件含有两列,第一列是 SNP 名称,第二列是需要计数的碱基。
举例如下
1 | plink --allow-extra-chr --chr-set 95 --file 1 --recode-allele count_allele.txt --recodeA --out 1_count |
其中的 count_allele.txt 文件内容如下
1 | 1_242598 C |
VCF 文件
使用 –recode vcf-iid 输出为 VCF 文件,使用 –recode vcf-iid bgz 输出为 vcf.gz 压缩文件。
这里 vcf-iid 含义是只使用个体号作为 VCF 文件中的样本号,如果仅仅是 --recode vcf,那么生成的VCF文件的样本号是家系号和个体号合并的结果(以_连接) 。
需要注意的是,plink 默认会将等位基因频率高的碱基 (A2) 设置为 ref ,将等位基因频率低的碱基 (A1) 设置为 alt 。如果你需要指定每个位点 的 ref 和 alt ,比如你可以设置与某一个 VCF 中的 ref 和 alt 保持一致,你可以添加一下选项:
1 | --a2-allele ref.vcf 4 3 '#' |
参数说明:
--a2-allele ref.vcf 4 3 '#'
:
--a2-allele
命令表示 指定A2基因(Major)
ref.vcf
是参考群体的vcf文件(必须是未压缩文件)
4
是参考群体A2基因所在列
3
是SNP名称列
#
是需要跳过的行
需要注意的地方如下
- ref.vcf 和 plink 文件中的染色体,snp名称,物理位置,两个碱基必须一致
- 需要转换的 plink 文件中的位点必须均在 ref.vcf 中。
位点顺序
经过 plink 处理之后,位点会按照物理位置排序。
样本顺序
按照特定样本顺序排序
除了merge操作外,plink 不会改变样本顺序,会保持原顺序。
如果想要改变样本顺序,那么可以使用 –indiv-sort <mode name> [filename] 选项,这个选项有 4 种模式,说明如下:
-
'‘none’/‘0’: 保持原顺序。这是 PLINK 1.9 中 除merge 外的所有操作的默认值。
-
‘natural’/‘n’: 按照家系号和个体号排序,这是PLINK 1.9 在merge操作中的默认值。举例如下
‘id2’ < ‘ID3’ < ‘id10’
-
‘ascii’/‘a’: 同样按照家系号和个体号排序,不过是按照 ASCII 码。举例如下
‘ID3’ < ‘id10’ < ‘id2’
-
‘file’/‘f’: 采用另一个文件中的顺序。这个文件前两列是家系号-个体号,空格或tab分隔。
如果想要使用 file 模式来改变样本顺序,先生成二进制文件,再生成正常的plink格式,分两步走。
1 | plink --allow-extra-chr --chr-set 95 --chr 1-18 --file raw --indiv-sort file keep_id.txt --make-bed --out sorted |
merge操作后样本顺序
merge 后得到的 plink 文件,其默认会按照家系号和个体号排序,即 –indiv-sort 中的 natural 模式。
质控
样本质控
除了常规的 call rate 之外,还有杂合子比例,性别检查,PCA分析剔除离群样本,IBS和IBD。
杂合子比例
使用选项 –het 。如果杂合子比例过高,可能是样本污染;如果杂合子比例过低,可能是近交造成的。
建议先根据 LD 质控(–indep-pairwise 50 5 0.2),只使用独立的位点计算样本杂合子比例。
结果文件后缀为 .het
,其内容举例如下,这里 O(HOM) 是样本的纯合子比例,E(HOM) 是期望的纯合子比例,N(NM) 是剔除缺失后的位点总数,根据 plink 官方文档 F 系数计算公式为 (<observed hom. count> - <expected count>) / (<total observations> - <expected count>)
FID | IID | O(HOM) | E(HOM) | N(NM) | F |
---|---|---|---|---|---|
0 | 1 | 30340 | 2.92E+04 | 45590 | 0.0693 |
0 | 2 | 28816 | 2.92E+04 | 45588 | -0.02362 |
因此每个样本的杂合子比例计算公式为: ,有人建议剔除群体杂合子比例均值的 3 倍标准差之外的样本。
性别检查
使用 –check-sex 选项, 使用 X 染色体上的位点检查性别,看 F 值,大于0.8为公,反之为母(一般小于0.2),0.2 到 0.8 之间的个体性别会标记为 “PROBLEM” 。
这个 F 值是根据 X 染色体上的近交系数(纯合子比例)得到的,所以需要事先剔除 X 染色体上的伪常染色体区域的位点。
1 | plink --chr-set 95 --allow-extra-chr --file all_sex_maf --check-sex --out all_sex_maf |
离群个体
使用 --pca
进行 PCA 分析,根据散点图手动剔除离群个体。
IBS 和 IBD
IBS
IBS 就是样本间同态相同的比例,实际计算两个样本间相同碱基的比例。
首先我们先看一个位点的情况,下面的 IBS Count 就是两个个体在这个位点中的共同的碱基数目。
个体1 | 个体2 | IBS State |
---|---|---|
AA | AA | 2 |
AA | Aa | 1 |
AA | aa | 0 |
扩展到多个位点,IBS 计算公式如下:
plink 计算 IBS 利用 –distance ibs 选项,生成 .mibs
为后缀的结果文件,例如
1 | plink --allow-extra-chr --chr-set 95 --file hapmap1 --distance square ibs flat-missing |
其中 square 指输出文件为矩阵格式,flat-missing 就是对应地从分母剔除缺失位点。
IBD
IBS 就是样本间同源相同的比例,其不仅仅是具有相同的等位基因,而且要来自于同一祖先。
通常 IBD 无法观测得到,PLINK 通过隐马尔科夫模型来计算 IBD = 0,1,2 的概率,然后计算 PI_HAT 作为 IBD 估计值,其计算公式类似 IBS ,即 P(IBD=2) + 0.5 P(IBD=1) 。
建议先根据 LD 质控(–indep-pairwise 50 5 0.2),只使用独立的位点计算样本杂合子比例。
在 plink 中使用 –genome 选项来计算 IBD ,还可以通过 –min <minimum PI_HAT value> 和 —max <maximum PI_HAT value> 来移除 PI_HAT 小于或大于某个值的样本对的行(例如 –min 0.2)。
输出结果后缀为 .genome
,其各列内容解释如下,这里主要关注 Z0,Z1,Z2,PI_HAT,DST (这就是 IBS) 。
由于输出文件中包含所有样本两两配对的结果,因此文件会比较大。
1 |
|
三种亲缘关系的IBD测试
路径:/mnt/data/zhouziwen/test/isolate/2021/2021-03-25-plink_test/IBD
同时使用杜长大三个品种的数据(总共样本数目为 3150),其系谱均检验均正确,质控 LD 后剩余 5260 个位点,使用这些位点计算 IBD,查看亲子、全同胞及半同胞三种亲缘关系的样本间的情况 。
按照理论来说,这三种亲缘关系的 Z0,Z1,Z2 和 PI_HAT的期望为
Z0 | Z1 | Z2 | PI_HAT | |
---|---|---|---|---|
亲子 | 0 | 1 | 0 | 0.5 |
全同胞 | 0.25 | 0.5 | 0.25 | 0.5 |
半同胞 | 0.5 | 0.5 | 0 | 0.25 |
下面是亲子的直方图,所有样本 Z0 均接近于 0,Z1 较大,但是 Z2 并不等于 0 。
全同胞的结果如下
半同胞的结果如下
测试结论
-
三种关系中 ,Z2 值均偏大,远高于期望值,因此这里的 IBD 算法不够精确。
-
三种关系中 ,PI_HAT 值均偏大,远高于期望值,估计是收到了 Z2 值偏大的影响。
-
这里最好区分的是亲子关系,其 Z0 值接近于 0,Z1 值均大于 0.6 。但是这里半同胞和全同胞的 PI_HAT 的分布有交叉,无法完全分开。
-
总体来说,IBD 的计算结果只能说勉强符合要求,应该只能用于筛选样本,比如挑选或剔除亲缘关系比较近的样本对。
-
计算 IBD 前必须根据 LD 质控位点,不然得到的结果更加糟糕。
位点质控
杂合子比例
位点杂合子比例可以通过 --hardy
进行计算,输出文件后缀维 .hwe
,内容举例如下。
其中 O(HET) 便是位点的观察杂合度,GENO 的三个数值对应着 A1A1,A1A2,A2A2 三种基因型的数目。
CHR | SNP | TEST | A1 | A2 | GENO | O(HET) | E(HET) | P |
---|---|---|---|---|---|---|---|---|
0 | snp1 | ALL(NP) | T | C | 4337/533/26259 | 0.01712 | 0.252 | 0 |
0 | snp2 | ALL(NP) | G | T | 8035/1174/22550 | 0.03697 | 0.3956 | 0 |
剔除多等位基因位点
使用 –snps-only just-acgt 选项会剔除非 SNP 位点(例如 indel 位点),对于存在三等位基因的 SNP ,会将存在三等位基因的基因型替换为0。
举个例子,输入文件,第一个SNP是一等位基因,第二个SNP是三等位基因,第三个SNP是INDEL 位点。
1 1 0 0 0 0 A A A C I D
2 2 0 0 0 0 A A A T D D
3 3 0 0 0 0 A A T T D D
4 4 0 0 0 0 A A A T D D
输出文件,保留了第一个SNP(即会保留无变异的位点),第二个SNP将存在三等位基因的基因型替换为0,直接剔除了第三个SNP。
1 1 0 0 0 -9 A A 0 0
2 2 0 0 0 -9 A A A T
3 3 0 0 0 -9 A A T T
4 4 0 0 0 -9 A A A T
注意:–make-bed 默认就会执行 –snps-only,因此如果生成二进制文件过程中出现多等位基因的 SNP ,不会报错,只会给出下面的提示:
1 | Warning: Variant 1 triallelic; setting rarest alleles missing. |
LD 过滤
使用 –indep-pairwise <window size>[‘kb’] <step size (variant ct)> <r^2 threshold> 选项来过滤位点,常用参数组合为 50 5 0.5 或 50 5 0.2 ,一般不动前面两个参数,修改第三个的值。
这三个参数说明如下:
- 窗口大小,单位 Kb
- 步长,单位位点数目,表示每次计算后移动窗口的位点数目
- r2 阈值
其它
flip-翻转位点
合并失败,首先检查是不是存在正反链问题(查看bim文件),如果是正反链问题,用 flip 命令进行部分位点的翻转。如果不是,估计就是基因型文件分型出问题了。
1 | -- flip <SNP ID list> |
Given a file containing a list of SNPs with A/C/G/T alleles, –flip swaps A↔T and C↔G. A warning will be given if any alleles are not named A, C, G, or T.
merge 注意事项
使用 –merge 合并时有以下注意事项:
-
合并后的基因型文件位点是所有被合并的基因型文件位点的并集,不存在的位点基因型会设为缺失。
-
如果多个被合并文件中存在共同样本(最好事先避免这一点),在第一个文件中已经存在的标记数据,默认情况下不会被第二个文件覆盖。例如第一个文件某个体某snp为AA,第二个文件该个体该SNP为AG,最终合并文件中还是AA。覆盖情况可以通过–merge-mode调整。
-
不同map相同snp的名称必须保持一致,不然会被 plink 软件视为不同的位点。
-
两个芯片文件的snp编码方式要保持一致,比如不能一个ACGT,一个1234。
计算位点间的 r2
r2 和 D’ 定义
首先我们先看一下理论,首先我们假设有两个位点,其等位基因分别为 Aa, Bb
等位基因 | A | a | B | b |
等位基因频率 | ||||
配子类型 | AB | Ab | aB | ab |
连锁平衡时的理论频率 | ||||
实际频率 | r | s | t | u |
偏离值 | +D | -D | -D | +D |
首先我们可以先求出所有样本中所有位点的等位基因频率,就是这里的 pA qA pB qB 。
如果这两个位点之间不存在连锁,那么这两个位点的组合的频率应该正好等于位点等位基因之和,例如 P(AB) = P(A)P(B) (孟德尔的自由组合定律)。
这里我们计算出第一个统计量,D 值 (连锁不平衡系数,coefficient of linkage disequilibrium),其计算公式为(缺证明)
如果连锁平衡,则 。
r2 计算公式为
这里r 是两个位点之间的相关系数(缺证明)
D’ 定义如下,这里 是
可以取到的最大值。D’ 取值范围同样为 [0,1]。 当等位基因频率较低时,D’ 有偏,偏高。
计算 r2
使用 –r2 计算位点间碱基计数的相关系数的平方看,举例如下,会生成 .ld
为后缀的结果文件。
1 | plink --allow-extra-chr --chr-set 95 --file pick --r2 |
可以设置一下过滤的参数
-
–ld-window 10 :表示计算LD的区间(SNP数目,默认为10),表示距离小于这个值的标记对都要进行LD的计算。
-
–ld-window-kb 1000 :默认为1Mb,表示只对距离在1Mb之内的SNP位点进行分析。
-
–ld-window-r2 0.2 :这个参数只能和 --r2参数搭配使用,默认值为0.2对输出结果进行过滤,只输出r2大于该参数的r2值
如果按照默认参数,就是只算两个位点对距离小于1000kb,中间间隔的SNP数目小于9(10-1)的位点对的r2,并且只输出 r2 大于 0.2 的结果。
指定位点和其它位点的r2
举个例子,我只需要位点Chr11_11501307 前后 5kb 距离位点与这个位点的 r2 ,那么我们首先需要提取这个区域的位点基因型。
然后我们计算 r2 时,只需要加上参数 --ld-snp Chr11_11501307
,同时设置其它过滤参数,避免结果被过滤(比如 --ld-window-r2 0
就是不基于 r2 大小过滤结果)。
1 | plink --allow-extra-chr --chr-set 95 --bfile pick --r2 --ld-snp Chr11_11501307 --ld-window 10000 --ld-window-kb 1000 --ld-window-r2 0 |
输出的结果文件第一个位点便均是 Chr11_11501307 。
类似地选项还有 --ld-snps
(后面跟着多个位点名称), --ld-snp-list
(后面跟着一列位点名称的文件)