编译：刘娟，编辑：景行、江舜尧。
原创微文，欢迎转发转载。
导读
全基因组关联研究(GWAS)发现多个与血细胞性状相关的基因组位点，这些基因位点的生物学相关性有待研究。研究者整合中性粒细胞基因表达和全转录组关联研究(TWAS)(N=196)以及中性粒细胞计数的全基因组关联研究(N=173480)，共发现174个全转录组关联研究显著基因富集于调控中性粒细胞的主控转录因子靶基因。CRISPR/Cas9技术敲除CD34+造血和祖细胞(HSPCs)中位于5q13.2，TAF9染色体的全转录组关联研究候选基因，体外试验发现其显著影响中性粒细胞。研究者还发现89个仅对剪接结有意义的独特基因，凸显超出基因表达的替代剪接在粒系形成中的重要性。本研究结果突出基因编辑后全转录组关联研究的优势以确定全基因组关联研究位点与造血相关的功能。
原名：Functional annotation of genetic associations by transcriptome-wide association analysis provides insights into neutrophil development regulation
译名：全转录组关联分析发现中性粒细胞发育调控新机制
期刊：Communications Biology
IF：4.165
发表时间：2020年12月
通讯作者：饶书权
通讯作者单位：哈佛医学院
DOI号：10.1038/s42003-020-01527-7
1 全转录组关联分析发现与中性粒细胞计数相关的易感基因
研究者为找到与中性粒细胞计数相关的基因，使用CD16+中性粒细胞基因表达参考板(N=196)和全基因组关联分析汇总统计数据分析欧洲血统（英国Biobank和INTERVAL研究)的173480名参与者样本。使用FUSION软件进行基于汇总的全转录组关联分析。FUSION整合表达参考板(单核苷酸多态性-表达关联)、全基因组关联分析汇总统计(单核苷酸多态性-中性粒细胞计数关联)和连锁不平衡(LD)参考板(单核苷酸多态性-单核苷酸多态性关联)评估表达顺式遗传成分与表型(基因表达-中性粒细胞计数)之间的关联。在实际应用中，以表达参照板作为连锁不平衡参照板，利用1Mb侧翼窗内的局部单核苷酸多态性，采用稀疏混合线性模型估计每个基因的顺式单核苷酸多态性表达效应大小。全转录组关联分析共发现174个转录组全基因显著的中性粒细胞计数关联(P<0.05/6430)(图1)。在所有全转录组关联分析中，174个关联基因中只有26个是全基因组关联分析风险位点中最接近单核苷酸多态性的基因。多个全转录组关联分析在骨髓造血或中性粒细胞发育中表现出显著作用，包括IKZF1、HNRNPK和NCOR1等。为找到未被全基因组关联分析单独识别的调控中性粒细胞发育的新调控区域，研究者接下来检验174个全转录组关联分析与全基因组显著的中性粒细胞计数全转录组关联分析位点间的重叠。值得注意的是，174个全转录组关联分析中的56个位于34个独立的1Mb区域，其中全转录组关联分析位点的全基因组关联分析关联统计低于全基因组意义(lead SNPGWAS P>8.31E-09)(表1)。这些发现可能是由于全基因组关联分析通过单个基因对中性粒细胞分化产生部分独立的影响，也可能是由于全转录组关联分析(N=6430个基因模型)相对于全基因组关联分析(2950万个推测变异体)的检测负担有所减轻。这些证据表明全转录组关联分析hits可能代表中性粒细胞发育的因果基因，凸显来自全转录组关联分析中性状相关组织eQTL重要性。图1.中性粒细胞计数的全转录组关联研究曼哈顿图(N=173480)。NEUT#表示中性粒细胞计数。 表1. 33个独立转录组全显性基因位点距任何GWAS发现的变异至少500 kb。2 全转录组关联分析中受中性粒细胞谱系特异性主控转录调节因子(TFs)的调控中性粒细胞谱系由复杂TF调控网络相互作用决定的，涉及PU.1、C/EBPα、C/EBPε、RUNX1、Fli-1原癌基因和FLI1等。为验证中性粒细胞计数全转录组关联分析hits基因是否受中性粒细胞谱系决定TFs的调控，参考转录调节因子(TRs)检测工具Lisa（该工具旨在结合人和小鼠的表观遗传学数据来识别调控查询基因集的TRs）。在1316个TRs中，包括CEBPA、FLI1、MED12和RUNX1(所有P<7.01E-08)，多个已知参与中性粒细胞规范和扩张的TF被优先考虑(图2a)。这些TF在不同血细胞类型中的相对表达量见图2b。为进一步了解全转录组关联分析hit的生物学相关通路，研究者进行KEGG通路和GO分析。尽管没有通路富集的证据，但中性粒细胞计数的全转录组关联分析hits在多个GO术语中均显著富集，其中最显著的是选择性剪接(P=2.72E-05，FDR=0.04) (图2c)，这与先前发现的选择性剪接在血液系统谱系中起关键作用相一致。图2.全转录组关联分析hits的功能分析揭示原理。a利用Lisa(Silico缺失分析的表观遗传景观和MARGE的第二代后代)系统识别全转录组关联分析hits的转录调控因子(TRs)。b 14种人造血细胞TRs基因表达的箱图。Y轴表示每个样本中基因的log2TPM值(每千碱基百万转录本)。c GO富集分析全转录组关联分析hits(中性粒细胞计数)。d 15种人类造血细胞类型的全转录组关联分析hit(中性粒细胞计数)的热图。行表示全转录组关联分析命中率;列表示细胞类型。从BloodSpot网站(http://servers.binf.ku)下载各基因归一化表达数据。dk / bloodspot HemaExplorer数据集)。HSC代表造血干细胞;MPP代表多能祖细胞;CMP代表普通骨髓祖细胞;GMP代表粒细胞单核细胞祖细胞;CLP代表共同淋巴祖细胞;MEP代表巨核-红细胞祖细胞;EB代表成红血球细胞;MK代表巨核细胞;HPC代表造血祖细胞;mDC代表髓样树突状细胞;pDC代表树突状细胞。3 造血过程中全转录组关联分析的阶段特异性表达模式为更好地了解全转录组关联分析靶基因在造血系统尤其是中性粒细胞谱系中的潜在功能，利用(http://servers.binf.ku.dk/bloodspot/)研究全转录组关联分析靶基因在7种中性粒细胞谱系相关细胞细胞类型全转录组关联分析靶基因的表达动态。在174个基因中，可用表达数据的基因共119个，层次聚类分析将其分为6个亚组(图2d)。大多数全转录组关联分析靶基因(79/119)在HSPC(如U2AF1和SOX7)、共同髓系祖细胞(CMP)/GMP(如PSMD3和TIMM50 )或两者(如TAF9和ZBTB16)中均有高表达，表明这些基因参与早期髓系细胞分化和/或扩张。研究者还发现多个在髓系或中性粒细胞阶段特异高表达的基因，如IFITM3(图2d)，这些基因也可能通过调节中性粒细胞的存活来解释中性粒细胞计数的变化。此外，多个基因在中性粒细胞系相关和非相关细胞类型中特异性高表达，如在浆细胞样DC细胞中的SEC16A和SH3D19，这可能与这些基因的多方面功能或中性粒细胞与其他细胞类型之间共同的遗传结构有关。4 spTWAS优先排序数十个独立于全转录组关联分析的因果基因为确定中性粒细胞发育所必需的其他剪接事件，研究者接下来进行拼接全转录组关联分析(spTWAS)，并在Bonferroni校正后的165个独特基因中发现825个spTWAS与中性粒细胞计数的关联 (图3a)。超过半数的基因(89/165)没有表现出明显的基因水平的全转录组关联，因此代表独立的hits(图3b)。几个spTWAS hits参与髓系发育，包括CSF3R、NLRP3和SPI1等。研究者发现突出的CSF3R，一种控制粒细胞产生、分化和功能的细胞因子，CSF3R基因存在多个spTWAS关联，其中chr1p34.3位点(rs3917932 lead SNPGWAS P=2.06E-39)(图3c)，共有7个剪接QTL(spQTL)在不同位置调控CSF3R选择性剪接，其中3个位于CSF3R的内含子或外显子(图3d，e)，代表潜在的剪接顺式调控变异。值得注意的是，CSF3R的整体表达没有明显的eQTLs (图3f )，表明选择性剪接可以独立地占一个风险位点。进一步通过KEGG通路和GO分析发现，中性粒细胞计数spTWAS位点在磷蛋白(P=2.70E-08)、选择性剪接(P=2.34E-07)和天然免疫(P=1.11E-04)中显著富集，这些数据与前人提出的spQTL解释复杂性状遗传力相当部分的观点一致，甚至可能比eQTL更多。图3.结合点spTWAS优先排序数十个独立于全转录组关联分析关联的因果基因。a 中性粒细胞计数结合点spTWAS的曼哈顿图。b 全转录组关联分析和结合点spTWAS的中性粒细胞计数维恩图。c 中性粒细胞计数全基因组关联分析区域曼哈顿图,在染色体1p34.3上CSF3R基因位点周围150 kb内。通过质量控制措施的基因型和筛选单核苷酸多态性根据其染色体位置绘制P值。每个圆代表一个单核苷酸多态性。箭头所示为GWAS-lead SNP rs3917932和junction spTWAS-lead spQTL。每个圆的颜色表示rs12187268配对r2值的范围。基因位置和转录方向注释在下方。d CSF3R基因位点剪接连接和spQTL位置。e代表性spQTL基因型(左为rs3917925，右为rs3917981)在CD16+中性粒细胞中划分归一化连接的盒图。f染色体1p34.3 CSF3R单核苷酸多态性的区域顺式eQTL分析。图中显示x轴上单核苷酸多态性的位置(染色体长度，Mb)和y轴上单核苷酸多态性基因表达关联的−log10(P值)。全转录组关联分析发现每个位点的多个hit基因(图1)，原因是多个基因在同一位点的表达相互关联，这与全基因组关联分析中连锁不平衡对单核苷酸多态性水平关联的影响相当。研究者将全转录组关联分析的靶基因分组在1个兆碱基内，发现一些位点(44/87)有单个靶基因，另一些位点有2-6个候选基因。为验证这些信号是由于多个相关特征还是条件独立，研究者进行条件和联合分析。例如，ZNF76解释了第6染色体位点的变异(rs914547 lead SNPGWAS P=1.42E-08，ZNF76 lead SNPGWAS P=1)，同时发现SMIM8解释了0.908全基因组关联分析信号( rs7774080 lead SNPGWAS P=8.29E–17，SMIM8 lead SNPGWAS P=1.26E–02)。条件分析和联合分析在一定程度上解决基因优先排序问题，但基因水平关联的精细映射仍具有挑战性，因为(1)建模的差异和偏差，预测的表达只不过是不完美地捕获顺式表达，(2)因果基因可能在具有广泛LD的eQTL区域无法进行测试。6 实验验证支持TAF9作为染色体5q13.2的因果基因根据全转录组关联分析发现非因果基因的可能性，实验验证对于候选全转录组关联分析建立因果关系至关重要。5q13.2染色体上的全基因组关联分析位点发现多个全转录组关联分析hits，条件分析和联合分析确定因果基因。染色体5q13.2上的单核苷酸多态性rs11745591与中性粒细胞计数和中性粒细胞百分数在同一方向上都有显著的相关性。多个单核苷酸多态性中，单核苷酸多态性 rs11745591位点表现出较强的LD (全基因组关联分析P=5.89E-15)。该位点存在6个全转录组关联分析关联，包括TATA-box结合蛋白关联因子9 (TAF9，全转录组关联分析 P=1.48E-12)、腺苷酸激酶6(AK6，全转录组关联分析 P=1.44E-12 )和包含125(CCDC125，全转录组关联分析 P=1.80E-12)的卷曲结构域，以及全转录组关联分析P值较低但显著的GTF2H2C、OCLN和CDK7。接下来，我们关心是否有一个全转录组关联分析hit基因，尽管是非因果的，但由于与该区域的另一个因果基因相关，所以优先排序。有趣的是，TAF9和AK6无论是总的基因表达量(corr.=0.98)，还是预测的基因表达量(corr.=1)，都显示出完美的基因表达。在5q13.2染色体上，与CCDC125、GTF2H2C、ocn和CDK7相比，TAF9和AK6的SNPs对表达的影响具有高度一致性(图4c)。值得注意的是,该区域的多个单核苷酸多态性对靶基因的表达具有显著的顺式效应,可能是由于强LD与因果单核苷酸多态性有关。进一步研究发现TAF9和AK6尽管翻译自不同的阅读框架，编码不相关蛋白质,但可能具有相同的调控区域，从而解释它们之间表达相关性的高概率。为确定TAF9(或AK6)或该位点的任何其他基因是否导致这种因果关系，研究者进行体外功能分析。设计针对6个基因的单导RNA (sgRNA)，通过3xNLS-SpCas9:sgRNA RNP电穿孔编辑人CD34+HSPCs，然后用G-CSF培养中性粒细胞分化(图5a)。由于不同来源的CD34+HSPCs具有高度可变的固有基因表达，因此选择具有预测靶基因相对高表达基因型的CD34+HSPCs进行功能检测(TAF9和AK6的rs12187268，CCDC125的rs12513817, GTF2H2C的rs299085, OCLN的rs6862253, CDK7的rs6879078)(图5b),以最大化目标基因敲除对中性粒细胞成熟的影响。所有的sgRNAs都产生高效的基因编辑，87.0%－96.3%的片段插入(图5c)，它们可能破坏蛋白质编码序列或触发无义介导的mRNA衰变，从而敲除靶基因。由于CDK7是一种细胞必需基因，对人类细胞适应度必不可少，与对照组相比，CDK7的敲除导致细胞数量减少，因此排除CDK7进行进一步分析。除TAF9外，与中性位点靶向相比，编辑所有sgRNAs的mRNA水平均显著降低(图5d)。考虑到TAF9只有一个外显子，并且片段插入诱导的编码区框移不能触发无义介导的衰变，因此我们检测TAF9的蛋白水平。与预期的一样，靶向sgRNA的TAF9蛋白水平显著低于对照组(图5d)。TAF9靶向性,但不是任何其他候选基因的靶向性,导致中性粒细胞在体外成熟显著增加(图5e,f),与全基因组关联分析结果一致的是，rs12187268等位基因预测的TAF9表达量较低(β=−0.77,P = 2.51E−18，图4c)，而中性粒细胞计数较高(effect size=0.025,P=3.18E−12)。图4. 在染色体5q13.2处，全转录组关联分析发现中性粒细胞计数的多个位点。a 中性粒细胞计数与全转录组关联分析hits的区域相关性。每个图的顶部显示该区域的所有基因。边缘相关的全转录组关联分析基因用蓝色表示，共同显著的基因用绿色表示。底部的板显示在绿色基因的预测表达之前(灰色)和之后(蓝色)对每个SNP的全基因组关联分析曼哈顿区域图。全基因关联分析统计数据计算双边p值。b在染色体5q13.2处预测全转录组关联分析基因表达相关性(corr.)矩阵。CD16+中性粒细胞中按代表性eQTL基因型分层的归一化基因表达盒图。利用线性回归模型和加性基因型效应发现eQLTs。NEUT#代表中性粒细胞计数；Chr代表染色体。图5.在染色体5q13.2敲除TAF9，而非其他TWAS敲除导致体外中性粒细胞扩张增加。a体外实验验证示意图。通过SpCas9 RNP电穿孔(EP)靶向中性位点对人CD34+HSPCs进行编辑，在染色体5q13.2(TAF9、AK6、CCDC125、GTF2H2C、OCLN和CDK)中分别进行TWAS hits。在体外向中性粒细胞分化10天后，细胞培养具有免疫表型特征。b基因组DNA Sanger测序显示eQTLs基因型(箭头表示)与染色体5q13.2处的TWAS匹配。c 3xNLS-SpCas9后HSPC的编辑效率:用sgRNA进行sgRNA电穿孔。用Sanger测序法测定电穿孔后6天的基因编辑情况。中性粒细胞分化4天后RNP编辑CD34+HSPCs中TWAS表达水平(mRNA或蛋白)。RT-qPCR检测mRNA水平，western blot检测蛋白水平。e具有代表性的RNP编辑CD34+HSPCs细胞培养(第12天)的流式细胞术。CD34-CD117-CD16+细胞群，中性粒细胞。f电穿孔后12天CD34-CD117-CD16+细胞数量的扩增，用细胞总扩增乘以流式细胞术测定的分数计算。在造血系统中，GWAS研究发现数千种与各种血细胞特征相关的变异。后GWAS时代的挑战已从识别与血细胞性状相关基因区域转向寻找驱动每个信号的准确调控变体和靶基因以及它们作用的细胞类型，以最终发现健康和疾病中造血调控的机制。个别位点的研究已经明确造血系统的重要调控因子，如胎儿血红蛋白表达的关键调控因子BCL11A，但相关遗传变异与潜在表型的靶基因连接的低通量仍然是血液生物学和临床实践的难题。为解决这一问题，使用CD16+中性粒细胞的参考表达板对中性粒细胞计数开展全转录组关联分析，发现174个候选致病基因，这些基因来自中性粒细胞计数GWAS位点和34个GWAS中未报道的额外风险位点，凸显TWAS对候选致病基因优先排序的优势。多个TFs被确定为骨髓生成的必要调节因子，包括PU.1、RUNX1、C/ EBPs、IRF8、LEF1等。考虑到(1)髓系谱系中关键TFs失调可能导致骨髓形成的严重缺陷, (2)GWAS发现常见变体往往具有较小的个体效应,全转录组关联分析中没有发现这些TFs。目前全转录组关联分析采用的基因表达模型是基于顺式eQTLs而非反式eQTLs，这也可能导致全转录组关联分析相关性缺失。Võsa和同事发现转基因eQTLs通常会聚集在已知疾病病因学中起核心作用的转基因基因上，并且可能具有更多的信息。然而，根据我们的无偏倚TF预测模型对于几个关键TF所预测的基因，全转录组关联分析结果更加丰富。C/EBPα通过启动和激活HSPCs中的髓系相关基因来指导髓系分化，并与髓系红系祖细胞室中的其他TFs竞争基因组占用以促进中性粒细胞分化。C/EBPα的表达是由造血系统特异性增强因子决定的，该因子受多种TFs调控，包括FLI1和RUNX1。大型中介复合物的调节成分MED12，使一般RNA聚合酶II转录机制和增强子结合调节因子之间的接触，SMAD1对早期骨髓生成至关重要，这在模式生物中得到证实。总的来说，这些数据进一步强调这些谱系限制TFs在骨髓生成中的主导调控作用，为研究TFs在骨髓生成中的作用机制提供线索。与Gsev等的研究一致, 全转录组关联分析频繁发现多个基因/位点,含有2-6基因/位点的43个位点,可能由于(1)总表达相关,(2)相关的预测表达式在因果和非因果之间的基因,(3)共享全基因组关联分析hits, 和(4) 一个GWAS位点的多个因果变异。一对基因在顺式中协同调控越紧密，全基因组关联分析hits之间的LD越强，如AK6和TAF9，就越难以根据全基因组关联分析和表达数据区分因果关系。此外，不同基因组之间的序列相似性可能导致短读序列比对错误，导致共表达分析的假阳性。TAF9和AK6之间共享观察序列。排除共有外显子序列后，TAF9和AK6的总基因表达存在显著相关性(相关系数=0.65,P=2.2E−16),但其他基因对的表达相关性可能存在假阳性。预测表达相关性往往高于总表达相关性，即使总表达相关性较低，也可能发现非因果基因。CRISPR/Cas9基因编辑技术和体外中性粒细胞分化系统发现TAF9是唯一驱动染色体5q13.2全基因组关联分析信号的全转录组关联分析hits。TAF9的敲除导致中性粒细胞扩张倍数显著增加。TAF9编码一种TATA盒结合蛋白(TBP)相关因子，作为特定启动子所需的转录辅激活因子/适配器。虽然尚无直接证据表明TAF9与中性粒细胞发育有关，但之前的研究揭示基础生物学过程如何在造血过程中具有不同和意想不到的作用，如Diamond-Blackfan贫血(DBA)的核糖体蛋白。尽管CMP/GMP阶段观察到TAF9的高表达，研究者无法使用体外中性粒细胞分化系统来分析受影响的细胞类型，因为祖细胞只占整个细胞培养和GMP的一小部分并迅速分化进入下一阶段直到中性粒细胞成熟。TAF9敲除小鼠或人CD34+HSPCs异种移植小鼠模型有助于TAF9功能的进一步研究。尽管本研究具有显著优势，但由于其设计，仍存在一些固有的缺陷。首先，利用SNPs加权线性组合的基因表达预测模型可能会混淆全转录组关联相关性。因为包含SNPs而没有对附近的基因表达产生调节作用，统计数据被夸大。TAF9基因位点的总表达量与期望表达量之间的基因表达相关矩阵存在差异。第二，表达参考板(N=6430个基因模型)是基于相对较小的样本量。因此，可能有一些对中性粒细胞发育重要的特征是我们目前无法使用全转录组关联分析捕捉到的。第三，部分全基因组关联分析hits可能不会通过eQTLs影响靶基因的表达，而是在基因位点上影响蛋白质结构或选择性剪接，这样靶基因就不会被全转录组关联分析捕获，这部分解释为什么某些全基因组关联分析位点没有候选致病基因。值得注意的是，当弱eQTLs被归为全基因组关联分析位点时，尽管与性状相关，全基因组关联分析可能无法识别因果基因。第四，只纳入分析欧洲血统样本。用于功能检测的CD34+HSPCs来自具有特定顺式eQTL基因型(rs12187268 TT预测较TC或CC更高的TAF9表达)的欧洲血统样本，目的是最大化靶基因敲除对中性粒细胞发育的影响。这一结果能否推广到其他种族，还有待进一步研究。综上所述，全转录组关联分析是一种强有力的统计方法，可以识别中性粒细胞分化的大、小效应基因，对这些基因的进一步研究将为中性粒细胞发育的生物学和遗传学提供更多的见解。
更多推荐
1 科研 | Science：人类血细胞中蛋白质编码基因的全基因组转录组学分析
