距离上一篇文章的解读已经过了一年，我又有机会拿gnomAD的文章讲组会文献分享了（狗头）（其实这篇文章六一就写好了，但是当时要发的时候简书在系统维护。。。过两天就忘记提交了。。。突然今天又想起来了。。。）

本文不是对文章的直接翻译，而是个人阅读文章后的整理和解读，细节之处如有遗漏可能是我觉得不影响整体理解，但如有理解偏颇之处欢迎指出问题。

《Landscape of multi-nucleotide variants in 125,748 human exomes and 15,708 genomes》 Nature Communications 

标题直译：在125,748个人类外显子组和15,708个人类基因组数中的多核苷酸变异景观/图谱

文章链接：https://www.nature.com/articles/s41467-019-12438-5

一、背景及概念整理

1、什么是multi-nucleotide variants (MNVs)？

    单词直译可以理解为多核苷酸变异，文中的定义可译为：在一个个体中存在于同一单倍型上的两个或多个相近变异的群集。这里主要是和传统意义上的单核苷酸变异single nucleotide variants (SNVs) 进行了区别。在Fig1.a中我们可以看到，左图的两个单点突变分别位于两条染色体上，形如复合杂合子，这不是MNV；右图展示了MNV最重要的点，指两个点突变同时发生在同一单倍型/染色体上。另外，这里的图就是个示例，两个点突变之间的距离可以是≥2的，但是也不要太远，这个概念同时也是区别于GWAS关注的LD上的SNP的，本文主要讨论了距离1~2bp的MNVs。

【啊。。。因为文章具体分析了1-100bp的MNV的calling效率，所以MNV的距离限制和定义在组会上被讨论好久。。。老板觉得可能是作者被reviewer问了calling效率的评估，所以补了1-100bp的所有结果，但其实本质上研究分析的只有1-2bp的MNV。。。如果还有问题可以评论留言讨论或者私信讨论哦！】

Figure 1.a. Definition and an example of an MNV.

2、MNV造成的影响？

    这里主要举例了两点距离≤2bp的MNVs，因为由密码子的变化导致氨基酸的变化是最好理解的，也是最直接影响功能的。在Fig1.b中可以看到，如果一个MNV上的两个突变位置正好位于同一个密码子上，他们分别以SNV形式变异时会导致一个错义突变和同义突变，但是同时变异时，会形成一个终止密码子，导致可能影响蛋白质功能的无义突变。（类似这样排列组合就会有很多的氨基酸变异结果，在Result2中会给这些变化进行具体的分类和解释。）

Figure 1.b. Impact of MNVs in coding regions. 

    当然，除了上述改变之外，还会有很多影响功能的方式和可能性，但是这这篇文章中没有做过多的分析和解释。如有需要，各位研究者可以具体案例具体分析。

3、对MNV的认知程度和研究基础

    根据introduction中的内容可以进行简单总结（对应的ref研究）：

    MNV的发生率：MNV事件发生的频率远高于随机概率（在trio WGS研究中得到），并且在人群中普遍存在（在1000G和ExAC中总结的）；每一代的germline MNV发生率，约等于SNV发生率的1-3%（trio WGS研究）。

    MNV的发生机制：如聚合酶zeta在DNA复制中产生的错误(酵母实验)；受序列上下文的具体碱基影响（在1000G中分析得到）；取决于其位于的基因组功能域和两个突变位点之间的距离（在trio数据中的发现）。

    MNV的影响：比SNV更有可能造成危害（在疾病trio中的发现）。

4、关于MNV的关键科学问题：

① MNVs在全基因组范围下的分布和完整的频率如何？

    除了之前基于ExAC的分析外，没有对成千上万个深度测序个体的整个基因组（包括非编码区）的mnv进行分析。

② 如何在遗传疾病分析中识别和解释MNVs？

    几乎所有现有的临床变异注释工具都将遗漏MNV，这可能在遗传疾病的家族分析中导致漏诊和假阳性遗传诊断。

二、全文框架整理

    基于上述背景和关键科学问题，结合文章不同块面使用的不同数据集，我自己整理了如下框架便于理解：

    1、如何更好的辨识和筛选出MNVs？（Result1 + Methods）

    2、MNVs在外显子组范围内，对蛋白质功能/氨基酸变化的影响情况分布（Result2）

    3、MNVs在全基因组范围内，突变的类型、来源、频率等分布规律，及其和基因组区域的关系（Result3-5）

    4、MNVs在遗传病（罕见病）中的作用 （Discussion）

三、文章具体展开

1、如何辨识和筛选出MNVs？

    这其中涉及到3个步骤：① phasing，② calling，③ filtering

① Phasing 确定单倍型，这里可以用到3种方法：

    1) 只基于测序结果的read-based phasing，普遍使用于所有测序结果，无关乎样本量；

    2) 基于家系测序结果的family-based phasing，加入父母的信息能更准确的判断变异来源从而确定单倍型情况；

    3) 基于群体单倍型规律的population-based phasing，常见于GWAS研究的前序步骤（工具如SHAPEIT），常用于LD区域相关的phasing处理。

    本文使用read-based的GATK HaplotypeCaller作为phasing算法工具；为了评估其表现，以trio-based的算法对应结果作为金标准。

② Calling 识别MNV：

    使用Hail这个工具(https://github.com/hail-is/hail)，利用window_by_locus这个参数调整来识别不同距离的MNVs；此外，研究只分析了常染色体。

    结合上述的phasing和calling策略，文章先评估了phasing的效果和calling时适合的距离：

Result1: 当MNVs≤10bp时，read-based phasing (GATK HaplotypeCaller)表现良好。

    GATK HaplotypeCaller可以成功对>85%的临近杂合子对进行单倍型分类（trio-based只能成功区分>60%）；在比较两种方法的结果时，>99.8%的trio-based结果可以在HaplotypeCaller中被识别出来（Fig1.c）。在下图的右侧和sup-table1-2中可以看到在calling distance ≤10bp时，敏感性和精确性都可以保持较高，>10bp后出现断崖式下降（Sup-Fig1.d）。

Figure 1.c  Graphical overview of the analysis of phasing sensitivity and specificity using trio samples from our gnomAD callset.

    综上分析，本文后续的分析中即以如下标准进行MNVs的辨识和筛选：（其中的③Filtering以quality control为主要目的）

2、MNVs在外显子组范围内，对蛋白质功能/氨基酸变化的影响情况分布

    本段分析只用到了125,748个外显子组数据，和≤2bp的MNVs，旨在分析单核苷酸内的MNVs的功能变化情况。

    这里需要引入多个作者定义的概念 （supplementary information）：

Gained nonsense: Neither of the individual SNPs is a nonsense/stop-gained mutation, but the MNP is. 

获得性无义突变：原来两个SNP都不是无义突变，MNP组合起来变成了无义突变。

Rescued nonsense: at least one of the individual SNPs is a nonsense/stop-gained mutation, but the MNP is not.

被拯救的无义突变：原来两个SNP至少有一个导致无义突变，MNP组合起来就不是无义突变了。

Gained Missense: the individual SNPs are synonymous, but the MNP results in a missense variant.

获得性错突变：原来两个SNP都是同义突变，MNP组合起来变成了一个新的错义突变。

Lost Missense: at least one of the individual SNPs is a missense variant, but the MNP is synonymous.

错过的错义突变：原来两个SNP至少有一个导致错义突变，MNP组合起来就是同义突变。

Changed Missense: at least one of the individual SNPs is a missense variant and the MNP is a new missense variant with a different resulting amino acid.

被改变的错义突变：原来两个SNP至少有一个导致错义突变，MNP组合起来是全新的错义突变。

Partially Changed Missense: The MNP is composed of two different missense variants that when considered together have the same amino acid outcome as only one of the variants (e.g. Missense A + Missense B = Missense A).

部分被改变的错义突变：原来两个SNP分别导致不同的错义突变，MNP组合起来是两者其一的错义突变。

Unchanged: Either the outcome of the MNP is identical to that of the individual SNPs or one of the SNPs is a synonymous variant that does not change the outcome of an adjacent non-synonymous variant

不改变的突变：要么MNP的结果与单个SNP的结果相同；要么其中一个SNP是同义变异，它不会改变相邻非同义变异的结果。

Result2：Rescued nonsense 比 Gained nonsense在（高约束性）基因中富集得更多；LoF相关的MNV注释非常重要。

    经统计，有31,575个MNVs位于同一个密码子中，他们的突变类型分布如Fig2.a。

    从功能上来说，gained / rescued nonsenses更可能影响基因的功能，所以关注到这两类突变上来看：1633个基因携带gained / rescued nonsenses，其中包括41个与疾病相关的基因（Fig2.b)。

    有≈6%的gained nonsenses、≈20%的missense\unchanged、≈22%的rescued nonsenses，他们都是在LOEUF decile <20%的3941个基因中（约束性较高、较易和疾病相关的基因）；由此可见，没有MNV注释时，LoF (loss-of-function)注释错误显著增多（Fig2.c）。

Figure 2. Functional impact of MNVs.

3、MNVs在全基因组范围内，突变的类型、来源、频率等分布规律，及其和基因组区域的关系

    由于在所有≤100bp的MNVs中，有87.9%的MNV是bp=1的adjacent MNV，所以接下来的这部分的分析都是基于15,708个全基因组数据中call出来的adjacent MNVs进行的分析。

Result3-4：adjacent MNVs的全基因组分布符合三个主要突变起源的假设，每种假设分别解释了相当一部分MNVs。

    首先，对于adjacent MNVs的产生机制，之前的研究认为主要由3种原因导致：1) 由独立的单点突变组合而成（颠换or转换、CpG区域均会影响单点突变的发生率）；2) pol-zeta导致的特定类型的DNA复制错误（TC>AA,GC>AA & GA>TT,GC>TT；3) 重复序列处的聚合酶滑脱事件。(Fig3.a)

Figure 3.a. Three major categories of the mutational origin of MNVs. 

    最常见的MNV模式CA -> TG（C->T的CpG转换Ti和A -> G的转换Ti组合）；最不常见的MNV模式是TA -> GC（两个非CpG颠换Tv的组合）。（Fig3.b，Sup-Fig4 c.e）下图最右的图可以发现有些突变来源的类型和第一个来源假设（SNV组合）并不一致，接下来换个角度关注后两种假设。

        以两个位点分别的突变频率是否一致来评估他们是否是同时突变的one-step MNV，统计分析发现，90.5% 的 GA- > TT 和 80.5% 的 GC- > AA （poly-zeta error的典型突变型，紫色）都是one-step MNV，显著高于所有突变型的平均值39.9%。（重复区域的TA>AT，one-step MNV比例也很高）

    关注到重复序列区域，突变类型及来源分布可以看到聚合酶滑脱事件的主要类型（AA>TT,AT>TA,TA>AT）在重复序列区域上显著富集。

    不同的MNV模式所占比例相差很大；每种MNV模式对应的起源占比各不相同，但总的来说，83.2%的MNVs可以找到对应的起源机制之一。

Figure 4.a. Distribution of MNVs across genome.

    Result5：一个基因组区域的MNV密度高度依赖于周围序列的CpG甲基化状态

    接下来研究了MNV不同模式的分布在功能注释类别之间的差异，使用了编码序列、增强子、启动子等13种不同的功能注释，以及来自ENCODE的DNA甲基化注释。

    MNV密度：在每个区域中，WX>YZ突变序列数量占原始WX序列的比例。

    涉及CpG-Ti相关的突变模式中，MNV密度与甲基化水平呈正相关；相反地，Non-cpg Tv 相关的突变模式以及与pol-zeta滑移相关的突变模式，MNV密度与甲基化状态呈负相关。（Fig.4.b-c）

    在基因编码序列周围选择了7个主要的区域注释，计算了这些区域中可能由不同突变来源解释的mnv的比例。在所有区域中，发现MNV信号主要由CpG-Ti所主导；Pol-zeta error的主导性次之，除了在转录起始区域（TSS，在这七个注释中甲基化率最低的区域）；Non-cpg Tv和重复区域的聚合酶滑脱的比例始终低于(或几乎等于)总信号的5%。(Fig.4.d)

Figure 4.b-d. Distribution of MNVs across genomic regions.   

4、MNVs在遗传病（罕见病）中的作用 

    Discussion：MNVs对罕见疾病的诊断贡献很少（即时用了文章的pipeline，作用也很有限）

    虽然文章的背景中提到了MNV在疾病中的意义，比如在Genome Res. 2019中发现：在受影响的儿童中，de novo MNVs 在先前与发育障碍相关的基因中显著富集。

    但在本文的结果部分其实只有一小部分提到了：在6072个罕见病家庭外显子组中（包括4275例病例样本）发现了16个gained nonsense和110个changed missense（CADD > 20分高，gnomAD频率低(n < 10)）；在人工评判后，这些mnv都不是罕见疾病的原因。但是本研究中管道的应用可能会发现以前的漏诊。

四、take-home message

    对我而言，这篇文章的意义就是：

    1、了解了MNV的phasing、calling、filtering流程

    2、MNVs在遗传病（罕见病）的诊断中作用有限，如果要分析，优先级排后面。。。