Cell | 单细胞表观组和转录组整合分析造血分化过程中的调控机

文献链接：Integrated Single-Cell Analysis Maps the Continuous Regulatory Landscape of Human Hematopoietic Differentiation - ScienceDirect
发表期刊：Cell
影响因子：66.85
发表时间：2018年4月

文章发表时间比较早（scATAC-seq使用的还是低通量的fludigm C1的方法），作者开创性得设计了一系列scATAC-seq数据挖掘的算法、scATAC-seq联合scRNA-seq数据整合的算法，来确定早期人类造血过程中的调控异质性和连续分化轨迹。

关心下游数据挖掘分析的伙伴可以对该文献Method部分仔细研读和参考。

0. 研究背景

造血干细胞分化为探索多能细胞命运分化本质提供优质的模型。造血系统是由少量能够产生各种血细胞的造血干细胞活动维持的，多能干细胞能够通过多级分化获得各类谱系的血细胞。

0.1 实验设计

通过6位健康人供体取骨髓样本，流式分选8种不同的细胞群（跨越髓系、红系、淋巴系），分选策略如图1A-B，单细胞ATAC-seq流程如图1C（Fludigm C1：32张IFCs芯片，共计3072个单细胞CALs）

0.2 数据质控

• 质控标准： 阈值设置--8,268 fragments per cell；6,442 fragments in peaks per cell（图1E）
• 与bulk ATAC-seq数据对比：该方法获得的单细胞染色质可及性图谱与CD34+细胞的bulk ATAC-seq图谱非常相似（图1D）
• 可识别高变异的转录因子基序：ChromVAR识别出高变异的转录因子基序如GATA1、BATF和CEBPB，为已知的造血主调节因子（图1F）

图1. scATAC-seq揭示造血干细胞的染色质可及性：(A) 人类造血功能分化的示意图；(B) 对CD34+细胞的分选策略；(C) 在本研究中使用的单细胞ATAC-seq工作流程；(D) 沿着TET2位点的单细胞表观基因组图谱；(E) scATAC-seq 质控标准：横轴-峰的peak数目，纵轴-峰的fragment百分比，红色虚线表示阈值；(F) 单细胞表观基因组谱的TF基序变异性分析

1. scATAC-seq 数据分析

1.1 人造血干细胞大类分群

• 人造血干细胞分群与转录因子活性推断投影：使用TF Z-score评分对单细胞图谱进行分层聚类，根据其免疫表型定义的细胞类型身份对单细胞进行分类（图2A）；对单细胞图谱进行tsne分群和PCA三维分群，预测红系、髓系、淋巴细和pDC细胞群及分化轨迹（图2B-D）；与主谱系调节因子ID3、CEBPB和GATA1相关的基序在淋巴系、髓系和红系发育中显示出连续的活性梯度，而HSC和LMPP的细胞群显示出与HOX基序相关的更高的可及性（图2E-2H）

图2. 人造血干细胞谱系投影：图A. 单细胞表观基因组图谱：底部色条表示不同免疫细胞身份（通过表面蛋白marker定义），TF motif 可及性分数用Z-score（热图颜色）表示；图B. tsne分群图（分群依据：Z-score，不同cluster根据细胞类型着色，参考图A）；图C. PCA分群图（分群依据：Z-score，不同cluster根据细胞类型着色，参考图A）；图D. 拟时序分析；图E-H. 不同转录因子在各cluster中的可及性映射图

1.2 人造血干细胞亚群细分

作者通过五个主成分应用k-medoids聚类来从头定义了14个造血细胞cluster（图3A-B），这些cluster在很大程度上与基于细胞表面蛋白marker标记定义的造血细胞亚群重叠（图3C），同时与造血相关的TF基序可及性变化（图3D）亦相似。

• 不同祖细胞内存在异质性：不同祖细胞主要由两个或多个不同的cluster组成，如CMP（K2-K5）、MEP（K5-K7）、GMP（K9、K10）和pDCs（K12、K13）
• TF相关变异性分析可推断分化方向：通过TF Z score（高或低）对个体EIPP HSC进行分类，将高/低距离与置换的HSC EIPP曲线进行比较（图3E）；发现CTCF、核因子kB（NF-kB）（由RELA基序代表）和ETS基序在HSC中显著可变但与任何特定的分化方向不相关（图3F）。相反，GATA和MESP/ID基序（由GATA2和MESP1基序表示）TF Z score分别与红细胞和淋巴轨迹显著相关（图3G、3H）

图3. 不同cluster的分子特征：图A. 单细胞表观遗传PCA分群图（14个cluster通过不同颜色区分）；图B. 在PCA空间上描述的数据驱动质子的中心；图C. 不同cluster对应细胞类型的百分比；图D. 不同cluster的可及性marker热图；图E. EIPP HSC 细胞群中共享相似基序的转录因子（EIPP指：表观基因组和免疫表型纯（EIPP）簇）；图F-H. 造血干细胞群在 RELA GATA2 MESP1转录因子基序可及性分析：箭头表示信号偏倚的方向

1.3 拟时序分析HSC发育轨迹中的异质性

造血干细胞发育过程中主调控因子表达和相关染色质变化的时间动态，为进一步的功能研究和分析与分化相关的调控变化提供了资源

• 参考scRNA-seq方法的拟时序分析：红系（K1,K3,K5,K6,K7）、髓系（K1,K2,K8,K14）、pDC（K1,K2,K8,K12,K13）、淋巴系（K1,K2,K9,K10,K11），如图4A-D
• GMP早期和晚期的发育调节因子异质性：GMP的两个cluster（K9和K10）显示髓样发育轨迹上的髓样定义因子SPI1（PU.1）和CEBP相关基序之间的可接近性的显著差异；为了进一步划分这一群体，对来自CD123表达的三个不同区域的细胞进行scATAC-seq、bulk ATAC-seq和bulk RNA-seq（图4E）。bulk ATAC-seq和bulk RNA-seq揭示了GMP-A和GMP-C群体中大量的染色质可接近性和转录组学差异（图4F）：差异表达基因包括重要的发育调节因子，包括HSPC TFs GATA2和TAL1的下调以及GMP-C细胞群中髓样基因SPIB、IRF8、TLR7和MPEG1的上调（图4F）。
• 单细胞表观基因组数据定义细胞群的可行性：来自三个细胞组分的scATAC-seq数据显示髓样分化的早期（GMP-A）和晚期（GMP-C）阶段的强烈分离（图4G，4H）。该方法验证了GMP中的异质性，并且证明了从单细胞表观基因组数据定义细胞群的数据驱动方法的精密度高于免疫表型分类细胞的方法。

图4. 识别连续的分化轨迹：图A-D.（A)红系、(B)淋巴系、(C) pDC和(D)髓系发育轨迹进展：不同cluster用不同颜色着色；图E. 根据CD123表达差异定义不同的GMP祖细胞分类：CD123低-（GMP-A，浅灰色），CD123中-（GMP-B，灰色），CD123高-（GMP-C，深灰色）；图F. GMP-C和GMP-A的差异表达gene火山图；图G. 髓系单细胞轨迹特征及沿轨迹密度；图H. 免疫表型定义细胞类型的髓系进展密度评分

• 沿骨髓发育轨迹聚类TF基序：骨髓发育过程中通过TF Z score评分聚类为6个cluster（图5A）：在造血干细胞中，与调控因子HOXB8和GATA1（cluster 1）相关的TF基序的可及性很高，并通过向CMP的分化而降低
• 主调节因子表达时间动态异质性分析：GATA基序可接近性的丧失（由GATA1基序代表）在HSC区室内开始，而HOX基序可接近性（由HOXB8基序代表）在HSC向CMP分化的转变处丢失，表明GATA基序可接近性的丧失可能是HSC谱系早期事件（图5B）；还观察到骨髓相关TF基序的两种不同的激活模式：cluster 4 TF（CEBPD-和SPIB相关基序）在CMP早期中开始显示活性并增加，cluster 5 TF（STAT1-、IRF8-和BCL11A-相关基序）在GMP-A中活性急剧增加至GMP-C，表示CEBP家族的TF（由CEBPD基序代表）为 形成myeloid-erythroid状态的起始因子（图5C）

2. scATAC-seq与scRNA-seq联合分析

作者进一步对 HSC, CMP 和 GMPs进行单细胞转录组测序（基于10x genomics 的高通量单细胞测序）， 共计捕获7,818 cells (2,268, 4,454, 1,096），大类分群如图 5D

• 单细胞表观组和转录组整合分析定义骨髓发育的不同阶段：作者设计了一套scATAC-seq、scRNA-seq数据整合的算法，使用已经发表数据集验证了其有效性和一致性（图5E），并绘制了已知髓系细胞分化过程中主要调节因子CEBPD、GATA2的动态（图5F-G）。并且筛选转录因子表达&染色质可及性相关性较高的11个转录因子，定义了骨髓发育的不同阶段。

图5. 髓系分化过程中转录因子的动态变化：图A. K-medoids聚类的TF基元可及性（左）和PWM标志（右），用于描绘髓系发育过程中的动态TF基元轮廓；图B-C. HSC活性TFGATA1（蓝色）和HOXB8（绿色），以及单核细胞活性调节因子CEBPD（黄色）和BCL11A（红色）在髓系进展中TF基序可及性Z-score评分的平滑谱。误差条（灰色）表示95%的置信区间；图D. scRNA-seq数据的t-SNE分群图；图E. scATAC-seq和scRNA-seq的髓系拟时间密度评分（数据整合自--已发表文献）；图F-G. CEBPD和GATA2的平均表达谱：单个细胞按照分类着色，CD34+细胞显示为黑色（对应平均基因表达表现为下图红色平滑谱）；图H. 高相关性gene-motif pairs沿髓系拟时间方向动态变化

作者筛选了具有高片段计数和在分化轨迹排序的细胞中具有显著可变性的调控元件，确定了14,005个顺式调控元件

• 有限数量的TF motif可及性模式可以在各个调节元件处诱导染色质可接近性变异：（参考图6A-C）在髓系调节器CEBPD周围的调节元件中，CEBPD-1远端元件“快速激活”，而CEBPD-2“远端元件缓慢激活”（图6B）；推断远端调节元件的动态激活模式与附近表达基因之间的相关性可用于将增强子连接至靶基因（图6C），还发现围绕CEBPD的动态调节元件与CEBPD表达高度相关（图6D）。
• 近端调控元件与附近基因表达的相关性更强：通过计算注释转录起始位点10 Mb内的动态基因的相关性，发现近端调控元件（<100 kb）与附近基因表达的相关性显著高于远端元件（>100kb）（图6E）
• scRNA-seq联合scATAC-seq的统计分析将增强子与靶基因启动子功能连接：作者发现，在启动子捕获HiC（PCHiC）数据中，调控元件与靶基因的相关性随着环置信度的增加而得到改善（图6F）；进一步测试先前定义的顺式连锁表达数量性状基因座（cis-eQTL）是否与使用这些整合的单细胞数据鉴定的增强子-基因相互作用重叠，发现cis-eQTL强烈富集scATAC/scRNA-seq相关峰-基因对（图6G）。因此，单细胞染色质可及性和基因表达之间的统计联系可以作为一种将增强子与靶基因启动子功能连接起来的手段

图6. 调节元件动力学连接远端元件与基因：图A. 每个细胞中CEBPD远端元件的数量（沿着髓系分化拟时序轨迹）：点的颜色表示细胞类型，值被平滑处理（蓝色曲线），误差条（灰色）表示95%的置信区间；图B. 髓系调节因子CEBPD附近的四个调控元件的顺式调控和表达动态；图C. 染色质可及性和RNA表达动态（沿着髓系拟时间轨迹）：根据峰值强度排序；图D. CEBPD基因周围的调控谱：动态增强子用灰色突出显示，蓝色表示peak-gene pairs显著相关，灰色表示peak-gene pairs不显著相关；图E. 基因组距离（相对目标gene）的平均Pearson相关系数；图F. 环置信度分类的平均Pearson相关系数；图G. 顺式eqtl重叠动态增强子上富集峰基因相关或启动子捕获HiC的p值(另请参见Figure S6)

Hi-C技术：Hi-C文库一次可以获取全基因组范围内互作的染色质片段，可以从全基因组的高度来研究染色质的空间结构特征。在Hi-C图谱中，染色质互作频率通过两个bin之间junction reads的数量来表示，和测序深度的概念类似，只有达到一定量的测序深度时才能够认为其代表的染色质互作信息是可靠的。
Capture Hi-C技术：就是在传统Hi-C文库的基础上，新增了一个捕获的过程，捕获目的片段用于后续的测序。Hi-C和Capture Hi-C的关系就好比全基因组测序和全外显子测序，Hi-C可以得到更加全面的信息，但是代价高昂，而Capture Hi-C只针对目标区域进行研究，同样的测序成本可以达到更高的测序深度，信息更加可靠，更加经济适用。

3. 总结与展望

遗传变异与基因表达调控元件的整合分析可能为了解与疾病相关的细胞类型或分化阶段提供新的见解。该文献内建立的这一套方法体系，证明了在单细胞中自然发生的调节异质性的相关性可以用来将调节元件与目标基因配对。这种将调控元件与基因连接起来的单细胞推理方法，可能对于推断稀有细胞或流式细胞仪无法明确定义细胞类型状态下的增强子-基因相互作用特别有用。