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 文件内容如下 (chromosome, SNP, cM, base-position, allele 1, allele 2)（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 文件。

需要注意的是，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名称列

​ # 是需要跳过的行

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 质控位点，不然得到的结果更加糟糕。

基于基因组亲缘关系剔除样本

使用选项 --rel-cutoff [maximum] 剔除基因组亲缘关系（应该是 G阵）较近的个体间的一个，会输出 plink.rel.id ，即最终保存的个体的ID。

位点质控

杂合子比例

位点杂合子比例可以通过 --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

剔除非SNP位点

使用 –snps-only just-acgt 选项会剔除非 SNP 位点（例如 indel 位点），–snps-only 会剔除所有存在至少一个等位基因是多字符的位点（也就是所只会保留两个等位基因均是单字符的位点）。just-acgt 会进一步剔除等位基因不在 {'A', 'C', 'G', 'T', 'a', 'c', 'g', 't', <missing code>} 范围内容的位点

对于存在三等位基因的 SNP ，PLINK默认就会将存在三等位基因的基因型替换为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.

VCF文件剔除3等位基因位点

当输入为 VCF 文件时，可以使用下面的选项剔除多等位基因位点

--biallelic-only ['strict'] ['list']

使用–biallelic-only 会剔除基因型中实际存在多等位碱基的位点（3个或以上），添加 strict 会无差别删除次等位基因ALT是2个及2个以上的情况（哪怕基因型中只出现了1个ALT），添加 list 会将剔除的位点放在 plink-temporary.skip.3allele 文件中。举例如下

plink --allow-extra-chr --chr-set 95 --vcf test.vcf --const-fid --biallelic-only list --recode --out out4

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的计算（换句话说，两个位点之间的SNP数目大于等于 (10-1) 的情况不计算 r2）。

• –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 的结果。

plink 有时输出的结果不全，因为有些位点对的 r2 是 nan 就没有输出，原因是对于对于两个位点同时有分型的样本（剔除缺失后），至少有一个位点没有变异 (来自于 https://groups.google.com/g/plink2-users/c/zD2aME9Kb3o)。

The only SNP-pairs that still get filtered under “–ld-window-r2 0” are those where the r^2 value is “nan”, i.e. at least one of the SNPs exhibits no variation at all, across all samples genotyped at both SNPs. If you want that case to be reported as well, you can use e.g. “–ld-window-r2 -1”. (This detail is mentioned in the plink 2.0 documentation.)

所以推荐使用 plink 2.0 版本进行分析

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

plink 2.0 版本分析r2

官网文档：https://www.cog-genomics.org/plink/2.0/ld#r

上面计算的r2 其实是碱基计数的相关系数的平方，不是按照连锁不平衡的公式计算的 r2

在 plink 2.0 中，计算 r2 有两种方式，--r2-unphased 对应碱基计数的相关系数的平方，--r2-phased 对应教科书上的基于单倍型频率计算的 r2 。如果要输出 nan 的结果，则 --ld-window-r2 的值设为 - 1，举例如下 （这里不需要设定 --ld-window-kb，因为 plink2 的默认值是无穷大）

./plink2 --allow-extra-chr --chr-set 95 --pedmap check_snp --r2-phased  --ld-snp-list need_snp_test --ld-window-kb 17 --ld-window-r2 -1