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

正文

请到「今天看啥」查看全文

- 图片说明

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

Para_02

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

Para_03

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

Para_04

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

Rare noncoding pQTLs showed tissue-specific enrichment

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

Para_01

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

1. 作为敏感性分析，我们尝试缓解这样一个担忧，即我们的发现可能是由用于测量循环蛋白质的平台的技术人工产物所驱动的。
2. 在551种被确定为Olink（Explore 3072）和SomaScan-v.4之间高度一致的蛋白质中，有261种具有高质量的顺式全基因组测序覆盖。
3. 我们的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

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

Para_02

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

Para_03

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

Para_04

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

Para_05

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

Para_06

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

Para_07

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

Para_08

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

Para_09

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

Methods

Para_01

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

UKB and WGS

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

Para_01

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

Human protein expression levels

人类蛋白质表达水平

Para_01

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

Genetic data filtering

基因数据过滤

Para_01

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

Structural variants

结构变异

Para_01

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

Genetic variant annotation

基因变异注释

Para_01

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

Para_02

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

Para_03

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

Association analyses

关联分析

Para_01

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

Single-variant association testing

单变量关联测试

Para_01

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

Para_02

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

UKB WGS coverage calculation and filtering

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

Para_01

请到「今天看啥」查看全文