把 GWAS 的信号接到剪接上：精神分裂症特异的 sQTL 图谱定位新的风险基因
下划线 = 中文容易“变味”或丢信息的词句：红线带编号 = 15 条注释卡（悬停看、点击跳转）。这一页没有一只老鼠：全部是人的基因型、脑组织的 RNA，和几套统计算法。看它的正确方式不是“药物有没有效”，而是“一个碱基凭什么被判定为致病的那一个”。英文原文附于文末
摘要（全译）
可变剪接（alternative splicing，AS）在精神分裂症（schizophrenia，SCZ）的发病中起着重要作用。既往研究把全基因组关联研究（GWAS）得到的遗传信号与表达数量性状位点（eQTL）联系了起来，但它与其它遗传调控机制——尤其是剪接 QTL（sQTL）——之间的相互作用，仍不清楚。本研究构建了一份较完整的疾病特异 sQTL 图谱，用来提供那些在 SCZ 中通过 RNA 剪接改变基因活性的遗传变异。
我们分析了 539 例 SCZ 患者的数据，共鉴定出 24,810 个显著 sQTL（FDR < 0.05），涉及 7083 个特有基因的 AS 事件。把这份图谱与一项大规模 SCZ GWAS 结合，我们用孟德尔随机化（MR）和共定位分析，锁定了 27 个显著的风险基因——它们受遗传性的 AS 调控，可能在 SCZ 中有因果作用。对这批基因，再在 539 例病例和 754 例对照中做差异剪接分析，又发现 12 个显著基因，它们可能因 AS 失调而升高 SCZ 风险。
值得注意的是，有五个基因（DPYD、LACC1、CCDC122、ANAPC7、DGKZ）在 MR 分析和差异剪接分析中表现出方向一致的剪接调控效应。对差异剪接基因做通路富集分析，提示了与 SCZ 相关的潜在生物学通路，尤其是突触传递和微管运动。此外，单细胞 RNA-seq 分析显示，若干基因优先表达于特定的脑细胞类型，包括少突胶质细胞、小胶质细胞和兴奋性神经元。
总的结论是：本研究突出了若干可能通过 AS 调控参与 SCZ 风险的易感基因；对这些基因的进一步刻画，有望推进对 SCZ 机制的理解和治疗的发现。
引言（精译）
精神分裂症（SCZ）是一种严重而复杂的精神障碍，以认知和思维的异常为特征，全球终生患病率约 1%。由于高发病率与高死亡率，SCZ 给个人、家庭和社会带来巨大的经济和医疗负担。但 SCZ 的发病机制至今仍很不清楚，现有治疗的获益也有限。因此，迫切需要找到有效的、特异的生物标志物，以开发针对 SCZ 的治疗策略。
SCZ 的遗传度估计至少为 80%，说明遗传因素在发病中占主导。全基因组关联研究（GWAS）的出现，为剖析 SCZ 的遗传病因提供了前所未有的机会。过去十年里，GWAS 已鉴定出数百个与 SCZ 相关的风险位点。然而，把 GWAS 的发现解释成生物学洞见并用于临床，仍是一个巨大的挑战。
pre-mRNA 的可变剪接（AS）是转录后基因调控中的一个关键步骤：它把内含子序列去掉，再把外显子按不同方式连接起来。人类中超过 95% 的多外显子基因会发生 AS，这极大增加了转录组与蛋白质的多样性。最新研究显示，AS 广泛存在于神经与免疫系统中，而异常的 AS 与多种脑疾病相关，尤其是在 SCZ 中。因此，弄清 AS 的调控机制，对理解 SCZ 的发病是必要的。
已有证据表明，AS 的调控可以由遗传变异控制，sQTL 已被广泛用于探索人类疾病中由 AS 调控的遗传变异。但我们对影响 AS 的遗传变异及其在 SCZ 中的致病机制，了解仍然有限。因此，把 sQTL 信息与 GWAS 结果结合起来的整合方法出现了，并在探索“其剪接表达水平受风险变异影响”的潜在风险基因方面显示出前景。
孟德尔随机化（MR）是这类整合方法中的代表：它把与剪接定量表达相关的风险变异当作工具变量（instrumental variables，IVs），来推断某个暴露（即受 RNA 剪接调控的风险基因）对某个结局（即疾病）的因果影响。把 GWAS 的遗传学发现与 sQTL 数据整合，MR 可以推断出那些可能在 SCZ 中有因果作用的风险基因。此外，考虑到 RNA 剪接在脑内高度异质，基于“非目标疾病样本”的 sQTL 数据可能掩盖疾病特异 AS 调控的作用，从而漏掉重要的生物学洞见。因此，用疾病特异 sQTL 数据做 MR 整合分析，将为疾病特异的 AS 调控机制提供新的认识。
据作者所知，此前还没有针对 SCZ 的疾病特异 sQTL 资源。本研究先从两个公开联盟收集了同时具有基因型和转录组的 539 例 SCZ 样本，做了全基因组 sQTL 分析，找出影响 AS 的遗传变异，共鉴定出 24,810 个显著的 sQTL 单核苷酸多态性（SNP）。接着，作者把这些 SCZ 特异的 sQTL 数据与 SCZ 的 GWAS 整合，做了一项较全面的 MR 研究，提出 27 个“其遗传性 AS 调控可能在 SCZ 中有因果作用”的显著基因。再结合共定位与差异剪接分析的证据，最终给出 12 个有希望的风险基因。此外，单细胞转录组分析显示有 13 个基因富集于脑细胞类型，包括少突胶质细胞、抑制性神经元、兴奋性神经元、星形胶质细胞和小胶质细胞。作者认为，本研究提供了一份较全面的 SCZ 特异 sQTL 图谱资源，并给出了一组有强证据支持的、有希望的新药靶点。
先把五个词分清：GWAS / eQTL / sQTL / TWAS / MR
这一节不是翻译，是我给你补的地基。这篇论文的每一个数字，都建立在这几个词的区别上；分不清它们，你读到的就只是一串缩写。
GWAS 问的是“谁更常见”。全基因组关联研究拿几十万人的基因型，逐个位点比较“有病的人”和“没病的人”在这个位置上各个碱基的频率。它不测任何 RNA，也不测蛋白，只看统计关联。输出是一串风险位点（risk locus），每个位点里往往有几十上百个彼此高度相关的 SNP，论文里通常只报一个“索引 SNP”（index SNP）。
eQTL 问的是“产量”。表达数量性状位点需要两样东西：同一个人的基因型，和同一个人的 mRNA 表达量（芯片或 RNA-seq）。它在每个位点上做回归：这个位点的哪个碱基，和某个基因“表达得多还是少”相关？输出是一对一对的“SNP—基因表达”。
sQTL 问的是“版本”。剪接数量性状位点同样需要基因型，但另一个变量不是表达总量，而是剪接比例——同一个基因的几种异构体各占多少。输出是“SNP—剪接比例”。一句话记住：eQTL 管这个基因产出多少，sQTL 管这个基因产出的是哪一版。一个基因的 mRNA 总量可以完全不变，但版本换掉了，蛋白就可能是残废的。
TWAS 是把前两者当梯子用的工具。转录组范围关联研究先用一套参考 panel（比如 GTEx 的 eQTL 权重）训练一个模型：给定基因型，预测这个基因的表达；再用 GWAS 的汇总数据去检验“被预测出来的表达”与疾病是否相关。它输出的是一批基因而不是位点。这篇论文没有用 TWAS，但它出现在几乎每一篇同题材论文的讨论里，你得知道它和 MR 的区别：TWAS 检验的是“预测表达与性状的关联”，MR 检验的是“由遗传变异驱动的暴露对结局的因果效应”，前者更像关联，后者才在假设成立时谈因果。
| 缩写 | 全称 | 一句话在问什么 | 需要的数据 | 输出 | 本文用了吗 |
|---|---|---|---|---|---|
| GWAS | genome-wide association study | 哪些碱基在病人里更常见 | 基因型 + 病例/对照 | 风险位点 | 用（借 PGC 的汇总数据） |
| eQTL | expression QTL | 哪些碱基和 mRNA 产量相关 | 基因型 + 表达量 | SNP—表达 | 用（做组织间对照旁证） |
| sQTL | splicing QTL | 哪些碱基和剪接比例相关 | 基因型 + PSI | SNP—剪接比例 | 主角，自己算的 |
| TWAS | transcriptome-wide association study | 预测出来的表达与疾病相关吗 | 参考 panel 权重 + GWAS | 基因排序 | 没用 |
| MR | Mendelian randomization | 用遗传变异当工具，暴露对结局有因果吗 | sQTL/eQTL + GWAS | 因果效应量（OR） | 用了，是核心 |
再补两个你一定会碰到的词：LD 和共定位。
LD（连锁不平衡）说的是：同一条染色体上邻近的碱基，并不是随机组合的，它们倾向于“整块”一起遗传。所以一个真致病变异旁边那些无辜的邻居，也会在病例里显得更常见，P 值一样漂亮。GWAS 报的“342 个独立显著 SNP”就是把这种整块效应去掉之后剩下的数目；而原始显著信号可能有几万个。记住一件事：LD 让你无法从统计上分辨“这一块里哪一个才是真的”。
共定位（colocalization）就是为这个难题准备的。它同时看两个性状（比如基因剪接比例、疾病风险）在同一区域的关联信号，用贝叶斯方法比较五种情形：H0 两个性状与该区域都无关；H1 只与性状 1 有关；H2 只与性状 2 有关；H3 两个都有关但由不同的因果变异驱动；H4 两个都有关且由同一个因果变异驱动。H4 的后验概率（PPH4）越高，你越有底气说“这个 SNP 既改变剪接、又改变疾病风险，而且不是两个恰好挨着的信号”。本文的门槛是 PPH4 > 0.90。
MR（孟德尔随机化）则是把“基因型”当成一次天然的随机分组：一个人出生时拿到哪个等位基因，近似随机。于是如果你用一个只影响剪接、不影响其它路径的 SNP 当工具变量，就能推断“剪接比例变化 → 疾病风险”的方向。它有三个不能违反的假设：工具变量要真的影响暴露（相关性）、不能与混杂因素相关（独立性）、只能通过暴露影响结局（排他性）。这些假设无法被数据完全验证，只能靠敏感性分析和生物学常识去守。
一句话总结这一节：GWAS 给你“哪里可疑”，eQTL/sQTL 给你“可疑的那个碱基在动什么”，MR 和共定位给你“这个动法有没有可能是原因”。四步走完，才轮得到“这是个新靶点”这句话。
方法精译：从一亿条读数到一个 PSI
数据从哪来。作者用了两个队列的脑组织 SNP 基因型与 RNA-seq 数据：CommonMind Consortium（CMC，N = 328）和 Lieber Institute for Brain Development（LIBD，N = 211），取人脑背外侧前额叶（DLPFC）区域。最终进入分析的是 539 例 SCZ 患者。基因型用 PLINK v1.9 合并，共 27,332,850 个标记。RNA-seq 用 fastp v0.23.4 做质控，用 STAR v2.5.2a 比对到 GRCh37，并用 WASP 校正比对偏倚。
剪接怎么变成数字。作者用 LeafCutter 做无注释的 AS 定量：它不依赖现成的转录本注释，而是直接从 RNA-seq 里找“外显子—外显子”的连接，把共用同一对剪接位点的内含子归成一个簇；簇里每个内含子的使用率，就是 PSI（percent spliced in，剪接比例）。参数为 minclureads 50、mincluratio 0.001、maxintronlen 500000；只保留在 40% 以上样本中出现的内含子；再按 GENCODE v19 的外显子坐标把内含子簇映射到基因。研究者一共得到 282,570 个 AS 事件。
sQTL 怎么算。做cis-sQTL 分析，窗口是内含子簇上下游 1000 kb；用 FastQTL 的线性回归模型检验“PSI 值—SNP 基因型”的关联。为控制遗传、生物学与技术层面的混杂，作者把年龄、性别、RNA 完整性指数（RIN）、群体结构和测序文库批次作为协变量回归掉。显著性用置换检验（1000 次）加 Benjamini–Hochberg 校正，FDR < 0.05。
GWAS 数据。用 PGC 的 SCZ GWAS 汇总数据（Trubetskoy 等），跨族裔共 74,776 例病例、101,023 例对照，报告了 342 个独立显著 SNP、287 个区域。为与 sQTL 的欧洲人样本匹配，MR 只用了其中的欧洲人部分：53,386 例病例、77,258 例对照、7,659,767 个 SNP。
MR 怎么配工具变量。用 TwoSampleMR 0.5.6。第一步 harmonize：确保暴露与结局里的效应等位基因是同一个，否则效应方向会整体反号。第二步选工具变量：做LD clumping（窗口 5000 kb、r2 < 0.01），保证选出的 SNP 彼此独立；再要求次等位基因频率（MAF）> 0.01、F 统计量 > 10（避免弱工具偏倚）。最后每个剪接区只剩一个合格 SNP，因此全部用 Wald ratio 估计；多重检验用 Bonferroni 校正（5179 个独立检验），阈值取 P < 1.93×10−6。
共定位。用 Coloc v5.1.0，输入同一区域内两个性状的关联统计量，计算上述五个假设的后验概率，以 PPH4 > 0.90 作为“共享同一个因果变异”的门槛。
差异剪接、富集、单细胞。在 539 例病例和 754 例对照之间，用 Wilcoxon 秩和检验比较各剪接区的 PSI，取名义 P < 1×10−3。通路富集用 MR 结果里 P 最小的前 500 个基因（因为 Bonferroni 之后没有一个显著），工具是 clusterProfiler 4.10.0，做 GO 与 WikiPathways 富集。单细胞分析用 Mathys 等提供的 24 名认知正常者的前额叶数据，Seurat 5.0.3，剔除表达于 3 个以下细胞的基因与表达少于 200 个基因的细胞，得到 6 类脑细胞，Wilcoxon 检验加 FDR < 0.05。
和 GTEx 对账。为看这套 sQTL 是不是“只是把已知的再算一遍”，作者用 GTEx 的 DLPFC 数据（v10，N = 268）跑了同一套流程，阈值取 P < 2.66×10−6（0.01/3764），并把这套数据也拿去做一遍 MR。
结果一：24,810 个 sQTL 与 7,083 个基因（精译）
为研究 SCZ 中疾病特异的 RNA 剪接遗传调控，作者用 CMC 与 BrainSeq 数据集里 539 例 SCZ 样本的 PSI 矩阵和基因型，做了全基因组 cis-sQTL 分析。经过严格质控，共保留 282,570 个 AS 事件和 27,332,850 个 SNP。最终在 SCZ 样本的 PFC 脑组织中，鉴定出 24,810 个显著 sQTL SNP（FDR < 0.05），涉及 7083 个含 sQTL 的基因（sGenes）。
为看这些 sQTL SNP 在基因组上分布在哪，作者算了每个 sQTL SNP 到最近剪接连接点的距离：与既往研究一致，sQTL SNP 富集在剪接连接点附近。此外，约 38.2% 的 sQTL SNP 位于其对应 AS 事件所在基因的基因体内。
这两个数字值得你记下来。第一，7,083 个基因意味着这套图谱覆盖了相当一部分表达基因——它不是几个明星基因的个案研究，而是一张资源型图谱。第二，38.2% 落在基因体内，说明另外六成左右的 sQTL 落在基因间区或其它位置：这些位点不改变蛋白编码序列，却能改变剪接，这正是后面讨论“非编码变异怎么致病”的落点。
结果二：用 MR 把“关联”翻成“方向”（精译）
为找出“通过影响 RNA 剪接而在 SCZ 风险中有因果贡献”的易感基因，作者把 SCZ GWAS 汇总数据（53,386 例病例、77,258 例对照）与 SCZ 特异 sQTL 数据（N = 539）整合，做了 MR 分析。要特别注意的是：最终分析中只用 Wald ratio 估计，因为经过严格的工具变量筛选（MAF > 0.01、LD pruning r2 < 0.01、F 统计量 > 10）之后，每个剪接区只剩一个合格 SNP。这种情况本来就要求用 Wald ratio——它是只有一个工具变量时的统计最优方法，能在不需要 IVW 荟萃分析所要求的“工具同质性”假设的前提下，给出无偏的因果效应估计。因此，对多重检验做 Bonferroni 校正后，MR 关联以 P < 1.93×10−6 为显著。
作者发现，所有 MR 显著的结果都通过了共定位分析的检验（PPH4 > 0.90）。最终在 31 个内含子使用区中鉴定出 27 个与 SCZ 风险显著相关的基因。其中，chr1:98293752–98386440 内含子使用区对应的 DPYD 关联最强（P = 1.06×10−18）。其它显著的潜在易感基因包括 DGKZ、ANAPC7、FTSJ2、BCL2L12、IRF3、GPM6A、MPHOSPH9、LRRN3 和 IMMP2L。
有一个结果特别有意思：KANSL1 的剪接定量表达（即 PSI 值）被遗传性升高后，在 chr17:44230332–44248221 这个剪接位点上与 SCZ 风险降低相关（OR = 0.97，P = 1.56×10−7），而在 chr17:44172067–44248221 位点上却与风险升高相关（OR = 1.03，P = 1.56×10−7）。作者据此认为，同一个基因的不同内含子使用区在 SCZ 中具有不同的生物学功能。
读到这里请你自己核一下：上面两个 P 值一模一样（1.56×10−7）。两个不同的估计给出完全相同的 P 值，概率极低，多半是原文或表格的复制错误。这类地方不用纠结，但你要在自己的论文里避免。
再把 MR 结果与 PGC3 那项最大的 SCZ GWAS 对比：27 个基因里有 19 个落在已知的 SCZ 易感位点上，包括 DPYD、DGKZ、BCL2L12、IRF3、MPHOSPH9、LRRN3、IMMP2L、FXR1、SUGP1、KANSL1、ANKRD45、TBC1D5、NDUFAF7、PRKD3、FAM114A2、LACC1、CCDC122、AKT3 和 GPM6A。这些受 AS 影响的 MR 结果，有助于在每一个 GWAS 信号里 pinpoint 出潜在的目标基因。更重要的是，有 8 个带 AS 事件的基因不与已知的 SCZ 关联位点重叠：FTSJ2、CSPG4P12、CCDC92、ZNF664、CRELD2、FAM49B、APOBEC3C 和 APOBEC3D。这说明把 SCZ 特异的 sQTL 数据纳进来，能够帮助找出 GWAS 之外的新的目标基因。
结果三：差异剪接、通路富集与单细胞（精译）
差异剪接。剔除不合要求的样本后，539 例 SCZ 病例和 754 例健康对照进入差异剪接分析。在 31 个显著的内含子使用区里，有 10 个区（对应 12 个 sGenes）在病例与对照之间表现出差异剪接定量表达（名义 P < 1×10−3），说明这些重叠的内含子使用区是有希望的功能性遗传位点。
值得注意的是，MR 分析显示：有一个内含子使用区的剪接表达上调与 SCZ 风险升高相关（OR > 1.00），另有三个区的剪接表达上调与风险降低相关（OR < 1.00）。这四个区分别是 chr1:98293752–98386440（对应 DPYD）、chr13:44449063–44453767（对应 LACC1 与 CCDC122）、chr12:110814021–110819557（对应 ANAPC7）、chr11:46391100–46392863（对应 DGKZ）。与这些预测一致，差异剪接分析在上述四个剪接区都验证了 SCZ 病例中剪接表达的变化。合起来，这些表达数据为五个受 RNA 剪接影响、可能在 SCZ 中起主要作用的基因提供了方向一致的证据：DPYD、LACC1、CCDC122、ANAPC7 和 DGKZ。此外，另外七个基因在 MR 与差异剪接分析中的效应方向不一致，作者认为这可能来自 sQTL 与 GWAS 数据集所用组织的生物学异质性。
顺便说一处原文的不一致：这一节的小标题写的是“10 个内含子使用区中 16 个基因的剪接失调”，正文却写“10/31 区、对应 12 个 sGenes”，摘要和讨论又都说 12 个基因。三个地方三个数字，你可以确认是标题写错了，但这也提醒你：不要默认论文内部的数字必然自洽。
通路富集。为了解 MR 挑出的基因调控哪些生物学过程，作者取了 MR 结果里按 P 值排序的前 500 个基因（因为 Bonferroni 校正后没有一个显著），用 GO 和 WikiPathways 两套方法做富集。GO 富集显示，这些基因在突触组织（synapse organization）、基于微管的运动、认知、沿微管运输等生物学过程中显著富集；其中最显著的正是突触组织。这与既往“突触功能与 SCZ 病理存在强遗传关联”的报道一致，为“突触发育是 SCZ 风险中被破坏的关键过程”提供了独立支持。Wiki 通路富集进一步显示 ADHD 与 ASD 通路，以及与 ASD 相关的突触信号通路——考虑到越来越多证据表明 ADHD、ASD 与 SCZ 之间存在共享的遗传风险，这个结果并不意外。
单细胞。为看这些易感基因在各种脑细胞类型中的表达，作者用认知正常者 PFC 的单细胞 RNA-seq 数据做聚类。在 MR 分析鉴定出的 27 个 SCZ 易感基因中，有 13/27（DPYD、CCDC122、TBC1D5、GPM6A、DGKZ、CCDC92、LRRN3、MPHOSPH9、FAM49B、IMMP2L、AKT3、PRKD3、KANSL1）富集于一种或多种细胞类型，包括少突胶质细胞、抑制性神经元、兴奋性神经元、星形胶质细胞和小胶质细胞。其中三个基因显示出对特定细胞类型的富集证据：DPYD、CCDC122 和 KANSL1。例如 DPYD 与 CCDC122 在少突胶质细胞中高表达。
讨论与局限（精译）
作者开头就说清了缺口：GWAS 迄今已找到 200 多个 SCZ 风险位点，但它们如何带来风险仍不清楚；而 RNA 剪接已被报道在 SCZ 发生中起关键作用。考虑到 sQTL 在 SCZ 病例中一直没有被很好地刻画，作者用 539 例 SCZ 样本的基因型与 RNA-seq 数据系统做了全基因组 sQTL 分析，再用这套 sQTL 与最新的 SCZ GWAS 做 MR 整合，得到 31 个内含子使用区中的 27 个“可能经 AS 调控导致 SCZ 发病”的因果基因；并发现 12 个基因在病例与对照之间剪接表达异常。关键在于，其中 5 个（DPYD、LACC1、CCDC122、ANAPC7、DGKZ）在 MR 与差异剪接分析中方向一致——作者认为这强烈提示它们可能是 SCZ 的潜在新治疗靶点。
为什么盯着 DPYD。DPYD（位于 chr1:98293752–98386440 内含子使用区）在 27 个显著 MR 结果中关联最强（P = 1.06×10−18），其在 PFC 脑组织中由遗传决定的 AS 调控可能对 SCZ 有因果作用。MR 显示，遗传性升高的 DPYD 剪接表达与脑组织中 SCZ 风险升高相关（OR = 1.26）；与差异剪接分析对照，DPYD 的剪接表达水平在 SCZ 病例中相对正常对照上调最显著（P = 5.25×10−5）。这些一致的结果强烈提示 DPYD 是脑 SCZ 的潜在因果基因。此外，DPYD 落在一个已被报道的 GWAS 命中信号之内，其单细胞分析又把它指向少突胶质细胞。几条证据叠起来，作者认为 DPYD 是有希望的治疗靶点。
关于细胞类型。基因调控高度依赖情境，常表现出细胞类型和发育阶段的特异性；因此某些基因可能只在特定的细胞环境里才贡献 SCZ 遗传风险。13 个 MR 鉴定的易感基因富集于一种或多种细胞类型，其中 DPYD 与 CCDC122 在少突胶质细胞中特异高表达，说明这些易感基因的细胞类型特异表达模式与 SCZ 的病理生理密切相关。越来越多证据表明少突胶质细胞相关基因的异常表达会严重损害髓鞘形成或维持，而髓鞘异常会破坏突触传递的精度与同步、导致突触功能障碍，进而影响正常的神经环路功能，最终造成 SCZ 患者的认知、情感和行为症状。
关于“是不是只是把已知的再算一遍”。为评估 SCZ 与 GTEx DLPFC 组织之间 sQTL 的共享程度，作者统计了两套数据都显著的 sGenes 数目：7083 个 SCZ sGenes 中有 850 个（约 12%）与 GTEx 前额叶的 sGenes 重叠；剩下 6233 个不重叠的 sGenes 代表潜在的 SCZ 特异调控者。更进一步，考虑到病例特异的 sQTL 有可能来自反向因果，以及人群间 LD 模式差异会影响 MR，作者又用 GTEx 的 sQTL 数据配 SCZ GWAS 做了一遍 MR：用 GTEx 的 cis-sQTL 当工具变量，得到 20 个达到 MR 显著性的基因；27 个基因里有 10 个在两套 sQTL 数据中均显著，其余 17 个只在本研究的 SCZ 特异 MR 中显著。而在这 10 个重复的基因里，有 4 个在两种组织中的 MR 效应方向不同。例如 DGKZ：在 SCZ 脑组织中，遗传性升高的表达与 SCZ 风险降低相关（OR = 0.88；95% CI 0.86–0.91；P = 1.55×10−13），在 GTEx 健康脑组织中却与风险升高相关（OR = 1.16；95% CI 1.12–1.21；P = 1.32×10−13）。作者认为，这说明整合疾病特异的 sQTL 数据能提供关于 SCZ 机制的新认识。
局限。作者列了四条：①sQTL 数据与 GWAS 汇总统计主要来自欧洲人，结论外推到其它人群受限，因此有必要在其它族裔人群中扩展研究；②本研究的 sQTL 数据主要来自脑组织整体（bulk），应使用外部的单细胞 sQTL 数据来研究 SCZ 相关的因果基因；③MR 中只纳入 cis-sQTL SNP 作为工具变量，以维持假设（原文在此处被截断，指的是排他性假设）；④（讨论后文亦提到）需要实验验证。此外，病例特异的 sQTL 存在反向因果的可能，这也是作者自己点出的。
这篇在轨迹里的位置
把这一篇放到王俊阳的发表线上看，它是“换语言”的那一格。2024 年那两篇脑出血（艾司氯胺酮、鼻内 rIL-10）是湿实验：给药、行为学、WB、免疫荧光。2025 年的三篇突然全是人的数据：J Neurochem 上缺血性卒中患者的 TCR 免疫组库（12 人里第 3）、这篇 SCZ 特异的 sQTL（7 人里第 5）、BMC Medicine 上 rs1950834 与 LRFN5（14 人里第 8）。一年之内，从“老鼠走得好不好”跳到“这个碱基改变的是哪一版转录本”。
这一篇作者只有 7 人，通讯是夏俊峰（Junfeng Xia），做的是计算遗传学。全文没有一处实验验证，所有结论都由公开数据和统计算法得出——这恰好是“干实验”最纯粹的样子。全文里也没有附作者贡献声明（我核过全文，没有 Contributions 一节），所以我不会替他编“他负责了哪一块”。从 7 人里第 5 的位置看，比较可能的是数据整理与分析链条中的某一环；你可以把这句话当成一个待确认的假设，而不是结论。
再往下看一步：这篇和 03 页那篇 LRFN5 恰好是一对。这一篇用统计把“关联”推到“因果”（MR＋共定位＋差异剪接），但停在“13 个基因富集在少突胶质细胞”这个层面；03 页那篇则反过来，从一个 SNP（rs1950834）出发，用报告基因、EMSA、ChIP、CRISPR 单碱基编辑、小鼠行为学，把“关联”推到“功能”。两页并排读，你会看到同一个问题的两端：一端问“哪些基因值得追”，另一端问“这个碱基到底怎么起作用”。
2026 年他在这条线上又进了一步：Nature Communications 上 1p36.23 变异经 RERE 影响 SCZ 风险，16 人里排第 2、共同第一作者。从“第 5 作者的计算分析”到“共同一作的位点功能研究”，这是同一套语言练熟之后的自然结果。
给刚进组的你
先做一件最具体的事：找一张 LeafCutter 的输出图，把某个内含子簇的 PSI 自己手算一遍。PSI 就是“这个内含子被保留的比例”，分子是支持这条内含子的 reads，分母是支持这条内含子加上跨越它的那对剪接连接的 reads。你在组会上听到“这个外显子跳跃增加”，说的就是这个比值变了。会算它，你就不会被“剪接上调/下调”这种说法绕晕。
第二件事：把 MR 的三个假设抄在笔记本第一页——相关性、独立性、排他性。然后回头看这篇的 31 个区为什么最后每个只剩一个 SNP。你会明白为什么作者说“只能用 Wald ratio”：不是不会算 IVW，而是数据只允许一种算法。这种“方法被数据逼着走”的细节，比结论更值得学。
第三件事：如果你以后要做这类分析，记住两个坑。一是 LD：同一区块里几十个 SNP 的 P 值都很好，那不是几十个发现，是一个发现。二是反向因果：病例特异的 sQTL 可能是疾病本身改变了剪接机器，而不是剪接变化导致疾病。共定位能排除“两个独立信号碰巧挨着”，排不掉反向因果——这一点作者自己在讨论里承认了。
最后，可以带着这三个问题去问师兄师姐：①咱们组有没有自己的 eQTL/sQTL 数据，还是全用 GTEx？②差异剪接分析的阈值（这里名义 P < 1×10−3）为什么不用 FDR？③如果要做实验验证，第一个会选 DPYD 还是 KANSL1，为什么？这三个问题没有标准答案，但能问出来，说明你已经读进去了。
▲ 本部分翻译难点（对照学习用）
附录 · 英文原文
Integrating genetic regulation and schizophrenia-specific splicing quantitative expression with GWAS prioritizes novel risk genes for schizophrenia
Abstract
Alternative splicing (AS) plays a vital role in the pathogenesis of schizophrenia (SCZ). Previous studies have linked the genetic signals from genome-wide association studies (GWAS) with expression quantitative trait loci (eQTL), but the interplay with other genetic regulatory mechanisms, particularly splicing QTL (sQTL), remains unclear. Here, we constructed a comprehensive disease-specific sQTL map to provide genetic variants that could alter gene activity through RNA splicing in SCZ. We analyzed data from 539 SCZ patients, identifying a total of 24,810 significant sQTLs (FDR < 0.05) involving in AS events of 7083 unique genes. By combining this with a large-scale SCZ GWAS, we employed Mendelian randomization (MR) and colocalization analyses to pinpoint 27 significant risk genes with genetic AS regulation that may play a causal role in SCZ. Additional differential splicing analysis of these genes in 539 cases and 754 controls revealed 12 significant genes that may increase SCZ risk due to their AS dysregulation. Notably, five genes (DPYD, LACC1, CCDC122, ANAPC7, and DGKZ) showed consistent splicing regulation effects in both MR analysis and differential splicing analysis. Pathway enrichment analysis of differentially spliced genes revealed potential biologically pathways relevant to SCZ, particularly in synaptic transmission and microtubule movement. Furthermore, single-cell RNA-seq analysis revealed that several genes were preferentially expressed in specific brain cell types, including oligodendrocytes, microglia, and excitatory neurons. Overall, our findings highlight several susceptibility genes that may contribute to SCZ risk by AS regulation. Further characterization of these genes could advance mechanistic understanding and therapeutic discovery for SCZ.
Introduction
Schizophrenia (SCZ) is a severe and complex psychiatric disorder characterized by abnormalities in cognition and thought, with a worldwide lifetime prevalence of around 1% [[1]–[3]]. Due to the high morbidity and mortality, SCZ imposes enormous economic and medical burdens on individuals, families and societies [[4], [5]]. However, the pathogenic mechanism of SCZ development is still largely unclear, and existing therapeutic treatments have shown limited benefits [[6]]. Therefore, there is an urgent need to identify effective and specific biomarkers for the development of therapeutic strategies for SCZ. The heritability of SCZ is estimated as being at least 80%, indicating that genetic factors play a dominant role in the pathogenesis of SCZ [[7]]. The emergence of the genome-wide association study (GWAS) has created an unprecedented opportunity to dissect the genetic etiology s of SCZ. Over the past decade, GWAS has identified hundreds of risk loci associated with SCZ [[8]–[11]]. However, interpretation of the GWAS findings into biology insights and clinical applications remains a great challenge.
Alternative splicing (AS) of pre-mRNA is an essential step in the post-transcriptional gene regulation that removes intronic sequences and links exons specifically [[12]–[14]]. Over 95% of multi-exon genes in humans are subjected to AS, greatly increasing transcriptome and protein diversity [[12], [15]]. Most recent studies have shown that AS is widely present in the nervous and immune systems, and aberrant AS is associated with a variety of brain disorders, especially in SCZ [[16]–[21]]. Therefore, unraveling the regulation mechanisms of AS is essential to better understand the pathogenesis of SCZ. Existing evidence indicates that AS regulation can be controlled by genetic variants, and splicing quantitative trait loci (sQTL) has been widely used to explore genetic variants of AS regulation underlying human disease [[22]–[26]].
Our current understanding of genetic variants that affect AS and their underlying pathogenic mechanisms in SCZ is still limited. Therefore, integrative approaches that combine sQTL information with GWAS findings have emerged and shown promise in exploring the potential risk genes whose splicing expression levels are affected by the identified risk variants. Mendelian randomization (MR) is a representative integrative approach that uses risk variants associated with splicing quantitative expression as instrumental variables (IVs) to infer the causal influence of an exposure (i.e., risk gene affected by RNA splicing regulation) on an outcome (i.e., disease) [[27]–[29]]. By integrating the GWAS genetic findings and sQTL data, MR could infer risk genes that may have a causal role in SCZ. Furthermore, given that RNA splicing is highly heterogeneous in the brain [[30]], sQTL data based on non-target disease samples may obscure the role of disease-specific AS regulation, resulting in important biological insights that would be missed. Therefore, MR integrative analysis using disease-specific sQTL data will provide new insights into the disease-specific AS regulation mechanism.
To the best of our knowledge, there have been no disease-specific sQTL resources for SCZ. In this study, we first collected 539 SCZ samples with both genotype and transcriptome from two public consortiums and performed a genome-wide sQTL analysis to identify genetic variants that affect AS. In total, we identified 24,810 significant sQTL single-nucleotide polymorphisms (SNPs) using stringent filtering criteria. Furthermore, we performed a comprehensive MR study by integrating the identified SCZ-specific sQTL data with GWAS of SCZ, and proposed 27 significant genes whose genetically AS regulation may have a causal role in SCZ. By combining evidence of colocalization and differential splicing analysis, we identified 12 promising risk genes for SCZ. In addition, the single-cell transcriptomic analysis revealed that 13 genes are enriched in brain cell types, including oligodendrocytes, inhibitory neurons, excitatory neurons, astrocytes, and microglial cells. Collectively, our study offers a comprehensive resource of SCZ-specific sQTL map and provides a set of promising novel drug targets with strong evidence for SCZ.
Materials and methods
SNP genotyping and RNA sequencing (RNA-seq) data of SCZ participants
We used the SNP genotyping and RNA-seq data in brain tissue from two cohorts, including the CommonMind Consortium (CMC [[31]]) and the Lieber Institute for Brain Development (LIBD [[32], [33]]). Briefly, all quality-controlled DNA genotyping and raw RNA-seq files for the dorsolateral prefrontal cortex (DLPFC) region of the human brain were downloaded from the LIBD database (http://eqtl.brainseq.org) and the CMC portal (https://www.synapse.org/CMC). Eventually, a total of 539 SCZ cases from the CMC (N = 328) and LIBD (N = 211) datasets were included in this study, consisting mainly of European participants (Supplementary Table). For genotyping data, the downloaded DNA genotyping files from each SCZ individual was subsequently merged with PLINK v1.9 [[34]] and totaling 27,332,850 genotyped or imputed markers were used for the following sQTL analysis. For RNA-seq data of LIBD, the downloaded RNA-seq FASTQ data were cleaned using fastp v0.23.4 [[35]] and then aligned to the GRCh37 genome assembly by STAR v2.5.2a [[36]] (converted into BAM files). Furthermore, the obtained BAM files of LIBD were merged with the CMC downloaded BAM files, and the WASP [[37]] tool was then employed to remove reads with potential mapping bias. More detailed information about the sample collection, DNA and RNA extraction and sequencing, quality control and statistical analyses have been described in previous publications [[31]–[33]].
Splicing quantification
To determine splicing quantification on RNA-seq data from DLPFC brain regions of 539 SCZ participants, we chose to quantify AS events using the unannotated LeafCutter [[38]] algorithm. The intron usage rates (i.e., percentage spliced-in value (PSI) value) calculated by LeafCutter were used as the splicing quantification indicator in this study. Specifically, we converted the BAM files into an intron junction file by using the bam2junc.sh script from LeafCutter. Using the leafcutter_cluster.py script, intron clustering was then performed with the following parameters: “-minclureads 50, -mincluratio 0.001, and -maxintronlen 500000”. We mapped intron clusters to genes based on exon coordinates from GENCODE v.19 annotation and the introns present in more than 40% of all samples were selected for further analysis. The PSI values were computed using the prepare_phenotype_table.py script from LeafCutter for the qualifying introns.
Identification of SCZ sQTLs
To identify sQTL SNPs in SCZ, we conducted a cis-sQTL analysis (within 1000 kb up/downstream of the intron clusters) of the obtained PSI matrix and genotype data. With the linear regression models implemented in the FastQTL [[39]] software package, the association between PSI values of AS events and SNP genotypes (i.e., sQTLs) was examined. To control for potential confounders (such as genetic, biological, and technical factors), we performed covariance correction analysis to regress out relevant covariates, including age, sex, RNA integrity number, population structure, and sequencing library batch effects. More details on the sQTL analyses were given in our previous work [[17]]. The intron level sQTL P-values were obtained by applying the permutation procedure (adaptively permute 1000 times) from FastQTL, Finally, sQTLs with Benjamini-Hochberg correction (FDR < 0.05) were considered statistically significant.
SCZ GWAS
To maximize statistical power, summary genetic association data from the largest available GWAS of SCZ [[40]] were used as outcome data in further MR analysis. Summary level data for 7,659,767 SNPs were obtained from the Psychiatric Genomics Consortium (PGC) data portal. Briefly, Trubetskoy et al. performed a large-scale trans-ancestry SCZ GWAS consisting of 74,776 SCZ cases and 101023 control individuals, which included European, Asian, African American, and Latino ancestry populations. Ultimately, they reported a total of 342 genome-wide independent significant SNPs located in 287 distinct genomic regions. To avoid biases due to variations in LD and allele frequencies, only GWAS from European populations (53,386 cases and 77,258 controls) were considered in this study. Detailed information on sample collection, genotyping, quality control, and statistical analyses can be found in the original publication [[40]] and the PGC website (https://www.med.unc.edu/pgc).
MR analysis
Genetic variation associated with RNA splicing was used as an IV to assess the causal association between exposure (i.e., SCZ sQTL data) and outcome (i.e., SCZ GWAS data). MR analysis was conducted using the “TwoSampleMR” R package (version 0.5.6) [[41]]. Prior to the MR analysis, we harmonized the exposure and outcome data to ensure the same effect allele of the SNP was used in both the sQTL and GWAS datasets. IVs were then performed linkage disequilibrium (LD) clumping using a window of 5000 kb and a low LD (r2 < 0.01) between IVs to ensure that the IVs (i.e., SNPs) were independent. MR analyses employ the Wald ratio method when only one cis IV is considered, and the inverse variance weighted (IVW) method when two or more cis IV are considered. Specifically, IVW combines the Wald ratio estimates of each individual SNP into one causal estimate for each risk factor. As our work included only one IV, we did not undertake any sensitivity analyses [[42], [43]]. To account for multiple testing, a Bonferroni correction was applied to adjust for 5179 independent tests (0.01/5179, P = 1.93 × 10−6, 5179 is the number of effective splicing sites valid AS events used for MR analyses). More details of the MR analysis can be found in the original papers [[44], [45]].
For the MR analysis, the same IV (i.e., SNP) should influence both exposure factor and outcome factor, rather than sharing IVs coincidentally due to LD. To assess the probability of the same IV being responsible between SCZ and sQTL, we further performed colocalization analysis for SCZ risk using the Bayesian approach implemented in the R package Coloc v.5.1.0 [[46]] (https://github.com/chr1swallace/coloc). Specifically, colocalization analysis was conducted to adjust such spurious results and posterior probabilities for five hypotheses (H0, H1, H2, H3, H4) were calculated. The correct hypothesis above is H4, and PPH4 (posterior probability for hypothesis 4) specifically quantifies the probability that both traits are driven by the same causal variant within a splicing region. To ensure the reliability of the MR results, we set a strict significant threshold for the posterior probability (two significant associations sharing a common causal variant) at PPH4 > 0.90 in colocalization analysis. Further details about the principle of the colocalization analysis have been published previously [[45]].
Differential splicing analysis
To investigate the RNA splicing level of the MR significant results in SCZ cases compared with controls, we obtained publicly available transcriptome RNA-seq data of 539 SCZ cases and 754 controls from CMC and LIBD datasets. The LeafCutter [[38]] was employed to generate PSI matrices with the same processing procedure as previously described [[17]]. Then, differential splicing analyses was carried out using the Wilcoxon rank sum test to compare PSI values (i.e., the RNA splicing level) between SCZ cases and controls, and the P-value < 1 × 10−3 was considered significant. More information on sample and RNA-seq handling protocols can be found in the original publication [[31]–[33]].
Functional enrichment analysis
There was no significant correlation the MR-identified genes after applying the Bonferroni correction, which is likely due to the assumption of independent tests that this correction requires. Thus, the top 500 MR results with the smallest P-value were selected for the functional enrichment analysis, allowing for more genes to be included. We used the clusterProfiler package in R language (version 4.10.0) for functional enrichment analysis, including Gene Ontology (GO) and Wiki Pathways (WP) gene sets pathway enrichment analysis. For the enrichment calculation of Biological Processes, Human (org.Hs.eg.db) gene annotation with Entrez Gene identifiers was used.
single-cell RNA-seq (scRNA-seq) analysis
To explore if MR-significant genes were specifically expressed in specific brain cell populations, we performed a single-cell expression analysis. First, we downloaded raw scRNA-seq data (i.e., FASTQ files) of 24 cognitively normal individuals in the PFC brain region from Mathys et al. study [[47]] (https://www.synapse.org/#!Synapse:syn18485175). Furthermore, Seurat [[48]] (version 5.0.3) workflow was applied to scRNA-seq data for data preprocessing and analysis, including gene and cell quality control, normalization and transformation, and cluster annotation. Specifically, genes expressed in fewer than 3 cells and cells expressing less than 200 genes were excluded. To reduce noise and improve interpretability, we performed principal component analysis (PCA) for all highly variable genes. We employed the FindAllMarkers function in Seurat to find marker genes for each cluster. Our study focused on 6 brain cell types provided by Mathys et al. [[47]], including oligodendrocytes, inhibitory neurons, excitatory neurons, astrocytes, microglial, and oligodendrocyte progenitor cells. To determine if potential SCZ-causal are highly expressed in one particular cell type, we explored the cell-type specific expression of these genes using the Wilcoxon rank sum test. To control the false discovery rate, FDR correction was applied to all the genes analyzed, and genes with FDR < 0.05 were considered significant.
MR analysis in non-SCZ study participants
We further examined whether our MR findings are informative for non-SCZ study populations using sQTL data of non-SCZ populations from the Genotype-Tissue Expression (GTEx) project. Briefly, the GTEx project characterized and released sQTLs in 54 tissues of over 900 healthy individuals. We have downloaded the latest sQTL data (Brain_Frontal_Cortex_BA9.v10) of PFC brain tissues (N = 268) from GTEx v10 and performed MR analysis using identical pipelines and parameters as in our SCZ MR study. The threshold for significant associations with MR evidence was set at P < 2.66 × 10–6 (i.e., Bonferroni corrected P-value cutoff of 0.01/3764 effective splicing sites).
Results
Identification and characterization of SCZ-specific sQTLs
To investigate the disease-specific genetic control of RNA splicing in SCZ, we performed a genome-wide cis-sQTL analysis using 539 SCZ samples with both PSI matrix of AS events and genotype from CMC and BrainSeq datasets (Fig.). We harvested a total of 282,570 AS events and 27,332,850 genotyped SNPs were retained for further sQTL analysis after stringent quality control. Eventually, we identified 24,810 significant sQTL SNPs (FDR < 0.05) involving 7083 unique sQTL-harboring genes (sGenes) in brain PFC regions of SCZ samples (Fig. and Supplementary Table ). To investigate the genomic distribution of sQTL SNPs, we examined the distance between a sQTL SNP and the corresponding nearest splicing junction. Consistent with previous findings [[49]–[51]], sQTL SNPs were enriched around the splice junction (Fig. ). In addition, we observed that roughly 38.2% of sQTL SNP were located within the body of the gene where the corresponding AS event occurred (Fig. ).Fig. 1Identification and characterization of sQTLs.A The pipeline of sQTL discovery is based on 539 SCZ brain tissues from CMC and LIBD datasets. B The Manhattan plot shows the distribution of these sQTL SNPs on different chromosomes. C Position of sQTL SNPs in relation to the splice junction. D Percentage (%) of index sQTL SNPs (the most significant variant within per intron usage region) located in or outside the corresponding sGene.
MR analysis using SCZ-specific sQTL data identified 27 candidate SCZ susceptibility genes
To identify susceptibility genes that causally contribute to SCZ risk by affecting RNA splicing, we performed a SCZ-specific MR study through integrating SCZ GWAS summary genetic data (53,386 cases and 77,258 controls) with SCZ-specific sQTL data (N = 539 SCZ samples). Notably, Wald ratio estimates were exclusively applied in our final MR analysis because only one qualified SNP per splicing region survived rigorous IV selection criteria (minor allele frequency > 0.01, LD pruning r² < 0.01, and F-statistic > 10). This scenario inherently necessitates the Wald ratio approach (the statistically optimal method when a single IV is available), as it provides unbiased causal effect estimates without requiring instrument homogeneity assumptions essential for IVW meta-analysis. Therefore, the MR associations with P-value < 1.93 × 10−6 were considered statistically significant after Bonferroni correction for multiple tests. We found that all MR significant results were robust to colocalization analyses (PPH4 > 0.90). Consequently, we identified 27 genes within 31 intron usage regions that demonstrated a significant association with SCZ risk, supported by robust evidence (Fig. and Supplementary Table). Among which, one significant gene (DPYD) at chr1:98293752–98386440 intron usage region showed the most significant association (P = 1.06 × 10−18). Other significant potential susceptibility genes include DGKZ, ANAPC7, FTSJ2, BCL2L12, IRF3, GPM6A, MPHOSPH9, LRRN3, and IMMP2L. Interestingly, genetically raised splicing quantitative expression (i.e., PSI values) of gene KANSL1 was associated with reduced SCZ risk in chr17:44230332–44248221 splicing site (OR = 0.97; P = 1.56 × 10–7) and increased SCZ risk in chr17:44172067–44248221 splicing site (OR = 1.03; P = 1.56 × 10–7) (Fig. and Supplementary Table), indicating that distinct intron usage regions of the same gene have different biological functions in SCZ.Fig. 2MR Analysis using SCZ-specific sQTL data identified 27 significant risk genes for SCZ.A The Manhattan plot shows all significant sGenes within intron usage regions in MR analysis. The red line indicates the P threshold for Bonferroni correction. B The forest plot shows the results that 31 intron usage regions reached significance in MR analysis. Additionally, we calculated 95% confidence interval for the odds ratio. 95% CI 95%: confidence interval, OR odds ratio, SCZ schizophrenia.
In addition, by comparing MR significant results with the largest SCZ GWAS from PGC3, we found that 19/27 significant genes identified by MR were located at known SCZ susceptibility loci, including DPYD, DGKZ, BCL2L12, IRF3, MPHOSPH9, LRRN3, IMMP2L, FXR1, SUGP1, KANSL1, ANKRD45, TBC1D5, NDUFAF7, PRKD3, FAM114A2, LACC1, CCDC122, AKT3, and GPM6A. These overlapping MR results affected by AS might help to pinpoint potential target genes in each GWAS signal. More importantly, we found that 8/27 genes with AS events did not overlap with known association loci of SCZ, including FTSJ2, CSPG4P12, CCDC92, ZNF664, CRELD2, FAM49B, APOBEC3C, and APOBEC3D. These results indicated that incorporating SCZ-specific sQTL data might facilitate the identification of novel target genes beyond GWAS findings.
Splicing dysregulation of 16 genes within 10 intron usage regions identified by MR analysis in SCZ cases
Significant intron usage regions predicted in MR analyses whose genetically AS regulation might have essential roles in SCZ. After excluding, a total of 539 SCZ cases and 754 healthy controls from CMC and BrainSeq datasets were included in the differential splicing analysis for SCZ. Among the 31 significant intron usage regions, we observed that 10/31 intron usage regions (corresponding to 12 sGenes) were differentially splicing quantitative expression (nominal P-value < 1 × 10−3) in SCZ cases compared with controls (Supplementary Table), suggesting that these overlapping intron usage regions represent promising functional genetic loci for SCZ. Notably, MR Analysis revealed that the upregulation of splicing expression in one intron usage region is associated with an increased risk of SCZ (OR > 1.00) and the upregulation of splicing expression in three intron usage regions is associated with a decreased risk of SCZ (OR < 1.00). These regions included chr1:98293752–98386440 (corresponding to sGene DPYD), chr13:44449063–44453767 (corresponding to sGene LACC1 and CCDC122), chr12:110814021–110819557 (corresponding to sGene ANAPC7) and chr11:46391100–46392863 (corresponding to sGene DGKZ). Coherently with these predictions, differential splicing analysis validated the changes in splicing expression for the above four splicing regions in SCZ cases. Collectively, these expression data provided consistent and convergent evidence in favor of five genes affected by RNA splicing that may have a major role in SCZ, including DPYD, LACC1, CCDC122, ANAPC7 and DGKZ (Fig.). In addition, we found that the remaining seven genes had different directions of effect between MR and differential splicing analysis, possibly due to the biological heterogeneity of different tissues used in sQTL and GWAS datasets.Fig. 3Differential splicing analysis of the four MR-identified significant intron usage regions.Violin plots show splicing dysregulation of four intron usage regions (A-D) in SCZ cases compared with controls. SCZ schizophrenia.
Functional enrichment analysis of the identified MR genes revealed related biological processes
To get insights into the biological processes regulated by the top 500 MR-derived genes (ranked according to P-value), we conducted functional enrichment analysis using two different programs: GO term and Wiki pathway enrichment analysis. Specifically, GO enrichment analysis revealed a strong enrichment of MR genes in biological processes like synapse organization, microtubule−based movement, cognition, and transport along microtubule (Fig.). Among these, the most significantly enriched pathways of the GO analysis were synaptic organization. This finding aligns with previous reports of strong genetic associations between synaptic function and the pathology of SCZ [[52]–[54]], providing independent support for synaptic development as a key process disrupted in SCZ risk. Wiki enrichment analysis further revealed that the ADHD (attention deficit hyperactivity disorder) and ASD (autism spectrum disorder) pathways and synaptic signaling associated with ASD (Fig.). This result is not unexpected, given the increasing body of evidence indicating a shared genetic risk among ADHD, ASD, and SCZ [[55]–[58]].Fig. 4GO and Wiki pathway enrichment analysis of MR-identified risk genes.A GO pathway enrichment analysis of MR-identified risk genes. B Wiki pathway enrichment analysis of MR-identified risk genes.
Cell-type specific expression of the potential SCZ-susceptibility genes
To investigate the expression of SCZ-susceptibility genes in various brain-relevant cell types, we analyzed the activity of these genes across different cell types using scRNA-seq data from the PFC brain region of cognitively normal individuals. Clustering analysis on scRNA-seq data was performed using the Seurat pipeline to identify cell types or subpopulations. We selected 2000 hypervariable genes and the top ten genes were flagged (Fig.). Furthermore, we identified six clusters of different brain cell types (Fig.). Among 27 SCZ-susceptibility genes identified in MR analysis, 13/27 (DPYD, CCDC122, TBC1D5, GPM6A, DGKZ, CCDC92, LRRN3, MPHOSPH9, FAM49B, IMMP2L, AKT3, PRKD3 and KANSL1) were enriched in one or more cell types, including oligodendrocytes, inhibitory neuron, excitatory neurons, astrocyte, and microglial (Fig.). Of note, three genes show evidence of enrichment in a particular cell type, including DPYD, CCDC122, and KANSL1. For instance, the genes DPYD and CCDC122 were highly expressed in oligodendrocytes, as shown in Fig..Fig. 5Single-cell sequencing analysis of the 27 MR-identified risk genes in the human brain.A Dot plot of top 2000 hypervariable genes for brain cell type. The top ten genes were flagged. B UMAP plot shows six cell clusters. C Enrichment of the MR-identified risk genes in different cell types of the brain.
Discussion
Hitherto, GWAS has identified more than 200 risk loci for SCZ, but it is unclear how they confer SCZ risk. More importantly, RNA splicing has been reported to play a key role in the development of SCZ. Considering that sQTL has not been well characterized in SCZ cases, we systematically undertook a genome-wide sQTL analysis using genotype and RNA-seq data derived from 539 SCZ samples. To identify the potential risk genes at SCZ risk loci, we further performed a MR integrative study using the obtained SCZ-specific sQTL data with the latest SCZ GWAS results. We identified 27 potential causal SCZ genes within 31 intron usage regions that act via AS regulation to contribute to SCZ pathogenesis (Supplementary Table). Moreover, we found that 12 genes have displayed aberrant splicing expression in SCZ cases compared with controls. Of note, five of these genes (DPYD, LACC1, CCDC122, ANAPC7, and DGKZ; Supplementary Table) showed the same directions of effect in both the MR and differential splicing analysis. This strongly indicates that these genes could be potential new treatment targets for SCZ.
One interesting finding in this study was regarding gene DPYD associated with the risk of SCZ. Among the 27 significant MR results, DPYD (located in chr1:98293752–98386440 intron usage region) showed a strong significant association (P-value = 1.06 × 10−18), whose genetically AS regulation in PFC brain tissues may have a causal role in SCZ. MR results have shown that genetically increased gene DPYD splicing expression was associated with increased SCZ risk in brain tissue (OR = 1.26). By comparing with differential splicing analysis findings, we found that the splicing expression level of DPYD was most significantly upregulated in SCZ cases versus normal controls (P-value = 5.25 × 10−5). These consistent results strongly suggested that DPYD represents a potential causal gene for brain SCZ. Moreover, we found that DPYD was located within a reported GWAS hit signal, indicating its implication in the latest SCZ GWAS. Furthermore, our scRNA-seq analysis showed that DPYD points to specific cell types that they likely act through to oligodendrocytes to SCZ. These lines of evidence support that DPYD may be a promising treatment target for SCZ.
It is well established that gene regulation is highly context-dependent, often exhibiting cell-type and developmental stage specificity. Consequently, certain genes may only contribute to schizophrenia genetic risk within specific cellular contexts. For the single-cell enrichment analysis, we found that 13 SCZ susceptibility genes identified by MR analysis were enriched in one or more cell types. For instance, genes DPYD and CCDC122 exhibit specific high expression in oligodendrocytes, indicating that the cell-type-specific expression patterns of these susceptibility genes are closely related to the pathophysiology of SCZ. Increasing evidence indicates abnormal expression of oligodendrocyte-related genes, which may severely impair myelin formation or maintenance [[59], [60]]. Myelin abnormalities can disrupt the precision and synchrony of synaptic transmission, leading to synaptic dysfunction. Furthermore, such impairment affects normal neural circuit function, ultimately resulting in cognitive, emotional, and behavioral symptoms in SCZ patients [[61]–[63]].
To evaluate the extent of sQTL sharing between SCZ and GTEx DLPFC tissues, we first counted the number of sharing sGenes that were significant in both SCZ and GTEx. We observed that 850 (approximately 12%) of the 7083 sGenes in SCZ overlapped those in GTEx PFC brain tissues. In addition, the 6233 non-overlapping sGenes represent potential SCZ-specific regulators (Supplementary fig.). Crucially, it is possible that case-specific sQTLs are a result of reverse causation. Considering that differences in LD patterns between populations tend to affect the MR results, we then performed MR analysis using sQTL data of GTEx and SCZ GWAS data and full significant MR results are shown in Supplementary Table. Using cis-sQTLs from the GTEx dataset as proposed instruments, we identified 20 genes that reached MR significance. Furthermore, we found that 10/27 replicated genes showed significance in both sQTL datasets (Supplementary Table), while the remaining 17 potential risk genes were uniquely significant in our SCZ-specific MR analysis. In addition, 4 of the 10 replicated genes had different directions of the MR effect between tissues. For example, genetically raised gene DGKZ expression was associated with reduced SCZ risk in SCZ brain tissues (OR = 0.88; 95% CI, 0.86–0.91; P-value = 1.55 × 10−13) and increased SCZ risk in GTEx healthy brain tissues (OR = 1.16; 95% CI, 1.12–1.21; P-value = 1.32 × 10–13). These results suggest that integrating SCZ-specific sQTL data may provide novel insights into the mechanism of SCZ.
It is important to acknowledge the potential limitations of our study. Firstly, the sQTL data and GWAS summary statistics were mainly from individuals of European ancestry, limiting the generalizability of our findings to other populations. Therefore, expanding research studies to evaluate the association of susceptibility genes with SCZ in other ethnic populations seems essential. Secondly, the sQTL dataset used in our MR integration study was primarily from brain tissue of the human brain. External single-cell sQTL datasets should be utilized to investigate the causal genes associated with SCZ. Thirdly, we only included cis-sQTL SNPs as IVs in the MR study to maintain the assumption that IVs must be strongly associated with the exposure. However, this may ignore the complex role of trans-sQTL in the genetic regulation mediated for SCZ. Lastly, the identified SCZ risk genes require further experiment validation to verify their biological function.
In summary, by integrating unique SCZ-specific sQTL data with the latest SCZ GWAS data, we performed MR analysis and identified 27 candidate susceptibility genes that contribute to SCZ risk through AS regulation. Differential splicing analyses further validated these findings, highlighting five potentially causal genes with the same direction of effect that may pose a risk of developing SCZ. We also identified key pathways and brain cell type specificity important in the pathogenesis of the SCZ. Our findings not only advanced our understanding of the pathological mechanisms of SCZ but also provided valuable targets and directions for developing effective treatments.
Supplementary information
Supplemental information Supplementary Tables
Acknowledgements
AcknowledgementsWe thank the participants of the Psychiatric Genomics Consortium Working Group, CommonMind Consortium, BrainSeq Consortium, ROSMAP Project and Mathys H et al. study for generating and made the summary statistics available for us and making this work possible.