PLINK語法體驗(yàn)
by張成龍 郵箱：yianquanwl@qq.com 該教程參考鄧飛博客與自己總結(jié)
plink軟件是GWAS分析中常用的軟件仑撞，它也是一個(gè)數(shù)據(jù)格式，plink里面有很多非常強(qiáng)大的功能妖滔，運(yùn)算速度很快隧哮，是我日常分析中常用的軟件之一。(這里軟件使用的版本為plink 1.9)
這里座舍，我將plink軟件分為三部分：
- 格式轉(zhuǎn)換
- 常用質(zhì)控
- 文件提取
1.格式轉(zhuǎn)換
第一種常用的格式：plink格式
- 正常格式
map和ped：比如a.ped沮翔，a.map。 - 二進(jìn)制文件
bim曲秉，bed采蚀，fam：比如a.bed, a.bim, a.fam。 - 第二種常用的格式：
vcf格式承二。 - 第三種常用的格式：
hapmap格式榆鼠。
1.1 plink正常格式轉(zhuǎn)二進(jìn)制格式
比如這里有plink格式的文件，前綴為a的plink文件：
$ ls
a.map a.ped
將其轉(zhuǎn)化為二進(jìn)制文件：b.bed, b.bim, b.fam
plink --file a --out b
結(jié)果
$ ls b*
b.bed b.bim b.fam b.log
注意：如果染色體超過23亥鸠，比如30對(duì)染色體璧眠，需要設(shè)定--chr-set 30
如果有非數(shù)字染色體，比如性染色體读虏，需要設(shè)定--allow-extra-chr
常用的動(dòng)物都有對(duì)應(yīng)的參數(shù)责静，直接設(shè)定相關(guān)動(dòng)物就行，比如牛的--cow盖桥，下面是其它動(dòng)植物的灾螃。
如果沒有對(duì)應(yīng)的物種，直接設(shè)置染色體的條數(shù)以及允許非數(shù)字染色體即可揩徊。
--cow
--dog
--horse
--mouse·
--rice
--sheep
1.2 plink二進(jìn)制格式轉(zhuǎn)為正常格式（map和ped）
這里有plink格式的文件腰鬼，前綴為b的plink二進(jìn)制文件：
$ ls b*
b.bed b.bim b.fam b.log
將其轉(zhuǎn)化文件：c.map, c.ped
plink --bfile b --recode --out c
注意：
- –bfile，因?yàn)檩斎胛募*為二進(jìn)制塑荒，所以用–bfile熄赡，如果是一般格式，用–file即可
- –recode齿税，要輸出正常格式彼硫，所以用–recode指定，如果不加這個(gè)參數(shù)，默認(rèn)是輸出二進(jìn)制文件
- –out拧篮，輸出文件的前綴
結(jié)果：
$ ls *c*
c.hh c.log c.map c.ped
1.3 正常plink文件轉(zhuǎn)為vcf文件
這里有plink格式的文件词渤，前綴為c的plink二進(jìn)制文件：
$ ls *c*
c.hh c.log c.map c.ped
將其轉(zhuǎn)化文件：d.vcf
plink --file c --recode vcf --out d
注意：
- –file，用–file指定正常plink格式的文件
- –recode vcf串绩，要輸出vcf文件格式
- –out缺虐，輸出文件的前綴
1.4 二進(jìn)制plink文件轉(zhuǎn)為vcf文件
和正常plink文件類似，除了--file 變?yōu)?-bfile即可礁凡。
現(xiàn)有文件：
$ ls b*
b.bed b.bim b.fam b.log
將二進(jìn)制文件轉(zhuǎn)化為vcf文件：
plink --bfile b --recode vcf --out e
1.5 vcf文件轉(zhuǎn)化為plink文件
轉(zhuǎn)化為正常plink文件：
現(xiàn)有文件：
$ ls e.vcf
e.vcf
plink --vcf e.vcf --recode --out f
注意：
- –vcf 需要文件名完整高氮，不能只寫前綴，所以這里要寫--vcf e.vcf
- –recode 保存plink文件
保存為二進(jìn)制文件：
plink --vcf e.vcf --out g
$ ls g*
g.bed g.bim g.fam g.log
2.常用質(zhì)控
2.1 SNP缺失質(zhì)控
無論是測序還是芯片顷牌，得到的基因型數(shù)據(jù)要進(jìn)行質(zhì)控纫溃，而對(duì)缺失數(shù)據(jù)進(jìn)行篩選，可以去掉低質(zhì)量的數(shù)據(jù)韧掩。如果一個(gè)個(gè)體紊浩，共有50萬SNP數(shù)據(jù)，發(fā)現(xiàn)20%的SNP數(shù)據(jù)（10萬）都缺失疗锐，那這個(gè)個(gè)體我們認(rèn)為質(zhì)量不合格坊谁，如果加入分析中可能會(huì)對(duì)結(jié)果產(chǎn)生負(fù)面的影響，所以我們可以把它刪除滑臊。同樣的道理口芍，如果某個(gè)SNP，在500個(gè)樣本中雇卷，缺失率為20%（即該SNP在100個(gè)個(gè)體中都沒有分型結(jié)果）鬓椭，我們也可以認(rèn)為該SNP質(zhì)量較差，將去刪除关划。當(dāng)然小染，這里的20%是過濾標(biāo)準(zhǔn)，可以改變質(zhì)控標(biāo)準(zhǔn)贮折。
現(xiàn)有文件：
$ ls a*
a.map a.ped
某個(gè)SNP在樣本中缺失大于10%裤翩，刪除該SNP：--geno
plink --file a --geno 0.1 --recode --out re
某個(gè)在某個(gè)樣本中，SNP缺失大于10%调榄，刪除該樣本：--mind
plink --file a --mind 0.1 --recode --out re
2.2 最小等位基因頻率過濾
最小等位基因頻率怎么計(jì)算踊赠？比如一個(gè)位點(diǎn)有AA或者AT或者TT，那么就可以計(jì)算A的基因頻率和T的基因頻率每庆，qA + qT = 1筐带，這里誰比較小，誰就是最小等位基因頻率缤灵，比如qA = 0.3, qT = 0.7伦籍， 那么這個(gè)位點(diǎn)的MAF為0.3. 之所以用這個(gè)過濾標(biāo)準(zhǔn)蓝晒，是因?yàn)镸AF如果非常小，比如低于0.02鸽斟，那么意味著大部分位點(diǎn)都是相同的基因型拔创，這些位點(diǎn)貢獻(xiàn)的信息非常少利诺，增加假陽性富蓄。更有甚者M(jìn)AF為0慢逾，那就是所有位點(diǎn)只有一種基因型立倍，這些位點(diǎn)沒有貢獻(xiàn)信息琢歇，放在計(jì)算中增加計(jì)算量岁疼，沒有意義掺涛，所以要根據(jù)MAF進(jìn)行過濾健芭。
現(xiàn)有文件：
$ ls a*
a.map a.ped
某個(gè)SNP在的MAF小于0.01扣癣，那么該SNP刪掉：--maf 0.01
plink --file a --maf 0.01 --recode --out re
2.3 哈溫平衡過濾
「卡方適合性檢驗(yàn)八千！」 忘分，一個(gè)群體是否符合這種狀況缸剪，即達(dá)到了遺傳平衡策添，也就是一對(duì)等位基因的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這對(duì)基因在群體中達(dá)此狀態(tài)唯竹，就是達(dá)到了遺傳平衡乐导。如果沒有達(dá)到這個(gè)狀態(tài)，就是一個(gè)遺傳不平衡的群體浸颓。但隨著群體中的隨機(jī)交配物臂，將會(huì)保持這個(gè)基因頻率和基因型分布比例，而較易達(dá)到遺傳平衡狀態(tài)产上。應(yīng)用Hardy-Weinberg遺傳平衡吻合度檢驗(yàn)方法棵磷，把計(jì)算得到的基因頻率代入，計(jì)算基因型平衡頻率晋涣，再乘以總?cè)藬?shù)泽本，求得預(yù)期值（e）。把觀察數(shù)（O）與預(yù)期值（e）作比較姻僧，進(jìn)行χ2檢驗(yàn)规丽。病例組和對(duì)照組的基因型分布的觀察值和預(yù)期值差異無顯著性（P>0.05），符合遺傳平衡定律.
現(xiàn)有文件：
$ ls a*
a.map a.ped
某個(gè)SNP在哈溫平衡檢驗(yàn)中p值小于1e-5撇贺，那么該SNP刪掉：--maf 0.01
plink --file a --hwe 1e-5 --recode --out re
1. 文件提取
文件提取赌莺，可以提取plink個(gè)數(shù)中的樣本信息，也可以提取特定的SNP位點(diǎn)信息松嘶。
3.1 樣本提取--keep和-- remove
–keep艘狭， 提取樣本ID
–remove，刪除樣本ID
提取樣本文件的格式：
- 第一列：FID，家系ID
- 第二列：IID巢音，個(gè)體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
3.2 SNP提取--extract和-- exclude
- –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個(gè)SNP
SNP刪除
plink --file a --exclude id_snp.txt --recode --out re
SNP合并
一傲绣、合并文件
plink合并文件需要用到merge參數(shù)
如果是ped和map格式文件掠哥，則用以下命令：
plink --file data1 --merge data2.ped data2.map --recode --out merge
如果是二進(jìn)制文件和ped,map格式文件，則用以下命令：
plink --bfile data1 --merge data2.ped data2.map --make-bed --out merge
如果都是二進(jìn)制文件秃诵，則用以下命令：
plink --bfile data1 --bmerge data2.bed data2.bim data2.fam --make-bed --out merge
如果是合并多個(gè)文件续搀，則用以下命令：
/plink-1.07-x86_64/plink --noweb --bfile file --merge-list batch.txt --make-bed --out batch
batch.txt的文件格式如下：
file1.bed file1.bim file1.fam
file2.bed file2.bim file2.fam
更新SNP位置
假設(shè)更新 rs10002到位置580000，如下所示：
原始文件：
...
rs10001 500000
rs10002 580000
rs10003 540000
rs10004 560000
...
新的文件：
...
rs10001 500000
rs10003 540000
rs10004 560000
rs10002 580000
...
更新SNP位置可以采用plink的--update-name和--update-chr參數(shù)
具體命令如下：
./plink --bfile mydata --update-map rsID.lst --update-name --make-bed --out mydata2
或者
./plink --bfile mydata --update-map chr-codes.txt --update-chr --make-bed --out mydata2
rsID.lst的輸入格式如下：
SNP_A-1919191 rs123456
SNP_A-64646464 rs222222
...
chr-codes.txt的輸入格式如下：
rs123456 1
rs987654 18
rs678678 X
..
后續(xù)出現(xiàn)的小問題
plink合并不了菠净，有兩個(gè)方面的問題
- map文件問題
- ped文件問題
1.map文件要統(tǒng)一禁舷，包括pos名稱，統(tǒng)一修改毅往，建議采用·plink提取的方法牵咙，方法見前面。
2.ped第六列要一致煞抬，至于怎么一致有兩種方法霜大。第一正則表達(dá)式，具體方法參考正則表達(dá)式章節(jié)革答。
第二種方法战坤，先將plink文件轉(zhuǎn)換為vcf格式，方法見前面残拐。然后將vcf轉(zhuǎn)為plink途茫，此時(shí)轉(zhuǎn)換語句為
plink --vcf input.vcf --recode --out output --const-fid family_id
通過--const-fid將family id設(shè)置成一個(gè)常量，默認(rèn)值是0.這樣我們得到的文件就可以合并了溪食。
