分类： 07 LD score regression

popcorn 背景介绍

回顾（在单一族裔中估计遗传相关的方法）：通过Bivariate LD Score regression估计遗传相关性 genetic correlation

随着GWAS等遗传变异关联研究在不同人群中不断开展，相关研究数量不断增加，这为我们提供了宝贵的机会以研究遗传结构在不同人种中的差异，有助于解决这个在医疗研究与群体遗传学中十分重要的问题。

该文章基于Brielin等的研究，简单介绍一种估计跨族裔遗传相关的方法，popcorn，用于估计在不同群体中常见SNP的因果效应量的相关性。该方法纳入了研究中所有SNP的关联信息（区别于其他需要事先进行LD-pruning的方法），而且仅使用GWAS的概括数据，减少了计算负担，并避免了无法获取基因型数据的问题。

通常情况下，在单一群体中，两个表性之间的遗传相关的定义为对于两个表型的SNP效应量大小之间的相关系数。而在多群体的情况时，因为不同群体内等位基因频率有显著差异，定义并不统一。因为某个变异可能在某个特定群体中有较高的效应量但较低的频率，所以本方法不仅考虑了效应量大小的相关，同时也考虑遗传冲击（allelic impact）的相关性。

该文中，跨族裔遗传效应相关的定义为SNP效应量的相关系数，跨族裔遗传冲击相关的定义为特定群体的等位方差正态化（allele-variance-normalized）的SNP效应量之间的相关系数。

直觉上，遗传效应相关性（genetic-effect correlation）衡量的是同一个SNP造成表型改变的程度差异，而遗传冲击相关性（genetic-impact correlation）则是分别在各个群体中，相比于的稀有allele，给予common allele更高的权重。举个例子，如果只看效应相关性，计算得出的相关系数是不完全的，没有反应出两个群体间的效应的完整信息，因为即使效应量的遗传相关是1，但其在两个群体中，对其各自表型分布的冲击也可能由于allele frequency的不同而不同。

所以跨族裔遗传相关的两个核心问题就被引出：

1 同一个变异在两个群体中效应量大小的差异？

2 同一个变异的效应量在两个群体中的差异是否被不同群体的allele frequency所稀释或加重？

为了估计遗传相关，popcorn采用了贝叶斯方法，假设基因型是分别从各个群体内抽取的，并且效应量服从正态的先验分布（无穷小模型）。无穷小模型的假设会产生一个观测统计量（Z scores）的多元正态分布，其中该分布的协方差是遗传力与遗传相关的一个方程。这使得我们可以直接将Z score的膨胀纳入模型中，而不需要进行LD pruning。然后popcorn会最大化近似加权似然方程，来估计遗传力与遗传相关性。

方程推导可参考原文章，这里我们要估计的是以下两个值：

遗传效应相关：pge=Cor(β1，β2)

遗传冲击相关：pgi=Cor(δ1β1，δ2β2)

其中δ为allele的标准差

使用方法：

该文章作者的github仓库： https://github.com/brielin/Popcorn

要求：

popcorn只能使用python2，建议用conda新建一个环境
需要以下python包：
numpy 1.14.2
scipy 1.0.1
pandas 0.22.0
pysnptools 0.3.9
bottleneck 1.0.0
statsmodels 0.8.0
(to use --plot_likelihood) matplotlib 1.5.1

从github clone之后，安装：

cd Popcorn
python setup.py install

可以下载使用千人基因组项目提前计算好的LD分数 （用于EUR-EAS）：

https://www.dropbox.com/sh/37n7drt7q4sjrzn/AAAa1HFeeRAE5M3YWG9Ac2Bta .

popcorn 方法分为两步，

第一步为计算LD分数（使用下载好的数据可以跳过这一步，作者提供了EUR-EAS的LD分数）

第二步为计算遗传相关性

第一步，计算LD分数（可选）

python -m Popcorn compute -v 2 --bfile1 Popcorn/test/EUR_ref --bfile2 Popcorn/test/EAS_ref --gen_effect Popcorn/test/EUR_EAS_test_ge.cscore

-v 表示 输出log的详细程度

—bfile表示 输入的plink bfile文件前缀

—gen_effec 表示计算遗传效应相关 （计算遗传冲击相关时不加这个flag）

最后是输出的路径和文件名

第一步log文件

Popcorn version 0.9.6
(C) 2015-2016 Brielin C Brown
University of California, Berkeley
GNU General Public License v3

Invoking command: python /home/brielin/CoHeritability/Popcorn/__main__.py compute -v 2 --bfile1 Popcorn/test/EUR_ref --bfile2 Popcorn/test/EAS_ref Popcorn/test/EUR_EAS_test.cscore
Beginning analysis at DATE
50000 SNPs in file 1
50000 SNPs in file 2
49855 SNPs in file 1 after MAF filter
44021 SNPs in file 2 after MAF filter
43898 SNPs common in both populations
T1
T2
TC
Analysis finished at DATE

结果的格式如下所示：

1	14930	rs75454623	A	G	0.520874751491	0.58630952381	1.4252539451	1.78489131664	1.49214493091
1	15211	rs78601809	T	G	0.731610337972	0.503968253968	1.20243142615	2.3876817626	1.29403313474
1	15820	rs200482301	T	G	0.728628230616	0.605158730159	1.56556914093	1.42934198599	1.40055911225
1	55545	rs28396308	T	C	0.813121272366	0.804563492063	2.59950434496	2.1809181648	2.06610287872

各列分别代表Chromosome, Base Position, SNP ID, Allele 1, Allele 2, Frequency of A1 in POP1, Frequency of A1 in POP2, LD score in POP1, LD score in POP2, and Cross-covariance score

第二步，计算相关

python -m Popcorn fit -v 1 --use_mle --cfile Popcorn/test/EUR_EAS_test_ge.cscore --gen_effect --sfile1 Popcorn/test/EUR_test.txt --sfile2 Popcorn/test/EAS_test.txt Popcorn/test/EAS_EUR_test_ge_result.txt

–sfile1与–sfile2为输入的两个群体的GWAS文件，注意要与LD分数中POP1和POP2的顺序相一致

输入的GWAS数据至少包含一下列（列名应与下面的介绍保持一致）：

1 rsid 或者 SNP （与LD分数文件中的id保持一致）

注： 可以使用 chr 以及 pos 来替代

2 a1 和 a2

3 N 每个SNP的样本大小

4 beta和SE 或者 OR和p-value 或者 Z分数 （不要求对应a1为effect allele，但要保持所有snp一致即可）

5 AF （可选）即a2的allele frequency，用于A/T ， G/C的SNP的flip

输入文件例如：

rsid	a1	a2	af	N	beta	SE
rs10458597	C	A	0.59421	50000	-0.0152782557183	0.00912151085515
rs3131972	A	C	0.83102	50000	-0.000638379272685	0.0119663325893
rs3131969	A	C	0.8654	50000	-0.000737906302463	0.0131325861467
rs3131967	A	C	0.87044	50000	-0.000189441590945	0.0133506462369

执行代码后的，第二步log文件：

Popcorn version 0.9.6
(C) 2015-2016 Brielin C Brown
University of California, Berkeley
GNU General Public License v3

Invoking command: python /home/brielin/CoHeritability/Popcorn/__main__.py fit -v 1 --cfile Popcorn/test/EUR_EAS_test_ge.cscore --gen_effect --sfile1 Popcorn/test/EUR_test.txt --sfile2 Popcorn/test/EAS_test.txt Popcorn/test/EAS_EUR_test_ge_result.txt
Beginning analysis at DATE
Analyzing a pair of traits in two populations
49999 49999 SNPs detected in input.
43897 SNPs remaining after merging with score file.
43897 SNPs remaining after filtering on AF and self-compliment SNPs.
Initial estimate:
       Val (obs)  SE
h1^2   0.172133 NaN
h2^2   0.086708 NaN
pg     0.302372 NaN
Performing jackknife estimation of the standard error using 200 blocks of size 219 .
Jackknife iter: 200
      Val (obs)        SE          Z     P (Z)
h1^2   0.172133  0.012870  13.374388      0
h2^2   0.086708  0.008555  10.135759      0
pge    0.302372  0.053021  13.157553      0
Analysis finished at DATE
Total time elapsed: T

pge 则为跨族裔遗传相关的估计值 (加–gen_effect时的结果)

pgi 为跨族裔遗传冲击的估计值（不加–gen_effect时的结果）

这里的p值的原假设是估计值等于1，p<0.05则意味着估计值显著小于1

参考：

https://pubmed.ncbi.nlm.nih.gov/27321947/

1 什么是遗传结构？

人类群体遗传学研究中，“遗传结构”（ ‘genetic architecture ）是指 对于某一表型，影响广义表型遗传力（ broad-sense phenotypic heritability ）的遗传变异的总和特征，这是一个整体上的概念。

具体来说，遗传结构包括了 1. 影响某一表性的变异的数量(#)，2. 对表性影响效应的大小（Beta / OR），3. 群体中这些变异的频率（MAF），以及 4. 这些变异之间或是与环境的相互作用 （上位效应 与 G X E）。

2 遗传结构的种类？

遗传结构通常被描述为 单基因的（monogenic）、寡基因的（oligogenic）、多基因的（polygenic）、或是全基因的（omnigenic），表示一个，若干，大量，几乎全部遗传变异都对表型变异性有贡献。

例如，与二型糖尿病相关的变异均为效应量较小的常见变异，但与一型糖尿病相关的变异则既包括效应量较小的常见变异(common variants)，还包括部分有着较大效应量的低平率变异和稀有变异（low-frequency and rare variants）。

但无论某个性状或疾病，有单基因的、寡基因的、多基因的、还是全基因的遗传结构，这些遗传结构对表型的遗传贡献还是存在本质上的变异性（ variability ）。这些变异性可能是表型测量上的不同或不足与真实生物学异质性的一个方程。因此，对于不同疾病，遗传变异的数量或是其他相关属性都会有很大差异。（注意遗传结构的定义。）

3 了解疾病遗传构建的重要性？

遗传结构在疾病的筛查与诊断，增进生物学理解，药物开发，孟德尔随机化以及基因作图方面都有重要的作用。

4 如何估计某一表型的遗传结构：

目前主要的方法包括GWAS与gene-based tests。

5 什么因素会影响遗传结构？

表型（Phenotype.）：不同的表型与遗传变异的相关联的方式不同，与环境的相互作用不同，对不同表型测量的质量的也不尽相同。这些因素都会导致观测到的遗传结构出现差异。

选择（Selection）：选择是指遗传变异的频率在局部环境中向着有利的方向变化的进化过程。 遗传结构可能会受到以下影响：性状本身的性质，解释方差的变异的相对年龄与效应大小，以及一些群体的特征。

去渠道化（Decanalization.）：Canalization 渠道化是指，虽然某一物种保有大量的遗传变异，并且所处环境也有各种变化，但即使在这种情况下，表型也能保持稳定不变的现象，也就是说表型对于遗传与环境的变化是稳健的。去渠道化则是这个渠道化的逆过程，在某种环境的变化，或者某些效应量极大的遗传变异影响下，上述稳健的系统被破坏，这种长期稳态被破坏时会导致疾病的产生。

参考：

1.Timpson, N. J., Greenwood, C. M. T., Soranzo, N., Lawson, D. J. & Richards, J. B. Genetic architecture: The shape of the genetic contribution to human traits and disease. Nature Reviews Genetics 19, 110–124 (2018).

在复杂性状的GWAS研究中，大部分遗传力由并没有达到全基因组显著性水平的SNP关联贡献。然而当前许多利用功能信息与GWAS数据来研究疾病的方法仅仅使用达到显著性水品的loci里的SNP，并且只假设每个loci里只有一个因果SNP，或是完全不考虑LD。基于这些不足，本方法的作者们期望使用所有SNP的信息，并明确的将LD纳入模型，来估计每个功能分类的SNP遗传力，以提升power。

本方法是在基础的LD分数回归上的延伸，可以参考：

• 连锁不平衡分数回归 LD score regression
• 遗传相关 跨性状的LD分数回归（预留链接）

在此之前分割SNP遗传力的方法借助REML（ restricted maximum likelihood ），例如GCTA，但需要个体的基因型的原始数据，并且需要很大的计算资源。所以作者们开发了分层LD分数回归，只需要GWAS的 summary statistics ，以及从与目标群体相对应的参考群体中计算得到的LD信息。

 该方法的核心基于以下的一个事实，  GWAS中某个SNP的χ²  检验统计量包含了所有被该SNP标记的SNP的效应。因此对于一个多基因性状，一个有高LD分数的SNP，相比于低 LD分数的SNP ，总体上也会有较高的 χ²  检验统计量 。原因主要是这些SNP可能标记了单个有很大效应量的SNP，或是多个较弱效果的SNP。（也就是 “以一敌百” 或者 “人多力量大”）

所以我们将所有SNP分为具有不同遗传力的功能分类，那么与某个高遗传力分类存在LD的SNP就会有更高的 χ²  检验统计量 。如果与某个功能分类存在高LD的SNP有较高的 χ²  检验统计量 ，那么就定义这个功能分类有遗传力聚集。

在多基因模型下，某个SNP j 的 χ²  检验统计量 的期望值为

• N为样本大小
• C则是功能分类的索引
• l(j,C)则是SNP j 对于分类C的LD分数

• a则衡量了混淆因素的大小
• τC则表示了每个SNP对于功能分类C遗传力的贡献。

该方程使我们可以估计 τC 的大小（也就是所谓的分割的遗传力）。定义某个分类的聚集为 该分类SNP遗传力的比例除以该分类SNP总数的比例。

该方法的作者们基于多个公开的注释数据库，构建了不针对任何细胞类型的全基线模型 ‘full baseline model’ ， 包括 coding, UTR, promoter and intronic regions the histone marks monomethylation (H3K4me1) and trimethylation (H3K4me3) of histone H3 at lysine 等等。除此之外，还基于 全基线模型 ，构建了多个针对特定细胞类型的模型，包含针对细胞类型的注释等。

尽管Stratified LD score regression提供了一种便捷有效的分析 GWAS的 summary statistics 的方法，但我们也要同时注意该方法的不足之处：

• 为了达到足够的power，需要较大的样本量，或是较大的SNP遗传力，而且性状必须是多基因的
• 该方法需要针对研究群体的LD参考数据
• 该方法目前不支持自定义的array
• 该方法基于加性模型，没考虑上位或非加性的效应
• 该方法依赖于可用的功能注释，如果没有相应注释则无法检测
• 等等

使用方法：

下载，安装，配置环境，数据清理详见

https://alkesgroup.broadinstitute.org/LDSCORE/

从以上链接中我们需要下载以下内容（以欧洲群体为例）：

• 基线模型LD分数 ：baseline.* in 1000G_Phase1_baseline_ldscores.tgz
• 频率 ：1000G.mac5eur.* in 1000G_Phase1_frq.tgz
• 权重 ：weights.* in weights_hm3_no_hla.tgz

解压后即可使用：

python ldsc.py 
	--h2 BMI.sumstats.gz\
	--ref-ld-chr baseline.\ 
	--w-ld-chr weights.\
	--overlap-annot\
	--frqfile-chr 1000G.mac5eur.\
	--out BMI_baseline

• –h2 ： 计算分割的遗传力，参数为之前处理好的gwas summary statistics
• --ref-ld-chr ：下载的参考LD分数文件
• --w-ld-chr： 权重文件
• --frqfile-chr： SNP频率文件
• --overlap-annot： 表示基线模型中功能分类有重叠
• –out：指定出输出文件的前缀

参考：

http://www.github.com/bulik/ldsc

https://github.com/bulik/ldsc/wiki/Partitioned-Heritability

Partitioning heritability by functional annotation using genome-wide association summary statistics https://www.nature.com/articles/ng.3404

关键词：population structure，cryptic relatedness，spurious association，LMM，GRM

早期的GWAS研究使用的模型为固定效应的线性模型，但通常会受两方面混淆因素的影响：

1. 群体结构/分层 （population structure/ stratification）：研究群体中存在有不同祖先（ancestry）的亚群体（subgroup）
2. 隐性关联（cryptic relatedness）：研究样本之间存在未知的亲缘关系

如果在模型中故意忽略掉这些混淆变量（confounding variable），就很可能导致结果出现假阳性（false positive）或是虚假关联（spurious association），所以这是我们在进行GWAS研究时，必须要考虑的问题。

Population structure

首先来看群体结构/分层

考虑一个case-control研究，如下图所示，红色的群体在整体样本中占了case的大多数，那么一些对疾病并没有影响，但在红蓝两群体之间等位基因频率相差很大的genetic marker，就有可能会造成虚假关联（spurious association）。

为了解决这个问题，我们通常会采用genomic control的方法，或是通过PCA来矫正。详见：

• Genomic control (基因组控制与基因组膨胀系数λ Genomic control λGC)
• PCA （群体分层与主成分分析 Population structure & PCA）

Cryptic relatedness

由于将主成分作为协变量纳入模型以矫正群体分层的方法，适用于样本中不存在亲缘关系的情况下使用，而样本中存在亲缘关系时，就不一定有效，但通常我们的研究中都会存在有亲缘关系的样本。

当样本中存在隐性关联 / 错误认定的关联 Cryptic and / or misspecified relatedness时，也可能会造成上述的虚假关联。由于我们通常不能掌握所有研究对象的家谱，所以无法完全去相关个体，例如下面的情况，样本中看似不相关的个体间实际上存在隐性关联：

Linear Mixed Model

为了解决以上问题，相比早期的线性模型，我们就采用一种更为灵活的模型，线性混合模型来进行检验。

固定效应 fixed effects：

• W为 n x （w + 1）的矩阵，包含了截距，以及协变量
• β 则是 ( w + 1 ) × 1 的向量，表示协变量的效应量
• Gs为 n × 1 的向量，某个位点的基因型，每一项的值通常为 0，1，2（等位基因allele的拷贝数）
• γ则是标量，目标位点基因型的效应量

随机效应 random effects：

• g为长度为n的随机向量，表示多基因效应

• δ2g为加性遗传方差
• Ψ是成对遗传相关的矩阵
• e是长度为n的随机向量，服从以下的分布，其中δ2e为非遗传效应造成的方差，我们认为这项在个体间是相互独立的

通常成对遗传相关的矩阵Ψ是未知的，在计算中我们使用empirical GRM（genetic relatedness matrix，通过纳入研究的SNP计算得出）

目前已有多种GWAS检验方法基于LMM模型，包括：

• EMMAX
• GEMMA
• GMMAT
• Bolt-LMM
• fastGWAS
• SAIGE
• 等等

参考：

Lecture 6: GWAS in Samples with Structure ，Summer Institute in Statistical Genetics 2015