plink软件备忘录

这是我个人使用的plink软件备忘录，将我所有会用到但是又没有完全记住的选项和功能记录在这里，方便查看。

windows中使用plink

1. 进入plink官网，下载 plink

2. 将 plink.exe 加入到环境变量 （我的做法是建立快捷方式，命名为 plink ，放到 D:\shortcuts 中）

输入输出

输入

compound格式

使用 –compound-genotypes 选项，读取 compound格式的plink文件（一个SNP的两个碱基之间没有空格）

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。规则如下：

00  Homozygote "1"/"1" （first allele in .bim file）
 01  Heterozygote
 11  Homozygote "2"/"2"
 10  Missing genotype

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开始）。

|-magic number--| |-mode-| |--genotype data---------| 

01101100 00011011 00000001 11011100 00001111 11100111

|--genotype data-cont'd--| 

00001111 01101011 00000001 

这里其相应的 bim 文件内容如下（A1 为等位基因频率较小的碱基 ）

1       snp1    0       1       G       A
1       snp2    0       2       1       2
1       snp3    0       3       A       C

fam 文件如下

1 1 0 0 1 0
1 2 0 0 1 0
1 3 1 2 1 2
2 1 0 0 1 0
2 2 0 0 1 2
2 3 1 2 1 2

这里还有一个小问题是，在读取每个字节数据的时候是从后向前读的。如果把一个字节的8个位置从后到前标记伪 A 到 H，举例如下

01101100
HGFEDCBA

第一个SNP的前4个基因型则为 AB, CD,EF,GH 四个位置对应的分型，如下：

01101100
HGFEDCBA

      AB   00  -- homozygote (first)
    CD     11  -- other homozygote (second)
  EF       01  -- heterozygote (third)
GH         10  -- missing genotype (fourth)

上面例子中的 3 个位点的解析如下，第一个SNP的 A1 和 A2 碱基为 G 和 A ，因此 00 对应 GG，01 对应 GA ，11 对应 GG ，10 对应缺失。由于这里总共是 6 个样本，因此一个 SNP 需要占用 2 个字节，第二个字节多余的位置使用 00 填充。

因此第一个SNP的第一个样本是GG，之后是AA,00, AA, AA, AA。

           Genotype    Person    SNP
11011100 

      00   G/G         1 1       snp1
    11     A/A         1 2       snp1
  10       0/0         1 3       snp1
11         A/A         2 1       snp1


00001111 

      11   A/A         2 2       snp1
    11     A/A         2 3       snp1
  00       (null)
00         (null)


11100111
      
      11   2/2         1 1       snp2
    10     0/0         1 2       snp2
  01       1/2         1 3       snp2
11         2/2         2 1       snp2


00001111 
  
      11   2/2         2 2       snp2
    11     2/2         2 3       snp2
  00       (null) 
00         (null)


01101011

      11   C/C         1 1       snp3
    01     A/C         1 2       snp3
  01       A/C         1 3       snp3
10         0/0         2 1       snp3


00000001

      10   0/0         2 2       snp3
    00     A/A         2 3       snp3
  00       (null)
00         (null)

VCF 文件

联合使用 –vcf filename 和 –const-fid 来读取 vcf 文件，这里 –const-fid 是将所有个体的家系号均设为 0 。

注意 VCF 文件作为 plink 的输入文件，其位点必须严格按照物理位置顺序排序，不然就会遇到下面的报错。

Error: .bim file has a split chromosome.  Use --make-bed by itself to
remedy this.

你需要先用 –make-bed 生成标准的二进制文件（排序好的），再进行其他处理。

使用 –bcf 选项可以读取 BCF2 文件。

输出

二进制格式

使用 --make-bed 输出二进制格式，使用--make-just-bim 和 --make-just-fam 只输出 bim 或 fam 文件。

compound格式

使用 –recode compound-genotypes 选项，生成 compound格式的plink文件（一个SNP的两个碱基之间没有空格）

plink --allow-extra-chr --chr-set 95 --file plink59 --recode compound-genotypes --out plink59_compound

012格式

使用 --recodeA 选项可以输出 012 格式的基因型文件，默认是按照 A1 (最小等位基因) 进行计数。

可以通过添加 --recode-allele 选项指定每一个 SNP 是计数哪一个碱基，其后面接着输入文件名称，这个输入文件含有两列，第一列是 SNP 名称，第二列是需要计数的碱基。

举例如下

plink --allow-extra-chr --chr-set 95 --file 1 --recode-allele count_allele.txt --recodeA --out 1_count

其中的 count_allele.txt 文件内容如下

1_242598	C
1_10673082	T
1_10723065	T
1_11407894	A
1_11426075	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 保持一致，你可以添加一下选项：

--a2-allele ref.vcf 4 3 '#'

参数说明：

--a2-allele ref.vcf 4 3 '#' :

​ --a2-allele 命令表示 指定A2基因（Major）

​ ref.vcf 是参考群体的vcf文件（必须是未压缩文件）

​ 4 是参考群体A2基因所在列

​ 3 是SNP名称列

​ # 是需要跳过的行

需要注意的地方如下

1. ref.vcf 和 plink 文件中的染色体，snp名称，物理位置，两个碱基必须一致
2. 需要转换的 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格式，分两步走。

plink --allow-extra-chr --chr-set 95 --chr 1-18 --file raw --indiv-sort file keep_id.txt --make-bed --out sorted

## --indiv-sort 只能生成二进制文件，但是顺序果然是一样的。

plink --allow-extra-chr --chr-set 95 --bfile sorted --recode --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 染色体上的伪常染色体区域的位点。

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 为后缀的结果文件，例如

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) 。

由于输出文件中包含所有样本两两配对的结果，因此文件会比较大。

FID1	Family ID for first sample
IID1	Individual ID for first sample
FID2	Family ID for second sample
IID2	Individual ID for second sample
RT	Relationship type inferred from .fam/.ped file
EZ	IBD sharing expected value, based on just .fam/.ped relationship
Z0	P(IBD=0)
Z1	P(IBD=1)
Z2	P(IBD=2)
PI_HAT	Proportion IBD, i.e. P(IBD=2) + 0.5*P(IBD=1)
PHE	Pairwise phenotypic code (1, 0, -1 = AA, AU, and UU pairs, respectively)
DST	IBS distance, i.e. (IBS2 + 0.5*IBS1) / (IBS0 + IBS1 + IBS2)
PPC	IBS binomial test
RATIO	HETHET : IBS0 SNP ratio (expected value 2)

三种亲缘关系的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 。

全同胞的结果如下

半同胞的结果如下

测试结论

1. 三种关系中 ，Z2 值均偏大，远高于期望值，因此这里的 IBD 算法不够精确。

2. 三种关系中 ，PI_HAT 值均偏大，远高于期望值，估计是收到了 Z2 值偏大的影响。

3. 这里最好区分的是亲子关系，其 Z0 值接近于 0，Z1 值均大于 0.6 。但是这里半同胞和全同胞的 PI_HAT 的分布有交叉，无法完全分开。

4. 总体来说，IBD 的计算结果只能说勉强符合要求，应该只能用于筛选样本，比如挑选或剔除亲缘关系比较近的样本对。

5. 计算 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 ，不会报错，只会给出下面的提示：

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 命令进行部分位点的翻转。如果不是，估计就是基因型文件分型出问题了。

-- 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 合并时有以下注意事项：

1. 合并后的基因型文件位点是所有被合并的基因型文件位点的并集，不存在的位点基因型会设为缺失。

2. 如果多个被合并文件中存在共同样本（最好事先避免这一点），在第一个文件中已经存在的标记数据，默认情况下不会被第二个文件覆盖。例如第一个文件某个体某snp为AA，第二个文件该个体该SNP为AG，最终合并文件中还是AA。覆盖情况可以通过–merge-mode调整。

3. 不同map相同snp的名称必须保持一致，不然会被 plink 软件视为不同的位点。

4. 两个芯片文件的snp编码方式要保持一致，比如不能一个ACGT，一个1234。

计算位点间的 r2

r2 和 D’ 定义

首先我们先看一下理论，首先我们假设有两个位点，其等位基因分别为 Aa, Bb

等位基因	A	a	B	b
等位基因频率
配子类型	AB	Ab	aB	ab
连锁平衡时的理论频率
实际频率	r	s	t	u
偏离值	+D	-D	-D	+D

首先我们可以先求出所有样本中所有位点的等位基因频率，就是这里的 p_A q_A p_B q_B 。

如果这两个位点之间不存在连锁，那么这两个位点的组合的频率应该正好等于位点等位基因之和，例如 P(AB) = P(A)P(B) （孟德尔的自由组合定律）。

这里我们计算出第一个统计量，D 值 (连锁不平衡系数，coefficient of linkage disequilibrium），其计算公式为（缺证明）

如果连锁平衡，则 。

r2 计算公式为

这里r 是两个位点之间的相关系数（缺证明）

D’ 定义如下，这里 是 可以取到的最大值。D’ 取值范围同样为 [0,1]。 当等位基因频率较低时，D’ 有偏，偏高。

计算 r2

使用 –r2 计算位点间碱基计数的相关系数的平方看，举例如下，会生成 .ld 为后缀的结果文件。

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 大小过滤结果）。

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 (后面跟着一列位点名称的文件)