16p染色体上常见和罕见遗传对自闭症影响的统计和功能收敛 – 自然遗传学

　　我们确认本研究已获得马萨诸塞州布莱根将军机构审查委员会 (IRB) 的审查和批准。研究名称是认知和行为变异的分子研究（IRB：2015P002376）。首席研究员是 EBR iPSYCH 研究得到了丹麦数据保护局和丹麦科学伦理委员会的批准。 PGS

　　的生成

　　对于自闭症三人组中的自闭症 PGS 分析，我们使用了来自丹麦 iPSYCH 集合的 GWAS，因为与自闭症三人组样本没有样本重叠（19,870 名自闭症患者，39,078 名对照；补充表 2）。对于所有其他自闭症 PGS 分析（例如，PGS 表达分析），我们使用了对相同 iPSYCH 自闭症样本以及来自 PGC 的自闭症样本的荟萃分析（合并：26,067 名自闭症患者，46,455 名对照）。对于 ADHD 的分析，我们使用了不重叠的仅 iPSYCH 的 ADHD GWAS（25,895 名患有 ADHD 的人，37,148 名对照者）。

　　为了生成 PGS 权重，我们首先将 LDpred v.1.0.11 应用于 GWASs48 的边际效应大小。我们在无穷小的遗传结构模型下使用 LDpred，LD 参考来自 Hapmap 3 SNP（n = 503 个欧洲血统样本）。所有 PGS 都是使用 PLINK 1.9 中的 –score 函数计算的（参考文献 49）。由于 LDpred 从 GWAS 边际效应大小估计后验因果效应大小，我们将所有 SNP 包括在 PGS 分析中，包括在为 S-pTDT 构建分层 PGS 时。

　　自闭症家庭队列

　　SSC 和西蒙斯基金会支持自闭症研究 (SFARI) 的收集、估算和质量控制之前已经描述（补充表 3）12,50。来自 PGC 自闭症组 (PGC) 的自闭症三人组如前所述 12，其中包含来自多重家庭的先证者的修改。我们通过使用 PLINK 1.9 生成祖先的主要成分并通过相对于 HapMap 参考群体的目视检查来定义 PGC 的欧洲祖先子集进行分析（补充图 1）。如果父母和先证者都通过主成分分析 (PCA) 具有欧洲血统，我们将一个家庭定义为欧洲血统（5,283 个三重奏中的 4,335 个，82%）。

　　全基因组 pTDT

　　我们进行全基因组 pTDT 以评估 SSC、SPARK 和 PGC 中 S-pTDT 分析的能力。我们估计了如前所述的多基因传播，除了 PGC 中病例/伪对照基因型的适应方法（补充说明）。三个群组中每一个的结果显示在补充图 2 中。

　　S-pTDT

　　S-pTDT 与 pTDT 相同，只是它不是测试由所有 SNP 构建的 PGS 的传输，而是测试 PGS 的传输由 SNP 的一个子集构成：

　　其中 PGSS,C 是子 C 的分层 PGS，PGSS,MP 是中间亲本分层 PGS（两个亲本的平均值）。 S-pTDT 是单样本双边学生 t 检验，用于检验 S-pTDT 分布是否具有不同于 0 的均值。给定队列的给定 PGS 的 S-pTDT 估计值等于平均 S-pTDT 值对于队列中的所有家庭。

　　我们通过将 SNP 分成大小相等的集合来创建分层的 PGS。我们以两种方式改变这种划分，以完成对区域传播的全面调查：首先，我们将 SNP 划分为不同大小（2,000、3,000、4,000、5,000 和 6,000 个 SNP）的分区，其次，我们从染色体的开始或结束。例如，为了创建 2,000 个 SNP 的基因组块，我们确定了 1 号染色体上的第一个 PGS SNP（最接近第一个碱基对的 SNP），计算了 2,000 个 PGS SNP，并将其定义为第一个分区。然后，我们计算了 1 号染色体上接下来的 2,000 个 PGS SNP，将其定义为下一个分区，直到 1 号染色体上剩余的 SNP 少于 2,000 个。接下来，我们在 2 号染色体和所有剩余的染色体上重复相同的程序。通过改变分区中 SNP 的数量以及 SNP 计数是从染色体的开头还是结尾开始，我们产生了 2,006 个（通常是重叠的）分区（1,003 个从染色体的开头开始，1,003 个从染色体的末端开始）。在分区之前，我们将 PGS SNP 与所有三个自闭症三人组（SSC、SPARK 和 PGC）中存在的 SNP 进行子集化，以避免跨分区 SNP 缺失的偏差。然后，我们使用 PLINK 1.9 中的线性评分 (–score) 估计三个队列中每个分区的分层 PGS，并如上所述在每个分区上执行 S-pTDT。

　　分区长度和 SNP 计数可预测过度传输（补充图 4）。我们使用线性模型——S-pTDT ≈（PGS SNP 数量）+（碱基对中的分区长度）回归出预期的过度传输——并通过模型残差分布的 sd 对模型残差进行归一化。这个过程为每个分区产生一个残差 z 分数，它估计了分区相对于预期传输过度或传输不足的 sd 数量。如果分区包括相邻 SNP 之间的间隙 > 1 Mb，我们将该间隙的贡献调整到 1 Mb，这解释了 LD 的衰减，但避免了在上述过度传输模型中不适当地校正 S-pTDT 信号。对于 33-Mb 分区的分析，S-pTDT z-score 仅回归 SNP 数，因为所有分区的碱基对长度相同。

　　S-pTDT 关联的补充分析

　　首先，我们证实通过 GWAS 的自闭症相关基因座在 S-pTDT 分布中富集。我们将自闭症相关基因座定义为最近发表的自闭症 GWAS 中的五个基因座，仅通过自闭症分析就达到了全基因组意义（索引 SNP：rs910805、rs10099100、rs201910565、rs71190156 和 rs111931861）11。

　　接下来，我们分析了自闭症相关 CNV 在 S-pTDT 分布中的分布。我们从 SFARI Gene (https://gene.sfari.org/database/cnv) 上的集合中识别出自闭症相关的 CNV，然后识别出边界内至少有一个这些 CNV 的 S-pTDT 分区 (16p11.2, 16p13.11、16p13.3、16p12.2、2p16.3、15q13.3、7q11.23、17q12、3q29、1q21.1、17p11.2、8p23.1、17q11.2、2q11.2、22q11。 2, 22q13.3 和 5q35；补充图 10)。我们还估计了 33-Mb 分区的节段重复内容和 S-pTDT 之间的关联：我们通过计算每个分区中与加州大学至少一个节段重复重叠的核苷酸分数来注释每个分区的节段重复率，圣克鲁斯 (UCSC) 基因组浏览器51。使用 BEDTools v.2.30.0（补充图 11）52 进行覆盖率计算。

　　为了排除 16p CNV 载体驱动 S-pTDT 信号的贡献，我们在去除先证者在 16p 携带遗传或新发神经发育障碍相关 CNV 的三重奏后，在 SSC + SPARK 中重复 S-pTDT 分析（我们无法执行在 PGC 中进行此分析，因为我们没有对该队列进行外显子组测序）。我们从最近的自闭症测序研究中采用了基于文献的神经发育障碍相关 CNV 的定义。在 SSC + SPARK 分析中的 5,048 个三重奏中，我们删除了 51 个（1.0%）具有合格 CNV 并重复 S-pTDT 分析（补充图 9）。

　　我们检验了 S-pTDT 与发育中的人脑中可接近染色质密度相关的假设，该密度由 H3K27ac 组蛋白标记标记。我们分析了已发表的染色质免疫沉淀 (ChIP) 资源——对受孕后 7 周的人类皮层的两个生物学重复进行测序分析33。在每个复制中，我们将每个定义的 33-Mb 分区内的 H3K27ac 峰计数相加，将计数缩放为均值 = 0 和 sd = 1，然后平均两个复制之间的 z 分数。

　　为了评估 S-pTDT 在 16p 上发现的特异性，我们在 ADHD 中进行了类似的分析。我们使用了来自 PGC 的 1,634 名欧洲血统 ADHD 三人组和来自 iPSYCH 联盟的外部 ADHD GWAS，其中 25,895 名患有 ADHD 的人和 37,148 名对照者 53。如上所述，我们从染色体的开头将基因组划分为 2,000、3,000、4,000、5,000 和 6,000 个 SNP 的块，并估计每个分区的 S-pTDT。然后我们估计了一个残余 z 分数，回归了自闭症分析中的 SNP 数量和分区大小（补充图 12）。

　　我们接下来评估了 16p 的多基因信号是否可以通过该区域内的特定基因座来解释。为了执行此分析，我们将 16p 划分为 25 个独立于 LD 的块，每个块的大小约为 1.5 Mb，如前所述 31。然后，我们使用 SSC 和 SPARK 中仅 iPSYCH 的自闭症总结统计来估计这 25 个块中的每一个的 S-pTDT（补充图 8）。为了评估单个基因座的贡献，我们估计了 16p S-pTDT 信号的衰减，作为去除最过度传输的剩余 S-pTDT 块的函数。具体来说，我们 (1) 估计每个块的传输，(2) 将块从过度传输最多到最少排序，(3) 使用来自所有块的 SNP 估计过度传输，(4) 使用来自所有块的 SNP 估计过度传输从最相关的剩余块中减去SNP，并且（5）重复步骤4直到只剩下一个块（图1e）。例如，第一个块（“从 PGS 中删除的 16p 分区数 = 0”）包括 16p PGS 中的所有 7,658 个 SNP。下一个块（“从 PGS 中删除的 16p 分区数 = 1”）从 16p 中最相关的块中减去 287 个 SNP，从而使这个新块具有 7,371 个 SNP。

　　我们接下来评估了 16p 的区域多基因信号相对于基因组中相同大小的比较分区。由于 16p 跨越大约 33 Mb 的基因组，我们通过从染色体的开头开始并定义相邻的 33-Mb 块来构建 33 Mb 的控制分区（补充表 4）。我们通过 (1) 1000 Genomes Phase 3 EUR 中的第一个 SNP 和 (2) gnomAD54,55 中第一个基因的起始位置的最小值来定义染色体的起始坐标。类似地，我们通过 (1) 1000 Genomes Phase 3 EUR 中的最后一个 SNP 和 (2) 以及 gnomAD 中最后一个基因的末端位置的最大值来定义染色体的末端坐标。这种方法产生了 74 个分区，包括 16p。我们使用这些边界通过从给定分区内的所有 SNP 构建分层 PGS 来执行 S-pTDT。

　　基因密度

　　我们首先编制了一份用于基因密度分析的共有基因列表。我们将此共识列表定义为 (1) 具有唯一基因名称的常染色体基因和来自 gnomAD 的非缺失 pLI 约束估计值和 (2) 在基因型组织表达 (GTEx) 皮层 ('Brain_Cortex') 中具有估计特异性表达的基因的交集56。该共识列表包括 17,909 个基因。在最近的一项分析中，我们通过外显子组测序进一步用与自闭症有关的 102 个基因对这个列表进行了注释32。然后，如果基因体中点位于边界内，我们将基因映射到上述定义的 33-Mb 边界。我们建立了线性模型，根据所有基因的密度预测皮层中的特定表达（特定表达 t 统计量的前 10%），并从回归的残余 z 分数计算两侧 P 值。

　　GO 分析

　　我们进行了 GO 分析以评估 16p 上注释的生物途径中基因的富集（//geneontology.org）。我们使用来自基因密度分析的相同的 17,909 个基因作为参考基因。我们在三类注释中测试了 16p（中点 <33,000,000 bp，n = 433 个基因）上所有基因的富集：生物过程、分子功能和细胞成分。 Bonferroni 在每一类注释中的显着富集结果报告在补充表 1 中。

　　自闭症差异表达分析

　　我们分析了 16p 上的基因是否在被诊断患有自闭症的个体和对照之间的差异表达基因中过度表达。我们使用了先前发表的人脑 RNA-seq 差异表达数据集（n = 51 名自闭症患者，n = 936 名对照）34。我们将基因定义为在 Bonferroni 显着水平上差异表达的基因，校正了数据集与上述共有基因列表之间重叠的基因数量（n = 15,288 个基因，n = 83 个差异表达基因，n = 16p 上的 383 个基因）。我们对 16p 上差异表达基因的富集进行了 χ2 检验。

　　CRISPR-Cas9 介导的 ASD 相关位点

　　缺失 CRISPR 介导的 16p11.2 位点缺失的设计和差异表达分析在已发表的资源中进行了描述。 15q13.3 CNV 删除系的设计遵循相同的方案，删除定义为边界 Ch15 30,787,764-32,804,328 (GRCh37)。我们生成了 n = 11 个杂合缺失系，另外还有一个 n = 6 个对照暴露于 CRISPR 构建体，但不引导 RNA。在估计缺失对该区域的影响时，我们排除了缺失窗口±100 kb的基因，因为这些基因的顺式调节区域可能已被缺失本身的产生所扰乱。

　　PGS 表达分析

　　HBTRC/NIH NeuroBioBank (snRNA-seq)

　　我们从 HBTRC/NIH NeuroBioBank 的死后脑组织的背外侧前额叶皮层 (DLPFC) 生成配对基因型和单核表达谱。本报告作者即将发表的手稿中将详细描述表达谱的生成。简而言之，我们开发并优化了一个工作流程，用于创建和分析从每个池中 20 个不同供体的脑组织（DLPFC、BA46）采样的细胞核池。在这个工作流程中，我们首先从每个供体解剖确定数量的组织，从每个样本中获得相似质量的组织，同时小心地代表所有皮质层。然后立即合并冷冻组织样本，同时分离它们的细胞核；所有后续处理步骤，包括核分离、液滴封装和 snRNA-seq 文库的制备，都涉及所有供体。这种“dropulation”工作流程使我们能够最大限度地减少实验变异性，包括对信使 RNA 确定的任何技术影响以及无细胞环境 RNA 的任何影响。然后使用数百个转录 SNP 的组合将这些实验中的每个细胞核重新分配给其来源供体。尽管单个 SNP 等位基因在许多供体之间共享，但许多 SNP 的组合对于队列中的每个供体都是独一无二的。通过全局聚类和鉴定每个簇中表达的标记基因，将细胞核分配到七个主要细胞类别（星形胶质细胞、内皮细胞、γ-氨基丁酸 (GABA) 能神经元、谷氨酸能神经元、小胶质细胞、少突胶质细胞和多突胶质细胞）。中位细胞类型比例为：谷氨酸能神经元 47.9%、GABA 能神经元 18.8%、星形胶质细胞 13.5%、少突胶质细胞 8.0%、多突胶质细胞 5.2%、小胶质细胞 1.5% 和内皮细胞 1.0%。所有下游分析都使用来自谷氨酸能神经元的表达数据。用 VST 归一化处理细胞类型特异性基因逐个供体表达矩阵。

　　我们执行了一些预关联质量控制 (QC) 步骤。大多数基因分型样本是欧洲血统（1,707/1,770, 96%），我们使用 PCA 确定了这些样本用于下游分析（补充图 18）。接下来，我们将任何样本识别为来自群组平均值（3/125 样本）的平均表达＞3 sd 的表达异常值。这产生了 122 个样本的最终 EUR 子集。我们在计数矩阵中的 122 个样本中鉴定出表达高于中位数的基因，其中唯一的处理步骤是对所有样本的表达计数和进行标准化处理。我们在所有后续的区域 PGS 表达分析中使用了这些基因（n = 8,878）。

　　CommonMind (bulk cortical RNA-seq)

　　我们接下来分析了来自 CommonMind 联盟中捐赠者的配对基因型和大量 DLPFC 表达数据。 CommonMind 出版物中描述了表达式计数矩阵的生成。在 CommonMind 中，我们将分析限制在来自国家心理健康研究所 (NIMH) 人脑收集核心 (HBCC) 和匹兹堡大学 (PITT) 生物库的捐赠者，因为之前的分析表明与上述 snRNA-seq 资源的一致性增加。我们使用 Seurat 中的 SCTransform 包分别在 HBCC 和 PITT 中对计数矩阵进行了方差稳定，目的是与单核资源的方法（参数：do.scale = FALSE、do.center = FALSE、return.只有.var.genes = FALSE，seed.use = NULL，n_genes = NULL)58。

　　CommonMind 的收藏在祖先上是异质的：最大的两个群体是非洲和欧洲的祖先捐赠者。因此，我们使用 PCA 来识别非洲血统（n = 193）和欧洲血统（n = 229）的供体，随后分别进行分析（补充图 19）。为了与单核资源保持一致，我们将分析限制在被诊断患有精神分裂症的供体或对照组。

　　每队列关联和荟萃分析

　　对于所有样本，我们使用最大的自闭症 GWAS（iPSYCH + PGC；见补充表 2）计算区域自闭症 PGC，使用 Plink 1.9 评分和基因型 QC（SNP 缺失 <1%，次要等位基因频率 > 0.1%，插补信息 > 95%）。

　　我们进行了两类局部 PGS-基因表达关联。第一类是每个基因的关联，如图 3a 所示。第二个是平均基因关联，如图3b，c所示。为了在数据集之间保持一致，我们将所有分析限制在 HBTRC 数据中谷氨酸能神经元中表达最多的基因的一半（n = 8,878 个基因）。图 3a 中的关联是对三个队列进行元分析的每个基因关联。我们将三个队列的个体水平表达和基因型 PGS 结合起来；在连接个体队列分析中使用的 PGS 和表达矩阵之前，我们在队列内将每个基因表达和每个分区 PGS 缩放为均值 = 0 和 sd = 1。每个基因关联遵循线性模型：基因表达 ≈ 区域PGS + 精神分裂症诊断状态 + 血统（二进制表示是/否非洲血统）+ 单细胞（二进制是 / 否）。关联 t 统计量来自区域 PGS 协变量。为了最大程度地检测平均效应，我们使用排列评估了平均 PGS-表达关联的显着性。具体来说，我们计算了 16p 中的均值（t 统计量），然后对每个队列中的 PGS-供体 ID 进行洗牌，执行关联并计算均值（t 统计量），重复 1,000 次。排列 P 值是观察到的 PGS 比排列后的 PGS 更负的次数。

　　对于第二类关联，我们首先对每个分区的基因表达进行平均，然后进行关联。对于每个队列关联，我们使用线性模型：基因的平均表达≈区域 PGS + 精神分裂症诊断状态。对于综合分析，我们使用线性模型：平均基因表达 ≈ 区域 PGS + 精神分裂症诊断状态 + 血统（二元表示是/否非洲血统）+ 单细胞（二元是 / 否）。我们使用补充图 21 中的主要成分对遗传祖先进行了敏感性分析。Hi

　　-C 分析

　　LCL 资源

　　每染色体 Hi-C 计数矩阵从 http://hic.umassmed.edu 以 1-Mb 分辨率下载GM06990 LCL43。由于计数矩阵是在 hg18 中构建的，我们使用国家生物技术信息中心 (NCBI) 的基因组重映射服务 (www.ncbi.nlm.nih.gov/genome/tools/remap) 将 33-Mb 分区从 hg19 转换为 hg18 .我们将 33-Mb 分区的边界与计数矩阵中最接近的 1-Mb 截止值相匹配。对于这个分析，我们没有分析跨越着丝粒的分区，产生 56 个分区进行分析。对于每个分区，我们将原始分区内接触频率估计为 Hi-C 计数矩阵的非对角元素的平均值。

　　妊娠中期皮质板资源

　　我们从来自三个捐赠者的妊娠中期皮质板样本资源中下载了 0.1-Mb 分辨率、Hi-C 接触矩阵，来自 NCBI 的基因表达综合 (GEO)。如上所述，我们将 33-Mb 分区的边界映射到 Hi-C 矩阵的 0.1-Mb 边界。与 LCL 相比，对角线元素被置零；因此，估计的原始分区内接触频率被估计为矩阵所有元素的平均值。我们分析了与 LCL 分析相同的 56 个分区。

　　由于我们的假设涉及基因组大区域的平均接触行为，而不是对拓扑相关域或基因增强子相互作用的更细粒度的分析，我们分析了每个队列中可用的最大 bin 窗口以增加信噪比 59 .

　　接触模型

　　原始分区内接触频率随基因密度和节段重复接触而变化（补充图23）。在 LCL 中，这种协方差可能是由于映射到具有增加的片段重复内容的区域的 Hi-C 读取数量增加。在皮质线中，接触矩阵中有大块归零的元素，其比率与节段重复内容密切相关，可能是由于在由于节段重复内容而难以映射的区域中故意将元素归零。在调节节段重复内容后，基因密度仍然是接触频率的重要预测因子，促使我们也调节基因密度并从以下模型中提取归一化残差：接触频率≈基因计数 + 节段重复内容。我们在图 4b 中的主要分析报告了每个分区的这些归一化残差。

　　16p11.2 CNV 的端粒区域分析

　　我们接下来分析了 16p11.2 基因座与 16 号染色体远端基因密集起始点之间的染色质接触。我们在妊娠中期皮质板数据中进行了这项分析，只是因为 1-Mb bin 分辨率LCL 资源没有足够的分辨率。我们根据 16 号染色体开始处的基因密集片段定义了端粒区域，从 0 Mb 到 16p11.2 删除数据集中该窗口中最终大脑表达基因的终点后最近的 100 kb 片段。 5.2 兆字节）。我们将 16p11.2 基因座定义为 Ch16：29.5–30.2 Mb。为了定义控制接触区域，我们首先计算了由 0-5.2 Mb（端粒区域）和 29.5-30.2 Mb（16p11.2 基因座）定义的接触矩阵跨越的最小 (24.3 Mb) 和最大 (30.2 Mb) 距离。然后，我们将 16p 上的控制接触定义为距离大于 24.3 或 <30.2 且不位于上述端粒-CNV 接触范围内的所有接触。这些结果对于包含具有“0”的接触矩阵的元素是稳健的，这可能反映了节段重复丰富的区域（端粒-CNV 与对照 P < 1 × 10-10 对于两种方法）。

　　报告摘要

　　有关研究设计的更多信息，请参阅与本文链接的自然研究报告摘要。