利用染色质配置的计算分析，将阿尔茨海默氏病变异体映射到其目标基因

摘要

我们提出了一个协议，用于识别使用三维染色质相互作用的全基因组关联研究（GWAS）识别的非编码变异的功能影响。

摘要

全基因组关联研究（GWAS）已经成功地确定了数百个与人类特征和疾病相关的基因组位点。然而，由于大多数全基因组显著（GWS）位点都属于非编码基因组，许多功能影响仍不得而知。由Hi-C或其衍生物识别的三维染色质相互作用可以通过将非编码变异与其可操作基因联系起来，为这些位点进行说明提供有用的工具。在这里，我们概述了一个协议，以映射GWAS非编码变体，以映射到他们的假定基因使用阿尔茨海默氏病（AD）GWAS和Hi-C数据集从人类成人脑组织。通过应用精细映射算法来识别假定因果单核苷酸多态性 （SNPs）。然后，使用基于 Hi-C 的增强剂-促进者相互作用，将 SNPs 映射到其假定目标基因。生成的基因集代表 AD 风险基因，因为它们可能受 AD 风险变异的调控。为了进一步深入了解AD背后的分子机制，我们使用发育性大脑表达数据和大脑单细胞表达图谱来描述AD风险基因。该协议可以扩展到任何GWAS和Hi-C数据集，以识别各种人类特征和疾病背后的假定靶基因和分子机制。

引言

全基因组关联研究（GWAS）在揭示一系列人类特征和疾病的遗传基础方面发挥了关键作用。这种大规模的基因分型已经发现了数千种与表型相关的基因组变异，范围从身高到精神分裂症风险。然而，尽管GWAS在识别疾病和特征相关位点方面取得了巨大成功，但机械学对这些变异如何促成表型的理解一直具有挑战性，因为大多数表型相关变异都存在于非编码中人类基因组的一小部分。由于这些变异通常与预测的调控元素重叠，它们可能会改变附近基因的转录控制。然而，非编码位点可以影响基因的转录，线性距离超过一兆碱基，使受每个变异影响的基因难以识别。三维（3D）染色质结构在调解远距离调控位点和基因启动子之间的连接方面起着重要作用，可用于识别受表型相关单核苷酸多态性（SNPs）影响的基因。

基因调控是由一个复杂的过程，其中涉及增强剂激活和染色质循环形成，物理连接增强剂的基因启动子，转录机制可以定向^1，2，3。由于染色质环通常跨越几百千基（kb），需要3D染色质结构的详细地图来破译基因调控机制。已发明多种染色质构象捕获技术，以识别3D染色质架构^4。在这些技术中，Hi-C 提供了最全面的架构，因为它捕获了全基因组的 3D 染色质相互作用配置文件。Hi-C数据集已经迅速被调整，以解释非编码全基因组显著（GWS）位点5，6，7，8，9，10，11，12，13，因为它可以将非编码变异链接到其假定的目标基因基于染色质相互作用配置文件。

在本文中，我们概述了使用染色质相互作用谱计算预测GWAS风险变异的假定目标基因的协议。我们应用这个协议来映射ADGWS位点¹⁴到他们的目标基因使用Hi-C数据集在成人人脑^9。由此产生的AD风险基因由其他功能基因组数据集为特征，包括单细胞转录组和发育表达谱。

研究方案

1. 工作站设置

1. 安装 R（3.5.0 版）和 RStudio 桌面。打开 RStudio。
2. 通过在 RStudio 中的控制台窗口中键入以下代码，在 R 中安装以下库。
   如果（！BiocManager" % in% 行名（已安装.包（）））
   安装.软件包（"生物管理"，存储库="https://cran.r-project.org"）
   生物管理器::安装（"基因组学"）
   生物经理:安装（"生物母体"）
   生物经理:安装（"WGCNA"）
   安装.软件包（"重塑"）
   安装.软件包（"ggplot2"）
   安装.软件包（"程序"）
   安装.软件包（"gProfileR"）
   安装.软件包（"整洁"）
   安装.软件包（"ggpubr"）
3. 下载文件。
   注: 在此协议中，所有文件都需要下载到 */工作目录。
   1. 通过单击材料表中提供的链接下载以下文件。
      1. 下载用于 AD 的精细映射可信 SNPs（Jansen 等人的第¹⁴号补充表 8 ）。
         注: 分析之前，在 41588_2018_311_MOESM3_ESM.xlsx 中打开工作表 8，删除前三行，并将工作表另存为 Supplementary_Table_8_Jansen.txt，格式为 tab 分隔格式。
      2. 从心理代码下载成人大脑中 10 kb 分辨率的 Hi-C 交互配置文件（如下所述anchored_chromatin_loops.
         注: 此文件具有以下格式:染色体、TSS_start、TSS_end、Enhancer_start和Enhancer_end。在使用其他 Hi-C 数据集的情况下，此协议要求以高分辨率 （5~20 kb） 处理 Hi-C 数据集。
      3. 从 PsychENCODE 下载单个单元格表达式数据集。
         注:这些来自神经典型控制样本。
      4. 从 BrainSpan 下载开发表达式数据集（如下所述为devExpr.rda）。
         注:267666527是一个压缩文件，因此解压缩267666527以提取"columns_metadata.csv"、"expression_matrix.csv"和"rows_metadata.csv"以生成devExpr.rda（见第3节）。
   2. 从Gencode版本19下载外音坐标（参见补充文件，描述为Gencode19_exon.bed和Gencode19_promoter.bed）。
      注: 启动器定义为转录起始站点 （TSS） 上游的 2 kb。这些文件具有以下格式:染色体、开始、结束和基因。
   3. 从生物市场下载基因注释文件（参见补充文件，如下所述为基因Anno.rda）。
      注:此文件可用于匹配基于Ensembl基因ID和HUGO基因命名委员会（HGNC）符号的基因。

2. 生成用于可信 SNPs 的 GRANGES 对象

1. 通过在 RStudio 中的控制台窗口中键入以下代码在 R 中设置。
   库（基因组学）
   选项（字符串 AsFactors = F）
   setwd（"*/工作"） = 这是通往工作目录的路径。
   criSNP = read.delim（"Supplementary_Table_8_Jansen.txt"，标题=T）
   信元SNP = 信元SNP_credSNP$可信.因Causal="是"，]
2. 通过在 RStudio 中的控制台窗口中键入以下代码来创建 GRanges 对象。
   信任 = GRanges（克里德SNP$Chr，IRanges（克里德SNP$bp，credSNP$bp），rsid_credSNP$SNP，P_credSNP$P）
   保存（信任，文件="AD_credibleSNP.rda"）

3. 位置映射

注: 对于每个步骤，在 RStudio 中的控制台窗口中键入相应的代码。

1. 在 R 中设置。
   选项（字符串AsFactors_F）
   库（基因组学）
   负载（"AD_credibleSNP.rda"） = （参见2）
2. 启动子/外置SNPs与基因的位置映射
   1. 加载启动子和外音区域，并生成 GRange 对象。
      外扰子 = 读取表（"Gencode19_exon.bed"）
      外向范围 = GRanges（exon_，1]，IRanges（exon_，2]，exon_，3_），基因[exon]，4]）
      启动器 = 读取表（"Gencode19_promoter.bed"）
      促进者 = GRanges（促进者，1]，IRanges（促进者，2]，促进者[，3]），基因[促进者]，4]）
   2. 将可信的 SN 与外音区域重叠。
      olap = 查找重叠（信任范围、外子范围）
      信元 = 信元[查询Hits（olap）]
      mcols（克瑞克森） = cbind（mcols（克瑞克森），mcols（外兆[主题Hits（olap）]）
   3. 将可信的 SN 与促进者区域重叠。
      olap = 查找重叠（信任范围、促进器）
      信任发起人 = 信任范围[查询Hits（olap）]
      mcols（信任者） = cbind（mcols（信任发起人），mcols（促进者[主题Hits（olap）]）
3. 使用染色质相互作用将SNPs与其假定目标基因联系起来。
   1. 加载 Hi-C 数据集并生成 GRange 对象。
      hic = read.table（"发起人anchored_chromatin_loops.bed"，跳过=1）
      名字（hic）= c（"chr"，"TSS_start"、"TSS_end"、"Enhancer_start"、"Enhancer_end"）
      hicranges = GRanges（hic$chr，IRanges（hic$TSS_start，hic$TSS_end），增强器=hic$Enhancer_start）
      olap = 查找重叠（奇数、启动子）
      hic促进者 = hicranges_查询Hits（olap）]
      mcols（活动发起人） = cbind（mcols（活动发起人），mcols（促进器[主题Hits（olap）]）
      hic增强剂 = GRanges（seqnames（hicpromoter），IRanges（hicpromoter$增强剂，hic促进器$增强剂=10000），基因=hic促进者$基因）
   2. 将可信的 SN 与 Hi-C GRange 对象重叠。
      olap = 查找重叠（信任范围，希克增强器）
      信元 = 信任范围[查询Hits（olap）]
      mcols（信任）= cbind（mcols（信任），mcols（希克增强剂[主题Hits（olap）]）
4. 编译由位置映射和染色质相互作用曲线定义的 AD 候选基因。
   • AD 的候选基因:
   ADgenes = 减少（联合，列表（信不自他基因，信不自生自新$基因，可信基因））
   • 将恩森布尔基因 ID 转换为 HGNC 符号
   负载（"基因Anno.rda"）
   ADhgnc = 基因Anno1_匹配（ADgenes，基因Anno1$ensembl_gene_id），"hgnc_symbol"]
   ADhgnc = ADhgnc[ADhgnc！=""=
   保存（ADgenes、ADhgnc、文件="ADgenes.rda"）
   write.table（ADhgnc，文件="ADgenes.txt"，行.names_F，col.names_F，报价=F，sep="\t"）

4. 发展表达轨迹

注: 对于每个步骤，在 RStudio 中的控制台窗口中键入相应的代码。

1. 在 R 中设置。
   库（重塑）;库（ggplot2）;库（基因组学）;图书馆（生物硕士）
   图书馆（"WGCNA"）
   选项（字符串AsFactors_F）
2. 处理表达式和元数据。
   datExpr = read.csv（"expression_matrix.csv"，标题 = FALSE）
   datExpr = datExpr_，-1]
   达特梅塔 = read.csv（"columns_metadata.csv"）
   datProbes = read.csv（"rows_metadata.csv"）
   datExpr = datExpr_datProbes$ensembl_gene_id！""，""
   datProbes = datProbes_datProbes$ensembl_gene_id！""，""，*
   datExpr.cr= 折叠行（datExpr，行组 = datProbes$ensembl_gene_id，行 ID = 行名（datExpr））
   datExpr = datExpr.cr$datET 折叠
   gename = 数据.帧（datExpr.cr$group2row）
   行名（datExpr） = gename$组
   1. 指定发展阶段。
      达特梅塔 $单位 = "产后"
      idx = grep（"pcw"，达特梅塔美元年龄）
      达特梅塔$单位[idx] = "产前"
      idx = grep（"年"，达特梅塔$年龄）
      达特梅塔$单位[idx] = "产后"
      达特梅塔元单位 = 系数（达美元$单位，水平=c（"产前"，"产后"））
   2. 选择皮质区域。
      达特梅塔$区域 = "子CTX"
      r = c（"A1C"、"STC"、"ITC"、"TCx"、"OFC"、"DFC"、"VFC"、"MFC"、"M1C"、"S1C"、"IPC"、"M1C-S1C"、"PCx"、"V1C"、"Ocx"）
      datMeta$区域_datMeta$structure_acronym %in% r= = "CTX"
      datExpr = datExpr_，其中（达特梅塔$Region="CTX"）*
      达特梅塔 = 达特梅塔[其中（达特梅塔$Region="CTX"），*
      保存（datExpr、datMeta、文件="devExpr.rda"）
3. 提取AD风险基因的发育表达图谱。
   负载（"ADgenes.rda"）
   exprdat = 应用（datExpr_匹配（ADgenes、行名（datExpr））），2，均值，na.rm_T）
   dat = 数据.帧（区域_datMeta$Region，单位_datMeta$单位，Expr_exprdat）
4. 比较AD风险基因的产前和产后表达水平。
   pdf（文件="developmental_expression.pdf"）
   ggplot（dat，aes（x=单位，y=Expr，填充=单位，alpha=单位））=ylab（"标准化表达式"）=geom_boxplot（异常值.大小= NA）= ggtitle（"大脑表达式"） = xlab（""）= scale_alpha_manual（值=c（0.2，1）=theme_classic（）= 主题（图例.)
   dev.off（）

5. 单元格类型的表达式配置文件

注: 对于每个步骤，在 RStudio 中的控制台窗口中键入相应的代码。

1. 在 R 中设置。
   选项（字符串AsFactors_F）
   负载（"ADgenes.rda"）
   负载（"基因Anno.rda"）
   目标名称 = "AD"
   目标基因 = ADhgnc
   单元格exp = read.表（"DER-20_Single_cell_expression_processed_TPM_backup.tsv"，标题=T，填充=T）
   单元格exp[1121，1] = 细胞exp[1120，1]
   细胞外号 = 细胞exp_-1120，*
   行名（单元格） = 单元格exp_，1]
   细胞外号 = 细胞exp_，-1]
   datExpr = 刻度（cellexp，center=T，刻度=F）
   datExpr = datExpr_，789:ncol（datExpr）]
2. 提取AD风险基因的细胞表达谱。
   exprdat = 应用（datExpr_匹配（目标基因、行名（datExpr））），2，均值，na.rm_T）
   dat = 数据.帧（组=目标名称，单元格名称（exprdat），Expr_exprdat）
   dat $celltype = 未列表（l 应用（strsplit（dat$cell，拆分=[.]），''''，1））
   dat = dat_-grep（"Ex_在"，dat$细胞类型），*
   dat $细胞型 = gsub（"Dev"，"胎儿"，dat$细胞型）
   dat$细胞类型 = 因子（dat$细胞型，水平=c（"神经元"，"星形细胞"，"微胶"，"内皮"，
   "奥托根德罗茨"，"OPC"，"胎儿"））
   pdf（文件="singlecell_expression_ADgenes.pdf"）
   ggplot（dat，aes（x=细胞类型，y=Expr，填充=细胞类型）））*
   ylab（"标准化表达式"） = xlab（"） = geom_violin（） = 主题（轴.text.x=element_text（角度 = 90，hjust =1））） = 主题（图例"无"） |
   ggtitle（粘贴0（"AD风险基因的细胞表达特征"））
   dev.off（）

6. AD风险基因的基因注释扩充分析

1. 通过在终端中键入以下命令下载并配置 HOMER。
   姆克迪尔·荷马
   cd 主
   wget http://homer.ucsd.edu/homer/configureHomer.pl
   perl ./configureHomer.pl -安装
   perl ./configureHomer.pl -安装人-p
   perl ./configureHomer.pl -安装人-o
2. 通过在终端中键入以下命令来运行 HOMER。
   导出 PATH=$PATH:*/工作/主站/箱
   findMotifs.pl [/工作/ADgenes.txt 人类 /工作/
3. 通过在 RStudio 中的控制台窗口中键入以下代码来绘制丰富术语。
   库（格普布）
   选项（字符串AsFactors_F）
   pdf（"GO_enrichment.pdf"，宽度=15，高度=8）
   plot_barplot = 函数（名称、名称、颜色）*
   输入 = read.delim（粘贴0（dbname，".txt"），标题=T）
   输入 = 输入，c（-1，-10，-11）*
   输入 = 唯一（输入）
   输入$FDR = p.调整（exp（输入$logP））
   input_sig = 输入 [输入 $FDR < 0.1，]
   input_sig$FDR = -log10（input_sig$FDR）
   input_sig = input_sig[订单（input_sig$FDR），]
   p = ggbarplot（input_sig，x = "术语"，y = "FDR"，填充 = 颜色，颜色 = "白色"， sort.val = "asc"，ylab = 表达式（-log=10]（斜体（FDR）），xlab = 粘贴0（名称，术语），旋转 = TRUE，标签 = 粘贴0（input_sig$Target.Genes.in.term，"/"input_sig$Genes.in.term），字体.label = 列表（颜色 = "白色"，大小 = 9），lab.vjust = 0.5，lab.5
   p = p= geom_hline（yintercept = -log10（0.05），线型 = 2，颜色 ="浅灰色"）
   返回（p）
   }
   p1 = plot_barplot（"biological_process"，"GO生物过程"，"#00AFBB"）
   p2 = plot_barplot（"kegg"，"KEGG"，"#E7B800"）
   p3 = plot_barplot（"反应"，"反应"，"#FC4E07"）
   gg排列（p1，p2，p3，标签 = c（"A"，"B"，"C"），ncol = 2，nrow = 2）
   dev.off（）

结果

此处描述的过程适用于原始研究¹⁴定义的一组800个可信的SNPs。位置图显示，103个SNPs与启动子重叠（43个唯一基因），42个SNPs与外显子重叠（27个唯一基因）。位置映射后，84% （669） SNPs 仍未被指示。在成人大脑中使用Hi-C数据集，我们能够根据物理接近性将另外208个SNPs与64个基因联系起来。我们总共将 284 个 AD 可信 SNPs 映射到 112 个 AD 风险基因（图 1A）。AD风险基因与淀粉样蛋白前体蛋白、淀粉样蛋白-β形成和免疫反应有关，反映了^{AD15、16、17、18}的已知生物学（图1B-D）。AD风险基因的发育表达特征显示产后显著富集，表明与年龄相关的AD风险升高（图2A）。最后，AD风险基因在微胶质、脑原免疫细胞中高度表达（图2B）。这与复发性发现一致，即AD具有很强的免疫基础，微胶质是AD发病机制^14、19、20的中心作用者。

图1:定义ADGWS位点的假定目标基因。（A） 从前 29 个 AD 位点派生的可信 SNPs 被分类为启动子SNPs、外音SNPs和未指定非编码的SNPs。发起人和外音SNPs通过位置映射直接分配给其目标基因，而成人大脑中的色母相互作用曲线则另外用于根据物理相互作用绘制SNPs。（B-D）AD风险基因中的GO（ B）、KEGG（C）和反应素（D）项的富集使用协议第6节所述。x 轴表示已更正的错误发现率 （FDR） -log10（P 值）。绘制了 FDR < 0.1 的丰富术语。灰色垂直线表示 FDR = 0.05。APP淀粉样蛋白前体蛋白。分子，每个术语中表示的AD风险基因的数量;分母，每个术语中的基因数量。请点击此处查看此图的较大版本。

图2:AD风险基因的表征。（A） 与产前皮层相比， AD 风险基因在产后皮层中高度表达.（B） 小提琴图描绘了不同细胞类型中不同细胞类型的基因表达值（规范化表达）的分布。这些结果表明，AD风险基因在微胶质中高度表达，与以前的研究¹⁴一致。请点击此处查看此图的较大版本。

补充文件1。请点击此处查看此文件（右键单击下载）。

补充文件2。请点击此处查看此文件（右键单击下载）。

补充文件3.请点击此处查看此文件（右键单击下载）。

讨论

在这里，我们描述了一个分析框架，它可用于基于位置映射和染色质相互作用对GWS位点进行功能性说明。此过程涉及多个步骤（有关更多详细信息，请参阅此评论¹³）。首先，鉴于染色质相互作用的轮廓是高度细胞类型的特异性，Hi-C数据从适当的细胞/组织类型获得，最好地捕获疾病的基础生物学需要使用。鉴于AD是一种神经退行性疾病，我们使用成人大脑Hi-C数据⁹来对GWS位点进行分文。其次，由于联动不平衡 （LD），每个 GWS 位点通常具有多达数百个与特征相关的 SNPs，因此，通过使用精细映射算法^21、22或使用高通量方法（如大规模并行报告器 （MPRA）²³或自转录活动区域排序）计算预测因果关系，获得假定因果（"可信"）SNPs 非常重要。STARR-seq）²⁴.对于这里描述的工作，我们使用了Jansen等人¹⁴中报道的可信的SNPs。第三，启动子和外音SNPs是基于位置映射进行加法。我们使用一个简单的位置映射策略，其中 SNPs 在与启动器（定义为转录起始位上游的 2 kb）或外源子重叠时映射到基因。但是，可以通过评估外音 SNPs 的功能后果来进一步阐述此方法，例如 SNP 是否导致无意义的中介衰变、误解变异或无意义变异。第四，来自相应组织/细胞类型的染色质相互作用曲线可用于根据物理接近性为其假定目标基因分配SNPs。我们使用锚定到启动器的交互配置文件，但我们可以通过将增强剂活动（由组蛋白 H3 K27 乙酰化或染色质可访问性引导）或外音相互作用考虑进一步优化或扩展交互配置文件。在这个过程中的一个重要考虑因素是使用一致的人类基因组构建。例如，如果汇总统计的基因组位置不基于 hg19（即 hg18 或 hg38），则应获得参考基因组的适当版本，或者需要使用提升²⁵将汇总统计数据转换为 hg19。

我们应用这个框架来识别AD GWAS的假定目标基因，将284个SNPs分配给112个AD风险基因。使用发育表达谱²⁶和细胞类型特定表达谱^9，我们随后证明，该基因集与已知的AD病理学一致，揭示了细胞类型（微胶质）、生物功能（免疫反应和淀粉样蛋白β），以及年龄风险升高。

虽然我们提出了一个框架，其中标定了AD的潜在目标基因及其基础生物学，但值得注意的是，基于Hi-C的注释可以展开，以注释任何非编码变异。随着更多的全基因组测序数据可用，我们对非编码罕见变异的理解也越来越多，Hi-C将为解释疾病相关基因变异提供关键资源。因此，从多种组织和细胞类型中获得的Hi-C资源简编对于促进广泛应用这一框架以收集对各种人类特征和疾病的生物见解至关重要。

材料

Name	Company	Catalog Number	Comments
10 kb resolution Hi-C interaction profiles in the adult brain from psychencode			http://adult.psychencode.org/
Developmental expression datasets			http://www.brainspan.org/
Fine-mapped credible SNPs for AD (Supplementary Table 8 from Jansen et al.14)			https://static-content.springer.com/
HOMER			http://homer.ucsd.edu/
R (version 3.5.0)			https://www.r-project.org/
RStudio Desktop			https://www.rstudio.com/
Single cell expression datasets			http://adult.psychencode.org/