Plink提取tagSNPs

1）識別haplotype block

在提取tagSNPs之前需要先識別haplotype block，這一步的代碼很多教程里都給了出來：

plink --file mydata --blocks

這里的mydata是由vcf文件經(jīng)過plink轉(zhuǎn)換來的，

如果是ped/map格式的文件(mydata.ped/mydata.map)，就用–file

如果是bed/fam格式的文件(mydata.bed/mydata.fam)，就用–bfile

當(dāng)我自己運行的時候出現(xiàn)了小問題（好像是因為表型的問題），文件不能正確的輸出，根據(jù)軟件給出的的提示加入?yún)?shù)：no-pheno-req

修改后的代碼為：

plink --file mydata --blocks no-pheno-req

這一步輸出兩個文件，分別是plink.blocks和plink.blocks.det

plink.blocks文件：

plink.blocks.det文件：

2）提取tagSNPs

網(wǎng)上的教程和官方文檔也給出了相應(yīng)的代碼：

plink --bfile mydata --show-tags mysnps.txt

這里的mysnps.txt是從plink.blocks文件中提取出來的snp編號（一行一個）：

但是運行結(jié)束后只得到了一個plink.tag文件，以為是哪里出錯了，于是查閱官網(wǎng)，官方的解釋：

plink.tags

that lists all the SNPs in the dataset that tag the SNPs in mysnps.txt (including the SNPs in the original file). A message is also written to the LOG file that indicates how many new SNPs were added

If the option

–list-all

is also added, then an additional file is generated that gives some more details for each target SNP (i.e. each SNP listed in mysnps.txt, in the above example) regarding how many and which tags were set for it. The file is named

plink.tags.list

就是說如果用上述代碼，最后只會輸出一個plink.tag文件，如果再加上–list-all，則會再輸出一個plink.tags.list文件，修改后代碼為：

plink --file mydata --list-all --show-tags mysnps.txt

plink.tag文件和mysnp.txt文件類似，plink.tag,list文件結(jié)果如下：