使用GWAS概括性数据进行条件分析 GCTA – COJO

GCTA-COJO: multi-SNP-based conditional & joint association analysis using GWAS summary data

本文内容:

背景介绍
GCTA-COJO使用方法
参考样本的选择
参考

背景介绍:

通常情况下GWAS或是meta分析检验的是基于单个SNP模型的相关,所估计的效应量为被检验SNP的边际效应,但相比于多SNP模型的联合效应,单个SNP模型的边际效应并没有考虑该SNP与周围SNP的LD

这会带来两个问题:

  1. 如果两个SNP成负相关,那么这两个SNP的效应都会被减弱。
  2. 如果两个SNP都达到了显著的阈值,事后很难通过LD来确定这两个SNP的相关的程度。

如果我们有原始样本层面的基因型,以及单SNP检验后的概括性数据,就可以在没有表型数据的情况下,将边缘效应转化为联合效应。

但一般情况下,对于一个包含很多队列的meta分析,我们无法获取个体的基因型数据,所以不能计算X′X(转化为联合效应时必须的步骤),但X′X本质上是一个SNP基因型的协方差矩阵,我们可以通过1.meta-analysis的基因型频率,与2.从可以获得基因型的参考样本计算而得的SNP之间的LD,来进行估计。这样的基于meta分析的基因频率与参考队列的LD的估计过程也可以用在相应的条件分析之中。(详细推导过程见原文)

使用方法:

COJO内置于GCTA软件中:

输入文件为

  1. 来自meta分析的概括性数据 summary-level statistics from a meta-analysis GWAS
  2. 参考样本的基因组文件 plink格式

-cojo-file [test.ma](<http://test.ma/>) 用于指定概括性数据:

概括性数据输入文件格式如下所示,

SNP A1 A2 freq b se p N 
rs1001 A G 0.8493 0.0024 0.0055 0.6653 129850 
rs1002 C G 0.0306 0.0034 0.0115 0.7659 129799 
rs1003 A C 0.5128 0.0045 0.0038 0.2319 129830
...

需要包含,各列需要按次顺序排列,列名则没有要求

1.SNP的ID,
2.效应等位A1,      #注意 A1,A2不能颠倒
3.参考等位A2,
4.效应等位的频率,   #也就是A1的频率
5.效应的大小,
6.标准差,
7.p值,
8.以及样本大小 

注意

1.对于病例对照研究,效应大小应当是 log(odds ratio),以及其相对应的标准差。

  1. 即使你的研究对象只是一部分SNP,在准备输入文件时,也要准备所有SNP的概括性数据,主要原因是GCTA-COJO需要所有SNP的概括性数据来计算表型方差。可以使用--extract选项来限制COJO分析的区域。

COJO分析中可用的选项:

--cojo-slct

通过逐步法筛选出独立关联的SNP。结果将保存在*.jma文件中,还会有一个额外的文件*.jma.ldr保存SNP两两间的LD。

--cojo-top-SNPs 10

通过逐步法筛选出指定数量的独立关联的SNP,不设P值得阈值,找到指定数量的SNP为止。输出文件同--cojo-slct

--cojo-joint

不进行模型筛选,直接估计所有纳入的SNP的联合效应。输出文件同上。

--cojo-cond cond.snplist

校正给定列表中的SNP后进行关联分析。结果将保存在*.cma文件中。

其中cond.snplist的格式如下所示:

rs1001
rs1002
...

--cojo-p 5e-8

基因组范围显著的阈值,默认为 5e-8 。当使用–cojo-slct时,该选项才有效。

-cojo-wind 10000

当两个SNP的距离超过指定值时,认为其为连锁平衡的。默认值为10000 Kb (i.e. 10 Mb)

-cojo-collinear 0.9

模型选择时,程序会检查SNP之间的共线性问题,如果某一SNP在回归中的R2大于阈值,那么该SNP将不会被选择。默认值为0.9 。

-diff-freq 0.2

检查GWAS概括性数据中频率与参考样本中频率的差值,差值过大的SNP会被去除。默认值为0.2.

-cojo-gc

选择后p值会进行基因组控制。 e.g. –cojo-gc 1.05.

COJO分析的使用例:

# Select multiple associated SNPs through a stepwise selection procedure
gcta64  --bfile test  --chr 1 --maf 0.01 --cojo-file test.ma --cojo-slct --out test_chr1

# Select a fixed number of of top associated SNPs through a stepwise selection procedure
gcta64  --bfile test  --chr 1 --maf 0.01 --cojo-file test.ma --cojo-top-SNPs 10 --out test_chr1

# Estimate the joint effects of a subset of SNPs (given in the file test.snplist) without model selection
gcta64  --bfile test  --chr 1 --extract test.snplist  --cojo-file test.ma --cojo-joint --out test_chr1

# Perform single-SNP association analyses conditional on a set of SNPs (given in the file cond.snplist) without model selection
gcta64  --bfile test  --chr 1 --maf 0.01 --cojo-file test.ma --cojo-cond cond.snplist --out test_chr1

输出文件格式

结尾为.jma的文件 (使用-cojo-slct 或 –cojo-joint时的输出)

Chr SNP bp freq refA b se p n freq_geno bJ bJ_se pJ LD_r
1 rs2001 172585028 0.6105 A 0.0377 0.0042 6.38e-19 121056 0.614 0.0379 0.0042 1.74e-19 -0.345
1 rs2002 174763990 0.4294 C 0.0287 0.0041 3.65e-12 124061 0.418 0.0289 0.0041 1.58e-12 0.012
1 rs2003 196696685 0.5863 T 0.0237 0.0042 1.38e-08 116314 0.589 0.0237 0.0042 1.67e-08 0.0 
...

输出文件的列依次是染色体号,SNP ID,位置,效应等位的频率,效应量大小,标准差,原始P值,估计的有效样本量,参考样本中的频率,joint analysis中该SNP的效应量大小,标准差,P值,以及列表中第i个SNP与第i+1个SNP的LD。

LD矩阵(使用-cojo-slct 或 –cojo-joint时的输出)

存储于后缀为jma.ldr的文件中

SNP rs2001 rs2002 rs2003 ...
rs2001 1 0.0525 -0.0672 ...
rs2002 0.0525 1 0.0045 ...
rs2003 -0.0672 0.0045 1 ...
...

结尾为.cma的文件 (使用–cojo-cond选项会生成)

Chr SNP bp freq refA b se p n freq_geno bC bC_se pC
1 rs2001 172585028 0.6105 A 0.0377 0.0042 6.38e-19 121056 0.614 0.0379 0.0042 1.74e-19
1 rs2002 174763990 0.4294 C 0.0287 0.0041 3.65e-12 124061 0.418 0.0289 0.0041 1.58e-12
1 rs2003 196696685 0.5863 T 0.0237 0.0042 1.38e-08 116314 0.589 0.0237 0.0042 1.67e-08
...

GCTA-COJO 分析中参考样本的选择(重要!)

  1. 如果你的概括性数据来自单一的GWAS,那么最好的参考样本就是这个GWAS中的样本。
  2. 当概括性数据来自Meta分析,个体的基因型无法获取时,可以使用一个样本量较大队列的数据。例如,Wood et al. 2014 Nat Genet 使用了ARIC cohort (data available from dbGaP)
  3. 建议使用样本量大于4000的参考样本
  4. 不建议使用 HapMap 或 1000G 的参考样本,因为其样本量太小。

实际使用中可以使用填补后的基因型数据来估计LD,例如可以从BBJ中随机抽取20000个无关联的个体,提取其填补后的基因型作为EAS群体的参考面板。

参考:

Conditional and joint analysis method: Yang et al. (2012) Conditional and joint multiple-SNP analysis of GWAS summary statistics identifies additional variants influencing complex traits. Nat Genet 44(4):369-375. [PubMed ID: 22426310]

GCTA software: Yang J, Lee SH, Goddard ME and Visscher PM. GCTA: a tool for Genome-wide Complex Trait Analysis. Am J Hum Genet. 2011 Jan 88(1): 76-82. [PubMed ID: 21167468]

使用GCTA (GREML)来估计SNP-遗传力 SNP Heritability

GWAS研究中发现的显著SNP只能解释人类群体中复杂性状很小一部分的遗传变异。那么剩余的遗传力在哪里? 很多时候这部分遗传力并没有丢失,而是由于部分snp效应太小以至于无法达到显著水平而没有被检测到。

与单SNP关联检验相对,GCTA中的GREML(genome-based restricted maximum likelihood)方法使用线性混合模型( linear mixed models , LMMs),将全部SNP的效应作为随机效应拟合。

模型如下:

y是 一个n × 1 的表型向量,n是样本大小; β 是固定效应的向量,诸如性别、年龄、以及一个或多个主成分(PC); u 是SNP效应的向量,服从如下分布:

I 是一个 n × n 的单位矩阵, ɛ 是残余效应的向量,服从以下分布:

W是标准化的基因型矩阵,其中第ijth 个元素为: 

其中 xij 是第j个个体的第i个SNP的参考allele的拷贝数, pi该 allele 的频率。如果我们定义  A = WW/N 并且定义 σ2g 为全部SNP解释的方差 , 即 σ2g=Nσ2u, (N为SNP的数量),那么方程 1 等价于:

其中g是一个n x 1的向量,表示各个个体的全部的遗传效应,服从 g∼N(0,Aσ2g) 。A 可以理解为个体两两间遗传关系组成的遗传关系矩阵 genetic relationship matrix (GRM) 。于是基于这个线性混合效应模型我们就可以通过估计的GRM,利用REML( restricted maximum likelihood )算法来无偏地估计全部SNP解释的方差 σ2g (继而估计SNP遗传力)。

基于以上定义,个体j与个体k的遗传关联可以由以下方程估计:

实际操作中我们首先需要排除掉有亲缘关系的个体,主要原因是该模型的目的是估计所有SNP解释的方差,如果纳入有亲缘关系的个体,会由于表型关联( phenotypic correlations )而造成偏差,例如由于共同的环境因素。即使没有上述的偏差,那么对它的解读也会有别于 “无关联的”个体:一个基于家系的估计值捕捉的是所有因果变异的贡献,而该方法捕捉的则是与被基因分型的SNP处于LD的因果变异的贡献


下面简单介绍如何利用GCTA软件估计snp-遗传力:

GCTA下载链接:https://cnsgenomics.com/software/gcta/#Download

GCTA语法上类似PLINK,文件格式也几乎通用,使用起来十分便捷:

第一步: 估计GRM

gcta64 --bfile test --autosome --maf 0.01 --make-grm --out test --thread-num 10

输入为plink格式的bed/bim/fam文件,

  • –autosome 只是用常染色体上的SNP
  • –maf 0.01 过滤掉maf小于0.01的snp
  • –make-grm 计算GRM
  • –out 输出文件的前缀
  • –thread-num 10 要使用的线程数

计算完成后GRM会存储在以下文件中: test.grm.bin, test.grm.N.bin and test.grm.id

我们也可以去除掉有亲缘关系的个体:

gcta64 --grm test --grm-cutoff 0.025 --make-grm --out test_rm025

第二步:利用GCTA-GREML估计SNP遗传力

gcta64 --grm test --pheno test.phen --reml --out test --thread-num 10
  • –pheno 表型文件
  • –reml 利用reml算法计算snp遗传力

计算结果储存在 test.hsq的文件中

当然我们以可以加入协变量,例如前10个主成分 PC1-10:

gcta64 --grm test --pheno test.phen --reml --qcovar test_10PCs.txt --out test --thread-num 10

对于较大的数据量,我们还可以分染色体进行计算,详见:https://cnsgenomics.com/software/gcta/#Tutorial

参考:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3014363/

Yang, J. et al. Common SNPs explain a large proportion of the heritability for human height. Nature Genetics 42, 565–569 (2010).

https://cnsgenomics.com/software/gcta/#Overview

https://cnsgenomics.com/software/gcta/#Tutorial