GWAS入门 – 综述推荐与导读

前言

受推上业内大佬启发,本文将总结对于初学GWAS有较大帮助的综述文章,这些文章多由领域内的leading scientist执笔,引用上千,有较大影响力。对于想快速了解十几年来GWAS发展的同学来说,是不可错过的文章。本文基于Abdel Abdellaoui的推文以及作者个人经验。如有其他推荐,欢迎补充。

综述推荐与导读

第一篇

Hirschhorn, J. N., & Daly, M. J. (2005). Genome-wide association studies for common diseases and complex traits. Nature reviews genetics6(2), 95-108.

https://www.nature.com/articles/nrg1521

首先介绍最早的关于GWAS介绍的Review之一,于GWAS刚刚萌芽的2005年发表,那时人类基因组测序刚刚完成,dbSNP开始建立,Hapmap项目也开始启动,这些项目奠定了GWAS研究发展的基础。这篇综述该介绍了GWAS相比于传统遗传学方法的优缺点,当时可用的测序高通量测序方法,以及GWAS研究中需要注意的核心问题等。可以说是将传统遗传学与现代基因组学衔接的一篇开山之作之一,值得一读。

个人推荐指数:10

第二篇

Balding, D. J. (2006). A tutorial on statistical methods for population association studies. Nature reviews genetics7(10), 781-791.

https://www.nature.com/articles/nrg1916

该综述介绍了早期GWAS研究中可用的的统计学工具。简要的介绍了GWAS研究核心的遗传学与统计学原理 ,并简要梳理了GWAS各个环节上会用到的基础的统计学原理与工具。该文章对于初学者理解GWAS的检验原理有很大帮助,后续的GWAS检验方法基本是基于这些基本原理的扩展与补充,万变不离其宗。

个人推荐指数:8

第三篇

McCarthy, M. I., Abecasis, G. R., Cardon, L. R., Goldstein, D. B., Little, J., Ioannidis, J., & Hirschhorn, J. N. (2008). Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nature reviews genetics9(5), 356-369.

https://www.nature.com/articles/nrg2344

该文发表于2008年,正是第一波GWAS的热潮结果发表后的时期,文中基于第一波GWAS的文章,总结了当时GWAS的研究现状,着重梳理了当时GWAS研究的不足与挑战,为接下来的GWAS研究指出了方向。

个人推荐指数:9

第四篇

Manolio, T. A., Collins, F. S., Cox, N. J., Goldstein, D. B., Hindorff, L. A., Hunter, D. J., … & Visscher, P. M. (2009). Finding the missing heritability of complex diseases. Nature461 (7265), 747-753.

https://www.nature.com/articles/nature08494

寻找复杂疾病的“丢失的遗传力“自始至终都是GWAS研究中的一个热门话题,该文总结了丢失的遗传力可能的来源并给出了可能研究的方法。是一篇较有影响力的文章。

个人推荐指数:8

第五篇

Ioannidis, J., Thomas, G., & Daly, M. J. (2009). Validating, augmenting and refining genome-wide association signals. Nature Reviews Genetics10(5), 318-329.

https://www.nature.com/articles/nrg2544

该文发表于2009,当时GWAS研究已经发现大量的与疾病关联的位点,但这些位点大多都只是真正引起功能改变的因果变异的marker,如何确定在大量的关联中找出因果变异变成了一个不可回避的问题。该文章总结了可以提高GWAS结果可靠性与寻找因果变异的早期方法。

个人推荐指数:7

第六篇

Marchini, J., & Howie, B. (2010). Genotype imputation for genome-wide association studies. Nature Reviews Genetics11(7), 499-511.

https://www.nature.com/articles/nrg2796

该文对基因型插补(genotype imputation)方法进行了总结,介绍了相关的基本概念与常用指标,该文对于理解基因型插补的理论基础有较大帮助。

个人推荐指数:6

第七篇

Price, A. L., Zaitlen, N. A., Reich, D., & Patterson, N. (2010). New approaches to population stratification in genome-wide association studies. Nature reviews genetics11(7), 459-463.

https://www.nature.com/articles/nrg2813

群体分层一直是GWAS研究中一个必须要妥善应对的问题,该文总结了对于群体分层的处理方法。文章较短,但梳理得简洁明了,对于理解 λgc,PCA,线性混合模型等帮助很大。推荐阅读。

个人推荐指数:9

第八篇

Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., McCarthy, M. I., Brown, M. A., & Yang, J. (2017). 10 years of GWAS discovery: biology, function, and translation. The American Journal of Human Genetics101(1), 5-22.

https://www.sciencedirect.com/science/article/pii/S0002929717302409?via%3Dihub

对GWAS问世10年来GWAS研究发展与成果的总结,该文介绍了GWAS的科学基础,并基于大量GWAS研究总结出了一些普遍的结论,同时举出了三个被广泛研究的复杂疾病的典型例子。

个人推荐指数:8

第九篇

Pasaniuc, B., & Price, A. L. (2017). Dissecting the genetics of complex traits using summary association statistics. Nature reviews genetics18(2), 117-127.

https://www.nature.com/articles/nrg.2016.142

该文发表于2017年,随着GWAS的summary statistics不断积累,使用summary statistics的下游分析方法也如雨后春笋般出现,该文总结了使用GWAS summary statistics来对疾病分析的post-GWAS方法,例如gene-based analysis,fine-mapping,以及PRS等。

个人推荐指数:8

第十篇

Tam, V., Patel, N., Turcotte, M., Bossé, Y., Paré, G., & Meyre, D. (2019). Benefits and limitations of genome-wide association studies. Nature Reviews Genetics20(8), 467-484.

https://www.nature.com/articles/s41576-019-0127-1

该文总结了GWAS研究的优势与不足,对于加深对GWAS的理解与了解未来发展方向有较大帮助。

个人推荐指数:7

第十一篇

Uffelmann, E., Huang, Q. Q., Munung, N. S., De Vries, J., Okada, Y., Martin, A. R., … & Posthuma, D. (2021). Genome-wide association studies. Nature Reviews Methods Primers1(1), 1-21.

https://www.nature.com/articles/s43586-021-00056-9

发表于2021年,目前最新的GWAS完整流程讲解,总结较为全面,可以查漏补缺,值得一读。

个人推荐指数:9

参考

RaPID高效检测IBD片段 (IBD segments detection)

本文内容:

RaPID简介
RaPID算法
RaPID教程
注意事项
参考

RaPID简介

RaPID是一款用于在基因组数据中检测IBD(identity-by-descent segments)片段的工具,其特点是应用了PBWT(Positional Burrows-Wheeler Transform)算法,能够快速高效地完成检测。BWT被广泛应用于文件压缩与序列比对等领域,序列比对的BWA软件也利用的是这个算法。

IBD回顾:GWASLab:状态同源 / 血缘同源 IBS/ IBD

RaPID算法示意图

RaPID算法简介

第一步,通过从每个窗口(wi)内随机取一个变异位点,对输入的定像后的基因型进行多重随机投影。

第二步,对每个投影后的面板使用PBWT(Positional Burrows-Wheeler Transform)算法,找出截断值L以上长度的精确匹配序列。

第三步,收集上一步的精确匹配序列,当某一区域精确匹配序列多于一定数量时,就判定为IBD片段。

图中例子各参数为:长度的截断值为10,窗口大小为5,重复次数为4,成功判定需要的次数为2,不足2次的精确匹配被舍弃


RaPID教程

下载地址:

https://github.com/ZhiGroup/RaPID

RaPID输入文件:

  1. 遗产图谱文件格式(tab间隔):

<site_number> <genetic_location>

每一行包含位点索引以及位点的遗传位置,位点的顺序要与输入VCF文件中位点顺序相同。

该文件的位点顺序应该是单调上升的,且不能有缺失。

RaPID也提供了两个python脚本可以用来对异常图谱进行过滤与插值填充:

对 genetic mapping file 进行过滤

python filter_mapping_file.py <genetic_map> <filtered_map>

对遗传位置进行插值填充, 使位点与输入VCF中位点一一匹配

python interpolate_loci.py <filtered_map> <vcf_input_gzip> <output_map_file>

2. 定相过的VCF文件

参考:GWASLab:Eagle2单倍型定相工具 Haplotype phasing

RaPID_v.1.7 版本 具体使用方法

./RaPID_v.1.7 \
	-i <input_file_vcf_compressed> \
	-g <genetic_mapping_file> \
	-d <min_length_in_cM> \
	-o <output_folder> \
	-w <window_size> \
	-r <#runs> \
	-s <#success>

使用作者提供的演示数据:

./RaPID_v.1.7 \
	-i 4k_1e7.vcf.gz \
	-g 4k_1e7_e0.001.g \
	-d 5 \
	-w 250 \
	-r 10 \
	-s 2 \
	-o output_folder

RaPID会将在样本两两之间检测到的所有IBD片段以以下格式输出:

<chr_name> <sample_id1> <sample_id2> <hap_id1> <hap_id2> <starting_pos_genomic> <ending_pos_genomic> <genetic_length> <starting_site> <ending_site>

执行后会得到名为 results.max.gz 的结果

Create sub-samples..
Done!

zcat results.max.gz | head
chr1	0	256	0	0	189	9999752	9.99956	0	94991
chr1	6	1510	1	1	3282417	9616920	6.3345	31000	91249
chr1	18	1801	0	1	189	5186009	5.18582	0	49249
chr1	33	435	0	0	3523968	8873991	5.35002	33250	84249
chr1	42	1854	0	0	1002623	6026274	5.02365	9500	57249
chr1	44	794	1	1	3282417	9797410	6.51499	31000	92999
chr1	51	1394	1	0	3337555	9999752	6.6622	31500	94991
chr1	73	88	1	1	149874	6820265	6.67039	1500	64749
chr1	80	199	0	1	1831845	6974084	5.14224	17500	66249
chr1	82	1849	0	0	3445503	8974078	5.52858	32500	85249

注意事项

如果输入文件的变异密度较低的话,需要相应的减小窗口大小 -w 的值。 例如,对于UKBB,位点密度比示例数据中低了80倍,那么窗口小就应相应的从-w 250减小到-w 3.

参考:

  1. Naseri, Ardalan, Xiaoming Liu, Kecong Tang, Shaojie Zhang, and Degui Zhi. “RaPID: ultra-fast, powerful, and accurate detection of segments identical by descent (IBD) in biobank-scale cohorts.” Genome biology 20, no. 1 (2019): 143;
  2. Ultra-fast Identity by Descent Detection in Biobank-Scale Cohorts using Positional Burrows-Wheeler Transform Ardalan Naseri, Xiaoming Liu, Shaojie Zhang, Degui Zhi bioRxiv 103325;

major/minor/reference/alternative/risk/effect allele 概念解析

这些名词很容易混淆而引起不必要的错误或误解,早期的遗传统计学软件,例如plink并没有很重视allele概念上的明确区分,但近年新出的软件或旧软件的新版本为保证统一性已经开始注意此问题。

本文内容

第一组 频率上的 major 与 minor allele 
第二组 参考基因组的 reference (ref) 与 alternative (alt) allele
第三组 关联检验的 reference (non-risk 或者 non-effect)与 risk/effect allele

首先第一组概念 major 与 minor allele

major allele 与 minor allele 通常针对某一大小确定的特定群体而言,频率最高的allele为该群体的major allele, 频率次高的为 minor allele,对于最常见的bi-allelic SNP来说,两个allele频率一高一低,就是这个群体中这个snp的major和minor allele,对于tri- 或者quad-allelic SNP (位点有三种或四种碱基的SNP)而言,minor allele则是频率第二高的那个allele

注意点:

区分major与minor的依据是 某一大小确定特定群体的 allele 频率

plink1.9目前采用的是major与minor allele的概念,软件会自动计算频率,对原始数据进行操作时会自动改变allele的排序,如果你使用plink1.9 的—frq选项计算频率,你会发现输出的文件中是MAF ,minor allele frequency,不会高于0.5

PLINK1.9中,A1为minor,A2为major allele,所以这里MAF是指A1(minor allele)的频率。。

CHR    SNP    A1   A2          MAF  NCHROBS
1      SNP1    T    C       0.1258    10000
1      SNP2    A    G       0.1258    10000


第二组 reference (ref) 与 alternative (alt) allele

reference allele 在这里是指某一参考基因组上该位点的allele,该位点上其他的allele则称为alternative allele。注意,这里reference 与 alternative allele与频率无关,唯一的决定因素是所选的参考基因组。参考基因组上的allele多为major allele,但这只是巧合,不能以此为依据将major和 reference allele划上等号,也有部分reference allele在该群体中为minor allele。

与plink1.9不同,plink2使用的概念则是reference 与 alternative allele,进行操作时不会自动依据频率而改变ref与alt的排序,使用plink2 的—frq选项计算频率,你会发现输出的文件中是alternative allele frequency (不是MAF),取值范围为[0,1]

PLINK2中则明确区分了reference 与 alternative allele的概念,例如上述的两个SNP,根据参考基因组对齐后,SNP1在参考基因组中的ref为T,那么alt就为C,这里计算的alt的频率为0.8742,按概念来说在该群体中,SNP1的T为ref allele,但却又是minor allele , 而C为alt,却又是major。 对于SNP2来说ref 则为 major,alt 为minor。

#CHROM	ID	REF	ALT	ALT_FREQS	OBS_CT
1	SNP1    	T	C	0.8742	10000
1	SNP2    	G	A	0.1258	10000

小窍门:使用plink2可以将自己手头数据的ref与alt allele与对应参考基因组对齐,示例代码如下:

plink2 \
       --bfile testfile \
       --ref-from-fa -fa hg19.fasta \  从参考基因组的fasta文件来决定plink文件中的ref
       --make-bed \
       --out testfile_fa


第三组 reference 与 risk/effect allele

在这里的概念再次改变,同样的reference allele,在与 risk/effect allele并列时,则指的是GWAS关联检测中的reference allele (non-risk 或者 non-effect),也就是效应量beta(或odds ratio)估计时的参考,概念上与上述ref与alt的组合无关,但为了保持统一性,近年来研究中关联检验的reference 也会与 reference genome保持一致,以避免混淆等。(注意:早期多以minor allele为关联检验的ref allele,这也是容易产生混淆的点)

risk allele 则很好理解,就是对疾病发生有贡献的那个allele,在复杂疾病的研究中,一般情况下risk allele经常为minor allele,但也会有例外。effect allele的概念也类似,就是我们想要研究其对疾病或表型效应的allele,所以通常是对表型或疾病有贡献的allele,关联检验结果中effect一栏指的就是effect allele的效应。

理解了以上概念后,我们在分辨allele时就能得心应手了。

GWAS QC – 估计杂合度并去除异常值 Heterozygosity rate

为了在GWAS分析中得到可靠的结果,我们需要严格进行质量控制,其中很重要的一环就是筛掉杂合性异常的样本,因为通常情况下过高或过低的杂合度一般意味着样本存在污染,或是存在近亲相交。

plink提供了—het的选项,可以通过距估计来估计F系数:

这里的F系数为近交系数:亲缘系数/ 近交系数/血缘系数 coefficient of kinship / inbreeding / relationship 

使用的具体的计算公式是:

注意:

1.当样本量较小时,由于估计的期望偏差可能很大,应当通过 -read-freq 选项来额外提供期望数。

2.该计算方法没有考虑SNP之间的LD,所以在进行估计前最好进行LD-pruning

一般情况下给出的合理范围是平均值正负3个标准差,但实际操作上如果数据来源可靠(例如来自某个标准化biobank,样本制备与测序流程可靠时,可以尝试4个标准差,或是略过此项QC以保留更多样本)

也可以使用—ibc选项(来源于GCTA软件),计算三种基于不同方法的F系数:

第一种,基于加性基因型值的方差:

第二种,基于偏离期望值的纯合性的大小,与上述—het选项相同。

第三种,基于联合配子之间的相关:

下面利用PLINK的-het简单进行演示,示例数据来自https://www.nature.com/articles/nprot.2010.182#MOESM365

# 首先进行LD-pruning
plink \
  --bfile gwa \
  --indep-pairwise 50 5 0.2 \
  --out gwa

#使用LD-prune后的SNP来计算杂合度
plink \
	--bfile gwa \
	--extract gwa.prune.in \
	--het \
	--out gwa

输出文件为 gwa.het

使用head命令查看:

head gwa.het

FID     IID       O(HOM)       E(HOM)        N(NM)            F
   0   A2001        54037    5.433e+04        82915     -0.01037
   1   A2002        54362     5.43e+04        82865     0.002085
   2   A2003        54120    5.435e+04        82932    -0.008133
   3   A2004        54538    5.436e+04        82955     0.006211
   4   A2005        54462    5.435e+04        82932     0.003898
   5   A2006        54387    5.437e+04        82970     0.000653
   6   A2007        54093    5.435e+04        82940    -0.009004
   7   A2008        54239    5.432e+04        82910    -0.002818
   8   A2009        54103    5.437e+04        82982    -0.009274

O(HOM)为观测到的纯合基因型数,E(HOM)为期望值,N为总数。

F = (O-E) / (N-E) 也就是我们要计算的F系数估计值。

首先画图确认分布,这里使用python

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

het = pd.read_csv("gwa.het","\\s+")
plt.figure(figsize=(10,5))
plt.hist(het["F"],bins=100)
plt.title("het F")

然后要剔除平均值3个标准差以外的异常样本,这里还是使用python:

#python
#抽出3个标准差以外的样本的ID
toRemove = het.loc[((np.mean(het["F"])-3*np.std(het["F"]))>het["F"])|(het["F"]>(np.mean(het["F"])+3*np.std(het["F"]))),["FID","IID"]]
#输出要剔除样本的FID与IID
toRemove.to_csv("to_remove.id"," ",header=None,index=None)

在终端中查看是否输出成功:

#在终端中查看
head to_remove.id

207 A2208
483 A2484
752 A2753
913 A2914
1049 A3050
1165 A3166
1356 A3357
1891 A3892
269 A270
393 A394

大功告成,to_remove.id 可以通过—remove选项来剔除出后续分析

参考:

Clarke, G., Anderson, C., Pettersson, F. et al. Basic statistical analysis in genetic case-control studies. Nat Protoc 6, 121–133 (2011). https://doi.org/10.1038/nprot.2010.182

https://www.cog-genomics.org/plink/1.9/basic_stats#ibc

Am J Hum Genet. 2011 Jan 7; 88(1): 76–82.

为什么在PCA或估计GRM时要去除长LD区域 Remove long-LD region

本文内容

1.背景介绍
2.为什么在PCA等时要去除LongLD区域?
3.长LD区域起始位置的列表
4.使用PLINK去除长LD区域里的SNP:

1.背景介绍:

LD :连锁不平衡 linkage disequilibrium LD

PCA :群体分层与主成分分析 Population structure & PCA

2.为什么在QC时要去除LongLD区域?

人类基因组中存在若干长LD的区域,这些区域多位于染色体的着丝粒附近,还有一些位于HLA等区域。如下图所示:

这些区域跨度很长(长度超过2Mb),单次LD-pruning无法完全去除互相成LD的SNP,在进行诸如PCA,或是计算GRM,进行基于LMM模型的GWAS分析时,我们应当去除掉这些区域。

长LD区域的形成并不一定是因为选择,其他原因诸如倒位多态性(inversion polymorphism)也可能造成长LD区域的存在。在进行研究时,应当谨慎区分这些区域形成的原因。如果在计算模型中没有对这些长LD区域进行处理,就可能影响群体遗传结构中对于本地群体的估计,造成系统性的偏倚。

3.长LD区域起始位置的列表(hg38,hg19与hg18参考基因组版本)

hg38版本

Chr	Start	Stop
chr1	47761740	51761740
chr1	125169943	125170022
chr1	144106678	144106709
chr1	181955019	181955047
chr2	85919365	100517106
chr2	87416141	87416186
chr2	87417804	87417863
chr2	87418924	87418981
chr2	89917298	89917322
chr2	135275091	135275210
chr2	182427027	189427029
chr2	207609786	207609808
chr3	47483505	49987563
chr3	83368158	86868160
chr5	44464140	51168409
chr5	129636407	132636409
chr6	25391792	33424245
chr6	26726947	26726981
chr6	57788603	58453888
chr6	61109122	61357029
chr6	61424410	61424451
chr6	139637169	142137170
chr7	54964812	66897578
chr7	62182500	62277073
chr8	8105067	12105082
chr8	43025699	48924888
chr8	47303500	47317337
chr8	110918594	113918595
chr9	40365644	40365693
chr9	64198500	64200392
chr9	88958735	88959017
chr10	36671065	43184546
chr10	41693521	41885273
chr11	88127183	91127184
chr12	32955798	41319931
chr12	34639034	34639084
chr14	87391719	87391996
chr14	94658026	94658080
chr17	43159541	43159574
chr20	4031884	4032441
chr20	33948532	36438183
chr22	30060084	30060162
chr22	42980497	42980522

hg19版本

Chr	Start	Stop ID
1 48000000 52000000 1
2 86000000 100500000 2
2 134500000 138000000 3
2 183000000 190000000 4
3 47500000 50000000 5
3 83500000 87000000 6
3 89000000 97500000 7
5 44500000 50500000 8
5 98000000 100500000 9
5 129000000 132000000 10
5 135500000 138500000 11
6 25000000 35000000 12
6 57000000 64000000 13
6 140000000 142500000 14
7 55000000 66000000 15
8 7000000 13000000 16
8 43000000 50000000 17
8 112000000 115000000 18
10 37000000 43000000 19
11 46000000 57000000 20
11 87500000 90500000 21
12 33000000 40000000 22
12 109500000 112000000 23
20 32000000 34500000 24

hg18版本

Chr	Start	Stop	ID
1	48060567	52060567	hild1
2	85941853	100407914	hild
2	134382738	137882738	hild3
2	182882739	189882739	hild4
3	47500000	50000000	hild5
3	83500000	87000000	hild6
3	89000000	97500000	hild7
5	44500000	50500000	hild8
5	98000000	100500000	hild9
5	129000000	132000000	hild10
5	135500000	138500000	hild11
6	25500000	33500000	hild12
6	57000000	64000000	hild13
6	140000000	142500000	hild14
7	55193285	66193285	hild15
8	8000000	12000000	hild16
8	43000000	50000000	hild17
8	112000000	115000000	hild18
10	37000000	43000000	hild19
11	46000000	57000000	hild20
11	87500000	90500000	hild21
12	33000000	40000000	hild22
12	109521663	112021663	hild23
20	32000000	34500000	hild24
X	14150264	16650264	hild25
X	25650264	28650264	hild26
X	33150264	35650264	hild27
X	55133704	60500000	hild28
X	65133704	67633704	hild29
X	71633704	77580511	hild30
X	80080511	86080511	hild31
X	100580511	103080511	hild32
X	125602146	128102146	hild33
X	129102146	131602146	hild34

4.使用PLINK去除长LD区域里的SNP:

我们可以使用PLINK来去除长LD区域里的SNP,主要分为两步:

1.将上一节中的列表拷贝进high-ld.txt文件中(使用时记得去掉header),使用--make-set选项提取区域中的SNP

2.在分析时利用--exclude去除掉上一步所生成列表中的SNP

plink --file mydata --make-set high-ld.txt --write-set --out hild
plink --file mydata --exclude hild.set --recode --out mydatatrimmed

参考:

https://genome.sph.umich.edu/wiki/Regions_of_high_linkage_disequilibrium_(LD)

Price et al. (2008) Long-Range LD Can Confound Genome Scans in Admixed Populations. Am. J. Hum. Genet. 86, 127-147

更新:

20220905 修改表述错误,更新PCA链接,并增加hg38版本

亲缘系数/ 近交系数/血缘系数 coefficient of kinship / inbreeding / relationship

遗传学中关于亲缘关系的常见的几个系数辨析:

亲缘系数 coefficient of kinship ,有时也称共祖系数 coefficient of coancestry, 是对个体间亲缘关系的直接衡量,定义为从两个个体随机抽取一个同源等位基因,所抽取的两个等位基因是血缘同源(IBD)的概率(即这两个等位基因相同,且来自同一个祖先)。

常用的用来衡量亲缘关系的,且容易混淆的系数还有如下两个,近交系数与血缘系数。

近交系数 coefficient of inbreeding ,由怀特( Wright, Sewall )最早定义,是指 一个个体的某个基因座上的两个等位基因为血缘同源 (IBD) 的概率,衡量的是这个个体父母间亲缘程度的大小,反映的是近亲交配的程度。有时也称为固定系数 fixation index ,通常用 F 表示。

血缘系数 coefficient of relationship ,该系数也由怀特定义,衍生自他对近交系数的定义 是指有共同祖先的两个个体间,基因型一致的概率,通常用 r 来表示。数值上是 亲缘系数 的两倍。

对于常见亲缘关系这三个系数的理论值如下图所示:

与个体的关系亲缘系数 近交系数 F 血缘系数 r
自己,同卵双胞胎1/21/21
亲兄弟姐妹1/41/41/2
母父,儿女*1/41/41/2
祖父母,孙子孙女1/81/81/4
舅舅,舅妈,侄女,侄子1/81/81/4
表兄弟1/161/161/8
同父异母或同母异父的兄弟姐妹1/81/81/4

*注意,虽然近交系数与亲缘系数值相等,当他们的概念并不相同。

亲缘系数是指两个个体间的亲缘关系,如 ( 母父,儿女 )这一项,亲缘系数为1/4的意思是从一个个体与其亲生父母或孩子个随机抽取一个同源等位基因,这两个等位基因相同,且来自同一个祖先 的概率是1/4.

而近交系数在这里的值为虽然也为1/4,但其表达的意思是如果这个个体与其 亲生父母或孩子近交后,子代的某个基因座上的两个等位基因为血缘同源的概率 为1/4.

参考:

Wright, Sewall (1922). “Coefficients of inbreeding and relationship”. American Naturalist. 56 (645): 330–338. doi:10.1086/279872.

群体分化系数 Fst Fixation index

群体分化系数(Fst,Fixation index)是用来衡量两群体间遗传距离的指标,多基于群体的SNP数据来估计。

目前主流的定义有两种,分别基于等位基因频率,或是血缘同源(IBD)。

如果 \bar{p} 是某个等位基因在整个群体里的频率, \sigma _{S}^{2}是等位基因在不同亚群体之间的被群体大小加权后的频率的方差(组间方差),\sigma _{T}^{2}是整个群体的等位基因频率的方差。那么Fst可以被定义为:

F_{{ST}}={\frac  {\sigma _{S}^{2}}{\sigma _{T}^{2}}}={\frac  {\sigma _{S}^{2}}{{\bar  {p}}(1-{\bar  {p}})}}

Wright的定义表示Fst衡量了群体结构可以解释的遗传变异的量。换句话说,衡量的是不属于亚群内多样性的多样性(组间多样性)所占总体多样性的比值,其中多样性通过两个随机抽取的等位基因是不同的概率估计,也就是2p(1-p)。

如果在第i个群体的等位基因频率为pi,相对大小为ci,那么Fst可以表示为:

F_{{ST}}={\frac  {{\bar  {p}}(1-{\bar  {p}})-\sum c_{i}p_{i}(1-p_{i})}{{\bar  {p}}(1-{\bar  {p}})}}={\frac  {{\bar  {p}}(1-{\bar  {p}})-\overline {p(1-p)}}{{\bar  {p}}(1-{\bar  {p}})}}

或者我们可以将Fst表示为:

F_{{ST}}={\frac  {f_{0}-{\bar  {f}}}{1-{\bar  {f}}}}

其中f0是给定两个来自同一亚群体的个体,这两个个体血缘同源(IBD)的概率,

{\bar {f}}则是  给定两个来自总体的个体,这两个个体血缘同源(IBD)的概率。

通过这样的定义,Fst也可以被理解为相比于整体,两个个体在亚群体中相似性的高低。


实践中,Fst定义中所需要的数据一般都很难直接测量,所以通常我们都采用估算的方法。对于DNA序列数据,一个最简单的估计值就是:

F_{{ST}}={\frac  {\pi _{{\text{Between}}}-\pi _{{\text{Within}}}}{\pi _{{\text{Between}}}}}

其中\pi _{{\text{Between}}}\pi _{{\text{Within}}}分别代表两个不同亚群或相同亚群的个体之间,成对等位基因之间不同的平均值(average number of pairwise differences)。 

参考:

https://en.wikipedia.org/wiki/Fixation_index

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

解释复杂疾病的四种主流模型 CDCV/RAME/infinitesimal/Broad-sense-heritability

一,常见疾病-常见变异假说 Common Disease Common Variant Hypothesis

该模型主要是指常见病的易感性主要由中等数量的常见变异引起。此假设有两个关键点,一是常见变异的效应相比于稀有变异应该较小,二是常见变异只有较小效应而常见病又有遗传力的话,那么一定有多个常见变异共同影响疾病易感性。早期的GWAS都是基于这一简单地假设。但这一假设有明显的不足,CDCV无法解释消失的遗传力的问题( missing heritability ),即基于CDCV的GWAS所发现的常见变异只能解释很小一部分推测的遗传力。所以其他模型也开始得到重视以解决这一问题。

二,主要效应的稀有变异模型 The rare alleles of major effect (RAME) model

RAME 模型主要是指常见病病因其实异质性非常高,也就是说,少量MAF<0.01的稀有变异可以促进疾病的发展。主要效应的原位变异可以解释个位数百分比的遗传力。该模型主要聚焦于由于单倍剂量不足或是功能获得性突变而引起的显性效应,这类效应能使得风险上升两倍或者更多。

三,无穷小模型 The infinitesimal model

近年来GWAS研究中无穷小模型逐渐流行。该模型认为复杂疾病的遗传变异是由于大量的,效应很弱(相对风险低于1.2)的变异引起。该模型解释了丢失的遗传力其实大部分是被隐藏了,由于大量对疾病有较弱效应的变异无法在检验中达到预设的显著阈值。目前很多GWAS关联检验方法都基于这一模型。

四,广义遗传力模型:非加性 GxG 与 GxE 相互作用,以及表观遗传效应

Broad sense heritability model: non-additive G×G and G×E interactions and epigenetic effects

广义遗传力模型认为常见变异的效应与稀有变异的效应不足以解释丢失的遗传力。该模型的主要支持依据就是目前在模式生物数量遗传研究中发现的基因型-基因型相互作用(G×G interactions; 也叫上位效应epistasis),以及基因型-环境相互作用(G×E interactions)。除此以外还包括研究越来越多的表观遗传效应,最为显著的就是亲源效应的遗传贡献,与DNA甲基化的继承。

图1,基于4种不同疾病模型的GWAS的特征。

X轴是SNP在染色体上的位置,Y轴是单个SNP所解释的疾病易感性方差的百分数(注意:一般曼哈顿图中y轴是-log10(P))。1,在CDCV模型中,少量的中等效应的loci会产生很强的信号。2,在RAME模型中,少数的稀有变异的因果效应在个别个体中有较强的效应,但不足以解释大部分方差。3,无穷小模型有少数显著的loci,若干在LD区块中的SNP也会显著。4,如果关联只出现在某种环境中(绿色与橙色的信号),那么在一个同时有两种环境混合的群体中,整体效应就会减弱(如箭头所示),只有很少的关联能够被检测,降低了所能解释的方差。

参考:

https://www.nature.com/articles/nrg3118

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

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

易感性尺度遗传力与观测尺度遗传力 Liability scale heritability & observed scale heritability

什么是易感性阈值模型?

参考:易感性-阈值模型 Genetic liability, Threshold model

什么是易感性尺度遗传力,以及观测尺度遗传力?

遗传力通常可以基于不能直接观测的易感性尺度(hl2),或是基于可观测的二分型表型(ho2)来计算。

通常遗传学家更为关注基于易感性尺度的遗传力(Liability scale heritability),主要原因是如果使用观测尺度遗传力(observed scale heritability),容易产生较大的偏差,原因有以下几点:

1 对于数量性状( quantitative traits)来说,测量的尺度与遗传力的表现尺度是一致的(都是连续的),对于二分类性状来说(binary traits), 病例或对照的状态 (case-control status) 是基于0 -1 尺度的,但遗传力在易感性尺度上才是容易解释的。

2 病例的确定。 在病例对照研究中,通常病例所占的比例都要(远)高于人群中的流行率。但病例的确定通常会对遗传力的估计造成偏差。

如下图所示,A为数量变量的易感性分布,B二分类的形状在人群中(即随机)的易感性分布 ,C为病例对照试验中的 易感性分布 ,刻意地过多纳入病例,会引起遗传力估计上的偏差。

3 SNP的质量控制(QC)。相比于数量性状的研究,病例对照实验的QC更为重要。对于数量性状来说,试验或测量上的误差一般不会与表型值相关,但 病例对照实验 中 病例组与对照组一般都是独立采集的,所以试验误差容易造成病例与其他病例更为相似,对照也与其它对照更为相似。这样人工造成的误差会影响基于基因组相似性计算的遗传力。


如何转换?

详细推导过程可以参考下面这篇文章:

Lee SH, Wray NR, Goddard ME, Visscher PM. Estimating missing heritability for disease from genome-wide association studies. Am J Hum Genet. 2011;88(3):294-305. doi:10.1016/j.ajhg.2011.02.002

Lee SH 等人的文章介绍了转换的方法,当case和control并不是从人群中随机抽取的时候,转换公式如下:

hl2为 易感性尺度的遗传力 ,ho2为 观测尺度遗传力

K为人群中病例的比例(流行率)

P为抽取的样本中病例的比例

z则为正态分布的密度函数的在所取阈值t处的值

使用R进行转换

#K = pop prevalence 人群中流行率
#P = proportion of cases in study 样本中病例比例
#h2 = Heritability estimate (on observed scale) 观测尺度遗传力
#T = liability threshold 易感性阈值
#zv = 正态分布的密度函数的在所取阈值t处的值

K=0.0659
P=0.0659
h2=0.0365
zv <- dnorm(qnorm(K))

h2_liab <- h2 * K^2 * ( 1 - K)^2 / P / (1-P) / zv^2
h2_liab

参考文献:

https://www.pnas.org/content/111/49/E5272

Estimating Missing Heritability for Disease from Genome-wide Association Studies

易感性-阈值模型 Genetic liability, Threshold model

易感性-阈值模型是遗传流行病学中重要的理论模型之一。

Liability 易感性

易感性是 遗传因素 环境因素 对某多因子疾病的效应 的总称。因为易感性是一个隐变量(latent variable),实际上很难直接测量某个特定个体的易感性,但我们可以通过某群体中发病个体数量来估计该群体对于疾病的易感性。

Threshold Model 阈值模型

易感性-阈值模型(liability -threshold model)是我们用来分析非孟德尔遗传的(non-Mendelian)分类(categorical)表型(例如二分类表型 binary traits)的基本模型之一,如下图所示,通常我们认为易感性服从一个N(0,1)的正态分布。当某个个体的集聚的易感性超过了阈值T时就会发病,阴影区域的面积表述了在这个群体中该疾病的流行率(prevalence)

对于某个疾病可以有多个阈值,分别对应疾病的不同严重程度。

参考:

https://onlinelibrary.wiley.com/doi/full/10.1002/0470011815.b2a05036