plink软件是我平时工作中最常用的软件之一,它的特点有两个:
快,真的是快,我用perl或者Python编写的代码运行需要50s,plink不到1s完成,在C语言面前,我掌握的语言是苍白的。所以,好好利用plink软件,对于速度的提升非常显著。
功能强大,我在使用plink的过程中,它逐渐给我惊喜,仔细研究说明文档非常有必要。下面我的笔记中,这些功能是比较常用的,后面我也会更新关于plink做GWAS的内容,欢迎继续关注。1
准备写一系列plink软件常用的命令,最近在数据分析时,需要将基因型的数据转化为0-1-2的形式,编程实现效果太差,100万的数据,plink十几秒完成,真的是厉害,非常值得学习,所以,开始搞起!
格式说明链接: http://zzz.bwh.harvard.edu/plink/data.shtml#map
❝ map格式的文件, 主要是图谱文件信息, 主要包括染色体名称, 所在的染色体和所在染色体的坐标. ❞
1, map文件没有行头
2, map文件包括四列: 染色体, SNP名称, 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 snp1 0 1
1 snp2 0 2
1 snp3 0 3
格式说明链接:http://zzz.bwh.harvard.edu/plink/data.shtml#ped
❝ bed格式的文件, 主要包括SNP的信息, 包括个体ID, 系谱信息, 表型和SNP的分型信息. ❞
1, 数据没有行头, 空格或者tab隔开的文件 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 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
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.
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 0
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 0
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 0
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 0
plink软件是GWAS
分析中常用的软件,它也是一个数据格式,plink里面有很多非常强大的功能,运算速度很快,是我日常分析中常用的软件之一。
之前写了一系列的GWAS
教程,点击这里查看,这里继续进行。看到我的学习笔记帮助了一些同学,我也由衷的感到高兴。
格式转换
「第一种常用的格式:plink
格式」
map
和ped
:比如a.ped,a.mapbim
,bed
,fam
:比如a.bed, a.bim, a.fam「第二种常用的格式:vcf
格式」
「第三种常用的格式:hapmap
格式」
比如这里有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
「注意:」
--chr-set 30
--allow-extra-chr
--cow
,下面是其它动植物的。如果没有对应的物种,直接设置染色体的条数以及允许非数字染色体即可。--cow
--dog
--horse
--mouse·
--rice
--sheep
这里有plink格式的文件,前缀为b的plink二进制文件:
$ ls b*
b.bed b.bim b.fam b.log
将其转化文件:c.map, c.ped
plink --bfile b --recode --out c
「注意:」
结果:
$ ls *c*
c.hh c.log c.map c.ped
这里有plink格式的文件,前缀为c的plink二进制文件:
$ ls *c*
c.hh c.log c.map c.ped
将其转化文件:d.vcf
plink --file c --recode vcf --out d
「注意:」
和正常plink文件类似,除了--file
变为--bfile
即可。
现有文件:
$ ls b*
b.bed b.bim b.fam b.log
将二进制文件转化为vcf文件:
plink --bfile b --recode vcf --out e
结果预览:
「转化为正常plink文件:」
现有文件:
$ ls e.vcf
e.vcf
plink --vcf e.vcf --recode --out f
「注意:」
--vcf e.vcf
保存为二进制文件:
plink --vcf e.vcf --out g
结果:
$ ls g*
g.bed g.bim g.fam g.log
❝ 无论是测序还是芯片,得到的基因型数据要进行质控,而对缺失数据进行筛选,可以去掉低质量的数据。如果一个个体,共有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 re
❝ 次等位基因频率怎么计算?比如一个位点有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 re
❝ 「卡方适合性检验!」 ,一个群体是否符合这种状况,即达到了遗传平衡,也就是一对等位基因的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 re
文件提取,可以提取plink个数中的样本信息,也可以提取特定的SNP位点信息。
「提取样本文件的格式:」
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
完成。
「提取样本文件的格式:」
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
完成。
这里,模拟一个plink文件的数据,8个样本,8个SNP位点,通过手动Excel计算样本杂合度和位点杂合度,比较plink计算杂合度的方法。
$ 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 CC
$ 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 487423
--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
这里:
其中F的计算方法:
可以这样认为,F值越小(包括负值),杂合度越高,F值越高,纯合度越高!
「Excel对比」将ped文件,copy到Excel中,手动计算纯合和杂合的个数,进行统计:
这里,使用参数--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
2/0/6
,第一个是次等位基因纯合个数,第二个是杂合个数,第三个是主等位基因纯合个数「对比Excel结果:」
$ 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
plink1.9 官网:https://www.cog-genomics.org/plink/1.9/
最近,有同学问我Illumina芯片数据中Top链,Forward链,以及与refSNP的区别,我查了一些资料,汇总如下。如果有错误,还请留言区批评指正,我也是现学现卖。
参考: 2
这里介绍一下芯片的几种格式:
这是大多数人应该研究的等位基因。它是指参考基因组中前链的等位基因。注意,不同版本的参考基因组往往在具有共同SNP的特定位置上可能有所不同;例如,如果GRCh37在参考链的SNP中有一个次等位基因,GRCh38则倾向于将该等位基因变回主等位基因,因此理想情况下,正向等位基因应该始终与基因组结合构建以唯一标识SNP。
前链(Forward)和后链(Reverse)来源于dbSNP数据库。
这是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.
如果两个多态性中的一个是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.
主流的还是用Forward链多一点,当然如果你之前的数据是Top链,那还是要用Top才可以合并。
注意:Top链和Forward不是对应的!
Top链的位点分型,有时候和Forward是一致的,有时候是不一样的。因为Top链的规则是“如果多态性是A/T或C/G,那么穿过环绕序列(SNP的上游或下游的两个核苷酸)找到一对明确的核苷酸,然后应用类似的规则:如果A或T在SNP的5’侧,那么它是Top链,否则就是BOT链”,它会根据SNP的上下游确定。
比如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。
好了,就到这里。
脚注: