bioinfo05-GWAS学习

• 参考：
  • PLINK 1.9 (cog-genomics.org)[1]

前言

发现plink2 和plink 差别还是挺大的，没什么plink2 教程，还是用老版。

1-基本配置

直接下载plink 即可。（2,1 版本的参数差异还挺大的，这里我使用 1.9版本）

ps：本来想尝试一下python 写的hail，但发现软件老是报错。Hail | Install Hail on GNU/Linux[2]

2-输入数据类型

参考：PLINK 1.9 / 2 基本使用方法 （未完工） – GWAS实验室 – GWASLab[3]

常用格式如下。

ped+map

• ped

pedigree file。ped文件一般包含6+2N列，第一至六列分别为 1. Family ID 2. Individual ID 3.Mother ID 4.Father ID 5.Sex 6.Phenotype。第六列以后为各个SNP的等位基因，两列一组，可以使用具体的碱基，也可以使用拷贝数(0,1)。

• map

map，与ped文件相伴随的文件，主要包含ped文件中SNP的位置信息。一般包含4列。分别是1. 染色体号 2.SNP ID 3.遗传图距（单位为摩根或厘摩，通常分析不需要这一列，使用哑值（dummy value） 0 填充） 4.碱基对坐标。每行一个SNP，顺序与ped文件中的SNP相对应。

因为纯文本格式占用大量储存空间，实际操作中尽量使用二进制格式，一组ped/map文件可转换成一组bed/bim/fam文件。

bed+bim+fam

• bed

不同于在基因组比对时，使用的记录位置信息的bed 文件，这里为二进制格式，存储基因型，可以想象成ped文件中除去前6列，剩下基因型数据组成的矩阵。

• fam

纯文本，存储pedigree （ped文件的前六列），也就是样本的一些家系信息。

• bim

纯文本，存储SNP索引 （map文件 + 两列allele信息 ）。

plink 输入

输入主要是上述的ped/map 或 bed/bim/fam文件。

此外，还可以提供phenotype/covariates file (optional)，表型与协变量文件，纯文本，该文件非必须，但表型与协变量通常使用单独的纯文本文件提供（为了准备与使用上的便捷）。该文件通常包含1.family ID 2.Individual ID, 第三列及以后为表型或协变量（常用的如年龄，性别，主成分等）。

我们可以通过-bfile 作为输入文件给plink。

3-质控SNP文件

检查missing比例

这里是如何判断为丢失的呢？

查看每个样本或SNP 各自的missing数目。

plink2 --bfile HapMap_3_r3_1 --missing   

会生成两个文件：

$ ls plink2*
plink2.log  plink2.smiss  plink2.vmiss

其中：

$ head -2 plink2.smiss
#FID    IID     PHENO1  MISSING_CT      OBS_CT  F_MISS
1328    NA06989 N       4203    1457897 0.00288292

$ head -2 plink2.vmiss
#CHROM  ID      MISSING_CT      OBS_CT  F_MISS
1       rs2185539       0       165     0

smiss 对应每个individual，vmiss 对应每个snp。

进行质控

# Delete SNPs with missingness >0.2.
plink2 --bfile HapMap_3_r3_1 --geno 0.2 --make-bed --out HapMap_3_r3_2

# Delete individuals with missingness >0.2.
plink2 --bfile HapMap_3_r3_2 --mind 0.2 --make-bed --out HapMap_3_r3_3

检查sex 分布：

plink --bfile HapMap_3_r3_3 --check-sex 

这个选项plink 2 版本没有了。

# The following two scripts can be used to deal with individuals with a sex discrepancy.
# Note, please use one of the two options below to generate the bfile hapmap_r23a_6, this file we will use in the next step of this tutorial.

# 1) Delete individuals with sex discrepancy.
grep "PROBLEM" plink.sexcheck| awk '{print$1,$2}'> sex_discrepancy.txt
# This command generates a list of individuals with the status “PROBLEM”.
plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6 
# This command removes the list of individuals with the status “PROBLEM”.

# 2) impute-sex.
#plink --bfile HapMap_3_r3_5 --impute-sex --make-bed --out HapMap_3_r3_6
# This imputes the sex based on the genotype information into your data set.

4-关联分析

这里--assoc 即是进行最简单的关联分析的方式：

./plink --bfile test --assoc --pheno pheno_test.txt \
--out result

详细的内容可以参考：(13条消息) 笔记 GWAS 操作流程3：plink关联分析_育种数据分析之放飞自我的博客-CSDN博客_plink关联分析[4]

5-转换文件

将bed+bim+fam 组合格式转换为文本的格式：

./plink --bfile test --recode --out test_new

将记录的 CC, AT 等转换为0,1,2：

./plink --file test_new --recodeA test_0407

就可以将结果转换了：

详细的页面说明参考：Data management - PLINK 1.9 (cog-genomics.org)[5]

参考资料

[1]

PLINK 1.9 (cog-genomics.org): https://www.cog-genomics.org/plink/1.9/

[2]

Hail | Install Hail on GNU/Linux: https://hail.is/docs/0.2/install/linux.html

[3]

PLINK 1.9 / 2 基本使用方法 （未完工） – GWAS实验室 – GWASLab: https://gwaslab.com/2021/03/28/plink-1-9-2-%E4%BD%BF%E7%94%A8%E6%96%B9%E6%B3%95/

[4]

(13条消息) 笔记 GWAS 操作流程3：plink关联分析_育种数据分析之放飞自我的博客-CSDN博客_plink关联分析: https://blog.csdn.net/yijiaobani/article/details/105294261

[5]

Data management - PLINK 1.9 (cog-genomics.org): https://www.cog-genomics.org/plink/1.9/data#recode