MR-MEGA 跨族裔GWAS荟萃分析 Meta-analysis

MR-MEGA 简介

MR-MEGA (Meta-Regression of Multi-Ethnic Genetic Association 多族裔遗传相关的荟萃回归) 是一款通过跨族裔荟萃回归来检测并精准定位复杂表型关联信号的工具。 该工具利用不同群体的全基因组概括性数据通过MDS方法推导遗传变异的轴。SNP等位的效应则根据其标准差的大小进行加权,加权后的效应就可以纳入线性回归的框架中,同时纳入遗传变异的轴作为协变量。该方法的灵活性使得我们可以将异质性分解为来自祖先以及残差这两个成分,也因此可以提高因果SNP精准定位的精度。

背景

此方法之前跨族裔荟萃分析中常使用的软件为同一研究组开发的MANTRA,相比于传统的固定以及随机效应荟萃分析,MANTRA在存在异质性时能够提高检验的power并提高跨族裔精确定位的精度。但由于MANTRA利用马尔科夫蒙特卡洛法来估计模型参数的后验分布,对于研究数较多的荟萃分析进行计算时,需要大量的计算资源。所以MANTRA的作者们又开发了MR-MEGA的方法,与MANTRA相比,MR-MEGA在保持power与精度的基础上,大幅削减了所需的计算资源。

方法核心原理

考虑对某一表型的 K 个GWAS研究,将第k个GWAS中的第j个SNP的参考等位频率表示为pkj,首先构建一个GWAS两两间欧几里得距离的矩阵,其中各项为由常染色体SNP计算得到,表示为D=[dkk’]。其中:

这个式子里,Ij为一个指示变量,表示第j个变异是否纳入距离计算之中。

接下来对这个距离矩阵D利用MDS(multi-dimensional scaling)方法进行降维,以得到T个遗传变异的轴,第k个GWAS的轴记为Xk。注意,这里遗传变异的轴的个数的选择应当基于GWAS中群体多样性,限制为T≤K-2。

对于第k个GWAS中的第j个变异,分别用bkj与vkj来表示参考等位的效应大小与其方差。于是我们就可以建立一个关于参考等位的效应的线性回归方程,

其中αj是截距,βtj是第j个SNP在第t个轴上的效应。第k个GWAS贡献则通过第j个SNP的倒方差来计算权重,记为v-1kj。于是我们可以将截距α视为各个轴上0所表示的群体中第j个SNP的期望效应值。

以此我们就能得到汇聚后的效应值,以及异质性检验的结果。


使用方法:

下载页面: https://genomics.ut.ee/en/tools/mr-mega

linux系统中下载完成后解压,

unzip MRMEGA_v*.zip #解压
make #编译
./MR-MEGA #测试是否安装成功

make后将MR-MEGA添加进环境(可选)

输入文件:

输入文件 1. 首先是各个GWAS的概括性数据,格式要求如下

每个GWAS的概括性数据都必须有以下的列,并由空格或分隔符分隔:

1) MARKERNAME – snp name
2) EA – effect allele
3) NEA – non effect allele
4) OR - odds ratio
5) OR_95L - lower confidence interval of OR
6) OR_95U - upper confidence interval of OR
7) EAF – effect allele frequency
8) N - sample size
9) CHROMOSOME  - chromosome of marker
10) POSITION - position of marker

当表型是数量型表型时,使用beta:

  1. BETA – beta
  2. SE – std. error

以下的列可选,默认是所有的SNP都在正链上:

  1. STRAND – marker strand (if the column is missing then program expects all markers being on positive strand)

例如:

MARKERNAME STRAND CHROMOSOME POSITION IMP EA NEA EAF N BETA SE
rs12565286 + 1 761153 0 G C 0.3 1200 -0.02 0.0403
rs2977670 + 1 763754 0 C G 0.23 1200 -0.01 0.40612
rs12138618 + 1 790098 0 G A 0.97 1200 -0.07 0.37
rs3094315 + 1 792429 0 G A 0.01 1199 0.0258 0.1012
rs3131968 + 1 794055 0 G A 0.27 1200 -0.373 0.0101
rs2519016 + 1 805811 0 T C 0.04 1200 0.26 0.3472
rs12562034 + 1 808311 0 G A 0.65 1200 0.0092 0.2

输入文件 2, 除此之外我们还要准备一个存储GWAS概括性数据文件名的列表,(默认名称为 mr-mega.in)

Pop1.txt.gz
Pop2.txt.gz
Pop3.txt.gz
Pop4.txt.gz
Pop5.txt.gz
Pop6.txt.gz
Pop7.txt.gz
Pop8.txt.gz

MR-MEGA基本语法

./MR-MEGA  [--name_pos <string>] ...  [--name_chr <string>] ... 

              [--name_n <string>] ...  [--name_strand <string>] ... 

              [--name_or_95u <string>] ...  [--name_or_95l <string>] ... 

              [--name_or <string>] ...  [--name_se <string>] ... 

              [--name_beta <string>] ...  [--name_eaf <string>] ... 

              [--name_nea <string>] ...  [--name_ea <string>] ... 

              [--name_marker <string>] ...  [-f <string>] ...  [--pc <int>]

              [-t <double>] [--no_std_names] [--debug] [--qt] [--gco]

              [--gc] [--no_alleles] [-m <string>] [-o <string>] [-i

              <string>] [--] [--version] [-h]

—name_xxx都是用来制定自定义列名

-f 可以对列进行过滤 例如 INFO>0.4

-pc <int>指定纳入回归的pc数量,默认为4,但要注意,PC数必须小于研究数-2,如果有5个研究,那么能够纳入PC的最大数量是2.

-qt 数量表型

-o <string> 指定输出前缀

-i <string> 制定输入的文件列表

例:

./MR-MEGA -i mr-mega.in -o mr-mega-out —pc 2

结果解读

MR-MEGA成功运行后会生成两个后缀分别为*.result 和 *.log的文件:

result文件中包含以下列:

#基本信息
MarkerName - unique marker identification across input files
Chromosome - chromosome of marker
Position - physical position in chromosome of marker
EA - allele, which effect was measured across input files
NEA - other allele
EAF - average effect allele frequency (weighted by the samplesize of each input file)
Nsample - total number of samples
Ncohort - total number of cohorts, where the marker was present

#效应
Effects - effect direction across cohorts (+ if the effect allele effect was positive, - if negative, 0 if the effect was zero, ? if marker was not available in cohort)
beta_0 - effect of first PC of meta-regression
se_0 - stderr of the effect of first PC of meta-regression
(beta_1)
(se_1)
(...)

#统计量与p值
chisq_association - chisq value of the association
ndf_association - number of degrees of freedom of the association
P-value_association - p-value of the association

#异质性与p值
chisq_ancestry_het - chisq value of the heterogeneity due to different ancestry
ndf_ancestry_het - ndf of the heterogeneity due to different ancestry
P-value_ancestry_het - p-value of the heterogeneity due to different ancestry
chisq_residual_het - chisq value of the residual heterogeneity
ndf_residual_het - ndf of the residual heterogeneity
P-value_residual_het - p-value of the residual heterogeneity

lnBF - log of Bayes factor
Comments - reason why marker was not analysed

参考:

https://genomics.ut.ee/en/tools/mr-mega

Reedik Mägi, Momoko Horikoshi, Tamar Sofer, Anubha Mahajan, Hidetoshi Kitajima, Nora Franceschini, Mark I. McCarthy, COGENT-Kidney Consortium, T2D-GENES Consortium, Andrew P. Morris, Trans-ethnic meta-regression of genome-wide association studies accounting for ancestry increases power for discovery and improves fine-mapping resolution, Human Molecular Genetics, Volume 26, Issue 18, 15 September 2017, Pages 3639–3650, https://doi.org/10.1093/hmg/ddx280

通过Bivariate LD Score regression估计遗传相关性 genetic correlation

背景回顾:

LD score 回归: 连锁不平衡分数回归 LD score regression

双变量/跨性状LD分数回归(Bivariate/Cross-trait LD Score regression)

简介与原理:

该方法使用GWAS的概括性数据来估计一对表型的遗传相关。该方法基于以下的事实:GWAS检验中,对一个SNP效应量的估计通常也会包含与该SNP成LD的其他SNP的效应,也就是说一个与其他SNP成高LD的SNP,通常也会有更高的卡方检验量。而当我们把卡方检验量替换为,来自两个相关表型的GWAS的Z分数的乘积时,这个事实依然存在。所以基于此,我们可以将单表型的LD分数回归扩展成如下所示的 双变量/跨性状LD分数回归

详细推导见:

Bulik-Sullivan, B., Finucane, H., Anttila, V. et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 47, 1236–1241 (2015). https://doi.org/10.1038/ng.3406


下面简单介绍使用方法:

数据下载,数据清理等步骤与 单表型的LD分数回归基本相同,详见 连锁不平衡分数回归 LD score regression

我们只需要将 单表型的LD分数 中 -h2 的选项改为 -rg, 再加上要估计遗传相关的两个表型的数据清理后的文件即可:

ldsc.py \
--rg scz.sumstats.gz,bip.sumstats.gz \
--ref-ld-chr eur_w_ld_chr/ \
--w-ld-chr eur_w_ld_chr/ \
--out scz_bip
--rg :计算遗传相关性 (genetic correlation),参数为上一步处理好的数据文件名
--ref-ld-chr :使用的按染色体号分类的LD score文件名,参数为LD score文件所在文件夹的路径。默认情况下ldsc会在文件名末尾添加染色体号,例如 --ref-ld-chr eur_w_ld_chr/ 的意思就是使用 eur_w_ld_chr/1.l2.ldscore ... eur_w_ld_chr/22.l2.ldscore 这些文件。如果你的染色体号在其他位置,也可以使用@来告诉LDSC,例如 --ref-ld-chr ld/chr@ 。当然也可以用 --ref-ld 来指定一个整体的LD score文件。
--w-ld-chr:指定ldsc回归权重所用的LDscore文件,理论上对于SNPj的LD分数,应当包含这个SNP与其他所有SNP的R2之和,但实际操作中,LD score回归对于计算权重的SNP的选择并不敏感,所以一般情况下我们可以使用与 --ref-ld 与相同的文件。.
--out : 输出文件的路径与前缀

运行后我们就能得到如下结果:

Genetic Covariance
------------------
Total Observed scale gencov: 0.3644 (0.0368)
Mean z1*z2: 0.1226
Intercept: 0.0037 (0.0071)

Genetic Correlation
-------------------

Genetic Correlation: 0.6561 (0.0605)
Z-score: 10.8503
P: 1.9880e-27

Summary of Genetic Correlation Results
              p1               p2     rg    se       z          p  h2_obs  h2_obs_se  h2_int  h2_int_se  gcov_int  gcov_int_se
 scz.sumstats.gz  bip.sumstats.gz  0.656  0.06   10.85  1.988e-27   0.522      0.053   1.001      0.009     0.004        0.007

Genetic Correlation: 0.6561 (0.0605) 即为我们需要估计的遗传相关

除此之外,LDSC还会输出一个Summary of Genetic Correlation Results,当有多组遗传相关需要估计时,查看该表格会更加方便。

常见问题:

为什么rg会在[-1,1]的范围以外: (https://github.com/bulik/ldsc/issues/78)

ldsc的算法并没有限定[-1,1]的范围,当rg接近一时,加上估计误差有时就会导致rg在范围外

参考:

Bulik-Sullivan, B., Finucane, H., Anttila, V. et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 47, 1236–1241 (2015). https://doi.org/10.1038/ng.3406

https://github.com/bulik/ldsc/wiki/Heritability-and-Genetic-Correlation

相关链接:

基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression 

遗传结构 Genetic architecture

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 如何估计某一表型的遗传结构:

目前主要的方法包括GWASgene-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).

Variant Normalization 变异的标准化

背景

VCF文件是一种灵活的,可以表示包括SNPs, indels 和CNV等多种变异类型的文件格式,但是VCF文件中变异的表示方式并不是唯一的,对同一变异的表示方式上的混淆会对后续分析造成影响。所以我们需要对这样的变异进行标准化,使其拥有唯一的(unique)表示方式。

变异标准化的定义

变异表示的标准化包含两个部分:

1,节俭原则 2,左对齐。

这两部分分别对应变异的长度与位置。


节俭原则 Parsimony

在表示一个变异(variant )时,节俭原则(parsimony )是指,在不使任何等位长度为0的情况下,用最少的核苷酸数来表示这个变异。

所以节俭原则意味着两点:

1 尽可能的用最少的核苷酸数来表示

2 但等位(allele)的长度至少为1不能为0

当一个变异的两个等位最左侧的核苷酸相同时,且删去这个相同的核苷酸后不会造成等位长度为零时,我们称这个变异在左侧是有冗余核苷酸的。

我们来看下面这个例子,前三组多核苷酸多态(Multi Nucleotide Polymorphism ,MNP)都有冗余的核苷酸,最后一组蓝色的变异的则负荷上述的节俭原则。

左对齐 Left alignment

左对齐的意思是,将变异起始位置向左移动到不能在移动为止。这是一个与插入和缺失变异相关的概念,目的是明确变异的位置(相对于明确长度的节俭原则)。

左对齐的定义如下:

当且仅当在保持其等位长度一定的情况下,无法再向左移动起始位置时,这个变异是左对齐的。

我们来看下面这个例子(一个未标准化的短重复序列):

红色:替代等位为空,不符合节俭原则

绿色:没有左对齐,但节俭

橙色:左对齐但不右节俭

蓝色:左对齐但不左节俭

紫色:左对齐,而且节俭 √

如何进行变异标准化:

目前有多种软件支持对VCF文件进行标准化,我们只用提供参考基因组以及输入文件就可以轻松的进行标准化:

例如,使用常用的bcftools进行标准化,命令如下所示:

bcftools norm -f ref.fa in.vcf -O z > out.vcf.gz

参考:

https://genome.sph.umich.edu/wiki/Variant_Normalization#Left_alignment

使用ANNOVAR 对Variants进行功能注释 Annotation POST-GWAS analysis

在GWAS分析中,或是进行针对rare variants的gene-based test时,对SNP进行功能注释是必不可少的一个步骤,本文将简单介绍一款最为常用的SNP注释软件 ANNOVAR url: https://annovar.openbioinformatics.org/en/latest/

ANNOVAR是一款对基因变异进行功能注释的软件,可以对多种生物的变异进行注释(包括 human genome hg18, hg19, hg38, 以及 mouse, worm, fly, yeast 等)。 给定一个包含染色体,起点,终点,参考核苷酸与观测核苷酸, ANNOVAR可以进行如下的功能注释:

  1. 基于基因的注释 Gene-based annotation:主要针对SNP或CNV是否引起蛋白编码改变进行注释,可以灵活选用 RefSeq genes, UCSC genes, ENSEMBL genes, GENCODE genes, AceView genes等多种不同来源的基因定义系统。
  2. 基于区域的注释 Region-based annotation:针对基因组某一特定区域的变异进行注释,例如44个物种的保守区域,预测的转录因子结合位点,GWAS hit, ENCODE H3K4Me1/H3K4Me3/H3K27Ac/CTCF sites,ChIP-Seq peaks, RNA-Seq peaks等
  3. 基于筛选的注释 Filter-based annotation:使用某一特定的数据库进行筛选注释,例如注释变异的rs id,1000基因组项目中的MAF,或是ExAC、gnomAD等,再例如SIFT/ PolyPhen/ LRT/ MutationTaster/ MutationAssessor/ FATHMM/ MetaSVM/ MetaLR 分数等。

实战操作

ANNOVAR是由perl编写的程序,首先通过下载页面,填写注册表格(你还可以给作者留言,作者貌似指着这个留言写基金申请,靠用户感动评审~),学术用途可以免费试用,几分钟后邮件会发来下载链接。

下载并解压(将路径添加到环境中),

我们还需要下载注释用的参考数据库:

目前可用的参考库:https://annovar.openbioinformatics.org/en/latest/user-guide/download/

注意:annovar只能使用官方提供的注释数据库

对于初学者来说,使用ANNOVAR最简单的方法就是使用 table_annovar.pl 程序。 我们通过以下例子来解释annovar的用法,

# 1.下载注释用数据库
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/
annotate_variation.pl -buildver hg19 -downdb cytoBand humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar exac03 humandb/ 
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar avsnp147 humandb/ 
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar dbnsfp30a humandb/

# 2.使用table_annovar.pl对 example/ex1.avinput输入文件进行注释
table_annovar.pl example/ex1.avinput humandb/ -buildver hg19 -out myanno -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f -nastring . -csvout -polish -xref example/gene_xref.txt

如上所示,第一部分代码,下载了对应hg19的各种数据库到humandb/ 文件夹(可以自定义)里,数据库内容可以参考以下链接:https://annovar.openbioinformatics.org/en/latest/user-guide/filter/

包含了对各个数据库详细的介绍:

第二段代码使用table_annovar.pl主程序对example/ex1.avinput输入文件同时进行多数据库的注释,

  • -buildver hg19 : 参考基因组使用hg19
  • -out myanno : 输出前缀为 myanno
  • -remove : 注释完成后删除缓存文件
  • -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a : 所用的数据库包括 ExAC version 0.3 (exac03) dbNFSP version 3.0a (dbnsfp30a), dbSNP version 147 with left-normalization (avsnp147) 数据库
  • -operation gx,r,f,f,f :指定针对protocol的操作(与-protocol一一对应),可选的操作包括g (gene-based), gx (gene-based with cross-reference annotation from -xref argument), r (region-based) 以及 f(filter-based).
  • -nastring “.” : 缺失的注释用 “.”替代
  • -csvout 输出为csv格式
  • -xref example/gene_xref.txt 交叉引用文件,例如 已知的由这个基因变异引起的疾病

其中annovar的标准输入格式.avinput的具体内容如下:

第一至五列为必须,分别是染色体号,起始位点,结束位点,参考等位基因 reference allele 以及 替代等位基因 alternative allele,第五列之后可自由添加所需要的信息

1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
1 5935162 5935162 A T comments: rs1287637, a splice site variant in NPHP4
1 162736463 162736463 C T comments: rs1000050, a SNP in Illumina SNP arrays
1 84875173 84875173 C T comments: rs6576700 or SNP_A-1780419, a SNP in Affymetrix SNP arrays
1 13211293 13211294 TC - comments: rs59770105, a 2-bp deletion
1 11403596 11403596 - AT comments: rs35561142, a 2-bp insertion
1 105492231 105492231 A ATAAA comments: rs10552169, a block substitution
1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
16 50745926 50745926 C T comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
16 50756540 50756540 G C comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss

当然,除了自有的avinput格式外,ANNOVAR还支持VCF等多种常用格式输入文件(-vcfinput)。

table_annovar.pl example/ex2.vcf humandb/ -buildver hg19 -out myanno -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation g,r,f,f,f -nastring . -vcfinput -polish

另外,我们也可以利用ANNOVAR的核心程序 annotate_variation.pl,快速简便的完成单一类型的注释

# 基于基因
annotate_variation.pl -geneanno -dbtype refGene -buildver hg19 example/ex1.avinput humandb/
#基于区域
annotate_variation.pl -regionanno -dbtype cytoBand -buildver hg19 example/ex1.avinput humandb/ 
#基于筛选
annotate_variation.pl -filter -dbtype exac03 -buildver hg19 example/ex1.avinput humandb/

以基于基因的注释为例(用法三者类似),

第一步下载数据库 refGene

annotate_variation.pl -downdb -buildver hg19 -webfrom annovar refGene humandb/
NOTICE: Web-based checking to see whether ANNOVAR new version is available ... Done
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg19_refGene.txt.gz ... OK
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg19_refLink.txt.gz ... OK
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg19_refGeneMrna.fa.gz ... OK
NOTICE: Uncompressing downloaded files
NOTICE: Finished downloading annotation files for hg19 build version, with files saved at the 'humandb' directory

第二步 注释 (这里没有指定注释的数据库,因为 annotate_variation.pl 默认的参数是 –geneanno -dbtype refGene)

annotate_variation.pl -out ex1 -build hg19 example/ex1.avinput humandb/
NOTICE: The --geneanno operation is set to ON by default
NOTICE: Reading gene annotation from humandb/hg19_refGene.txt ... Done with 48660 transcripts (including 10375 without coding sequence annotation) for 25588 unique genes
NOTICE: Reading FASTA sequences from humandb/hg19_refGeneMrna.fa ... Done with 14 sequences
WARNING: A total of 333 sequences will be ignored due to lack of correct ORF annotation
NOTICE: Finished gene-based annotation on 15 genetic variants in example/ex1.avinput
NOTICE: Output files were written to ex1.variant_function, ex1.exonic_variant_function

第三步 查看结果

除了log文件外,还有三个后缀分别为.variant_function ,.exonic_variant_function 以及 .invalid_input 的文件生成:

文件1 .variant_function

cat ex1.variant_function 
UTR5 ISG15(NM_005101:c.-33T>C) 1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
UTR3 ATAD3C(NM_001039211:c.*91G>T) 1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
splicing NPHP4(NM_001291593:exon19:c.1279-2T>A,NM_001291594:exon18:c.1282-2T>A,NM_015102:exon22:c.2818-2T>A) 1 5935162 5935162 A T comments: rs1287637, a splice site variant in NPHP4
intronic DDR2 1 162736463 162736463 C T comments: rs1000050, a SNP in Illumina SNP arrays
intronic DNASE2B 1 84875173 84875173 C T comments: rs6576700 or SNP_A-1780419, a SNP in Affymetrix SNP arrays
intergenic LOC645354(dist=11566),LOC391003(dist=116902) 1 13211293 13211294 TC - comments: rs59770105, a 2-bp deletion
intergenic UBIAD1(dist=55105),PTCHD2(dist=135699) 1 11403596 11403596 - AT comments: rs35561142, a 2-bp insertion
intergenic LOC100129138(dist=872538),NONE(dist=NONE) 1 105492231 105492231 A ATAAA comments: rs10552169, a block substitution
exonic IL23R 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
exonic ATG16L1 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
exonic NOD2 16 50745926 50745926 C T comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
exonic NOD2 16 50756540 50756540 G C comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
exonic NOD2 16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
exonic GJB2 13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
exonic CRYL1,GJB6 13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss

第一列就表示这个变异的功能注释,是位于外显子还是内含子等等,可能的注释结果如下图所示:

第二个文件,.exonic_variant_function:包含了外显子上的变异所引起的具体的氨基酸变化。如下所示,

cat ex1.exonic_variant_function 
line9 nonsynonymous SNV IL23R:NM_144701:exon9:c.G1142A:p.R381Q, 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
line10 nonsynonymous SNV ATG16L1:NM_001190267:exon9:c.A550G:p.T184A,ATG16L1:NM_017974:exon8:c.A841G:p.T281A,ATG16L1:NM_001190266:exon9:c.A646G:p.T216A,ATG16L1:NM_030803:exon9:c.A898G:p.T300A,ATG16L1:NM_198890:exon5:c.A409G:p.T137A, 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
line11 nonsynonymous SNV NOD2:NM_022162:exon4:c.C2104T:p.R702W,NOD2:NM_001293557:exon3:c.C2023T:p.R675W, 16 50745926 50745926 C comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
line12 nonsynonymous SNV NOD2:NM_022162:exon8:c.G2722C:p.G908R,NOD2:NM_001293557:exon7:c.G2641C:p.G881R, 16 50756540 50756540 G comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
line13 frameshift insertion NOD2:NM_022162:exon11:c.3017dupC:p.A1006fs,NOD2:NM_001293557:exon10:c.2936dupC:p.A979fs, 16 50763778 5076377comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
line14 frameshift deletion GJB2:NM_004004:exon2:c.35delG:p.G12fs, 13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
line15 frameshift deletion GJB6:NM_001110221:wholegene,GJB6:NM_001110220:wholegene,GJB6:NM_001110219:wholegene,CRYL1:NM_015974:wholegene,GJB6:NM_006783:wholegene, 13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss

可能的取值如下所示,如果我们需要进行gene-based test来找到rare variants,那我们就可以依据下表选中我们想要纳入研究的variants,制作相应的group file。

第三个文件, .invalid_input 则是注释失败的输入文件中的变异。

参考:

https://www.nature.com/articles/nprot.2015.105

https://annovar.openbioinformatics.org/en/latest/user-guide/gene/

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

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

血缘系数 coefficient of kinship / kinship coefficient,有些地方称为近亲系数,有时也称共祖系数 coefficient of coancestry, (目前缺少权威的中文翻译),是对个体间血缘关系的直接衡量,定义为从两个个体随机抽取一个同源等位基因,所抽取的两个等位基因是血缘同源(IBD)的概率(即这两个等位基因相同,且来自同一个祖先)。

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

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

亲缘系数 coefficient of relationship / coefficient of relatedness ,该系数也由怀特定义,衍生自他对近交系数的定义 是指有共同祖先的两个个体间,基因型一致的概率,通常用 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.

SNP的LD剪枝与聚集 LD pruning & clumping

GWAS相关的研究中,很多时候我们需要从总的SNP数据中,基于SNP两两之间的LD,来抽取出一个不含互相关联SNP的子集,目前主要的两种方法分别是 LD pruning 与 clumping。

例如,

在进行主成分分析(PCA)时,我们需要事先对SNP进行LD pruning 以去除互相关联(处于LD)的SNP,以防止高LD区域过高的方差对结果的影响。

在计算风险分数PRS时,我们需要从显著的loci中选取具有代表性的SNP来计算分线分数,这时就需要进行clumping,基于LD的r2,以及GWAS所得到的p值,来筛选出这个LD区域中的代表SNP(重要性最高),这样我们可以获得更准确的风险分数。

LD pruning 与 clumping 方法的异同如下所示:

根据保留主要用途
PruningLD的R2处于LD的一对SNP中MAF最高的PCA
ClumpingLD的R2 与 SNP的P值 处于 LD的一对SNP中P值最显著的 PRS
Pruning 与 clumping 的主要区别

具体算法上,可以简单理解为:

Pruning:选取第一个SNP,然后计算这个SNP与窗口区间里第二个,第三个,等等的r2,当检测到高的相关性时,就会从这一对SNP中去除MAF较低的那个,保留 MAF 高的,也就是说这个过程中可能会去除掉我们选的第一个SNP。完成后下一步就是选取下一个SNP,重复这个过程。

Clumping:首先会依据从GWAS得到的p值对SNP的重要性进行排序,然后选取排序后的第一个SNP, 计算这个SNP与 窗口区间里 其他SNP的r2, 当检测到高的相关性时,就会从这一对SNP中去除重要性低的那个, 这个过程中我们选的第一个SNP一定会得到保留。 完成后下一步就是选取 p值 排序后的下一个 SNP,重复这个过程。


PLINK中提供了 Pruning 和 Clumping 的功能:

Pruning:

我们主要是用–indep-pairwise选项,也就是根据SNP两两之间的LD来pruning。

--indep-pairwise  <window size>['kb']  <step size (variant ct)>  <r^2 threshold>

例 --indep-pairwise 500 50 0.2
这三个参数代表的意思分别是: 窗口大小,每一步移动窗口的距离,以及判定关联的r2阈值
plink -bfile input --indep-pairwise 500 50 0.2 --out input_pruned

输出两个文件
input_pruned.prune.in    #pruning后保留的互不相关的SNP
input_pruned.prune.out  #去除掉的SNP

Clumping:

PLINK提供了多种参数选项,具体可以参考:https://www.cog-genomics.org/plink/1.9/postproc

参考:

https://www.cog-genomics.org/plink/1.9/postproc

https://www.cog-genomics.org/plink/1.9/ld

https://www.biostars.org/p/343818/

LiftOver 基因组坐标变换

本文内容:

  1. liftover 简介与网页版工具
  2. liftover 命令行工具使用方法
  3. 基于0的坐标系统与基于1的坐标系统之间的转换

人类参考基因组不同版本之间基因坐标不同,进行研究时我们需要统一基因组坐标系,从一个版本的坐标系向另一个坐标系转换的过程就称之为liftover,UCSC,ensembl等为我们提供了便捷的工具,本文以 UCSC 的liftOver工具为例:

首先UCSC liftOver工具提供了即开即用的网页版:http://genome.ucsc.edu/cgi-bin/hgLiftOver

选择目标物种,新旧基因组版本,粘贴或上传原文件就可以开始liftover。


但更多时候我们需要使用命令行的工具,以下,本文主要介绍liftOver的命令行版本,下载地址:http://hgdownload.soe.ucsc.edu/downloads.html

liftover的url: http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver

下载完成后添加到环境路径中,并通过chmod +x 添加执行权限,然后在命令行中输入liftOver验证是否安装成功:(如果成功会有如下界面)

liftOver主要语法如下:

liftOver <输入文件> <chain文件> <输出文件> <unmapped文件>

input/output 可以使用bed格式文件chain file 则需要我们根据转化前后的数据库在ucsc上下载,无法转换的条目则输出到unmapped文件里。

这里以hg19->hg38为例,进行liftover,首先下载hg19tohg38的chain文件,下载地址为: http://hgdownload.soe.ucsc.edu/downloads.html#source_downloads

点击hg19的liftover files后下载对应文件:

下载完成后不需解压即可开始使用:

liftOver input.bed hg19ToHg38.over.chain output.bed unmapped.txt

input.bed

就转换成了: output.bed


在进行LiftOver时,有一点需要我们注意,那就是文件中变异位点起始的编号,是基于0还是基于1的。

因为 UCSC 使用基于0的坐标系统,而 Ensembl 等使用基于1的坐标系统 ,不同工具切换时应该注意这一不同。

除此以外,一些文件格式是基于1的(GFF, SAM , VCF),而另一些是基于0的(BED,BAM),不同文件间转换时也需要注意。

基于0和基于1的坐标系统示可以理解为:

基于1的坐标系:对核苷酸直接编号。

基于0的坐标系:对两个核苷酸之间的间隙编号。

enter image description here

在表示单核苷酸或多个核苷酸变异时,

基于1的坐标系:直接使用变异位点的编号。

基于0的坐标系:变异两边的位置作为起止。

enter image description here

表示insert或deletion时,

基于1的坐标系:Deletion直接使用相应位点编号,insertion则是插入位置两边的核苷酸编号。

基于0的坐标系: Deletion 是插入位置两边的间隙编号表示, Insertions 则直接由插入间隙的编号表示。

基于1的坐标系 与 基于0的坐标系 互相转换时,伪代码如下:

从 基于0的坐标系 向 基于1的坐标系 转换:

if (type=SNV){start=start+1; end=end;}
if (type=DEL){start=start+1; end=end;}
if (type=INS){start=start; end=end+1;}

从 基于1的坐标系 向 基于0的坐标系 转换:

if (type=SNV){start=start-1; end=end;}
if (type=DEL){start=start-1; end=end;}
if (type=INS){start=start; end=end-1;}

参考:

https://www.biostars.org/p/84686/

https://genome.ucsc.edu/cgi-bin/hgLiftOver

http://hgdownload.soe.ucsc.edu/downloads.html#liftover

群体分化系数 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/

gnomAD 数据库: The Genome Aggregation Database

一,数据库简介

gnomAD是目前收录范围最广的基因组变异数据库之一,包含了全世界各人种的变异数据。gnomAD 与有较长历史的dbSNP的主要不同点在于,dbSNP包括了通过各种各样研究方法不同的项目而发现的基因组变异,dbSNP对这些变异加以整理,给予ID,但 gnomAD 为了能够正确的算出等位的频率,对所纳入样本的二代测序数据进行了统一标准的解析,这是 gnomAD 的一大特点。另外,对于50bp以上的基因组结构变异, gnomAD 也有着较高质量的数据。

目前版本 gnomAD v3.1 所包含数据汇总如下:

PopulationDescriptionGenomes
afrAfrican/African-American20,744
amiAmish456
amrLatino/Admixed American7,647
asjAshkenazi Jewish1,736
easEast Asian2,604
finFinnish5,316
nfeNon-Finnish European34,029
midMiddle Eastern158
sasSouth Asian2,419
othOther (population not assigned)1,047

gnomAD 的前身是ExAC (Exome Aggregation Consortium), ExAC 只包含外显子组数据,目前已经被 gnomAD 取代。 gnomAD 的主要资助者是 the Broad Institute。

gnomAD的主页:

二,搜索基因

以ALDH1A1为例,页面最上端显示了基因的基本信息,右边Dataset处可以选择不同的subset(例如 non-cancer, non-neuro等等),Constraint处显示了synonymous,missense 以及 pLoF(基因丧失功能的可能性)的统计值,接下来的图标显示了每个部分的测序深度。

点击 show transcript 或者 show tissue也可以看到不同的transcript 或是在不同组织中的表达。

接下来是ClinVAr中收录的变异位点,与gnomAD中收录的位点的位置信息。

最后则是对变异位点的功能注释:

三,搜索变异

以ALDH2上的rs671为例,

我们可以找到该变异在数据库中的基本信息,包括频数,频率等,之后紧跟着对该变异的注释:

之后是该变异在ClinVar数据库中的信息:

然后是该变异在各群体中的详细的频率信息,以及年龄分布信息。

gnomAD还提供了该变异的质量信息。

最后还提供了便捷的浏览器,可以直观地在基因组中浏览该变异。

gnomAD内所有的信息均提供免费的下载服务,如果有需要也可以按需下载。

参考:

https://gnomad.broadinstitute.org/about