plink软件cookbook
plink软件是我平时工作中最常用的软件之一，它的特点有两个：
- 快
- 功能强大
快，真的是快，我用perl或者Python编写的代码运行需要50s，plink不到1s完成，在C语言面前，我掌握的语言是苍白的。所以，好好利用plink软件，对于速度的提升非常显著。
功能强大，我在使用plink的过程中，它逐渐给我惊喜，仔细研究说明文档非常有必要。下面我的笔记中，这些功能是比较常用的，后面我也会更新关于plink做GWAS的内容，欢迎继续关注。1
(PART) Part I 软件介绍
1 plink 软件介绍
准备写一系列plink软件常用的命令，最近在数据分析时，需要将基因型的数据转化为0-1-2的形式，编程实现效果太差，100万的数据，plink十几秒完成，真的是厉害，非常值得学习，所以，开始搞起！
1.1 map格式
格式说明链接: http://zzz.bwh.harvard.edu/plink/data.shtml#map
❝ map格式的文件, 主要是图谱文件信息, 主要包括染色体名称, 所在的染色体和所在染色体的坐标. ❞
1, map文件没有行头
2, map文件包括四列: 染色体, SNP名称, SNP位置, 碱基对坐标
- 染色体编号为数字, 未知为0
- SNP名称为字符或数字, 如果不重要, 可以从1编号, 注意要和ped文件SNP列一一对应
- 染色体的摩尔未知(可选项, 可以用0)
- SNP物理坐标
3, 如果只有SNP名称, 可以手动构建map文件, 第二列为SNP名称, 其它三列为0即可.
「英文格式说明」
By default, each line of the MAP file describes a single marker and must contain exactly 4 columns:
chromosome (1-22, X, Y or 0 if unplaced)
rs# or snp identifier
Genetic distance (morgans)
Base-pair position (bp units)
Genetic distance can be specified in centimorgans with the --cm flag. Alternatively, you can use a MAP file with the genetic distance excluded by adding the flag --map3, i.e.
plink --file mydata --map3
In this case, the three columns are expected to be
chromosome (1-22, X, Y or 0 if unplaced)
rs# or snp identifier
Base-pair position (bp units)
Base-pair positions are expected to correspond to positive integers within the range of typical human chromosome sizes.1.1.1 Example
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3- 这里有3个SNP, 分别名为snp1, snp3, snp3 「(第二列)」
- 这三个SNP在第一个染色体上 「(第一列)」
- 第三列为0
- 第四列为SNP所在染色体的坐标
1.2 ped格式
格式说明链接:http://zzz.bwh.harvard.edu/plink/data.shtml#ped
❝ bed格式的文件, 主要包括SNP的信息, 包括个体ID, 系谱信息, 表型和SNP的分型信息. ❞
1, 数据没有行头, 空格或者tab隔开的文件 2, 必须要有六列, 包括系谱信息, 表型信息
- 第一列: Family ID # 如果没有, 可以用个体ID代替
- 第二列: Individual ID # 个体ID编号
- 第三列: Paternal ID # 父本编号
- 第四列: Maternal ID # 母本编号
- 第五列: Sex (1=male; 2=female; other=unknown) # 性别, 如果未知, 用0表示
- 第六列: Phenotype # 表型数据, 如果未知, 用0表示
- 第七列以后: 为SNP分型数据, 可以是AT CG或11 12, 或者A T C G或1 1 2 2
3, 上面六列, 必须要有, 如果没有相关数据, 用0表示.
「英文格式说明」
As well as the --file command described above, PED and MAP files can be specified separately, if they have different names:
plink --ped mydata.ped --map autosomal.map
Note Loading a large file (100K+ SNPs) can take a while (which is why we suggest converting to binary format). PLINK will give an error message in most circumstances when something has gone wrong.
The PED file is a white-space (space or tab) delimited file: the first six columns are mandatory:
Family ID
Individual ID
Paternal ID
Maternal ID
Sex (1=male; 2=female; other=unknown)
Phenotype
The IDs are alphanumeric: the combination of family and individual ID should uniquely identify a person. A PED file must have 1 and only 1 phenotype in the sixth column. The phenotype can be either a quantitative trait or an affection status column: PLINK will automatically detect which type (i.e. based on whether a value other than 0, 1, 2 or the missing genotype code is observed).1.2.1 Example
1 1 0 0 1 0 G G 2 2 C C
1 2 0 0 2 0 A A 0 0 A C
1 3 1 2 1 2 0 0 1 2 A C
2 1 0 0 1 0 A A 2 2 0 0
2 2 0 0 2 2 A A 2 2 0 0
2 3 1 2 1 2 A A 2 2 A A- 数据包括两个家系 「(第一列)」
- 每个家系有三个个体 「(第二列)」
- 第三列父本编号
- 第四列母本编号
- 第五列性别
- 第六列表型值
- 第七列, 第八列为一个基因型
- 第九列, 第十列为第二个基因型
- 第十一列, 第十二列为第三个基因型
2 plink操作练习
2.1 练习1
ped格式是1 1 2 2 的格式转化为0 1 2 这里1 1之间有空格
「test1.ped」
1 1 0 0 1 0 1 1 2 2 1 1
1 2 0 0 2 0 2 2 0 0 2 1
1 3 1 2 1 2 0 0 1 2 2 1
2 1 0 0 1 0 1 1 2 2 0 0
2 2 0 0 2 2 2 2 2 2 0 0
2 3 1 2 1 2 1 1 2 2 1 1「test1.map」
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3「plink命令行」
plink --file test2 --recodeA「结果:」
FID IID PAT MAT SEX PHENOTYPE snp1_2 snp2_1 snp3_2
1 1 0 0 1 -9 0 0 0
1 2 0 0 2 -9 2 NA 1
1 3 1 2 1 2 NA 1 1
2 1 0 0 1 -9 0 0 NA
2 2 0 0 2 2 2 0 NA
2 3 1 2 1 2 0 0 0可以看出, 表型值缺失的部分变为了“-9”, 基因型值缺失的部分变为了NA, snp的major变为了0, snp的minor变为了2, 杂合变为了1.
2.2 练习2
ped格式是11 22 的格式转化为0 1 2 这里11中间没有空格「test2.ped」
1 1 0 0 1 0 11 22 11
1 2 0 0 2 0 22 00 21
1 3 1 2 1 2 00 12 21
2 1 0 0 1 0 11 22 00
2 2 0 0 2 2 22 22 00
2 3 1 2 1 2 11 22 11「test2.map」
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3「命令:」结果:plink.raw
FID IID PAT MAT SEX PHENOTYPE snp1_2 snp2_1 snp3_2
1 1 0 0 1 -9 0 0 0
1 2 0 0 2 -9 2 NA 1
1 3 1 2 1 2 NA 1 1
2 1 0 0 1 -9 0 0 NA
2 2 0 0 2 2 2 0 NA
2 3 1 2 1 2 0 0 02.3 练习3
ped格式是A T C G 的格式转化为0 1 2 这里A A之间有空格
「test3.ped」
1 1 0 0 1 0 A A T T G G
1 2 0 0 2 0 G G 0 0 A G
1 3 1 2 1 2 0 0 A T A G
2 1 0 0 1 0 A A A T 0 0
2 2 0 0 2 2 G G T T 0 0
2 3 1 2 1 2 A A T T G G「test3.map」
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3「plink命令行」
plink --file test3 --recodeA「结果不变」
FID IID PAT MAT SEX PHENOTYPE snp1_2 snp2_1 snp3_2
1 1 0 0 1 -9 0 0 0
1 2 0 0 2 -9 2 NA 1
1 3 1 2 1 2 NA 1 1
2 1 0 0 1 -9 0 0 NA
2 2 0 0 2 2 2 0 NA
2 3 1 2 1 2 0 0 02.4 练习4
ped格式是AT CG 的格式转化为0 1 2 这里AT之间没有空格
「test3.ped」
1 1 0 0 1 0 AA TT GG
1 2 0 0 2 0 GG 00 AG
1 3 1 2 1 2 00 AT AG
2 1 0 0 1 0 AA AT 00
2 2 0 0 2 2 GG TT 00
2 3 1 2 1 2 AA TT GG「test3.map」
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3「plink命令行」
plink --file test4 --recodeA「结果不变」
FID IID PAT MAT SEX PHENOTYPE snp1_2 snp2_1 snp3_2
1 1 0 0 1 -9 0 0 0
1 2 0 0 2 -9 2 NA 1
1 3 1 2 1 2 NA 1 1
2 1 0 0 1 -9 0 0 NA
2 2 0 0 2 2 2 0 NA
2 3 1 2 1 2 0 0 02.5 练习5
ped格式是A T 1 2 的格式转化为0 1 2 这里AT和12都包含
「test3.ped」
1 1 0 0 1 0 G G 2 2 C C
1 2 0 0 2 0 A A 0 0 A C
1 3 1 2 1 2 0 0 1 2 A C
2 1 0 0 1 0 G G 2 2 0 0
2 2 0 0 2 2 A A 2 2 0 0
2 3 1 2 1 2 G G 2 2 C C「test3.map」
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3「plink命令行」
plink --file test4 --recodeA「结果不变」
FID IID PAT MAT SEX PHENOTYPE snp1_2 snp2_1 snp3_2
1 1 0 0 1 -9 0 0 0
1 2 0 0 2 -9 2 NA 1
1 3 1 2 1 2 NA 1 1
2 1 0 0 1 -9 0 0 NA
2 2 0 0 2 2 2 0 NA
2 3 1 2 1 2 0 0 02.6 结论
- 1, ped文件中, SNP的分型是1 1 2 2 还是11 22 还是AA TT 还是 AA 22不影响结果
- 2, ped文件中, SNP转化为012的标准是, 主等位基因为0, 杂合为1, 次等位基因为2
- 3, plink命令中, 如果使用–file name, 那么ped和map文件名为: name.ped 和 name.map
3 plink软件格式转化
plink软件是GWAS分析中常用的软件，它也是一个数据格式，plink里面有很多非常强大的功能，运算速度很快，是我日常分析中常用的软件之一。
之前写了一系列的GWAS教程，点击这里查看，这里继续进行。看到我的学习笔记帮助了一些同学，我也由衷的感到高兴。
格式转换
「第一种常用的格式：plink格式」
- 正常格式
map和ped：比如a.ped，a.map - 二进制文件
bim，bed，fam：比如a.bed, a.bim, a.fam
「第二种常用的格式：vcf格式」
「第三种常用的格式：hapmap格式」
3.1 plink正常格式转二进制格式
比如这里有plink格式的文件，前缀为a的plink文件：
$ ls
a.map a.ped将其转化为二进制文件：b.bed, b.bim, b.fam
plink --file a --out b结果：
$ ls b*
b.bed b.bim b.fam b.log「注意：」
- 如果染色体超过23，比如30对染色体，需要设定
--chr-set 30 - 如果有非数字染色体，比如性染色体，需要设定
--allow-extra-chr - 常用的动物都有对应的参数，直接设定相关动物就行，比如牛的
--cow，下面是其它动植物的。如果没有对应的物种，直接设置染色体的条数以及允许非数字染色体即可。
--cow
--dog
--horse
--mouse·
--rice
--sheep3.2 plink二进制格式转为正常格式
这里有plink格式的文件，前缀为b的plink二进制文件：
$ ls b*
b.bed b.bim b.fam b.log将其转化文件：c.map, c.ped
plink --bfile b --recode --out c「注意：」
- –bfile，因为输入文件b*为二进制，所以用–bfile，如果是一般格式，用–file即可
- –recode，要输出正常格式，所以用–recode指定，如果不加这个参数，默认是输出二进制文件
- –out，输出文件的前缀
结果：
$ ls *c*
c.hh c.log c.map c.ped3.3 正常plink文件转为vcf文件
这里有plink格式的文件，前缀为c的plink二进制文件：
$ ls *c*
c.hh c.log c.map c.ped将其转化文件：d.vcf
plink --file c --recode vcf --out d「注意：」
- –file，用–file指定正常plink格式的文件
- –recode vcf，要输出vcf文件格式
- –out，输出文件的前缀
3.4 二进制plink文件转为vcf文件
和正常plink文件类似，除了--file 变为--bfile即可。
现有文件：
$ ls b*
b.bed b.bim b.fam b.log将二进制文件转化为vcf文件：
plink --bfile b --recode vcf --out e结果预览：
3.5 vcf文件转化为plink文件
「转化为正常plink文件：」
现有文件：
$ ls e.vcf
e.vcf
plink --vcf e.vcf --recode --out f「注意：」
- –vcf 需要文件名完整，不能只写前缀，所以这里要写
--vcf e.vcf - –recode 保存plink文件
保存为二进制文件：
plink --vcf e.vcf --out g结果：
$ ls g*
g.bed g.bim g.fam g.log(PART) Part II 数据质控
4 plink常用质控
4.1 SNP缺失质控
❝ 无论是测序还是芯片，得到的基因型数据要进行质控，而对缺失数据进行筛选，可以去掉低质量的数据。如果一个个体，共有50万SNP数据，发现20%的SNP数据（10万）都缺失，那这个个体我们认为质量不合格，如果加入分析中可能会对结果产生负面的影响，所以我们可以把它删除。同样的道理，如果某个SNP，在500个样本中，缺失率为20%（即该SNP在100个个体中都没有分型结果），我们也可以认为该SNP质量较差，将去删除。当然，这里的20%是过滤标准，可以改变质控标准。 ❞
现有文件：
$ ls a*
a.map a.ped「某个SNP在样本中缺失大于10%，删除该SNP：--geno」
plink --file a --geno 0.1 --recode --out re「某个在某个样本中，SNP缺失大于10%，删除该样本：--mind」
plink --file a --mind 0.1 --recode --out re4.2 次等位基因频率过滤
❝ 次等位基因频率怎么计算？比如一个位点有AA或者AT或者TT，那么就可以计算A的基因频率和T的基因频率，qA + qT = 1，这里谁比较小，谁就是最小等位基因频率，比如qA = 0.3, qT = 0.7， 那么这个位点的MAF为0.3. 之所以用这个过滤标准，是因为MAF如果非常小，比如低于0.02，那么意味着大部分位点都是相同的基因型，这些位点贡献的信息非常少，增加假阳性。更有甚者MAF为0，那就是所有位点只有一种基因型，这些位点没有贡献信息，放在计算中增加计算量，没有意义，所以要根据MAF进行过滤。 ❞
现有文件：
$ ls a*
a.map a.ped「某个SNP在的MAF小于0.01，那么该SNP删掉：--maf 0.01」
plink --file a --maf 0.01 --recode --out re4.3 哈温平衡过滤
❝ 「卡方适合性检验！」 ，一个群体是否符合这种状况，即达到了遗传平衡，也就是一对等位基因的3种基因型的比例分布符合公式：p2+2pq+q2=1,p+q=1,(p+q)2=1.基因型MM的频率为p2,NN的频率为q2,MN的频率为2pq。MN:MN：NN＝P2：2pq：q2。MN这对基因在群体中达此状态，就是达到了遗传平衡。如果没有达到这个状态，就是一个遗传不平衡的群体。但随着群体中的随机交配，将会保持这个基因频率和基因型分布比例，而较易达到遗传平衡状态。应用Hardy-Weinberg遗传平衡吻合度检验方法，把计算得到的基因频率代入，计算基因型平衡频率，再乘以总人数，求得预期值（e）。把观察数（O）与预期值（e）作比较，进行χ2检验。病例组和对照组的基因型分布的观察值和预期值差异无显著性（P>0.05），符合遗传平衡定律. 现有文件： ❞
$ ls a*
a.map a.ped「某个SNP在哈温平衡检验中p值小于1e-5，那么该SNP删掉：--hwe 1e-5」
plink --file a --hwe 1e-5 --recode --out re5 plink文件提取
文件提取，可以提取plink个数中的样本信息，也可以提取特定的SNP位点信息。
5.1 样本提取
- –keep， 提取样本ID
- –remove，删除样本ID
「提取样本文件的格式：」
- 第一列：FID，家系ID
- 第二列：IID，个体ID
1328 NA06989
1377 NA11891
1349 NA11843
1330 NA12341
1344 NA10850
1328 NA06984
1463 NA12877
1418 NA12275
13291 NA06986
1418 NA12272「样本提取」
plink --file a --keep id_sample.txt --recode --out re完成。
$ wc -l re*
2 re.hh
32 re.log
1431211 re.map
10 re.ped「样本删除」
plink --file a --remove id_sample.txt --recode --out re完成。
5.2 SNP提取
- –extract， 提取SNP ID
- –exclude，删除SNP ID
「提取样本文件的格式：」
- 一列：SNP名称ID
rs2185539
rs11240767
rs3131972
rs3131969
rs1048488
rs12562034
rs12124819
rs4040617
rs2905036
rs4245756「SNP提取」
plink --file a --extract id_snp.txt --recode --out re完成。
$ wc -l re*
179 re.hh
30 re.log
10 re.map
164 re.ped可以看到，map共10行，共提取10个SNP
「SNP删除」
plink --file a --exclude id_snp.txt --recode --out re完成。
6 计算杂合度
这里，模拟一个plink文件的数据，8个样本，8个SNP位点，通过手动Excel计算样本杂合度和位点杂合度，比较plink计算杂合度的方法。
6.1 模拟数据
6.1.1 ped数据
$ cat a.ped
FAMILY1 ID1 0 0 0 -9 CC CC AA GG AG GG GG GC
FAMILY1 ID2 0 0 0 -9 CC GC AG GG GG AA AG CC
FAMILY1 ID3 0 0 0 -9 GG CC AG GG GG GA AG GC
FAMILY1 ID4 0 0 0 -9 GG CC GG GG GG AA GG GG
FAMILY1 ID5 0 0 0 -9 GG CC GG GG GG AA AG GC
FAMILY1 ID6 0 0 0 -9 GG CC GG GG GG AA AA CC
FAMILY1 ID7 0 0 0 -9 GG CC GG AG AA AA GG CC
FAMILY1 ID9 0 0 0 -9 GG CC GG AG AA AA GG CC6.1.2 map数据
$ cat a.map
1 snp1 0 55910
1 snp2 0 85204
1 snp3 0 122948
1 snp4 0 203750
1 snp5 0 312707
1 snp6 0 356863
1 snp7 0 400518
1 snp8 0 4874236.2 计算样本的杂合度--het
$ plink --file a --het「结果查看：」
$ cat plink.het
FID IID O(HOM) E(HOM) N(NM) F
FAMILY1 ID1 6 5.32 8 0.2536
FAMILY1 ID2 5 5.32 8 -0.1195
FAMILY1 ID3 4 5.32 8 -0.4927
FAMILY1 ID4 8 5.32 8 1
FAMILY1 ID5 6 5.32 8 0.2536
FAMILY1 ID6 8 5.32 8 1
FAMILY1 ID7 7 5.32 8 0.6268
FAMILY1 ID9 7 5.32 8 0.6268这里：
- FID，家系ID
- IID，个体ID
- O(HOM)：观察到的纯合个数
- E(HOM)：期望的纯合个数
- N(NM)：没有缺失的SNP个数
- F：计算的值
其中F的计算方法：
- O: O(HOM)
- E: E(HOM)
- N: N(NM)
可以这样认为，F值越小（包括负值），杂合度越高，F值越高，纯合度越高！
「Excel对比」将ped文件，copy到Excel中，手动计算纯合和杂合的个数，进行统计：
6.3 计算SNP位点杂合度
这里，使用参数--hardy，可以给出位点的纯合和杂合个数：
$ plink --file a --hardy结果：
$ cat plink.hwe
CHR SNP TEST A1 A2 GENO O(HET) E(HET) P
1 snp1 ALL(NP) C G 2/0/6 0 0.375 0.01538
1 snp2 ALL(NP) G C 0/1/7 0.125 0.1172 1
1 snp3 ALL(NP) A G 1/2/5 0.25 0.375 0.3846
1 snp4 ALL(NP) A G 0/2/6 0.25 0.2188 1
1 snp5 ALL(NP) A G 2/1/5 0.125 0.4297 0.07692
1 snp6 ALL(NP) G A 1/1/6 0.125 0.3047 0.2
1 snp7 ALL(NP) A G 1/3/4 0.375 0.4297 1
1 snp8 ALL(NP) G C 1/3/4 0.375 0.4297 1- GENO，
2/0/6，第一个是次等位基因纯合个数，第二个是杂合个数，第三个是主等位基因纯合个数 - O(HET)，是杂合所在的比值
「对比Excel结果：」
6.4 计算SNP位点的基因频率
$ plink --file a --freq「结果查看：」
$ cat plink.frq
CHR SNP A1 A2 MAF NCHROBS
1 snp1 C G 0.25 16
1 snp2 G C 0.0625 16
1 snp3 A G 0.25 16
1 snp4 A G 0.125 16
1 snp5 A G 0.3125 16
1 snp6 G A 0.1875 16
1 snp7 A G 0.3125 16
1 snp8 G C 0.3125 16(APPENDIX) 附录知识
7 参考文献
plink1.9 官网：https://www.cog-genomics.org/plink/1.9/
8 Top链与Forward链资料整理
最近，有同学问我Illumina芯片数据中Top链，Forward链，以及与refSNP的区别，我查了一些资料，汇总如下。如果有错误，还请留言区批评指正，我也是现学现卖。
参考： 2
这里介绍一下芯片的几种格式：
- Forward allele
- A/B allele
- TOP/BOT allele
8.1 Forward allelel
这是大多数人应该研究的等位基因。它是指参考基因组中前链的等位基因。注意，不同版本的参考基因组往往在具有共同SNP的特定位置上可能有所不同；例如，如果GRCh37在参考链的SNP中有一个次等位基因，GRCh38则倾向于将该等位基因变回主等位基因，因此理想情况下，正向等位基因应该始终与基因组结合构建以唯一标识SNP。
前链（Forward）和后链（Reverse）来源于dbSNP数据库。
8.2 Illumina’s A/B allele coding
这是Illumina定义的类型，等位基因编码方法解决了上述问题，即等位基因不依赖于特定的基因组组合，而是基于实际的多态性本身。
简而言之，如果两个多态性中的一个是A或T，另一个是C或G，那么A或T被称为A等位基因，C或G被称为B等位基因。
Illumina’s A/B allele coding, or TOP/BOT strand definition, is explained in here in detail by Illumina. The allele coding method solves the problem aformentioned, that is, the alleles are not dependent on the specific genome assembly, but are based on the actual polymorphism itself. Briefly, if one of the two polymorphism is A or T, and the other one is C or G, then the A or T is refered to as A allele, and the C or G is refered to as B allele, and the strand with A or T is refered to as TOP and BOT strand, respectively.
有时候，也会用1/2来表示A/B
Sometimes, people often use 1/2 to denote Illumina’s A/B allele, since numeric coding is more convenient in many scenarios and since some old association software only recognize numeric coded alleles.
8.3 TOP/BOT allele
如果两个多态性中的一个是A或T，另一个是C或G，带有A或T的链分别称为TOP和BOT链.
Briefly, if one of the two polymorphism is A or T, and the other one is C or G, then the A or T is refered to as A allele, and the C or G is refered to as B allele, and the strand with A or T is refered to as TOP and BOT strand, respectively.
注意：
如果多态性是A/T或C/G，那么穿过环绕序列（SNP的上游或下游的两个核苷酸）找到一对明确的核苷酸，然后应用类似的规则：如果A或T在SNP的5’侧，那么它是Top链，否则就是BOT链。对于Top链，A和B等位基因分别表示A和T（或C和G）；而对于BOT链，A和B等位基因分别表示T和A（或G和C）
If the polymorphism is A/T or C/G, then walk through the surrouding sequence (the two nucleotides up or downstream of the SNP) to find a pair of unambiguous nucleotides, and then a similar rule is applied: if A or T is on 5’ side of the SNP, then it’s a TOP strand otherwise it’s a BOT strand. For TOP strand, A and B allele denote A and T (or C and G), respectively; whereas for BOT strand, A and B allele denote T and A (or G and C), respectively.
Illumina的编码方案不依赖于前链(Forward)的定义（因此正确的基因组组装），因此它几乎总是确保基因组构建之间的一致性，并确保新测序基因组序列或未组装基因组序列的即时等位基因指定。
另外，在Illumina BeadStudio软件中，可以指定AB类型，或者ACGT类型（TOP链），或者Forward链类型。TOP alleles是TOP链，但不一定是Forward链，具体解释如下：
When exporting genotypes from the Illumina BeadStudio software, the user can choose AB genotypes, or ACGT genotypes (commonly refered to as “TOP alleles”), or forward strand genotype in newer version of the software. The TOP alleles is the allele on the TOP strand, which may or may not be the forward strand: see the example above, the “fwd/B” means that dbSNP’s forward strand corresponds to Illumina’s BOT strand, so the “TOP allele” is the opposite as the “forward strand allele”. Unfortunately many users simply do not know or understand what is “TOP allele”: they simply take for granted that “TOP” means “forward” and then complain that there are many discordances when merging two different data sets (one coded as forward strand and one exported from BeadStudio). The convert_bim_allele.pl program that I describe in this article will solve problems like this.
8.4 用Top链还是Forward链
主流的还是用Forward链多一点，当然如果你之前的数据是Top链，那还是要用Top才可以合并。
注意：Top链和Forward不是对应的！
- Top链与Bot链对应
- Forward链与Reverse链对应
Top链的位点分型，有时候和Forward是一致的，有时候是不一样的。因为Top链的规则是“如果多态性是A/T或C/G，那么穿过环绕序列（SNP的上游或下游的两个核苷酸）找到一对明确的核苷酸，然后应用类似的规则：如果A或T在SNP的5’侧，那么它是Top链，否则就是BOT链”，它会根据SNP的上下游确定。
8.5 dbSNP中T>C是什么意思？
比如rs1004491这个SNP，在dbSNP数据库中是T/C突变
rs1004491 [Homo sapiens]
Variant type:SNVAlleles:T>C[Hide Flanks]
AAAGCCTTCTGAACTGAGTGAAAATACAGCCAAGATCTTGGCAAAGCTTC
TCCCTCAGTATTTAGACCAGGTAAGAATTTCTTGACTCATCTCCAACATA
[T/C]
GTGTTTACTGTGGAAAACACACATTTTATTTTCTTGCTATTGCATGTTAT
TGCTGGCCGGGGACCCAATTGCAGTCTCTTTAAGCCTTCAACAGTTGGCT之所以是T>C，是因为平均而言，这个位点T为主等位基因（major），C为次等位基因（minor）
下图可以看到，整体而言（209010个样本），T的频率为0.701，C的频率为0.298，当然对于少数的群体（比如这里的Asian）中，T为0.482，C为0.518，但整体而言T>C。
好了，就到这里。
脚注：
- 话说用Rmarkdown写博客，真的是很爽。↩︎
- https://gengen.openbioinformatics.org/en/latest/tutorial/coding/#introduction↩︎