plink - 关于提取某一个特定的SNP
程序员文章站
2022-03-11 21:41:25
...
一、 查看某一个snp的基因型频率
plink --bfile file --snp rs10402893 --freq --out rs10402893 --noweb
如果要计算全部snp频率则去掉“–snp rs10402893”命令就好
二、 计算某一个snp是否符合哈迪温伯格检验检验:
跟上面计算基因型频率类似,不过是将 --freq 改成 --hardy
plink --bfile file --snp rs10402893 --hardy --out rs10402893_hw --noweb
同样,如果要计算全部snp则去掉“–snp rs10402893”命令就好
结果如图所示:
第一列 snp 所在染色体
第二列 snp 名称
第三列 test的名称
第四列 Minor allele code 次要等位基因频率
第五列 Major allele code 主要等位基因频率
第六列 具体数据 也就是 AA Aa aa 的个数
第七列 观察到的2pq 的值
第八列 期望的2pq的值
第九列 对这个数据进行卡方检验,看显不显著
三、 提取某一个特定的SNP的信息
plink --bfile file --extract rs10402893.txt --make-bed --out rs10402893
会得到关于这一个snp的.bed .fam 和 .bim文件
如果想提取很多snp,则把rs号写进“–extract rs10402893.txt”的txt中,格式为一列rs号(即每一个rs号占一行)
- 如果是要提取所有样本的指定snp,比如有500个样本,我们要看一下这500个样本的遗传数据里在这个特定snp位点的碱基组成详细信息,那么这个既有样本ID信息,又有对应等位基因信息的文件应该是ped了。
根据需求,我们要生成的也应该是ped文件。因此把–make-bed换成
plink --file ped --extract rs10402893.txt --recode --out rs10402893
注意:-file可以自动读取ped文件,因此像-bfile命令一样,不需要文件后缀,不然就报错了哦亲:>
四:plink --keep命令提取表型
举例,如果自己的bfile中有1000个样本,然而自己做关联分析只需要其中的980个,另外20个被剔除,这个时候就需要用keep命令提取我们所需要的数据。
plink中的介绍为:
- bfile为提取数据的源头文件
- keep为要保留的ID,这里的文件需要自己准备:
–keep accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column
说明这个文件必须是要有两行:第一行family IDs,第2行within-family IDs(与fam文件类似) - out定义输出文件名
五、使用plink实现bed、bim与ped、map之间的转换
- plink将bed文件转化为ped,map文件
plink --bfile file --recode --out file --noweb
- plink将ped,map转为二进制格式(bed, bim,fam)
plink --file test --make-bed --out test3
@小插曲 我竟然在macbook上搞了好久怎么写txt文件,而不是什么多信息文本??
打开文本编辑器 - 文件 - 新建 完了之后点击菜单栏上的 格式 - 制作纯文本 就可以了!>