专栏首页育种数据分析之放飞自我笔记 | GWAS 操作流程3:plink关联分析--完结篇

笔记 | GWAS 操作流程3:plink关联分析--完结篇

飞哥注:具体数据文件,查看之前的系列博文。这里使用的是清洗后的数据。

「文件」

HapMap_3_r3_11.bed  HapMap_3_r3_11.bim  HapMap_3_r3_11.fam  HapMap_3_r3_11.log

1. 查看数据

这里,文件是bed文件,二进制不方便查看,我们将其转化为ped文件和map文件

注意,这里我使用的是ped和map格式,如果ped文件中有表型数据(第六列),如果想指定表型数据,用--pheno,包括三列:家系,个体,表型值。

plink --bfile HapMap_3_r3_11 --recode --out test

这里的数据:

  • 基因型个体:110个
  • SNP个数:1073743

2. plink关联分析的类型

❝**参考:**https://www.jianshu.com/p/286050959dbd?utm_campaign=haruki ❞

2.1 阈值性状(1,2)

plink的语境叫“case and control”,即表型值数据是两类数据:1,2,其中0和-9都表示缺失。可以选择的方法有卡方检验和逻辑斯蒂回归(X2关联分析和logistic分析)。

  • 「--assoc」,不允许有协变量
  • 「--logistic」,允许有协变量,如果考虑协变量,速度变慢。比assoc速度慢。

2.2 连续性状(定量性状)

这里的性状是连续性状,也就是除了1,2,0,-9外还有其它数值,「--assoc」会进行T检验(Student's test),还可以用**--linear**进行分析。

  • 「--assoc」,不允许有协变量,速度快
  • 「--linear」,允许有协变量,速度慢

3. 控制假阳性

因为plink进行关联分析时常常面对的是大量的SNP数据,容易产生假阳性,因此需要矫正。

  • Bonferroni,使用0.05/n计算出矫正后的p值作为阈值,其中n为检测SNP的个数。缺点是SNP之间可能存在LD,而且这种方法过于保守,容易产生假阴性。
  • FDR,是一种最小化假阳性预测比例的方法

plink的解决方法是--adjust,生成多种类型的p值。

「文件说明:」

.*.adjusted (basic multiple-testing corrections)
Produced by --adjust.

A text file with a header line, and then one line per set or polymorphic variant with the following 8-11 fields:

CHR	Chromosome code. Not present with set tests.
'SNP'/'SET'	Variant/set identifier
UNADJ	Unadjusted p-value
GC	Devlin & Roeder (1999) genomic control corrected p-value. Requires an additive model.
QQ	P-value quantile. Only present with 'qq-plot' modifier.
BONF	Bonferroni correction
HOLM	Holm-Bonferroni (1979) adjusted p-value
SIDAK_SS	Šidák single-step adjusted p-value
SIDAK_SD	Šidák step-down adjusted p-value
FDR_BH	Benjamini & Hochberg (1995) step-up false discovery control
FDR_BY	Benjamini & Yekutieli (2001) step-up false discovery control
Variants/sets are sorted in p-value order. (As a result, if the QQ field is present, its values just increase linearly.)
  • UNADJ:原始p值
  • GC:基因组矫正P值(依赖加性模型)
  • QQ:P-value的QQ图
  • BONF:Bonferroni 矫正结果
  • HOLM:Holm-Bonferroni (1979) adjusted p-value,矫正结果
  • FDR_BY Benjamini & Yekutieli (2001) step-up false discovery control,FDR方法

4. 阈值性状关联分析

数据:观测值一列是1和2,可以用的方法有:--assoc--logistic

4.1 accoc关联分析未矫正

plink --file test --assoc --out result

「结果文件说明:」

.assoc, .assoc.fisher (case/control association allelic test report)
Produced by --assoc acting on a case/control phenotype.

A text file with a header line, and then one line per variant typically with the following 9-10 fields:

CHR	Chromosome code
SNP	Variant identifier
BP	Base-pair coordinate
A1	Allele 1 (usually minor)
F_A	Allele 1 frequency among cases
F_U	Allele 1 frequency among controls
A2	Allele 2
CHISQ	Allelic test chi-square statistic. Not present with 'fisher'/'fisher-midp' modifier.
P	Allelic test p-value
OR	odds(allele 1 | case) / odds(allele 1 | control)
If the 'counts' modifier is present, the 5th and 6th fields are replaced with:

C_A	Allele 1 count among cases
C_U	Allele 1 count among controls
If --ci 0.xy has also been specified, there are three additional fields at the end:

SE	Standard error of odds ratio estimate
Lxy	Bottom of xy% symmetric approx. confidence interval for odds ratio
Hxy	Top of xy% approx. confidence interval for odds ratio
  • CHR 染色体
  • SNP SNPID
  • BP 物理位置
  • 。。。
  • CHISQ 卡方值
  • P P值 结果文件:

4.2 assoc关联分析矫正p值

 plink --file test --assoc --adjust --out result_adjust

结果生成三个文件:

result_adjust.assoc           result_adjust.log
result_adjust.assoc.adjusted

查看矫正后的值:

4.3 logistic不考虑协变量

plink --file test --logistic --out result_logistic

结果生成:

result_logistic.assoc.logistic  result_logistic.log

4.4 logistic 不考虑协变量矫正p值

plink --file test --logistic --adjust --out result_logistic_adjust

结果:

result_logistic_adjust.assoc.logistic           result_logistic_adjust.log
result_logistic_adjust.assoc.logistic.adjusted

查看矫正值

如果考虑协变量,加上--covar即可。

5. 可视化

「把文件改名字,方便后面代码里作图,这样不用修改代码了。」

cp result.assoc assoc_results
cp result_logistic.assoc.logistic logistic_results.assoc.logistic

注意,logistic里面会有NA,所以我们这里讲NA的行去掉。

awk '!/'NA'/' logistic_results.assoc.logistic >logistic_results.assoc2.logistic

5.1 曼哈顿图

R代码:注意,如果没有安装qqman,就运行代码:install.packages("qqman"),安装即可。

#install.packages("qqman",repos="http://cran.cnr.berkeley.edu/",lib="~" ) # location of installation can be changed but has to correspond with the library location
#library("qqman",lib.loc="~")
library(qqman)
results_log <- read.table("logistic_results.assoc_2.logistic", head=TRUE)
jpeg("Logistic_manhattan.jpeg")
manhattan(results_log,chr="CHR",bp="BP",p="P",snp="SNP", main = "Manhattan plot: logistic")
dev.off()

results_as <- read.table("assoc_results.assoc", head=TRUE)
jpeg("assoc_manhattan.jpeg")
manhattan(results_as,chr="CHR",bp="BP",p="P",snp="SNP", main = "Manhattan plot: assoc")
dev.off()

「assoc的卡方检验结果:」

「--logistic逻辑斯蒂回归结果:」

可以看到两者结果类似,assoc检测到了显著性位点,logistic没有检测到显著性位点。

5.2 QQ plot

R代码:

#install.packages("qqman",repos="http://cran.cnr.berkeley.edu/",lib="~" ) # location of installation can be changed but has to correspond with the library location
#library("qqman",lib.loc="~")
library(qqman)
results_log <- read.table("logistic_results.assoc_2.logistic", head=TRUE)
jpeg("QQ-Plot_logistic.jpeg")
qq(results_log$P, main = "Q-Q plot of GWAS p-values : log")
dev.off()

results_as <- read.table("assoc_results.assoc", head=TRUE)
jpeg("QQ-Plot_assoc.jpeg")
qq(results_as$P, main = "Q-Q plot of GWAS p-values : log")
dev.off()

「--assoc结果:」

「--logistic结果:」

6. 注意

注意,这是阈值性状的结果,分类性状,可以使用assoc和logistic,连续性状的话,如果没有协变量就用assoc,如果有协变量,就用linear即可。

7. 总结

这是使用plink计算GWAS分析的流程,包括数据的清洗,以及建模,以及出结果,以及可视化。

但是,这只是一个开始,后面再研究一下GEMMA,R中的GAPIT的操作方法,最后能够通过编码进行GWAS的分析,才算是掌握了这个分析方法。

GWAS学习笔记系列:

笔记 | GWAS 操作流程1:下载数据

笔记 | GWAS 操作流程2-1:缺失质控

笔记 | GWAS 操作流程2-2:性别质控

笔记 | GWAS 操作流程2-3:MAF过滤

笔记 | GWAS 操作流程2-4:哈温平衡检验

笔记 | GWAS 操作流程2-5:杂合率检验

笔记 | GWAS 操作流程2-6:去掉亲缘关系近的个体

本文分享自微信公众号 - 育种数据分析之放飞自我(R-breeding),作者:邓飞2013

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-05-05

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 笔记 | GWAS 操作流程1:下载数据

    这里,总结一下GWAS的学习笔记,GWAS全称“全基因组关联分析”,使用统计模型找到与性状关联的位点,用于分子标记选择(MAS)或者基因定位,这次学习的教程是p...

    邓飞
  • 如何使用TASSEL l 做GWAS 说明文档

    之前写的Tassel说明文档,虽然我都是使用命令行相关的软件,但是我发现,Linux,命令行对大多数人还是可望而不可即,分享一篇我做的说明文档,用示例数据,一步...

    邓飞
  • 如何使用Tassel 做GWAS 说明文档

    之前写的Tassel说明文档,虽然我都是使用命令行相关的软件,但是我发现,Linux,命令行对大多数人还是可望而不可即,分享一篇我做的说明文档,用示例数据,一步...

    邓飞
  • Python版组合数计算方法优化思路和源码

    总体说明:本文的优化思路并不局限于Python,但C、C++、C#、Java等语言无法使用内置类型直接表示大整数,需要通过数组等特定形式并自己实现大整数乘除法才...

    Python小屋屋主
  • 分别用逻辑回归和决策树实现鸢尾花数据集分类

    Aidol
  • Windows数据库编程接口简介

    数据库是计算机中一种专门管理数据资源的系统,目前几乎所有软件都需要与数据库打交道(包括操作系统,比如Windows上的注册表其实也是一种数据库),有些软件更是以...

    Masimaro
  • Python:爬虫系列笔记(2) -- 基本了解及urllib的使用

    1.什么是爬虫 爬虫,即网络爬虫,大家可以理解为在网络上爬行的一直蜘蛛,互联网就比作一张大网,而爬虫便是在这张网上爬来爬去的蜘蛛咯,如果它遇到资源,那么它就会抓...

    昱良
  • AI 挑战赛 | 基于百度 ApolloSpace 数据集的自动驾驶挑战赛

    美国时间 3 月 8 日,百度方面宣布 Apollo 自动驾驶开放平台正式加入 DeepDrive 深度学习自动驾驶产业联盟,并发布了 Apollo 自动驾驶数...

    朱晓霞
  • pca

    混乱的数据中通常包含三种成分:噪音、旋转和冗余。在区分噪音的时候,可以使用信噪比或者方差来衡量,方差大的是主要信号或者主要分量;方差较小的则认为是噪音或者次要分...

    pydata
  • 30分钟摸透iOS中谓词NSPredicate的来龙去脉

        在现代汉语的解释中,谓词是用来描述或判断客体性质、特征或者客体之间关系的词项。通俗的说,它是描述事物属性的。在iOS开发Cocoa框架中,有提供NSPr...

    珲少

扫码关注云+社区

领取腾讯云代金券