GWAS 全基因组关联分析 – 第 5 页 – GWASLab

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.

GWAS QC – 性别检查与阈值选择 Sex check

性别错误 (Sex discrepancy)：有时数据录入会有错误，我们应基于X染色体的杂合性，确定是否有性别录入错误的个体。PLINK默认对于男性，X染色体的纯合性（homozygosity）估计值应高于0.8，女性应低于0.2，但实际情况更为复杂，本文简要对性别检查与阈值选择进行讲解。

本文使用数据来源：https://onlinelibrary.wiley.com/doi/full/10.1002/mpr.1608

PLINK中提供了check-sex选项，通过计算X染色体的纯合性来判断样本性别是否有错误：—check-sex [女性阈值] [男性阈值] 默认的参数为 0.2 0.8，但该标准过于严格，并不适用与实际情况，下面会讲解。首先介绍基本操作：

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

#估计X染色体的纯合性
plink \
        --bfile HapMap_3_r3_1 \
        --extract HapMap_3_r3_1.prune.in \
	--check-sex \
        --out HapMap_3_r3_1

注意：

要提前确认X染色体上PAR区域已经被分离（可以使用—split-x 选项来进行分离）
默认的阈值是，男性>0.8，女性<0.2

log文件如下：

PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to HapMap_3_r3_1.log.
Options in effect:
  --bfile HapMap_3_r3_1
  --check-sex
  --extract HapMap_3_r3_1.prune.in
  --out HapMap_3_r3_1

191875 MB RAM detected; reserving 95937 MB for main workspace.
Allocated 3037 MB successfully, after larger attempt(s) failed.
1457897 variants loaded from .bim file.
165 people (80 males, 85 females) loaded from .fam.
112 phenotype values loaded from .fam.
--extract: 166240 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 53 nonfounders present.
Calculating allele frequencies... done.
Warning: 35 het. haploid genotypes present (see HapMap_3_r3_1.hh ); many
commands treat these as missing.
Total genotyping rate is 0.996691.
166240 variants and 165 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
are missing.)
--check-sex: 3993 Xchr and 0 Ychr variant(s) scanned, 1 problem detected.
Report written to HapMap_3_r3_1.sexcheck .

–check-sex: 3993 Xchr and 0 Ychr variant(s) scanned, 1 problem detected.从日志文件中可以看出我们检测出了1个不一致的样本。上述代码运行结果如下：

head HapMap_3_r3_1.sexcheck

FID       IID       PEDSEX       SNPSEX       STATUS            F
   1328   NA06989            2            2           OK      0.02151
   1377   NA11891            1            1           OK            1
   1349   NA11843            1            1           OK            1
   1330   NA12341            2            2           OK     0.009813
   1444   NA12739            1            1           OK            1
   1344   NA10850            2            2           OK       0.0533
   1328   NA06984            1            1           OK       0.9989
   1463   NA12877            1            1           OK            1
   1418   NA12275            2            2           OK      -0.0747

前两列为FID与IID，三四列为PED文件中的性别，与通过X染色体估计的性别，第五列表示是否一致，不一致时会表示为PROBLEM，第六列是估计的F系数。可以提取出不一致的样本：

cat HapMap_3_r3_1.sexcheck | grep PROBLEM
1349   NA10854            2            1      PROBLEM       0.9933

#将有问题的样本FID与IID输出到sex_error_to_remove.id文件中
cat HapMap_3_r3_1.sexcheck | grep PROBLEM | awk '{print $1,$2}'> sex_error_to_remove.id

后续PLINK分析可以使用 —remove sex_error_to_remove.id 来剔除性别错误的样本.以上结果的F系数分布如图所示：

实际情况：上面的例子是一种很理想的分布，可以看到样本被清晰地分开，但实际情况要复杂很多，例如在千人基因组项目Phase1的样本中，一些女性样本的F系数也大于0.8，即使进行LD-pruning后，女性样本最大的F系数也超过了0.66。所以一般来说，我们会进行两次—check-sex，第一次用于观察分布，在看到一个右侧紧凑的峰（男性），与左侧平缓的峰后（女性），根据其分界线来决定判定阈值，（真实的数据的分布大多如下图所示）

这里常用的阈值为 -check-sex 0.7 0.8

#如果使用真实数据，通过目视决定阈值后，更改参数，重新判定性别
plink \
        --bfile HapMap_3_r3_1 \
        --extract HapMap_3_r3_1.prune.in \
	   --check-sex 0.7 0.8 \
        --out HapMap_3_r3_1

参考：

https://onlinelibrary.wiley.com/doi/full/10.1002/mpr.1608

Basic statistics

孟德尔随机化系列之一：基础概念 Mendelian randomization I

本文是MR系列的第一篇，孟德尔随机化的简介，该系列会介绍孟德尔随机化的基础概念，统计方法分类，常见误区与实践操作等内容。

1.背景与目的
- 1.1 明确因果关系
- 1.2 RCT是金标准，但缺点明显
- 1.3 孟德尔随机化的本质
2.孟德尔随机化的统计学方法 – 工具变量
3.核心假设
- 3.1 关联性假设
- 3.2 排他性限制
- 3.3 独立性假设
4.孟德尔随机化的优势

1 背景与目的

1.1目的是明确因果关系

在关联分析中我们时常面对的一个问题便是我们很难确定一个变量是否是真正的因果变量，而非有其他未观测的因素同时影响这个变量与结果，造成这个变量与结果相关联。在循证医学中，或是制定干预策略时，明确因果性是十分必要的。

这个问题实际上与内生性 endogeneity 相关，包括：反向因果关系 reverse causation, 忽略的混淆变量造成的偏倚 omitted variable bias due to confounding, 测量误差measurement error, 以及双向因果关系bidirectional causality等等问题。（这里的内生性在统计学上是指在回归分析中，解释变量（x）与误差项相关。）

1.2. RCT是金标准，但缺点明显

一般来说，明确因果关系的金标准时随机对照试验 RCT randomized control trial (RCT)，即对受试者随机分为对照组和实验组，以研究某个因素的影响。但现实中，要完成随机对照试验的难度非常高，需要大量的人力物力，有时因为伦理问题，对某个因素的研究几乎是不可能的。这时我们就要借助其他方法，而孟德尔随机化就是其中之一。

1.3. 孟德尔随机化与RCT的相似性

孟德尔随机化（MR，Mendelian randomization）便是为了解决以上问题而开发的方法，MR与RCT直接相关，两者有很高的相似性，如下图所示。

孟德尔随机化的核心其实是利用了孟德尔第二定律，也就是自由组合规律(law of independent assortment），当具有两对（或更多对）相对性状的亲本进行杂交，在子一代产生配子时，在等位基因分离的同时，非同源染色体上的基因表现为自由组合，这一过程类似于随机对照试验中的随机分组，所以我个人理解的孟德尔随机化就是基于孟德尔第二定律的随机对照试验。

2 孟德尔随机化的统计学方法 – 工具变量

孟德尔随机化在统计学上的本质实际是利用工具变量（Instrumental variables）来研究因果性，这一方法常用在经济学研究中。

工具变量简单来说就是，一个与X相关，但与被忽略的混淆因素以及Y不相关的变量。在经济学研究中工具变量可以是政策改革，自然灾害等等，而在遗传学中，这个变量就是基因。

如果一个基因变异Z 是某个暴露因素X的因果变量，并且对结果Y没有直接因果关系，那么这个基因变异Z与结果Y的关联，只能通过X对Y的因果关系而被观察到（X->Y）。

2.1 两阶段最小二乘法

通常我们可以用两阶段最小二乘法（2SLS，2 stage least squared method）来估计X对Y的效应：

考虑一种最简单的单样本的情况，有一个基因变异Z，与Z相关的因素X，以及与Z不相关的结果Y，

我们想探究X与Y之间的因果关系。

第一阶段，X对工具变量进行回归，

第二阶段，Y对第一阶段X的预测值进行回归，

合并后可以化为Y直接对工具变量进行回归。

我们所关心的系数β2SLS实际上也等同于两段协方差的比值

2.2 两样本MR

另一种常见的情况则是两样本MR，如果我们有一个与X相关联的工具变量，我们只有在X对Y有因果关系的情况下，才能观测到这个工具变量与Y的关联。

这意味着βiv,y = βiv,x 乘以 βx,y。也就是说，我们可以不用通过X与Y的回归来估计β，

而是可以简单地通过 βx,y = βiv,y / βiv,x 来计算 X对Y的效应量。

这就意味着与两阶段最小二乘法相对，我们可以利用两个独立的GWAS 的概括性统计量来计算这个比值。这种方法通常叫做两样本MR.

当然MR还有其他更复杂的统计模型方法，这里不做深究，有兴趣的朋友的可以看这篇文章：预留链接

核心假设：

当然，既然是统计模型那就要满足模型的基本假设，通常情况下MR建立在几点基本假设之上，

主要假设：

3.1 遗传变异必须与暴露因素X强相关。（关联性假设），例如：弱工具变量的使用会导致估计出现偏倚。

3.2 遗传变异不能与结果直接相关。(排他性限制)，例如：可能影响因素包括多效性等。

3.3 遗传变异不能与任何可能的混淆因素相关 (独立性假设)，例如：人群分层

其他假设：

3.4 不存在选型交配 No genetic assortative mating，例如：人们经常会与自己教育和经济水平相似的人结婚。

3.5 对所有个体，IV对于X的影响方向是相同的。例如：潜在的上位效应与GxE基因与环境的相互作用都可能会影响此假设。

总结来看，孟德尔随机化以基因型作为工具变量的优势是：

4.1 遗传相关中，因果关系的方向是确定的，遗传多样性导致了不同的表型，反之则不成立

4.2 一般情况下我们所测量的环境暴露因素都或多或少与行为，社会，心理等因素相关，造成偏倚。但遗传变异则不受这些混淆因素影响。

4.3 相对来说，遗传变异与其效应的测量误差较小。

4.4 并不一定要找到因果SNP，一个与因果SNP处于LD的SNP即可满足假设条件。

4.5.目前GWAS的数据相对容易获取。

参考：

Melinda C. Mills, Nicola Barban, and F. C. T. An Introduction to Statistical Genetic Data Analysis. (2020).

Curr Epidemiol Rep . 2017;4(4):330-345. doi: 10.1007/s40471-017-0128-6. Epub 2017 Nov 22.

为什么在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版本

使用LocusZoom绘制Regional plot详解 – GWAS可视化

Locuszoom是一款非常好用的绘制region plot的软件，是GWAS研究中的必备软件之一。

Regional plot示例： https://my.locuszoom.org/gwas/881707/region/?chrom=3&start=45634967&end=46034967

Regional plot：横轴为染色体位置，左边纵轴为-log10（p），marker的颜色表示与leading SNP 的LD的大小，右侧纵轴则表示重组率。下方还标有gwas catalog中已知的与疾病相关的hit，以及该区域中基因的范围。该图主要用来用于精确定位因果SNP。

locuszoom官网： http://locuszoom.org/

该软件提供有最初的网页版（底层还是命令行版本），近来新发布的网页交互版（基于JavaScript），以及单独可下载的linux命令行版本（基于R/Python）。

网页版方便直接，命令行自定义程度更高。

本文以命令行版本为基础简要介绍locuszoom的使用方法，同时与网页版一一对照来方便理解与使用。

本文目录：

1. 命令行版本简介与软件与参考数据准备（**只看网页版的话可以跳过 1**）
2. 输入文件格式 
3. 指定绘制的区间范围
4. 指定LD文件，族裔群体，参考基因组版本
	**-》以上内容对应的网页版本**
5. 批量模式 ！ Batch mode（命令行版本的优势所在）
	**-》批量模式对应的网页版本**
6. 自制LD文件
7 额外功能
	7.1指定多个参考SNP
	7.2 使用GWAS catalog进行注释
	7.3 精确定位的可信集
	7.4 对多个SNP的注释
	7.5 绘制额外的BED文件轨道
参考url

使用方法：

命令行版本简介与软件与参考数据准备（只看网页版的话可以跳过 1）

命令行版本总体上与网页版一一对应，但可自定义设置的程度更高，目前仅支持linux版本，windows和mac用户只能使用网页版或通过虚拟机。

下载地址：https://github.com/statgen/locuszoom-standalone

注意我们需要下载：程序 + 数据库 + LD信息

（*单独下载软件无法绘制图像，或是需要使用自制的LD信息）

目前的最新版本为1.4版，发布于2017-05-01，共有23G。

除此之外，我们还需要准备以下有版本要求的python和R：

1.Python 2.7+ (注意要使用PYTHON2！)

2.R 3.0+.

locuzoom下载完成后，解压并添加进环境后即可使用：

cd <directory where you want to place locuszoom> 
tar zxf /path/to/locuszoom.tgz
ln -s bin/locuszoom /usr/local/bin/locuszoom #根据自己的环境变量制定

locuszoom的文件结构主要包括的内容如下所示：

locuszoom/
	bin/
		locuszoom (this is the locuszoom "executable")
		locuszoom.R (the R script which is used by locuszoom for creating the plots)
		dbmeister.py (script for creating custom user databases)
		lzupdate.py (script for creating an updated copy of the provided locuszoom database)
	conf/ (configuration file located here)
	data/
		database/ (SQLite file located here)
		hapmap/ (hapmap genotype files)
		1000G/ (1000G genotype files)
	src/ (source code for locuszoom)

命令行locuszoom语法：

locuszoom [输入文件] [选项]

一个简单的例子：

locuszoom --metal Kathiresan_2009_HDL.txt --refgene FADS1

输入文件格式

locuszoom主要使用METAL 或者 EPACTS的文件格式，这里主要介绍METAL格式，因为其准备起来更加简单便捷。

METAL格式：文件必须包含以tab分隔的下两列 1.markers 与 2.p-values

如下所

示

locuzoom在输入时的选项：
--metal 指定metal格式的输入文件名
--delim （可选）可以指定分隔符，默认为tab
--markercol  指定输入文件表示marker的列名，默认为"MarkerName"
--pvalcol 指定输入文件表示p值的列名，默认为"P-value"
--no-transform 当p已经转换为-log10（p）时，可以使用此选项跳过log转换

对应网页区块：

首先上传输入文件 , marker与P value的列名，以及指定分隔符

指定绘制的区间范围

在命令行版本中我们可以通过多种方式指定区间，包括

#2.1指定SNP与两侧的范围
--refsnp <your snp> --flank 500kb

#2.2指定SNP与区间的起止位点
--refsnp <your snp> --chr # --start <base position> --end <base position>

#2.3指定SNP与基因（绘制基因所在的范围）
--refsnp <your snp> --refgene <your gene>

#2.4 指定基因与两侧范围
--refgene <your gene> --flank 250kb

#2.5 指定基因与区间的起止位点
--refgene <your gene> --chr # --start <base position> --end <base position>

#2.6 区间的起止位点
--chr # --start <base position> --end <base position>

--flank 选项在这里指从起始和终止位点向外侧计算的距离，单位为kb，在没有指定参考SNP时，locuzoom会自动选择最显著的SNP作为参考SNP。

网页版中指定绘制的区间范围所对应区块，填写方法与命令行版本一致：

4 指定LD文件，族裔群体，参考基因组版本

目前locuszoom内置了以下的组合，可以根据自己需要进行指定：

Genotype files available for: 
--source hapmap
  --build hg18
    --pop YRI
    --pop CEU
    --pop JPT+CHB

--source 1000G_March2012
  --build hg19
    --pop AMR
    --pop ASN
    --pop AFR
    --pop EUR

--source 1000G_June2010
  --build hg18
    --pop YRI
    --pop CEU
    --pop JPT+CHB

--source 1000G_Nov2014
  --build hg19
    --pop AMR
    --pop ASN
    --pop AFR
    --pop EUR
  --build hg38
    --pop SAS
    --pop EAS
    --pop AMR
    --pop AFR
    --pop EUR

一个简单的使用例：

--pop ASN --build hg19 --source 1000G_March2012
#指定1000G_March2012数据，以hg19为参考基因坐标，来计算ASN族裔的LD

以上内容对应的网页版：

在Genome Build/LD Population中选择所使用的参考基因组版本与LD文件

有了以上信息就可以进行绘制。命令行中运行命令或是在网页版中点击Plot Data按钮。根据所绘制SNP数量多少，运行时间会有不同。等待即可。

5 批量模式！ Batch mode（命令行版本的优势所在）

当我们有多个基因座，需要绘制多张regional plot，可以使用批量模式。

只需要通过 --hitspec 选项来指定一个包含多个区域信息的文件即可。

文件需要包含以下列：

snp 2.chr 3.start 4.end 5.flank 6.run 7.m2zargs

该文件格式如下：

选项1-5与之单次模式的内容相同，run 是指是否指定该行进行绘制，m2zargs 则是对该行信息绘制时额外的参数。

对应的网页版区块与内容：

6 自制LD文件（命令行版本）

可以通过--LD选项来指定自己的LD文件，文件格式如下：

dprime列可以为空，但Rsq列必须为有效的数据。该文件以空格分隔，且必须包含header。

7 额外功能（未完待续，持续更新）

7.1指定多个参考SNP

7.2 使用GWAS catalog进行注释

7.3 精确定位的可信集

7.4 对多个SNP的注释

7.5 绘制额外的BED轨道

参考：

http://locuszoom.org/

https://genome.sph.umich.edu/wiki/LocusZoom_Standalone#User-supplied_LD

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：

BETA – beta
SE – std. error

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

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

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

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

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可以进行如下的功能注释：

基于基因的注释 Gene-based annotation：主要针对SNP或CNV是否引起蛋白编码改变进行注释，可以灵活选用 RefSeq genes, UCSC genes, ENSEMBL genes, GENCODE genes, AceView genes等多种不同来源的基因定义系统。
基于区域的注释 Region-based annotation：针对基因组某一特定区域的变异进行注释，例如44个物种的保守区域，预测的转录因子结合位点，GWAS hit， ENCODE H3K4Me1/H3K4Me3/H3K27Ac/CTCF sites，ChIP-Seq peaks, RNA-Seq peaks等
基于筛选的注释 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/

--rg ：计算遗传相关性 （genetic correlation），参数为上一步处理好的数据文件名

--out ： 输出文件的路径与前缀

`--rg` ：计算遗传相关性（genetic correlation），参数为上一步处理好的数据文件名

`--out` ：输出文件的路径与前缀