plink学习笔记

2021-01-12 本文已影响0人斩毛毛

官方文档学习更详细，轻用力戳此处

1. 基本格式

1.1 ped, map (二者为一组)

pep格式包括基因型信息，如下所示

A01     A01     0       0       0       0       G       G       T       T
A02     A02     0       0       0       0       T       G        G       T

前六列固定
第一列 Family ID
第二列 Individual ID
第三列 Paternal ID
第四列 Maternal ID
第五列 Sex (1=male; 2=female; other=unknown)
第六列 Phenotype
再往后，就是基因型数据，每两列为一个样本基因型，比如第一个样本基因型为GG,第二个为TT

map 格式包括基因型信息，如下所示

0       Chr00:49209     0       49209
0       Chr00:49287     0       49287

表示SNP的信息
第一列代表染色体
第二列代表snp的名称
第三列代表摩尔距离，一般用不上
第四列代表物理距离，这个有用的

1.2 bed fam bin 为一组

首先明确bed是一个二进制的文件，与fam，bin文件互相对应。基因型用0,1表示，具体如下

00 ref 纯合
11 alt 纯合
01 缺失
10 杂合

bim文件

0       Chr00:49209     0       49209  C  T
0       Chr00:49287     0       49287  G  T

该文件是在map的结果上在添加两列SNP位点

A01 A01 0 0 0 -9
A02 A02 0 0 0 -9

第一列 Family ID
第二列 Individual ID
第三列 Paternal ID （0表示无）
第四列 Maternal ID（（0表示无））
第五列 Sex (1=male; 2=female; 0=unknown)
第六列 Phenotype（0/-9 表示无）

1.3 ped 和bed互相转化

ped 转bed

plink --allow-extra-chr --noweb -file test--make-bed --out test1
# -file： 输入为ped
# --noweb 不显示网页

bed 转ped

plink --bfile test1 --recode --out test
#--bfile: 输入为bed
#--recode 输出格式ped（默认bed）

2. 可接受的格式

除了上面提到的bed，ped格式以外，我们比较常用的就是vcf,or bcf格式

vcf 格式转为ped格式

plink --vcf SNP.vcf.gz--recode --out test--const-fid --allow-extra-chr

# --vcf, 压缩或者不压缩都可以接受

$\color{red}{--allow-extr-chr}$ :因为plink本身是针对人类进行开发的，所以遇到格式不对的染色体会不识别，加上该参数即可，变为0
$\color{red}{--id-delim}$ ：读取vcf文件时，plink 默认为下划线分隔vcf名称，前面为family，后面为samples
$\color{red}{--double-id}$ ：将familyID和sampleID命名为相同
$\color{red}{--const-fid}$ ：将familyID命名为0，samples-9

$\color{red}{注意}$ , 将vcf变为ped时，map文件的SNP编号那列没有值，可以自行编号添加即可

3. 数据过滤

plink同样可以进行数据过滤，可根据maf，缺失率，ld，哈德温伯格平衡等

plink --bfile Test --mind 0.1 \
  --maf 0.05 --geno 0.1 --hwe 0.01 \
  --make-bed --out clean

## bfile: bed的前缀

主要的参数如下
--mind 过滤掉缺失率为10%的 $\color{red}{样本}$
--geno 过滤掉缺失率为10%的 $\color{red}{SNP}$
--maf 过滤掉maf低于0.05的 $\color{red}{SNP}$
--hwe 基于A和a不符合哈德温伯格平衡

LD的过滤

plink --bfile Test --indep-pairwise 50 5 0.5

结果得到两个文件，一个是plink.prune.in，一个是plink.prune.out，
就是标记名称

下一步从原来的bed文件中调取该位点

plink --bfile Test --extract plink.prune.in --make-bed --out pruneddata

缺失率查看

plink --vcf in.vcf.gz --missing --allow-extra-chr
# 输入vcf或者bed都可以

得到2个文件，plink.imiss，plink.lmiss

plink.imisss是一个样本缺失率

样本缺失率

plink.lmiss 是SNP位点缺失率

SNP位点缺失率

奇怪，我的这个SNP为什么是没有。。。(在*map文件自行添加)

查看MAF频率

plink --vcf in.vcf.gz --freq --allow-extra-chr

得到一个plink.frq 文件

MAF频率

最后一列为：Non-missing allele count

4. 提取（删除）样本（SNP位点）

删除（提取）样本

plink --bfile Test --keep samples --make-bed --out test

# --keep 保留样本
# --remove 删除该样本
# --recode 01 t ranspose. 输出pep文件以及转制后的tped文件，基因型用01表示
# ---output-missing-genotype 9 缺失值用9表示
#  --set-missing-var-ids @:# 输出ID，染色体：坐标

其中samples 的格式为两列
第一列：familey
第二列：sampels
与pep文件的前两列保持一致即可

删除（提取）SNP位点

plink --bfile Test --extract SNP_ID --make-bed --out test

# --extract 提取SNP位点
# --exclude 删除该点

其中SNP_ID 为一列，即SNP ID，与map第二列一样

5、计算样本杂合度

plink --bfile ../final/chr1 --het --out chr1