GWAS实战教程之解读PLINK的summary结果

在这一期内容中，小陈会带大家简单认识一下PLINK软件输出的GWAS summary结果。相信之前关注公众号的伙伴肯定对GWAS summary数据不陌生，因为它是我们做孟德尔随机化研究的基础，由于不同GWAS分析软件输出的summary结果不太一致，这让很多朋友很是头疼，今天我就以PLINK的输出格式和大家讲解一下，希望借此使大家理解GWAS summary数据。

GWAS summary数据是指反映SNP对表型影响相关信息的一类文件，从这个数据中我们可以知道哪些SNP对该表型有显著影响。

library(data.table)
res =fread("C:/Users/86151/Downloads/plink/myWES_chr2.assoc.logistic",header=T)
colnames(res) ## 查看列名
head(res)

从上图可以看到这个数据有9列，第1列代表染色体信息，第2列代表SNP的ID（绝大多时候是rsid），第3列就是base position，第4列的A1就是effect allele，第5列TEST表示检测的是某一变量与表型的关系，第6列NMISS表示的是样本数（原始样本1619个，经过指控之后只有1617个合格样本），第7列是OR值（只有当表型是二分类变量时才使用，BETA = logOR），第8列是统计检验量，第9列就是P值。

看到这里，很多人就困惑了，为什么没有SE这一列呢？在进行孟德尔随机化分析时我们可以利用OR和P值换算出SE的，在后面我会和大家简单讲解一下。

接下来，我带大家重点解读一下TEST这一列：

table(res$TEST)

利用table()函数我们看到TEST有8个类别，这里的ADD就是指加性模型（additive model）下的SNP对表型的影响，age是指年龄对表型的影响，PC1到PC5是指前5个主成分对表型的影响，sex则是指性别对表型的影响。因此，我们不难看出ADD才是我们真正需要的信息，而其它的都是协变量信息。

这个可以用如下模型解释：

model <- glm(phenotype ~ SNP + age + sex + PC1 + PC2 + PC3 + PC4 + PC5, family = binomial(), data = data)

接下来，我们筛选显著的SNP：

mapping = fread('./mapping_exm.tsv',header=T) # 读取之前的mapping文件
res = merge(res,mapping, by.x='SNP', by.y ='ID') # 将结果文件和mapping文件进行合并
res <- res[order(res$P),] # 将结果的P值从小到大进行排序
res <- res[which(res$TEST=='ADD'&res$P<1.47e-5),] # 挑出Bonferroni矫正后仍然显著的SNP

这里我们看到有6个外显子在Bonferroni矫正后仍然显著与吸烟行为相关（p < 0.05/3398）。

关于GWAS summary数据的解读就简单讲到这儿，希望大家能更深入理解一下！