全基因组测序分析鉴定出与循环蛋白水平相关的罕见、大效应非编码变异和调控区域

Basic Information

英文标题：Whole-genome sequencing analysis identifies rare, large-effect noncoding variants and regulatory regions associated with circulating protein levels
中文标题：全基因组测序分析鉴定出与循环蛋白水平相关的罕见、大效应非编码变异和调控区域
发表日期：24 February 2025
文章类型：Article
所属期刊：Nature Genetics
文章作者：Gareth Hawkes | Robin N. Beaumont
文章链接：

Abstract

Para_01

罕见非编码遗传变异对常见表型的贡献在很大程度上是未知的，这是由于历史上缺乏全人群规模的全基因组测序数据，以及难以将非编码变异分类为功能相似的组别。
为了开始解决这些挑战，我们使用全基因组测序数据进行了顺式关联分析，该数据包括11亿个变异，1.23亿个基于非编码聚合的测试和约5万名英国生物银行参与者的2907种循环蛋白水平。
我们发现了604个与循环蛋白水平相关的独立罕见非编码单变异关联。
与蛋白质编码变异不同，罕见非编码遗传变异几乎同样可能增加或减少蛋白质水平。
罕见非编码聚合测试确定了357个条件独立相关区域。
其中，74个（占21%）仅通过单变异测试无法检测到。
我们的研究结果对识别与常见人类表型相关的罕见非编码遗传变异及其作用具有重要意义，包括测试非编码变异聚合的重要性。

Main

Para_01

人类基因组中非编码区域的稀有遗传变异可导致严重的罕见病，但它们在常见复杂性状中的作用仍然 largely 未知。
基于阵列的推算和全基因组关联研究 (GWAS) 已经发现了数以万计的人类疾病与常见变异的关联，其中大多数位于编码区域之外。
然而，利用外显子组测序数据识别与常见表型相关的稀有变异的努力主要局限于编码区域，这些研究结合了单变异测试（统计功效低但功能解释更清晰）和聚合测试（统计功效更高但难以解释）。
我们评估非编码基因组中稀有变异的能力直到最近还受到限制，这是因为缺乏大规模研究中的全基因组测序 (WGS) 数据，并且难以定义具有生物学意义的非编码调控基因组单元。

Para_02

利用全基因组测序（WGS）数据识别非编码调控元件，可以为基因调控提供重要的见解，这补充了从外显子组测序和基于阵列的研究中获得的知识。
根据群体遗传学指标（如约束性），功能性非编码DNA的数量估计是编码序列的4到5倍，且分别有10%的启动子和6%的增强子受到与编码区域相当的突变约束。
全基因组测序使我们能够研究内含子、近端和远端调控元件的作用，并覆盖人群中的整个等位基因频率谱，包括那些即使在非常大的样本中也仅被观察到一次或两次的变异比例。

Para_03

在常见表型的背景下，关于全基因组测序（WGS）数据的研究非常少。
最近来自 TOPMed 的例子考虑了血脂水平（样本量 N = 66,000）和血压（样本量 N = 51,456），以及英国生物样本库（UKB）中的一项关于身高的研究（发现样本量 N = 200,003）：这些研究发现的新信号数量有限，可能是因为样本量相对较小，不足以检测新的罕见变异。
此外，WGS 还被用于研究多达 3,000 名个体的蛋白质水平相关的遗传变异。

Para_04

英国生物银行（UKB）发布循环蛋白数据，并结合全基因组测序（WGS），为测试罕见非编码遗传变异对常见、生物学上接近且测量良好的人类表型的影响提供了前所未有的机会。
三篇最近发表的研究描述了这些数据在2023年的发布，涉及54,306名个体中的2,923种循环蛋白。
首先，Eldjarn等人通过单变量分析识别出30,062个蛋白数量性状位点（pQTLs），这些数据来自UKB的150,119例全基因组测序分析，包含2,931种测量的蛋白水平，并将结果与冰岛全基因组测序队列的蛋白质组学数据进行了比较。
其次，Dhindsa等人通过外显子组测序分析识别出5,433个pQTLs，并在编码区域内进行聚合测试，鉴定出1,962个基因-蛋白关联关系。
最后，Sun等人通过基于传统插补的全基因组关联研究（GWAS）识别出14,287个pQTL单变量。

Para_05

利用全基因组测序（WGS）数据和循环蛋白水平作为示例特征，我们测试了两个相关假设：（1）罕见的非编码单基因变异体，目前通过推算或外显子组测序无法检测到，对常见人类表型有与编码变异相似的影响；（2）基因组调控区域中的罕见非编码遗传变异集合（群体）与人类表型相关，类似于编码序列中的基因水平负担分析。
重要的是，与之前关于英国生物银行（UKB）蛋白质组数据的三篇论文不同，我们使用了所有具有循环蛋白测量值的参与者的短读长全基因组测序检测到的DNA序列变异的完整等位基因频率范围，提供了约11亿个变异的信息，但将每种表型的搜索限制在衍生该蛋白的蛋白编码基因周围的顺式调控区域。

Para_06

我们使用来自英国生物样本库（UKB）的46,362名推断为欧洲遗传背景个体的注释全基因组测序（WGS）数据，对2,903种测量的循环蛋白水平进行了初步发现关联分析。
此外，我们分别在899名和1,098名推断为南亚和非洲遗传背景的个体中进行了单变异分析。
我们的分析鉴定了1,651个高质量单变异，对比了编码和非编码变异的影响，通过识别出357个独立位点突显了罕见变异非编码聚合检测的能力，并展示了在基于全基因组测序的关联研究中考虑测序覆盖度的重要性。

Results

Para_01

排除了十二种蛋白质，这些蛋白质或者是融合蛋白，或者无法直接与 HGNC 基因符号匹配（补充表 1）。
对于每种测量的蛋白质，我们在编码该蛋白质的基因周围的顺式窗口中进行了单变异体（小等位基因计数 (MAC) ≥ 5）和基因组聚合关联测试（小等位基因频率 (MAF) < 0.1%），顺式窗口从 5′- 和 3′-UTR 扩展 1 Mb，基于该基因任何 GENCODE 转录本的最极端 5′- 和 3′-末端。
我们使用 1 Mb 作为最近确定的边界距离，当一个 pQTL 更可能是一个顺式而不是反式关联时即适用此距离。
在运行时对循环蛋白质水平测量值进行秩逆归一化，并将年龄、年龄平方、性别、招募中心、全基因组测序 (WGS) 中心、采血后时间、禁食时间、40 个遗传主成分 (PCs) 和 Olink 批次 ID 包括为协变量（方法）。
尽管仅考虑每个蛋白质的 2-Mbp 顺式窗口，我们仍测试了 128,434,590 个单变异体（包括单核苷酸变异体和小插入或缺失 (indels)）和结构变异体关联，其 MAC ≥ 5。
我们使用 GCTA-CoJo 的修改版本识别独立变异体关联（方法）。

Para_02

我们使用 Ensembl 的 Variant Effect Predictor v.110.1 (VEP) 对所有遗传变异进行了注释（方法），并利用其输出将变异分类为基因中心型（例如，编码、预测的内含子剪接、未指定内含子和近端调控）和基因间调控型（例如，Ensembl 调控区域、非编码 RNA 和未指定基因间）以进行基于聚合的关联测试。
如果非编码变异位于基因的 UTR 区域或内含子区域内，或者距离 UTR 不超过 5 kbp，则将其注释为基因中心型。
此外，我们在重叠（1-kb 重叠）的 2-kb 滑动窗口中对所有非编码变异进行了聚合测试。
我们还通过约束度量（JARVIS）、保守性（基因组进化速率分析 GERP）和/或预测的有害性（组合注释依赖缺失评分 CADD）进一步对部分聚合单元内的变异进行了细分。
为了识别独立的、罕见的、非编码聚合关联，我们对非编码聚合测试调整了常见 pQTLs（MAF > 0.01）以及与蛋白自身编码基因相关的所有编码变异，无论变异频率如何，该基因此后称为同源基因（图 1）。
总共，我们对 123,598,575 个聚合进行了关联测试，其中包括 182,136,116 个 MAF < 0.1% 的变异，直到包括单例为止。
统计显著性基于模拟研究定义（方法）：单个变异 P < 2.95 × 10−10，基因组聚合 P < 8.71 × 10−9。

Fig. 1: Study design for pQTL single and aggregate genetic variant analysis.

- 图片说明

◉ 展示研究设计的流程图。

We identified 13,457 candidate pQTLs across 2,891 proteins

我们在2891种蛋白质中鉴定出13457个候选pQTLs

Para_01

对于2,891种循环蛋白质，我们确定了13,457个统计上独立的顺式-pQTL关联（MAC ≥ 5），其中293个为结构变异（图2和图3以及补充表2）。
我们发现了2,036种蛋白质中至少存在一个顺式关联，每种循环蛋白质中位数为5个独立pQTL；429种蛋白质与超过10个独立pQTL相关联，其中CD177的最大值为49个。
这些结果与之前使用150,000个英国生物样本库基因组的研究一致，并利用相同的数据对剩余样本进行了推算。

Fig. 2: Manhattan results of single-variant and aggregate pQTL analyses.

- 图片说明

◉ 显示顺式变异与循环蛋白水平相关性的曼哈顿图，调整了相关的常见变异和同源基因的所有编码变异后。◉ a-e，x轴表示基因组位置，y轴显示我们所有蛋白质的顺式结果的-log10P双侧值，分为单个变异（a），编码聚合体（b），以基因为中心的调控（近端）聚合体（c），基因间和内含子调控聚合体（d）以及滑动窗口聚合体（e）。◉ 红线代表Bonferroni显著性阈值（单个变异为P ≤ 2.95 × 10−10，聚合测试为P ≤ 8.71 × 10−9）。◉ P值来自混合线性模型的双侧检验。

Fig. 3: Effect size distributions of rare, independently associated pQTL variants.

- 图片说明

◉ 图a至d展示了罕见pQTL变异的影响大小（β）与小等位基因频率（MAF）的关系，突出显示了编码基因（a）和非编码基因（b）的pQTLs中对应的（红圈）和非对应的（蓝方块）效应，并按预测后果对编码变体进行分层，包括对应基因中的编码变体（c）以及非编码对应基因中的变体（这些变体被注释为对蛋白质编码基因具有调控作用）（d）。◉ 非编码、非对应变体按注释分层的效果大小展示在扩展数据图3中。

Over half of independent pQTLs may be coverage artefacts

超过一半的独立pQTL可能是覆盖范围伪影

Para_01

由于我们的研究旨在评估非编码变异对蛋白质水平的作用，我们将所有全基因组显著的变异（包括编码和非编码）共同建模（方法），以避免单倍型效应带来的混淆。
然而，我们随后发现，具有最大数量统计独立 pQTL 的蛋白质在对应基因的低覆盖区域中显著富集（补充表 2 和 3）。
例如，CD177（MANE Select 转录本 ENST00000618265.5）有 49 个独立的 pQTL，但在外显子 5 或内含子 6 上几乎没有覆盖（扩展数据图 1）。
此外，那些其对应基因表达主要由重复多态性的拷贝数驱动的蛋白质也显示出富集现象。
例如，FCGR3B（43 个独立信号）、AMY2A（25 个独立信号）和 AMY2B（27 个独立信号）是由已知具有高度可变拷贝数的基因编码的，其中淀粉酶水平主要由后两个基因的拷贝数决定。
此外，脂蛋白(a) 水平（34 个独立信号）主要由 Kringle 的重复多态性解释。
对于这些基因中的每一个，独立信号数量较大的最可能解释是真正的因果变异尚未被充分校正，从而导致了过多的人工独立 pQTL 数量。
我们之前在表达 QTL 中已经展示了这一点。

Para_02

为了确保我们只包含真正独立的因果变异，我们将后续分析限制在那些在位点上覆盖良好的基因（排除的数量 = 1,247；补充表3和方法）。
我们还排除了在联合条件分析前后回归系数存在显著差异的蛋白质，因为这可能是未检测到的因果变异的标志（排除的数量 = 618；补充表4）。
后续分析显示，这些低覆盖率区域在之前报道过的含有片段重复的基因组区域中最为富集（方法），这些片段重复的比值比（OR）为1.10（1.07, 1.12），每重叠一个片段重复的P值为1.46 × 10−17（补充表3），这表明这些低覆盖率区域并非英国生物银行全基因组测序数据集所独有。
经过这些过滤步骤后，我们保留了1,026种蛋白质以进行进一步分析。

We identified 1,651 high-quality rare pQTLs for 599 proteins

我们鉴定了1,651个高质量的罕见蛋白质数量性状位点（pQTLs），涉及599种蛋白质

Para_01

在这些过滤之后，我们得到了5,076个顺式-pQTLs，其中包括98个结构变异。
对于给定蛋白质的所有独立顺式-pQTLs共同解释的平均方差（样本内）为4.10%（中位数为2.02%），与之前报道的估计值相似。
在确定的5,076个变异关联中，762个（15.0%）属于罕见频率范围（0.1% ≤ MAF < 1%），925个（18.2%）非常罕见（MAF ≤ 0.1%），其中包括36个结构变异。
我们将两个频率区间中的1,651个单核苷酸或短插入/缺失变异称为罕见pQTLs。

Para_02

我们还分别对1,184名和1,027名遗传学推断为南亚裔和非洲裔的个体进行了单变异测试，这些个体同时具有全基因组测序（WGS）和Olink蛋白质组学数据。
在考虑的三种遗传学推断的祖先来源中（欧洲（EUR）、南亚（SAS）和非洲（AFR）），我们在EUR和SAS个体之间观察到pQTLs效应量之间的强相关性（r = 0.797，P < 1 × 10−300），而在EUR和AFR个体之间以及AFR和SAS个体之间的相关性较弱（分别为0.656，P < 1.00 × 10−300和0.527，P = 6.79 × 10−250）。
尽管SAS和AFR分析的样本量小得多，我们仍分别识别出228个和396个独立的pQTLs。
其中，有77个（33.8%）和70个（17.7%）也在我们对EUR个体的分析中是主要的pQTLs。

Para_03

我们将我们的单变异 pQTL 结果与 Eldjarn 等人进行了比较，他们分析了从 150,119 例英国生物样本库全基因组测序数据中推算出的基因组数据，这些数据来自 54,306 名具有蛋白质组学数据的个体。
我们发现他们的 3,386 个顺式 pQTL（在我们测试的区域内；占比 76.1%）中，有 2,575 个直接映射到英国生物样本库 DRAGEN 全基因组测序调用结果，并且与我们针对相同循环蛋白的至少一个信号呈强连锁不平衡（LD）（r² ≥ 0.8）。
当仅考虑罕见变异（MAF < 1%）时，重叠比例更大（904 个中有 848 个；占比 93.8%）。
平均而言，每种蛋白质我们鉴定出 5.0 个（中位数为 4）独立的顺式 pQTL（包括罕见和常见变异），而 Eldjarn 等人鉴定出 10.1 个（中位数为 8）。

Most coding pQTLs reduce circulating cognate protein levels

大多数编码蛋白的数量性状位点（pQTLs）降低了循环对应蛋白的水平

Para_01

正如预期的那样，并且与最近的外显子组测序研究一致，96.3%被注释为高可信度功能丧失的 pQTLs（由 LOFTEE 定义）与循环蛋白水平降低相关联。（补充信息，补充表 7 和扩展数据图 2）。

Para_02

基于聚合的罕见编码变异测试确定了545个与539种循环蛋白水平相关的基因，这些结果在调整了我们单变异分析中常见的（>0.1%）pQTLs之后得出（补充表8）。
其中，有14个基因未被Dhindsa等人进行的外显子组测序分析发现（补充表9），这可能是由于测序覆盖范围的差异，从而展示了全基因组测序（WGS）发现外显子测序无法捕捉到的编码变异的能力。

Noncoding pQTLs are enriched in 5ʹ-UTRs

非编码pQTL在5ʹ-UTRs中富集

Para_01

我们发现了604个与369种蛋白质相关的独立罕见非编码单变异（图2和图3以及补充表10），与罕见编码变异（985个高置信度的失活或错义变异）进行了比较。
我们根据注释和最严重的后果将pQTLs分为同源组和非同源组，优先考虑相对于同源基因的注释。
在604个非编码pQTLs中，343个（56.8%）被注释为对同源基因具有调控作用（即位于UTRs或内含子5千碱基范围内；图3d）。
相比之下，972个（98.7%）罕见编码变异被注释为与同源基因相关。

Para_02

先前比较两种不同平台的研究表明，cis-pQTLs 的高比例可能是由表位结合伪影引起的。
尽管我们无法获得 SomaLogic 数据以与 Olink 进行直接比较，但我们注意到，我们的 604 个罕见非编码变异中没有任何一个与编码变异至少处于‘低’连锁不平衡 (r2 ≥ 0.1) 状态。
我们还发现，罕见非编码变异的影响大小与编码 pQTL 的最大 r2 值之间没有关联的证据（P = 0.801），即使在调整变异频率后仍然如此（P = 0.763）。
这表明，表位效应对我们所发现的非编码 pQTL 关联影响甚微。

Para_03

非编码变异的平均绝对效应为1.15个标准差（中位数为0.86个标准差），分别相当于罕见功能缺失和错义pQTLs平均绝对效应的57.1%和74.1%。
此外，罕见非编码pQTLs在增加循环蛋白效应和减少循环蛋白效应之间的分布更加均匀（65.2%为减少效应），而罕见编码pQTLs则有86.3%表现为减少效应（异质性P值=5.41×10^-13）。
在89.8%的情况下，关联最强的罕见pQTL位于或靠近其对应基因。
单独考虑时，罕见编码pQTLs和非编码pQTLs分别有97.9%和66.9%位于或靠近其对应基因。
罕见非编码pQTLs分布在顺式位点上，距离对应基因的最大距离为993 kb，接近所测试顺式区域的边界。

Para_04

接着，我们测试了在不同注释类别中，标记为调控其同源基因的稀有变异的富集情况。
基于每个测试变异最严重的预测后果（与同源基因相关的后果被优先考虑），我们在以下类别中观察到 pQTLs 的富集：5ʹ-非翻译区 (优势比 OR = 23.6, Fisher 精确检验双侧 P 值 = 4.46 × 10−48)，3ʹ-非翻译区 (OR = 3.08, P = 4.34 × 10−7)，预测的内含子剪接位点 (OR = 208.5, P = 1.08 × 10−29)，非编码外显子 (OR = 2.69, P = 9.02 × 10−5) 和上游变异 (OR = 3.00, P = 3.80 × 10−20; 图 4)。
当我们考虑标记为调控顺式窗口中另一基因的稀有非编码 pQTLs 时，没有观察到富集的证据 (扩展数据图 4)。

Fig. 4: Distribution of annotations of lead, rare, pQTL, noncoding variants compared with all variants tested.

- 图片说明

◉ 在先导变异集（n = 343；深蓝色）中与所有测试变异（N = 1,367,793；浅蓝色）相比，如果变异被注释到其同源基因的情况下的比例。◉ P 值来源于双侧 Fisher 精确检验。◉ whiskers 显示的是 95% 置信区间（CI）比例。◉ 加粗红色文本表示在 P < 0.05/9 的情况下显示出富集或耗竭证据的比较，其中 9 是统计检验的数量。◉ 非同源分布见扩展数据图 4。

Para_05

基于确定的1,651个稀有单变异，我们在P < 0.05/(1,651 × 2,907)的条件下，识别出与2,907种蛋白质之一相关的1,040个额外独特的远端pQTL关联（补充表11）。
需要注意的是，我们包含了之前因顺式全基因组测序数据质量而排除的蛋白质的远端效应，因为排除依据是顺式-WGS数据的质量，而非表型质量。

Aggregate testing identified 357 regulatory regions

聚合测试确定了357个调控区域

Para_01

通过基于聚合的关联性测试，我们发现了357个与循环蛋白水平相关的条件独立罕见变异非编码区域（补充表12和13）。
如果聚合单元被注释到对应的基因本身，或者聚合的起始位置位于对应基因单元内，我们将每个独立关联的聚合关联标记为映射到该对应基因。
所有其他聚合被标记为非对应基因，尽管其中许多可能比下一个最近的基因更接近对应基因。
与所有测试的对应基因聚合单元相比，我们观察到注释为5ʹ-非翻译区（5ʹ-UTR）的对应基因聚合体富集的证据（图5，优势比OR = 6.36，P值 = 5.11 × 10⁻¹¹），以及预测为内含子剪接位点的变异体聚合体的富集（OR = 29.6，P值 = 4.89 × 10⁻¹⁷）。

Fig. 5: Distribution of annotations of significant cognate aggregates of rare noncoding compared with all aggregates tested.

- 图片说明

◉ 在铅变异集合中聚集物的比例（n = 125；深蓝色）与所有测试的聚集物（N = 1,456,986；浅蓝色）相比。◉ P 值来源于双侧 Fisher 精确检验。◉ 95% 置信区间的比例由须状图表示。◉ 加粗红色文本表示在 P < 0.05/10 的情况下，存在富集或耗竭的比较结果，其中 10 是统计检验的数量。◉ 非同源比较结果展示在扩展数据图 5 中。

Para_02

超过90%的非编码聚合信号是在不限制变异达到我们规定的保守性或约束阈值的情况下检测到的。
在357个（占5.04%）条件独立的非编码聚合关联中，有18个仅在选择高度保守（GERP > 2）的变异时被识别，而16个（占4.48%）仅在选择高度受限（JARVIS > 0.99）的变异时被识别。
大多数聚合区域如果仅通过单变异分析将会被遗漏：259个（占72.6%）非编码聚合中没有主要的pQTL，74个（占20.7%）不含任何全基因组显著性的单变异。
基于357个独立的非编码聚合，我们还确定了45个独特的跨非编码聚合-蛋白关联（P < 0.05/(357 × 2907)）

Para_03

绝大多数罕见的非编码聚集体是通过允许双向关联的统计测试识别出来的，并且假设很大一部分变异是非因果的。
在357个罕见非编码聚集体pQTLs中，仅有10个（2.80%）在负担框架下最显著相关（假设所有罕见变异都会产生相同方向的影响），这与基于编码区的聚集体形成强烈对比，后者的比例为38.6%。
这一差异表明，非编码区域中的罕见变异很可能导致性状增加和减少的混合效应，并且并非聚集体中的所有变异都是因果性的。
这一观察结果表明，改进对非编码变异的注释将进一步提高检测非编码聚集信号的机会。

Para_04

我们发现滑动窗口聚合检验能够识别出那些仅通过限定在注释区域内的检验无法检测到的关联。
例如，有34种（3.31%）循环蛋白具有滑动窗口pQTL。
对于这34种蛋白，当包含滑动窗口时发现了109个非编码pQTL，而不包含滑动窗口时则只有81个。
对于三种蛋白，唯一的pQTL是滑动窗口。
例如，在ACAT-V框架下，LAMP3基因内含子区域chr3:183135000-183137000与循环LAMP3水平相关联（P = 1.15 × 10−13），但该区域不包含任何Ensembl调控区域。
我们的结果表明，相对于仅测试先前已映射的区域，滑动窗口具有额外的价值，尽管我们未包含的细胞和组织特异性数据中的注释可能可以检测到这些区域。

Rare noncoding pQTLs showed tissue-specific enrichment

罕见的非编码 pQTL 显示出组织特异性的富集

Para_01

接下来，我们旨在确定罕见的非编码pQTL在组织相关的非编码调控区域中的存在程度。
首先，我们测试了这样一个假设：如果相关蛋白是分泌蛋白或信号蛋白，则更有可能识别出罕见的非编码pQTL，因为循环蛋白相比非分泌蛋白更能代表整体蛋白丰度。
其次，我们基于以下假设测试了另一个假设：在血液和肝细胞的Ensembl调控元件中，罕见的非编码pQTL会比其他20种组织类型更丰富（方法），因为这些是最相关的组织类型。
我们发现，在所有调控区域中，单变异pQTL的关联在所有蛋白组中高度富集，无论其分泌状态如何（图6、扩展数据图6和7以及补充表15-20）。
在增强子和CTCF结合位点内，分泌蛋白或信号蛋白的关联比非分泌蛋白更高度富集。
在非编码聚合关联中，并且仅限于通过滑动窗口识别出的关联，分泌蛋白和信号蛋白在所有调控区域中的富集程度高于非分泌蛋白。
在滑动窗口内，关联在预测活跃于血管的调控区域中最高度富集，其次是肝脏（补充表20）。

Fig. 6: QQ plot for enrichment of loci within Ensembl-predicted active regions within tissue groups.

- 图片说明

◉ 经验性的单侧 P 值用于评估在 22 个组织群组中，Ensembl 预测的活跃区域内的信号富集情况，按蛋白质分泌状态分类。◉ a,d,g，单个变异体的富集情况。◉ b,e,h，基于 Ensembl 调控区域的聚合测试的富集情况。◉ c,f,i，基于滑动窗口的聚合测试的富集情况。◉ a-c，所有预测活跃区域内的富集情况。◉ d-f，启动子区域的富集情况。◉ g-i，增强子区域的富集情况。◉ 每个图中从上到下的黑色曲线分别表示在零假设下预期 P 值分布的 95% 置信区间。◉ 每个注释类别中的蛋白质可以在补充表 21 中找到。◉ P 值是根据采样变体与注释区域重叠程度的比例计算得出的。◉ 蓝色点表示分泌或信号蛋白，红色点表示非分泌蛋白组。◉ 肝细胞和血管细胞用圆圈标记，其他所有细胞类型用三角形标记。

Para_02

作为敏感性分析，我们尝试缓解这样一个担忧，即我们的发现可能是由用于测量循环蛋白质的平台的技术人工产物所驱动的。
在551种被确定为Olink（Explore 3072）和SomaScan-v.4之间高度一致的蛋白质中，有261种具有高质量的顺式全基因组测序覆盖。
我们的pQTL关联在这261种蛋白质中富集：尽管它们仅占我们检测的蛋白质的25.4%，但却涵盖了我们所有pQTL的31.1%（两比例检验，P = 2.85 × 10−4），稀有pQTL的32.7%（P = 6.05 × 10−5）以及基于非编码聚合的pQTL的37.5%（P = 1.29 × 10−5）。

Discussion

Para_01

以循环蛋白质水平为例，我们已经证明了全基因组测序（WGS）数据分析能够发现与常见表型相关的罕见非编码变异及其聚合体。
通过全基因组测序数据，利用罕见变异聚合测试，我们可以考虑的变异数量比单独进行单个变异测试时多出一倍以上。

Para_02

然而，我们的结果还表明，蛋白质的‘统计独立’pQTL数量与顺式测序覆盖率呈负相关。
与以往研究不同，我们因此尝试考虑那些未完全捕捉到的区域：仅约三分之一的被测量蛋白质在其位点上具有我们认为高质量的全基因组测序数据。
尽管我们在排除的区域中观察到之前报道为问题区域的富集现象，这表明这些问题并非英国生物样本库（UKB）独有，但要充分理解每个顺式区域中观察到的低覆盖率仍需进一步分析。
这是对其他将使用此规模测序数据的研究者的一个关键警示点。

Para_03

我们已经确定了数百个新的非编码罕见聚合和单变异体关联，这些关联与认知基因两侧 1 Mb 范围内的顺式窗口中测量的蛋白质水平相关。
我们发现非编码关联的效果大小有时与编码关联的大小相似，但在促进蛋白质增加和减少之间的平衡性更强。
由于连锁不平衡 (LD) 和单倍型效应的复杂性质，要确定是编码信号驱动非编码信号还是反之非常困难。
为避免这种解释，我们采取了保守的方法，并对认知基因的所有编码变异进行了条件控制。

Para_04

我们在单变异 pQTLs 数量方面观察到与 Eldjarn 等人的研究结果存在一些差异。
这些差异可能部分由条件分析方法驱动：他们的分析使用了逐步向前的条件分析来定义条件独立的 pQTLs，与此相关的问题之前已被详细讨论过（例如，参考文献33），而我们则在 GCTA-CoJo 中实施了向前和向后的条件分析步骤。
这两项研究中 pQTLs 的差异突显了解释同一位点上的多个独立 pQTLs 的困难。

Para_05

我们发现5ʹ-非翻译区和预测的内含子剪接受体或供体位点富含罕见的非编码pQTLs。
由于UTRs和内含子通常不会在靶向编码外显子的外显子组测序中被捕获，我们的结果突显了全基因组测序在寻找新的罕见基因相关变异中的重要性。

Para_06

我们还展示了在非编码区域进行聚合测试的力量，这类似于已在编码区域中功能相似变异聚合的成熟方法。
通过测试稀有的基因组非编码元素聚合体，这些聚合体按照与基因的距离、预测的调控活性或使用滑动窗口等方式分组，我们发现了另外464个条件独立的区域，而这些区域无法仅通过单变异测试识别。

Para_07

与基于聚合的编码关联相比，非编码基因组聚合 pQTLs 在允许双向效应和/或稀疏因果关系的测试中更为丰富。
这一观察结果与以下事实一致：在非编码区域中，变异效应和功能区域的预测精度较低。
然而，我们利用当前注释和数据识别出非编码关联这一事实表明，随着功能注释的改进和群体遗传数据的积累，常见表型中可能会有更多发现。

Para_08

我们的研究存在一些局限性。
首先，由于我们无法获取其他研究的类似数据，因此无法在独立数据集中重复我们的结果。
然而，我们的大部分关联达到了远超阈值的统计置信水平。
此外，尽管我们的主要分析因样本量限制而局限于欧洲血统个体，但效应大小在不同血统群体中是一致的。
其次，我们不能完全确定已经考虑到所有可能由与编码变异或常见变异的连锁不平衡（LD）导致的残余混杂因素，包括复杂或重复区域中的变异。
出于类似原因，我们在同一区域内识别出的部分蛋白质数量性状位点（pQTLs）可能并非真正独立。
然而，由于这些关联具有不同的特征，包括对蛋白质增加和减少效应的更均匀分布，因此它们不太可能受到编码变异残余混杂的显著影响。
我们也没有观察到独立非编码 pQTLs 的效应大小与编码变异的残余 LD 之间存在任何关系，这表明我们的结果不太可能受到表位效应的干扰。
第三，所有循环蛋白均在血液中测量。
虽然大部分组织特异性蛋白仅在特定组织中表达，但我们受限于只能考虑循环蛋白水平。
第四，我们无法考虑与测量蛋白水平所使用的技术相关的结合效应。
然而，在 Olink 和 SomaLogic 平台之间显示出强相关性的 551 种黄金标准蛋白的子集中，我们确实观察到了关联的富集现象。
最后，尽管我们指出滑动窗口聚合识别了潜在的新调控区域，但不可能全面检查所有公开可用的调控图谱。
此外，测试所有可能的组织特异性注释会带来巨大的假设检验负担。

Para_09

总之，以循环蛋白质测量的示例特征为例，我们发现了循环蛋白质水平与罕见非编码变异之间的多种新关联。
我们的结果表明，可能有许多具有重大影响的罕见非编码变异等待被发现，这些变异会影响复杂的表型。

Methods

Para_01

本研究遵守所有适用的伦理规定。
英国生物银行（UKB）研究已获得西北研究中心伦理委员会的伦理批准。

UKB and WGS

英国生物银行（UKB）和全基因组测序（WGS）

Para_01

用于 UKB 的全基因组测序（WGS）通过 Illumina NovaSeq 6000 测序仪完成，平均覆盖度为 32.5×，采用 150-bp 双端测序读长。
测序所使用的基因组构建为 GRCh38，单核苷酸多态性和短插入缺失变异是通过 DRAGEN 3.7.8 联合分析得出的。

Human protein expression levels

人类蛋白质表达水平

Para_01

使用 Olink 技术，对多达 54,219 名 UKB 参与者的 2,932 种蛋白质的水平进行了分析，具体方法参见参考文献 17，由 UKB Pharma 蛋白组学项目执行。
在数据提供给研究人员使用之前，进行了质量控制程序，包括去除异常值。
在数据发布前，蛋白质水平还经过了 log2 转换。
经过质量控制过滤后，共有 54,189 名具有蛋白质表达数据的个体被批准用于分析。
Sun 等人未发现批次或平板混杂效应的证据。

Genetic data filtering

基因数据过滤

Para_01

如果每个样本基因型的 sum(LAD) < 8（其中 LAD 是局部等位基因深度）或基因型质量 (GQ) < 10，我们将任何 DRAGEN WGS 基因型调用设置为缺失，这适用于 UKB 使用 bcftools v.1.2 提供的 154,430 个项目变异调用格式 (pVCFs)。
在这些额外的质量控制步骤之后，单例的传递率（理论上应精确为 0.5，假设大多数变异未受到强烈负向选择）为 0.497，而 UKB 最初提供的数据为 0.456。
随后，我们删除了缺失率 >10% 的任何变异。
应用多等位分裂程序后，每个变异被分配一个唯一 ID (CHR:BP:REF:ALT)，然后基于（访问日期：2024 年 3 月 30 日）提供的 1000G b38 参考序列，使用 bcftools 对所有每条染色体上的 VCF 文件进行合并，并对插入/缺失进行了规范化和左对齐处理。
每个合并后的 pVCF 随后被转换为 plink（v.2.0）p(gen/var/sam) 格式。

Structural variants

结构变异

Para_01

在分析时，DRAGEN版本的结构变异（变异≥50 bp）仅以单样本格式提供。
我们使用 SURVIVOR v.1.07 与默认设置，将结构变异调用合并并 harmonize 成包含 54,219 名具有循环蛋白测量值个体的 pVCF 文件。
然后，我们将该项目 (p)VCF 转换为 plink p(gen/var/sam) 格式，并与 DRAGEN SNV 或 indel 调用的 plink 文件合并。

Genetic variant annotation

基因变异注释

Para_01

我们使用 Ensembl VEP（v.110）、LOFTEE 和 UTRannotator 对所有遗传变异进行了注释。
在可能的情况下，我们将每个变异分配到以下三个类别之一：编码区、近端调控区或基因间调控区。
如果一个变异对任何转录本的编码序列有预测影响，则被归类为编码区变异；如果变异位于转录本非翻译区（UTRs）5千碱基对范围内的区域，并且在任何转录本中不已经是编码区变异，则归类为近端调控区变异；最后，如果它既不是编码区变异，又映射到一个与基因无关的调控元件（详细信息见下文），则归类为基因间或内含子调控区变异。
此外，我们还对大小为 2,000 碱基对的滑动窗口中的变异进行了测试，无论每个窗口中有多少个变异，但排除了编码区变异以尽量减少假设重叠。

Para_02

然后，我们根据预测的后果和位置，将每个变异分配到我们称为掩码的分组中。
我们使用了五种已发表的变异评分来按后果对变异进行分组：
（1）GERP：GERP 分数是衡量变异水平保守性的指标。
如果一个变异的 GERP 分数 > 2，则我们认为它是高度保守的。
（2）PhastCon 分数：PhastCon 是基于窗口的物种间保守性度量方法，可以仅针对哺乳动物（phastCon 30）或所有物种（phast_100）。
我们测试了非编码基因组窗口，即排除任何包含外显子的窗口，并选择了 PhastCon 分数处于第 99 百分位的窗口。
（3）约束分数：约束是根据每个窗口的局部可变性和观察到的突变率在 1 kbp 窗口大小内计算得出的。
我们测试了约束 z 分数 ≥ 4 的窗口。
（4）SpliceAI 分数：SpliceAI 分数是一种衡量预信使 RNA 区域内每个变异是否为剪接供体、受体或两者都不是的可能性的指标。
如果一个变异的 AI 分数 > 50，则将其分类为高置信度的剪接位点。
（5）CADD：CADD 分数用于预测一个变异可能有多大的破坏性。
我们仅将 CADD 分数应用于编码变异，并且只有被 VEP 标记为高置信度的失活功能变异才会被考虑。
错义变异中 CADD 分数 > 25 的会被分离出来在单独的掩码中进行测试。
（6）JARVIS 分数：JARVIS 分数是基于约束测量的机器学习模型推导出来的，旨在更好地优先考虑关联研究中的非编码遗传变异。

Para_03

每个基因组掩码由多个变异组成，这些变异具有不同的后果，基于它们的位置、上述得分之一和/或预测的编码后果。
例如，要将一个变异分类为错义 CADD > 25，它必须改变基因转录本外显子的密码子，并被预测为高度有害。
在补充表 22 中，我们列出了分配给每个掩码和分类的完整后果列表。

Association analyses

关联分析

Para_01

我们对UKB中测量的2,932种蛋白质的顺式位点内的单变异和聚合变异进行了测试。
为了定义顺式窗口，我们首先将每种蛋白质映射到一个编码基因（少数排除情况见补充表1），并确定了Ensembl记录的每个基因的最长转录本。
基于最长转录本，我们将顺式窗口定义为转录本基因5ʹ-和3ʹ-非翻译区（UTRs）两侧各1 Mb的范围（受染色体起始和终止位置限制），以及编码序列和内含子序列中的变异。
所有关联分析均校正了年龄、性别、年龄平方、UKB招募中心（作为地理的代理）、前40个遗传主成分、全基因组测序批次、Olink板、禁食时间和采血后时间。

Single-variant association testing

单变量关联测试

Para_01

为了识别与蛋白质水平相关的顺式单个变异，我们首先使用 regenie v.3.3（参考文献41）在顺式窗口中对所有遗传变异进行了关联测试，这些变异的最小等位基因计数（MAC）至少为五。
然后，我们使用修改版的 GCTA-CoJo19（diff-freq = 0.2, cojo P = 2.95 × 10−10）在条件联合分析中选择了主要变异，以英国生物银行全基因组测序数据为连锁不平衡（LD）参考面板，并限制为具有蛋白质组学数据的个体。
测试显示，如果某个变异所解释的方差大于任何独立变异所能解释的最小方差的 900 倍，则 GCTA-CoJo 会过滤掉该变异（‘sqrt(ldlt_B.vectorD().maxCoeff()/ldlt_B.vectorD().minCoeff()) > 30’：gcta/meta/joint_meta.cpp 的第 732 行，网址为，访问日期为2024年3月17日）。
我们理解此过滤器的存在是为了捕捉由共线性引起的统计混淆，例如，如果用于计算 LD 的参考基因组和遗传数据之间的相关性较差时就会出现这种情况。
然而，在我们共同考虑常见和罕见变异的目的下，由于我们使用的 LD 参考面板与我们的发现数据集完全匹配，我们发现此过滤器错误地去除了效应较大的 pQTLs。
因此，我们移除了该过滤器并重新编译了 GCTA-CoJo，其代码可在获取。

Para_02

GCTA-CoJo 在超过某个阈值后假设 P = 0。
在这种情况下，我们根据以下 R 脚本重新计算了 log10(P)：log10(P-计算值) = log10(exp(1)) × (pt(−abs(β/标准误差), 自由度 = N-2, log = TRUE) + log(2))。

UKB WGS coverage calculation and filtering

英国生物银行全基因组测序覆盖度计算与过滤

Para_01

我们根据质量控制标准，从 UKB 全基因组测序数据中计算了覆盖率，涉及998名推断为欧洲血统的个体，这些个体同时具有蛋白质组学数据和全基因组测序数据。
我们首先计算了每个碱基对上深度大于8的样本数量。
随后，我们计算了每个同源蛋白编码基因在5ʹ-非翻译区和3ʹ-非翻译区之间碱基数的比例，要求超过90%的个体在这些区域的深度大于8。
我们对完整的顺式窗口重复了这一分析。
然后，一个蛋白的基因水平覆盖率和区域水平覆盖率均需满足：超过99.5%的碱基和超过90%的个体深度大于8。

Para_02

为了考虑可能遗漏或错误基因分型的单一有害编码位点，如果任何与该蛋白质相关的 pQTL 的边际效应估计值（来自我们的 regenie 分析）和联合效应估计值（GCTA-CoJo）之间的最大差异超过第 90 个百分位数，我们还会额外过滤这些蛋白质。
扩展数据图 8 提供了描述完整过滤流程的流程图。

Para_03

然后，我们测试了我们的结果与因低覆盖率而被排除的基因组区域之间的重叠，这些区域来自一个有问题的基因组区域数据库，该数据库由 Genome in a Bottle 和 precisionFDA Truth Challenge v.2 创建（版本 v.3.5，2024 年 9 月 10 日访问）。
这些区域包括"串联重复序列、所有大于 6 个碱基对的均聚物、所有大于 10 个碱基对的不完美均聚物、所有难以映射的区域、所有片段重复、GC 含量小于 25% 或大于 65% 的区域、‘问题启动子’、chrX/Y XTR 和扩增子、卫星序列以及‘其他困难’区域（仅针对 GRCh38 包括 T2T-consortium 提供的区域）"。
对于每种蛋白质，我们确定了至少有 1 个碱基对与基因区域（从 5ʹ- 到 3ʹ-UTR）或 2-Mbp 窗口重叠的每种类型的问题区域的数量。
这些变量在逻辑回归模型中进行了测试，以二元变量为依据描述该蛋白质是否因低覆盖率而被过滤，并调整了 2-Mbp 窗口的总长度（该长度自然随基因大小而变化）。

Rare variant genomic aggregate testing

罕见变异基因组聚合检测

Para_01

为了识别单变量分析能力不足的基因组非编码潜在调控区域，我们随后进行了非编码罕见变异（MAF < 0.1%）的基因组聚合关联测试。
为了检验非编码罕见变异聚合信号是否由与常见变异和/或单变量信号的残余连锁不平衡（LD）及单倍型结构引起或混淆，我们对每个达到 Bonferroni 校正 P < 0.05 的罕见变异聚合测试结果进行了以下步骤：
（1）为了生成我们的主要非编码发现结果，我们调整了在联合（CoJo）分析中被确定为独立信号的常见主效变异（MAF > 0.1%，~MAC = 40），并调整了所有被注释为该基因任何转录本中编码序列的遗传变异（无论其 P 值如何），这些基因映射到感兴趣的蛋白质上；
（2）为了识别独立的非编码聚合关联，如果至少一个聚合通过了我们的显著性阈值，则我们进行逐步回归分析。
对于每种蛋白质，从最显著相关的非编码聚合（按 P 值）开始，我们对达到全基因组显著性的聚合进行了额外的非编码聚合测试，并调整了顶级信号中的所有变异。
此过程重复进行，每次加入更多来自下一个最强相关聚合的变异，直到没有聚合达到全基因组显著性为止。
（3）为了确定我们的主要聚合发现结果是否可能是由单一低频主效变异引起的，我们识别了包含正好一个主效遗传变异的聚合关联。
（4）作为敏感性分析步骤，为了评估这些结果可能在多大程度上受到连锁不平衡（LD）的混淆影响，我们进一步调整了所有已识别的 pQTL 单变量。

Para_02

对于最大等位基因频率阈值为 0.1% 的变异，使用 regenie 进行了基因单元测试，测试基于补充表 22 中指定的遗传单元。
Regenie 执行四种类型的基因组单元测试：(1) 标准 BURDEN 测试，假设给定基因单元掩码中的每个变异对表型的影响大小和符号大致相同；
(2) SKAT 测试，允许单元中每个变异的关联符号有所不同；
(3) ACAT 测试，其中单元中每个变异的关联符号可以不同，并且只需要掩码中的少量变异具有关联性；
(4) ACAT-O，这是一种综合测试，结合了 BURDEN、SKAT 和 ACAT 测试，旨在最大化这三种测试的统计功效。

Para_03

我们对每个至少有一个变异的基因单元所对应的每个掩码进行了上述四种统计检验中的每一种。
此外，还对每个单元中的所有单例变异（MAC=1）进行了关联性测试；regenie 还为每个基因组单元估算了一个‘全掩码’关联强度，这是个体掩码测试统计量的汇总结果。
为了确保这不会导致非编码和编码关联统计量的混杂，我们将每个基因转录本拆分为一个编码转录本，用于测试所有编码掩码，以及一个近端转录本，用于测试所有近端掩码。
调控基因组单元根据其 ENSR 分配、1kb 约束窗口或 phastCon 保守窗口的程度进行分类。
我们根据它们覆盖的相应染色体区域命名滑动窗口掩码。
执行主要分析所需的代码可在 Zenodo 上获取。

Statistical significance

统计显著性

Para_01

统计学显著性是根据对20个随机生成的正态分布连续特征进行全基因组测序（WGS）分析所观察到的最小P值来定义的。
单变异体和聚合关联分析的最小P值被视为独立：单变异体的P值为2.95×10^-10；聚合分析的P值为8.71×10^-9。

Ensembl regulatory region enrichment

Ensembl调控区域富集

Para_01

我们计算了单个变异和聚合区域与Ensemble31调控区域的重叠富集程度，这些调控区域涵盖了118种组织或细胞系，并与遗传背景进行了比较。
对于每种组织，Ensembl还提供了每个区域是否处于活跃状态（或不活跃、受抑制等）以及调控活动类型（启动子、增强子、CTCF结合位点、转录因子结合位点或开放染色质区域）的预测。
随后，我们仅考虑被预测为活跃的区域，排除了细胞系和癌症衍生的组织，并将剩余的组织类型分组为22个超组。
我们使用DAVID软件将蛋白质分为分泌型、信号型或膜结合型蛋白质。
然后，我们将蛋白质归类为分泌型或信号型蛋白质，如果它们被注释为属于任一组，并将所有其他蛋白质标记为非分泌型。

Para_02

为了确定统计上的富集程度，我们进行了10,000次模拟的引导抽样。
对于每次模拟，我们随机选择与我们从基因组顺式调控区域中确定的独立信号数量相等的稀有非编码变异或聚合体，这些信号由我们的稀有非编码单变异或聚合体测试集合确定。
然后，我们确定随机选择的变异或聚合体集合与任何调控区域的重叠情况（我们对每个面板和组织类型的分层重新执行此操作），并将任何模拟中的重叠数量分布与我们在独立关联中观察到的重叠数量进行比较。
随后，我们为观察到的重叠分配了一个经验P值。

Statistics and reproducibility

统计学与可重复性

Para_01

本研究纳入了所有在 UKB 中同时测量了蛋白质水平和全基因组序列的个体。
参与者在研究招募时年龄介于37至80岁之间。
对遗传数据进行主成分分析（PC分析），以根据遗传祖先定义同质群体。
参与者按照遗传特征被分为欧洲样、南亚样和非洲样的遗传祖先群体。
在每个群体内分别进行分析，并在混合模型中调整主成分，以考虑隐性相关性。
本研究为观察性研究，因此不适用随机化或绑定。

Reporting summary

报告摘要

Para_01

有关研究设计的更多信息可在此文章链接的《自然组合报告摘要》中获取。

Data availability

Para_01

由于 UKB 的数据可用性和数据返还政策，数据不能公开共享。
符合访问其数据集条件的研究人员可以从 UKB 获得数据（）。
汇总统计数据可通过 Zenodo 在 .5281/zenodo.14203628 上获取。

Code availability

Para_01

与分析相关的代码可通过 GitHub 获取（）。

Basic Information

英文标题：Whole-genome sequencing analysis identifies rare, large-effect noncoding variants and regulatory regions associated with circulating protein levels
中文标题：全基因组测序分析鉴定出与循环蛋白水平相关的罕见、大效应非编码变异和调控区域
发表日期：24 February 2025
文章类型：Article
所属期刊：Nature Genetics
文章作者：Gareth Hawkes | Robin N. Beaumont
文章链接：

Abstract

Para_01

罕见非编码遗传变异对常见表型的贡献在很大程度上是未知的，这是由于历史上缺乏全人群规模的全基因组测序数据，以及难以将非编码变异分类为功能相似的组别。
为了开始解决这些挑战，我们使用全基因组测序数据进行了顺式关联分析，该数据包括11亿个变异，1.23亿个基于非编码聚合的测试和约5万名英国生物银行参与者的2907种循环蛋白水平。
我们发现了604个与循环蛋白水平相关的独立罕见非编码单变异关联。
与蛋白质编码变异不同，罕见非编码遗传变异几乎同样可能增加或减少蛋白质水平。
罕见非编码聚合测试确定了357个条件独立相关区域。
其中，74个（占21%）仅通过单变异测试无法检测到。
我们的研究结果对识别与常见人类表型相关的罕见非编码遗传变异及其作用具有重要意义，包括测试非编码变异聚合的重要性。

Main

Para_01

人类基因组中非编码区域的稀有遗传变异可导致严重的罕见病，但它们在常见复杂性状中的作用仍然 largely 未知。
基于阵列的推算和全基因组关联研究 (GWAS) 已经发现了数以万计的人类疾病与常见变异的关联，其中大多数位于编码区域之外。
然而，利用外显子组测序数据识别与常见表型相关的稀有变异的努力主要局限于编码区域，这些研究结合了单变异测试（统计功效低但功能解释更清晰）和聚合测试（统计功效更高但难以解释）。
我们评估非编码基因组中稀有变异的能力直到最近还受到限制，这是因为缺乏大规模研究中的全基因组测序 (WGS) 数据，并且难以定义具有生物学意义的非编码调控基因组单元。

Para_02

利用全基因组测序（WGS）数据识别非编码调控元件，可以为基因调控提供重要的见解，这补充了从外显子组测序和基于阵列的研究中获得的知识。
根据群体遗传学指标（如约束性），功能性非编码DNA的数量估计是编码序列的4到5倍，且分别有10%的启动子和6%的增强子受到与编码区域相当的突变约束。
全基因组测序使我们能够研究内含子、近端和远端调控元件的作用，并覆盖人群中的整个等位基因频率谱，包括那些即使在非常大的样本中也仅被观察到一次或两次的变异比例。

Para_03

在常见表型的背景下，关于全基因组测序（WGS）数据的研究非常少。
最近来自 TOPMed 的例子考虑了血脂水平（样本量 N = 66,000）和血压（样本量 N = 51,456），以及英国生物样本库（UKB）中的一项关于身高的研究（发现样本量 N = 200,003）：这些研究发现的新信号数量有限，可能是因为样本量相对较小，不足以检测新的罕见变异。
此外，WGS 还被用于研究多达 3,000 名个体的蛋白质水平相关的遗传变异。

Para_04

英国生物银行（UKB）发布循环蛋白数据，并结合全基因组测序（WGS），为测试罕见非编码遗传变异对常见、生物学上接近且测量良好的人类表型的影响提供了前所未有的机会。
三篇最近发表的研究描述了这些数据在2023年的发布，涉及54,306名个体中的2,923种循环蛋白。
首先，Eldjarn等人通过单变量分析识别出30,062个蛋白数量性状位点（pQTLs），这些数据来自UKB的150,119例全基因组测序分析，包含2,931种测量的蛋白水平，并将结果与冰岛全基因组测序队列的蛋白质组学数据进行了比较。
其次，Dhindsa等人通过外显子组测序分析识别出5,433个pQTLs，并在编码区域内进行聚合测试，鉴定出1,962个基因-蛋白关联关系。
最后，Sun等人通过基于传统插补的全基因组关联研究（GWAS）识别出14,287个pQTL单变量。

Para_05

利用全基因组测序（WGS）数据和循环蛋白水平作为示例特征，我们测试了两个相关假设：（1）罕见的非编码单基因变异体，目前通过推算或外显子组测序无法检测到，对常见人类表型有与编码变异相似的影响；（2）基因组调控区域中的罕见非编码遗传变异集合（群体）与人类表型相关，类似于编码序列中的基因水平负担分析。
重要的是，与之前关于英国生物银行（UKB）蛋白质组数据的三篇论文不同，我们使用了所有具有循环蛋白测量值的参与者的短读长全基因组测序检测到的DNA序列变异的完整等位基因频率范围，提供了约11亿个变异的信息，但将每种表型的搜索限制在衍生该蛋白的蛋白编码基因周围的顺式调控区域。

Para_06

我们使用来自英国生物样本库（UKB）的46,362名推断为欧洲遗传背景个体的注释全基因组测序（WGS）数据，对2,903种测量的循环蛋白水平进行了初步发现关联分析。
此外，我们分别在899名和1,098名推断为南亚和非洲遗传背景的个体中进行了单变异分析。
我们的分析鉴定了1,651个高质量单变异，对比了编码和非编码变异的影响，通过识别出357个独立位点突显了罕见变异非编码聚合检测的能力，并展示了在基于全基因组测序的关联研究中考虑测序覆盖度的重要性。

Results

Para_01

排除了十二种蛋白质，这些蛋白质或者是融合蛋白，或者无法直接与 HGNC 基因符号匹配（补充表 1）。
对于每种测量的蛋白质，我们在编码该蛋白质的基因周围的顺式窗口中进行了单变异体（小等位基因计数 (MAC) ≥ 5）和基因组聚合关联测试（小等位基因频率 (MAF) < 0.1%），顺式窗口从 5′- 和 3′-UTR 扩展 1 Mb，基于该基因任何 GENCODE 转录本的最极端 5′- 和 3′-末端。
我们使用 1 Mb 作为最近确定的边界距离，当一个 pQTL 更可能是一个顺式而不是反式关联时即适用此距离。
在运行时对循环蛋白质水平测量值进行秩逆归一化，并将年龄、年龄平方、性别、招募中心、全基因组测序 (WGS) 中心、采血后时间、禁食时间、40 个遗传主成分 (PCs) 和 Olink 批次 ID 包括为协变量（方法）。
尽管仅考虑每个蛋白质的 2-Mbp 顺式窗口，我们仍测试了 128,434,590 个单变异体（包括单核苷酸变异体和小插入或缺失 (indels)）和结构变异体关联，其 MAC ≥ 5。
我们使用 GCTA-CoJo 的修改版本识别独立变异体关联（方法）。

Para_02

我们使用 Ensembl 的 Variant Effect Predictor v.110.1 (VEP) 对所有遗传变异进行了注释（方法），并利用其输出将变异分类为基因中心型（例如，编码、预测的内含子剪接、未指定内含子和近端调控）和基因间调控型（例如，Ensembl 调控区域、非编码 RNA 和未指定基因间）以进行基于聚合的关联测试。
如果非编码变异位于基因的 UTR 区域或内含子区域内，或者距离 UTR 不超过 5 kbp，则将其注释为基因中心型。
此外，我们在重叠（1-kb 重叠）的 2-kb 滑动窗口中对所有非编码变异进行了聚合测试。
我们还通过约束度量（JARVIS）、保守性（基因组进化速率分析 GERP）和/或预测的有害性（组合注释依赖缺失评分 CADD）进一步对部分聚合单元内的变异进行了细分。
为了识别独立的、罕见的、非编码聚合关联，我们对非编码聚合测试调整了常见 pQTLs（MAF > 0.01）以及与蛋白自身编码基因相关的所有编码变异，无论变异频率如何，该基因此后称为同源基因（图 1）。
总共，我们对 123,598,575 个聚合进行了关联测试，其中包括 182,136,116 个 MAF < 0.1% 的变异，直到包括单例为止。
统计显著性基于模拟研究定义（方法）：单个变异 P < 2.95 × 10−10，基因组聚合 P < 8.71 × 10−9。

Fig. 1: Study design for pQTL single and aggregate genetic variant analysis.

- 图片说明

◉ 展示研究设计的流程图。

We identified 13,457 candidate pQTLs across 2,891 proteins

我们在2891种蛋白质中鉴定出13457个候选pQTLs

Para_01

对于2,891种循环蛋白质，我们确定了13,457个统计上独立的顺式-pQTL关联（MAC ≥ 5），其中293个为结构变异（图2和图3以及补充表2）。
我们发现了2,036种蛋白质中至少存在一个顺式关联，每种循环蛋白质中位数为5个独立pQTL；429种蛋白质与超过10个独立pQTL相关联，其中CD177的最大值为49个。
这些结果与之前使用150,000个英国生物样本库基因组的研究一致，并利用相同的数据对剩余样本进行了推算。

Fig. 2: Manhattan results of single-variant and aggregate pQTL analyses.

- 图片说明

Fig. 3: Effect size distributions of rare, independently associated pQTL variants.

- 图片说明

Over half of independent pQTLs may be coverage artefacts

超过一半的独立pQTL可能是覆盖范围伪影

Para_01

由于我们的研究旨在评估非编码变异对蛋白质水平的作用，我们将所有全基因组显著的变异（包括编码和非编码）共同建模（方法），以避免单倍型效应带来的混淆。
然而，我们随后发现，具有最大数量统计独立 pQTL 的蛋白质在对应基因的低覆盖区域中显著富集（补充表 2 和 3）。
例如，CD177（MANE Select 转录本 ENST00000618265.5）有 49 个独立的 pQTL，但在外显子 5 或内含子 6 上几乎没有覆盖（扩展数据图 1）。
此外，那些其对应基因表达主要由重复多态性的拷贝数驱动的蛋白质也显示出富集现象。
例如，FCGR3B（43 个独立信号）、AMY2A（25 个独立信号）和 AMY2B（27 个独立信号）是由已知具有高度可变拷贝数的基因编码的，其中淀粉酶水平主要由后两个基因的拷贝数决定。
此外，脂蛋白(a) 水平（34 个独立信号）主要由 Kringle 的重复多态性解释。
对于这些基因中的每一个，独立信号数量较大的最可能解释是真正的因果变异尚未被充分校正，从而导致了过多的人工独立 pQTL 数量。
我们之前在表达 QTL 中已经展示了这一点。

Para_02

为了确保我们只包含真正独立的因果变异，我们将后续分析限制在那些在位点上覆盖良好的基因（排除的数量 = 1,247；补充表3和方法）。
我们还排除了在联合条件分析前后回归系数存在显著差异的蛋白质，因为这可能是未检测到的因果变异的标志（排除的数量 = 618；补充表4）。
后续分析显示，这些低覆盖率区域在之前报道过的含有片段重复的基因组区域中最为富集（方法），这些片段重复的比值比（OR）为1.10（1.07, 1.12），每重叠一个片段重复的P值为1.46 × 10−17（补充表3），这表明这些低覆盖率区域并非英国生物银行全基因组测序数据集所独有。
经过这些过滤步骤后，我们保留了1,026种蛋白质以进行进一步分析。

We identified 1,651 high-quality rare pQTLs for 599 proteins

我们鉴定了1,651个高质量的罕见蛋白质数量性状位点（pQTLs），涉及599种蛋白质

Para_01

在这些过滤之后，我们得到了5,076个顺式-pQTLs，其中包括98个结构变异。
对于给定蛋白质的所有独立顺式-pQTLs共同解释的平均方差（样本内）为4.10%（中位数为2.02%），与之前报道的估计值相似。
在确定的5,076个变异关联中，762个（15.0%）属于罕见频率范围（0.1% ≤ MAF < 1%），925个（18.2%）非常罕见（MAF ≤ 0.1%），其中包括36个结构变异。
我们将两个频率区间中的1,651个单核苷酸或短插入/缺失变异称为罕见pQTLs。

Para_02

我们还分别对1,184名和1,027名遗传学推断为南亚裔和非洲裔的个体进行了单变异测试，这些个体同时具有全基因组测序（WGS）和Olink蛋白质组学数据。
在考虑的三种遗传学推断的祖先来源中（欧洲（EUR）、南亚（SAS）和非洲（AFR）），我们在EUR和SAS个体之间观察到pQTLs效应量之间的强相关性（r = 0.797，P < 1 × 10−300），而在EUR和AFR个体之间以及AFR和SAS个体之间的相关性较弱（分别为0.656，P < 1.00 × 10−300和0.527，P = 6.79 × 10−250）。
尽管SAS和AFR分析的样本量小得多，我们仍分别识别出228个和396个独立的pQTLs。
其中，有77个（33.8%）和70个（17.7%）也在我们对EUR个体的分析中是主要的pQTLs。

Para_03

我们将我们的单变异 pQTL 结果与 Eldjarn 等人进行了比较，他们分析了从 150,119 例英国生物样本库全基因组测序数据中推算出的基因组数据，这些数据来自 54,306 名具有蛋白质组学数据的个体。
我们发现他们的 3,386 个顺式 pQTL（在我们测试的区域内；占比 76.1%）中，有 2,575 个直接映射到英国生物样本库 DRAGEN 全基因组测序调用结果，并且与我们针对相同循环蛋白的至少一个信号呈强连锁不平衡（LD）（r² ≥ 0.8）。
当仅考虑罕见变异（MAF < 1%）时，重叠比例更大（904 个中有 848 个；占比 93.8%）。
平均而言，每种蛋白质我们鉴定出 5.0 个（中位数为 4）独立的顺式 pQTL（包括罕见和常见变异），而 Eldjarn 等人鉴定出 10.1 个（中位数为 8）。

Most coding pQTLs reduce circulating cognate protein levels

大多数编码蛋白的数量性状位点（pQTLs）降低了循环对应蛋白的水平

Para_01

正如预期的那样，并且与最近的外显子组测序研究一致，96.3%被注释为高可信度功能丧失的 pQTLs（由 LOFTEE 定义）与循环蛋白水平降低相关联。（补充信息，补充表 7 和扩展数据图 2）。

Para_02

基于聚合的罕见编码变异测试确定了545个与539种循环蛋白水平相关的基因，这些结果在调整了我们单变异分析中常见的（>0.1%）pQTLs之后得出（补充表8）。
其中，有14个基因未被Dhindsa等人进行的外显子组测序分析发现（补充表9），这可能是由于测序覆盖范围的差异，从而展示了全基因组测序（WGS）发现外显子测序无法捕捉到的编码变异的能力。

Noncoding pQTLs are enriched in 5ʹ-UTRs

非编码pQTL在5ʹ-UTRs中富集

Para_01

我们发现了604个与369种蛋白质相关的独立罕见非编码单变异（图2和图3以及补充表10），与罕见编码变异（985个高置信度的失活或错义变异）进行了比较。
我们根据注释和最严重的后果将pQTLs分为同源组和非同源组，优先考虑相对于同源基因的注释。
在604个非编码pQTLs中，343个（56.8%）被注释为对同源基因具有调控作用（即位于UTRs或内含子5千碱基范围内；图3d）。
相比之下，972个（98.7%）罕见编码变异被注释为与同源基因相关。

Para_02

先前比较两种不同平台的研究表明，cis-pQTLs 的高比例可能是由表位结合伪影引起的。
尽管我们无法获得 SomaLogic 数据以与 Olink 进行直接比较，但我们注意到，我们的 604 个罕见非编码变异中没有任何一个与编码变异至少处于‘低’连锁不平衡 (r2 ≥ 0.1) 状态。
我们还发现，罕见非编码变异的影响大小与编码 pQTL 的最大 r2 值之间没有关联的证据（P = 0.801），即使在调整变异频率后仍然如此（P = 0.763）。
这表明，表位效应对我们所发现的非编码 pQTL 关联影响甚微。

Para_03

非编码变异的平均绝对效应为1.15个标准差（中位数为0.86个标准差），分别相当于罕见功能缺失和错义pQTLs平均绝对效应的57.1%和74.1%。
此外，罕见非编码pQTLs在增加循环蛋白效应和减少循环蛋白效应之间的分布更加均匀（65.2%为减少效应），而罕见编码pQTLs则有86.3%表现为减少效应（异质性P值=5.41×10^-13）。
在89.8%的情况下，关联最强的罕见pQTL位于或靠近其对应基因。
单独考虑时，罕见编码pQTLs和非编码pQTLs分别有97.9%和66.9%位于或靠近其对应基因。
罕见非编码pQTLs分布在顺式位点上，距离对应基因的最大距离为993 kb，接近所测试顺式区域的边界。

Para_04

接着，我们测试了在不同注释类别中，标记为调控其同源基因的稀有变异的富集情况。
基于每个测试变异最严重的预测后果（与同源基因相关的后果被优先考虑），我们在以下类别中观察到 pQTLs 的富集：5ʹ-非翻译区 (优势比 OR = 23.6, Fisher 精确检验双侧 P 值 = 4.46 × 10−48)，3ʹ-非翻译区 (OR = 3.08, P = 4.34 × 10−7)，预测的内含子剪接位点 (OR = 208.5, P = 1.08 × 10−29)，非编码外显子 (OR = 2.69, P = 9.02 × 10−5) 和上游变异 (OR = 3.00, P = 3.80 × 10−20; 图 4)。
当我们考虑标记为调控顺式窗口中另一基因的稀有非编码 pQTLs 时，没有观察到富集的证据 (扩展数据图 4)。

Fig. 4: Distribution of annotations of lead, rare, pQTL, noncoding variants compared with all variants tested.

- 图片说明

Para_05

基于确定的1,651个稀有单变异，我们在P < 0.05/(1,651 × 2,907)的条件下，识别出与2,907种蛋白质之一相关的1,040个额外独特的远端pQTL关联（补充表11）。
需要注意的是，我们包含了之前因顺式全基因组测序数据质量而排除的蛋白质的远端效应，因为排除依据是顺式-WGS数据的质量，而非表型质量。

Aggregate testing identified 357 regulatory regions

聚合测试确定了357个调控区域

Para_01

通过基于聚合的关联性测试，我们发现了357个与循环蛋白水平相关的条件独立罕见变异非编码区域（补充表12和13）。
如果聚合单元被注释到对应的基因本身，或者聚合的起始位置位于对应基因单元内，我们将每个独立关联的聚合关联标记为映射到该对应基因。
所有其他聚合被标记为非对应基因，尽管其中许多可能比下一个最近的基因更接近对应基因。
与所有测试的对应基因聚合单元相比，我们观察到注释为5ʹ-非翻译区（5ʹ-UTR）的对应基因聚合体富集的证据（图5，优势比OR = 6.36，P值 = 5.11 × 10⁻¹¹），以及预测为内含子剪接位点的变异体聚合体的富集（OR = 29.6，P值 = 4.89 × 10⁻¹⁷）。

Fig. 5: Distribution of annotations of significant cognate aggregates of rare noncoding compared with all aggregates tested.

- 图片说明

Para_02

超过90%的非编码聚合信号是在不限制变异达到我们规定的保守性或约束阈值的情况下检测到的。
在357个（占5.04%）条件独立的非编码聚合关联中，有18个仅在选择高度保守（GERP > 2）的变异时被识别，而16个（占4.48%）仅在选择高度受限（JARVIS > 0.99）的变异时被识别。
大多数聚合区域如果仅通过单变异分析将会被遗漏：259个（占72.6%）非编码聚合中没有主要的pQTL，74个（占20.7%）不含任何全基因组显著性的单变异。
基于357个独立的非编码聚合，我们还确定了45个独特的跨非编码聚合-蛋白关联（P < 0.05/(357 × 2907)）

Para_03

绝大多数罕见的非编码聚集体是通过允许双向关联的统计测试识别出来的，并且假设很大一部分变异是非因果的。
在357个罕见非编码聚集体pQTLs中，仅有10个（2.80%）在负担框架下最显著相关（假设所有罕见变异都会产生相同方向的影响），这与基于编码区的聚集体形成强烈对比，后者的比例为38.6%。
这一差异表明，非编码区域中的罕见变异很可能导致性状增加和减少的混合效应，并且并非聚集体中的所有变异都是因果性的。
这一观察结果表明，改进对非编码变异的注释将进一步提高检测非编码聚集信号的机会。

Para_04

我们发现滑动窗口聚合检验能够识别出那些仅通过限定在注释区域内的检验无法检测到的关联。
例如，有34种（3.31%）循环蛋白具有滑动窗口pQTL。
对于这34种蛋白，当包含滑动窗口时发现了109个非编码pQTL，而不包含滑动窗口时则只有81个。
对于三种蛋白，唯一的pQTL是滑动窗口。
例如，在ACAT-V框架下，LAMP3基因内含子区域chr3:183135000-183137000与循环LAMP3水平相关联（P = 1.15 × 10−13），但该区域不包含任何Ensembl调控区域。
我们的结果表明，相对于仅测试先前已映射的区域，滑动窗口具有额外的价值，尽管我们未包含的细胞和组织特异性数据中的注释可能可以检测到这些区域。

Rare noncoding pQTLs showed tissue-specific enrichment

罕见的非编码 pQTL 显示出组织特异性的富集

Para_01

接下来，我们旨在确定罕见的非编码pQTL在组织相关的非编码调控区域中的存在程度。
首先，我们测试了这样一个假设：如果相关蛋白是分泌蛋白或信号蛋白，则更有可能识别出罕见的非编码pQTL，因为循环蛋白相比非分泌蛋白更能代表整体蛋白丰度。
其次，我们基于以下假设测试了另一个假设：在血液和肝细胞的Ensembl调控元件中，罕见的非编码pQTL会比其他20种组织类型更丰富（方法），因为这些是最相关的组织类型。
我们发现，在所有调控区域中，单变异pQTL的关联在所有蛋白组中高度富集，无论其分泌状态如何（图6、扩展数据图6和7以及补充表15-20）。
在增强子和CTCF结合位点内，分泌蛋白或信号蛋白的关联比非分泌蛋白更高度富集。
在非编码聚合关联中，并且仅限于通过滑动窗口识别出的关联，分泌蛋白和信号蛋白在所有调控区域中的富集程度高于非分泌蛋白。
在滑动窗口内，关联在预测活跃于血管的调控区域中最高度富集，其次是肝脏（补充表20）。

Fig. 6: QQ plot for enrichment of loci within Ensembl-predicted active regions within tissue groups.

- 图片说明

Para_02

作为敏感性分析，我们尝试缓解这样一个担忧，即我们的发现可能是由用于测量循环蛋白质的平台的技术人工产物所驱动的。
在551种被确定为Olink（Explore 3072）和SomaScan-v.4之间高度一致的蛋白质中，有261种具有高质量的顺式全基因组测序覆盖。
我们的pQTL关联在这261种蛋白质中富集：尽管它们仅占我们检测的蛋白质的25.4%，但却涵盖了我们所有pQTL的31.1%（两比例检验，P = 2.85 × 10−4），稀有pQTL的32.7%（P = 6.05 × 10−5）以及基于非编码聚合的pQTL的37.5%（P = 1.29 × 10−5）。

Discussion

Para_01

以循环蛋白质水平为例，我们已经证明了全基因组测序（WGS）数据分析能够发现与常见表型相关的罕见非编码变异及其聚合体。
通过全基因组测序数据，利用罕见变异聚合测试，我们可以考虑的变异数量比单独进行单个变异测试时多出一倍以上。

Para_02

然而，我们的结果还表明，蛋白质的‘统计独立’pQTL数量与顺式测序覆盖率呈负相关。
与以往研究不同，我们因此尝试考虑那些未完全捕捉到的区域：仅约三分之一的被测量蛋白质在其位点上具有我们认为高质量的全基因组测序数据。
尽管我们在排除的区域中观察到之前报道为问题区域的富集现象，这表明这些问题并非英国生物样本库（UKB）独有，但要充分理解每个顺式区域中观察到的低覆盖率仍需进一步分析。
这是对其他将使用此规模测序数据的研究者的一个关键警示点。

Para_03

我们已经确定了数百个新的非编码罕见聚合和单变异体关联，这些关联与认知基因两侧 1 Mb 范围内的顺式窗口中测量的蛋白质水平相关。
我们发现非编码关联的效果大小有时与编码关联的大小相似，但在促进蛋白质增加和减少之间的平衡性更强。
由于连锁不平衡 (LD) 和单倍型效应的复杂性质，要确定是编码信号驱动非编码信号还是反之非常困难。
为避免这种解释，我们采取了保守的方法，并对认知基因的所有编码变异进行了条件控制。

Para_04

我们在单变异 pQTLs 数量方面观察到与 Eldjarn 等人的研究结果存在一些差异。
这些差异可能部分由条件分析方法驱动：他们的分析使用了逐步向前的条件分析来定义条件独立的 pQTLs，与此相关的问题之前已被详细讨论过（例如，参考文献33），而我们则在 GCTA-CoJo 中实施了向前和向后的条件分析步骤。
这两项研究中 pQTLs 的差异突显了解释同一位点上的多个独立 pQTLs 的困难。

Para_05

我们发现5ʹ-非翻译区和预测的内含子剪接受体或供体位点富含罕见的非编码pQTLs。
由于UTRs和内含子通常不会在靶向编码外显子的外显子组测序中被捕获，我们的结果突显了全基因组测序在寻找新的罕见基因相关变异中的重要性。

Para_06

我们还展示了在非编码区域进行聚合测试的力量，这类似于已在编码区域中功能相似变异聚合的成熟方法。
通过测试稀有的基因组非编码元素聚合体，这些聚合体按照与基因的距离、预测的调控活性或使用滑动窗口等方式分组，我们发现了另外464个条件独立的区域，而这些区域无法仅通过单变异测试识别。

Para_07

与基于聚合的编码关联相比，非编码基因组聚合 pQTLs 在允许双向效应和/或稀疏因果关系的测试中更为丰富。
这一观察结果与以下事实一致：在非编码区域中，变异效应和功能区域的预测精度较低。
然而，我们利用当前注释和数据识别出非编码关联这一事实表明，随着功能注释的改进和群体遗传数据的积累，常见表型中可能会有更多发现。

Para_08

我们的研究存在一些局限性。
首先，由于我们无法获取其他研究的类似数据，因此无法在独立数据集中重复我们的结果。
然而，我们的大部分关联达到了远超阈值的统计置信水平。
此外，尽管我们的主要分析因样本量限制而局限于欧洲血统个体，但效应大小在不同血统群体中是一致的。
其次，我们不能完全确定已经考虑到所有可能由与编码变异或常见变异的连锁不平衡（LD）导致的残余混杂因素，包括复杂或重复区域中的变异。
出于类似原因，我们在同一区域内识别出的部分蛋白质数量性状位点（pQTLs）可能并非真正独立。
然而，由于这些关联具有不同的特征，包括对蛋白质增加和减少效应的更均匀分布，因此它们不太可能受到编码变异残余混杂的显著影响。
我们也没有观察到独立非编码 pQTLs 的效应大小与编码变异的残余 LD 之间存在任何关系，这表明我们的结果不太可能受到表位效应的干扰。
第三，所有循环蛋白均在血液中测量。
虽然大部分组织特异性蛋白仅在特定组织中表达，但我们受限于只能考虑循环蛋白水平。
第四，我们无法考虑与测量蛋白水平所使用的技术相关的结合效应。
然而，在 Olink 和 SomaLogic 平台之间显示出强相关性的 551 种黄金标准蛋白的子集中，我们确实观察到了关联的富集现象。
最后，尽管我们指出滑动窗口聚合识别了潜在的新调控区域，但不可能全面检查所有公开可用的调控图谱。
此外，测试所有可能的组织特异性注释会带来巨大的假设检验负担。

Para_09

总之，以循环蛋白质测量的示例特征为例，我们发现了循环蛋白质水平与罕见非编码变异之间的多种新关联。
我们的结果表明，可能有许多具有重大影响的罕见非编码变异等待被发现，这些变异会影响复杂的表型。

Methods

Para_01

本研究遵守所有适用的伦理规定。
英国生物银行（UKB）研究已获得西北研究中心伦理委员会的伦理批准。

UKB and WGS

英国生物银行（UKB）和全基因组测序（WGS）

Para_01

用于 UKB 的全基因组测序（WGS）通过 Illumina NovaSeq 6000 测序仪完成，平均覆盖度为 32.5×，采用 150-bp 双端测序读长。
测序所使用的基因组构建为 GRCh38，单核苷酸多态性和短插入缺失变异是通过 DRAGEN 3.7.8 联合分析得出的。

Human protein expression levels

人类蛋白质表达水平

Para_01

使用 Olink 技术，对多达 54,219 名 UKB 参与者的 2,932 种蛋白质的水平进行了分析，具体方法参见参考文献 17，由 UKB Pharma 蛋白组学项目执行。
在数据提供给研究人员使用之前，进行了质量控制程序，包括去除异常值。
在数据发布前，蛋白质水平还经过了 log2 转换。
经过质量控制过滤后，共有 54,189 名具有蛋白质表达数据的个体被批准用于分析。
Sun 等人未发现批次或平板混杂效应的证据。

Genetic data filtering

基因数据过滤

Para_01

如果每个样本基因型的 sum(LAD) < 8（其中 LAD 是局部等位基因深度）或基因型质量 (GQ) < 10，我们将任何 DRAGEN WGS 基因型调用设置为缺失，这适用于 UKB 使用 bcftools v.1.2 提供的 154,430 个项目变异调用格式 (pVCFs)。
在这些额外的质量控制步骤之后，单例的传递率（理论上应精确为 0.5，假设大多数变异未受到强烈负向选择）为 0.497，而 UKB 最初提供的数据为 0.456。
随后，我们删除了缺失率 >10% 的任何变异。
应用多等位分裂程序后，每个变异被分配一个唯一 ID (CHR:BP:REF:ALT)，然后基于（访问日期：2024 年 3 月 30 日）提供的 1000G b38 参考序列，使用 bcftools 对所有每条染色体上的 VCF 文件进行合并，并对插入/缺失进行了规范化和左对齐处理。
每个合并后的 pVCF 随后被转换为 plink（v.2.0）p(gen/var/sam) 格式。

Structural variants

结构变异

Para_01

在分析时，DRAGEN版本的结构变异（变异≥50 bp）仅以单样本格式提供。
我们使用 SURVIVOR v.1.07 与默认设置，将结构变异调用合并并 harmonize 成包含 54,219 名具有循环蛋白测量值个体的 pVCF 文件。
然后，我们将该项目 (p)VCF 转换为 plink p(gen/var/sam) 格式，并与 DRAGEN SNV 或 indel 调用的 plink 文件合并。

Genetic variant annotation

基因变异注释

Para_01

我们使用 Ensembl VEP（v.110）、LOFTEE 和 UTRannotator 对所有遗传变异进行了注释。
在可能的情况下，我们将每个变异分配到以下三个类别之一：编码区、近端调控区或基因间调控区。
如果一个变异对任何转录本的编码序列有预测影响，则被归类为编码区变异；如果变异位于转录本非翻译区（UTRs）5千碱基对范围内的区域，并且在任何转录本中不已经是编码区变异，则归类为近端调控区变异；最后，如果它既不是编码区变异，又映射到一个与基因无关的调控元件（详细信息见下文），则归类为基因间或内含子调控区变异。
此外，我们还对大小为 2,000 碱基对的滑动窗口中的变异进行了测试，无论每个窗口中有多少个变异，但排除了编码区变异以尽量减少假设重叠。

Para_02

然后，我们根据预测的后果和位置，将每个变异分配到我们称为掩码的分组中。
我们使用了五种已发表的变异评分来按后果对变异进行分组：
（1）GERP：GERP 分数是衡量变异水平保守性的指标。
如果一个变异的 GERP 分数 > 2，则我们认为它是高度保守的。
（2）PhastCon 分数：PhastCon 是基于窗口的物种间保守性度量方法，可以仅针对哺乳动物（phastCon 30）或所有物种（phast_100）。
我们测试了非编码基因组窗口，即排除任何包含外显子的窗口，并选择了 PhastCon 分数处于第 99 百分位的窗口。
（3）约束分数：约束是根据每个窗口的局部可变性和观察到的突变率在 1 kbp 窗口大小内计算得出的。
我们测试了约束 z 分数 ≥ 4 的窗口。
（4）SpliceAI 分数：SpliceAI 分数是一种衡量预信使 RNA 区域内每个变异是否为剪接供体、受体或两者都不是的可能性的指标。
如果一个变异的 AI 分数 > 50，则将其分类为高置信度的剪接位点。
（5）CADD：CADD 分数用于预测一个变异可能有多大的破坏性。
我们仅将 CADD 分数应用于编码变异，并且只有被 VEP 标记为高置信度的失活功能变异才会被考虑。
错义变异中 CADD 分数 > 25 的会被分离出来在单独的掩码中进行测试。
（6）JARVIS 分数：JARVIS 分数是基于约束测量的机器学习模型推导出来的，旨在更好地优先考虑关联研究中的非编码遗传变异。

Para_03

每个基因组掩码由多个变异组成，这些变异具有不同的后果，基于它们的位置、上述得分之一和/或预测的编码后果。
例如，要将一个变异分类为错义 CADD > 25，它必须改变基因转录本外显子的密码子，并被预测为高度有害。
在补充表 22 中，我们列出了分配给每个掩码和分类的完整后果列表。

Association analyses

关联分析

Para_01

我们对UKB中测量的2,932种蛋白质的顺式位点内的单变异和聚合变异进行了测试。
为了定义顺式窗口，我们首先将每种蛋白质映射到一个编码基因（少数排除情况见补充表1），并确定了Ensembl记录的每个基因的最长转录本。
基于最长转录本，我们将顺式窗口定义为转录本基因5ʹ-和3ʹ-非翻译区（UTRs）两侧各1 Mb的范围（受染色体起始和终止位置限制），以及编码序列和内含子序列中的变异。
所有关联分析均校正了年龄、性别、年龄平方、UKB招募中心（作为地理的代理）、前40个遗传主成分、全基因组测序批次、Olink板、禁食时间和采血后时间。

Single-variant association testing

单变量关联测试

Para_01

为了识别与蛋白质水平相关的顺式单个变异，我们首先使用 regenie v.3.3（参考文献41）在顺式窗口中对所有遗传变异进行了关联测试，这些变异的最小等位基因计数（MAC）至少为五。
然后，我们使用修改版的 GCTA-CoJo19（diff-freq = 0.2, cojo P = 2.95 × 10−10）在条件联合分析中选择了主要变异，以英国生物银行全基因组测序数据为连锁不平衡（LD）参考面板，并限制为具有蛋白质组学数据的个体。
测试显示，如果某个变异所解释的方差大于任何独立变异所能解释的最小方差的 900 倍，则 GCTA-CoJo 会过滤掉该变异（‘sqrt(ldlt_B.vectorD().maxCoeff()/ldlt_B.vectorD().minCoeff()) > 30’：gcta/meta/joint_meta.cpp 的第 732 行，网址为，访问日期为2024年3月17日）。
我们理解此过滤器的存在是为了捕捉由共线性引起的统计混淆，例如，如果用于计算 LD 的参考基因组和遗传数据之间的相关性较差时就会出现这种情况。
然而，在我们共同考虑常见和罕见变异的目的下，由于我们使用的 LD 参考面板与我们的发现数据集完全匹配，我们发现此过滤器错误地去除了效应较大的 pQTLs。
因此，我们移除了该过滤器并重新编译了 GCTA-CoJo，其代码可在获取。

Para_02

GCTA-CoJo 在超过某个阈值后假设 P = 0。
在这种情况下，我们根据以下 R 脚本重新计算了 log10(P)：log10(P-计算值) = log10(exp(1)) × (pt(−abs(β/标准误差), 自由度 = N-2, log = TRUE) + log(2))。

UKB WGS coverage calculation and filtering

英国生物银行全基因组测序覆盖度计算与过滤

Para_01

我们根据质量控制标准，从 UKB 全基因组测序数据中计算了覆盖率，涉及998名推断为欧洲血统的个体，这些个体同时具有蛋白质组学数据和全基因组测序数据。
我们首先计算了每个碱基对上深度大于8的样本数量。
随后，我们计算了每个同源蛋白编码基因在5ʹ-非翻译区和3ʹ-非翻译区之间碱基数的比例，要求超过90%的个体在这些区域的深度大于8。
我们对完整的顺式窗口重复了这一分析。
然后，一个蛋白的基因水平覆盖率和区域水平覆盖率均需满足：超过99.5%的碱基和超过90%的个体深度大于8。

Para_02

为了考虑可能遗漏或错误基因分型的单一有害编码位点，如果任何与该蛋白质相关的 pQTL 的边际效应估计值（来自我们的 regenie 分析）和联合效应估计值（GCTA-CoJo）之间的最大差异超过第 90 个百分位数，我们还会额外过滤这些蛋白质。
扩展数据图 8 提供了描述完整过滤流程的流程图。

Para_03

然后，我们测试了我们的结果与因低覆盖率而被排除的基因组区域之间的重叠，这些区域来自一个有问题的基因组区域数据库，该数据库由 Genome in a Bottle 和 precisionFDA Truth Challenge v.2 创建（版本 v.3.5，2024 年 9 月 10 日访问）。
这些区域包括"串联重复序列、所有大于 6 个碱基对的均聚物、所有大于 10 个碱基对的不完美均聚物、所有难以映射的区域、所有片段重复、GC 含量小于 25% 或大于 65% 的区域、‘问题启动子’、chrX/Y XTR 和扩增子、卫星序列以及‘其他困难’区域（仅针对 GRCh38 包括 T2T-consortium 提供的区域）"。
对于每种蛋白质，我们确定了至少有 1 个碱基对与基因区域（从 5ʹ- 到 3ʹ-UTR）或 2-Mbp 窗口重叠的每种类型的问题区域的数量。
这些变量在逻辑回归模型中进行了测试，以二元变量为依据描述该蛋白质是否因低覆盖率而被过滤，并调整了 2-Mbp 窗口的总长度（该长度自然随基因大小而变化）。

Rare variant genomic aggregate testing

罕见变异基因组聚合检测

Para_01

为了识别单变量分析能力不足的基因组非编码潜在调控区域，我们随后进行了非编码罕见变异（MAF < 0.1%）的基因组聚合关联测试。
为了检验非编码罕见变异聚合信号是否由与常见变异和/或单变量信号的残余连锁不平衡（LD）及单倍型结构引起或混淆，我们对每个达到 Bonferroni 校正 P < 0.05 的罕见变异聚合测试结果进行了以下步骤：
（1）为了生成我们的主要非编码发现结果，我们调整了在联合（CoJo）分析中被确定为独立信号的常见主效变异（MAF > 0.1%，~MAC = 40），并调整了所有被注释为该基因任何转录本中编码序列的遗传变异（无论其 P 值如何），这些基因映射到感兴趣的蛋白质上；
（2）为了识别独立的非编码聚合关联，如果至少一个聚合通过了我们的显著性阈值，则我们进行逐步回归分析。
对于每种蛋白质，从最显著相关的非编码聚合（按 P 值）开始，我们对达到全基因组显著性的聚合进行了额外的非编码聚合测试，并调整了顶级信号中的所有变异。
此过程重复进行，每次加入更多来自下一个最强相关聚合的变异，直到没有聚合达到全基因组显著性为止。
（3）为了确定我们的主要聚合发现结果是否可能是由单一低频主效变异引起的，我们识别了包含正好一个主效遗传变异的聚合关联。
（4）作为敏感性分析步骤，为了评估这些结果可能在多大程度上受到连锁不平衡（LD）的混淆影响，我们进一步调整了所有已识别的 pQTL 单变量。

Para_02

对于最大等位基因频率阈值为 0.1% 的变异，使用 regenie 进行了基因单元测试，测试基于补充表 22 中指定的遗传单元。
Regenie 执行四种类型的基因组单元测试：(1) 标准 BURDEN 测试，假设给定基因单元掩码中的每个变异对表型的影响大小和符号大致相同；
(2) SKAT 测试，允许单元中每个变异的关联符号有所不同；
(3) ACAT 测试，其中单元中每个变异的关联符号可以不同，并且只需要掩码中的少量变异具有关联性；
(4) ACAT-O，这是一种综合测试，结合了 BURDEN、SKAT 和 ACAT 测试，旨在最大化这三种测试的统计功效。

Para_03

我们对每个至少有一个变异的基因单元所对应的每个掩码进行了上述四种统计检验中的每一种。
此外，还对每个单元中的所有单例变异（MAC=1）进行了关联性测试；regenie 还为每个基因组单元估算了一个‘全掩码’关联强度，这是个体掩码测试统计量的汇总结果。
为了确保这不会导致非编码和编码关联统计量的混杂，我们将每个基因转录本拆分为一个编码转录本，用于测试所有编码掩码，以及一个近端转录本，用于测试所有近端掩码。
调控基因组单元根据其 ENSR 分配、1kb 约束窗口或 phastCon 保守窗口的程度进行分类。
我们根据它们覆盖的相应染色体区域命名滑动窗口掩码。
执行主要分析所需的代码可在 Zenodo 上获取。

Statistical significance

统计显著性

Para_01

统计学显著性是根据对20个随机生成的正态分布连续特征进行全基因组测序（WGS）分析所观察到的最小P值来定义的。
单变异体和聚合关联分析的最小P值被视为独立：单变异体的P值为2.95×10^-10；聚合分析的P值为8.71×10^-9。

Ensembl regulatory region enrichment

Ensembl调控区域富集

Para_01

我们计算了单个变异和聚合区域与Ensemble31调控区域的重叠富集程度，这些调控区域涵盖了118种组织或细胞系，并与遗传背景进行了比较。
对于每种组织，Ensembl还提供了每个区域是否处于活跃状态（或不活跃、受抑制等）以及调控活动类型（启动子、增强子、CTCF结合位点、转录因子结合位点或开放染色质区域）的预测。
随后，我们仅考虑被预测为活跃的区域，排除了细胞系和癌症衍生的组织，并将剩余的组织类型分组为22个超组。
我们使用DAVID软件将蛋白质分为分泌型、信号型或膜结合型蛋白质。
然后，我们将蛋白质归类为分泌型或信号型蛋白质，如果它们被注释为属于任一组，并将所有其他蛋白质标记为非分泌型。

Para_02

为了确定统计上的富集程度，我们进行了10,000次模拟的引导抽样。
对于每次模拟，我们随机选择与我们从基因组顺式调控区域中确定的独立信号数量相等的稀有非编码变异或聚合体，这些信号由我们的稀有非编码单变异或聚合体测试集合确定。
然后，我们确定随机选择的变异或聚合体集合与任何调控区域的重叠情况（我们对每个面板和组织类型的分层重新执行此操作），并将任何模拟中的重叠数量分布与我们在独立关联中观察到的重叠数量进行比较。
随后，我们为观察到的重叠分配了一个经验P值。

Statistics and reproducibility

统计学与可重复性

Para_01

本研究纳入了所有在 UKB 中同时测量了蛋白质水平和全基因组序列的个体。
参与者在研究招募时年龄介于37至80岁之间。
对遗传数据进行主成分分析（PC分析），以根据遗传祖先定义同质群体。
参与者按照遗传特征被分为欧洲样、南亚样和非洲样的遗传祖先群体。
在每个群体内分别进行分析，并在混合模型中调整主成分，以考虑隐性相关性。
本研究为观察性研究，因此不适用随机化或绑定。

Reporting summary

报告摘要

Para_01

有关研究设计的更多信息可在此文章链接的《自然组合报告摘要》中获取。

Data availability

Para_01

由于 UKB 的数据可用性和数据返还政策，数据不能公开共享。
符合访问其数据集条件的研究人员可以从 UKB 获得数据（）。
汇总统计数据可通过 Zenodo 在 .5281/zenodo.14203628 上获取。

Code availability

Para_01

与分析相关的代码可通过 GitHub 获取（）。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-03，如有侵权请联系 cloudcommunity@tencent 删除数据统计变量编码测试

科技改变生活-雨落星辰 - 所有的伟大,都源于一个勇敢的开始