空间表观基因组学揭示了不同级别胰腺神经内分泌肿瘤中与基质微环境相关的调控程序

《iScience》：Spatial epigenomics reveals stromal niche-associated regulatory programs across pancreatic neuroendocrine tumor grades

亮点：
- snATAC-seq定义了NF-PanNETs中与肿瘤分级相关的调控程序
- 免疫细胞和成纤维细胞状态显示出不同的染色质调控程序
- 空间ATAC-seq映射了原发组织中与肿瘤分级相关的状态
- 肿瘤的可及性状态与富含myCAF的空间区域相关

摘要：
非功能性胰腺神经内分泌肿瘤（NF-PanNETs）占胰腺神经内分泌肿瘤的大多数，并表现出显著的临床和生物学异质性，但其表观遗传调控和空间结构仍知之甚少。在这里，我们使用单核转座酶可及染色质测序（snATAC-seq）和空间ATAC-seq对多个肿瘤分级的NF-PanNETs进行了综合分析。snATAC-seq识别出与肿瘤分级相关的染色质可及性程序，这些程序的特点是不同的转录因子活性：低级别肿瘤中具有HNF1A和FOXA2的活性，中级别肿瘤中具有ETS家族因子的活性，高级别肿瘤中具有SNAI相关的调控程序，以及与CEBPB和TOX相关的免疫状态。空间ATAC-seq揭示了肿瘤-基质之间的空间协调相互作用，其中靠近基质的区域富含MYC/FOX相关的增殖程序，而远端区域则富含SNAI相关的侵袭程序。总体而言，我们的研究提供了跨肿瘤分级的NF-PanNETs的空间解析表观基因组图谱，为理解表观遗传异质性以及提出关于调控程序和肿瘤-基质相互作用的假设提供了资源，为未来的功能性和临床研究奠定了基础。

引言：
胰腺神经内分泌肿瘤（PanNETs）是胰腺中第二常见的上皮肿瘤类型。PanNETs被定义为起源于神经内分泌前体细胞或胰腺上皮细胞神经内分泌转分化的异质性恶性肿瘤。根据是否存在临床表现的激素分泌综合征，它们被分为功能性和非功能性（NF-PanNETs）亚型。NF-PanNETs占PanNETs的60%–90%，由于分子成像技术的进步，其发病率和患病率正在增加。尽管某些NF-PanNETs可能具有相对缓慢的进程，但认识到它们的内在恶性潜能及其不同的侵袭能力至关重要。大约40%–80%的NF-PanNETs患者在诊断时已有远处转移，最常见的是肝脏（40%–93%），并且预后明显较差。目前，临床治疗决策主要基于肿瘤分级和分期。系统性疗法，如化疗（例如替莫唑胺）、靶向治疗（例如哺乳动物雷帕霉素靶点[mTOR]抑制剂）、激素治疗（例如生长抑素类似物）和肽受体放射性核素治疗，在无进展生存期方面显示出临床益处；然而，对总体生存期的影响有限。因此，全面了解驱动NF-PanNETs进展的分子机制对于推进临床管理和治疗策略至关重要。

基因组突变谱和NF-PanNETs中的失调信号通路已通过全基因组/外显子测序和全转录组测序得到表征。NF-PanNETs经常携带DAXX/ATRX、MEN1和mTOR通路基因的突变。然而，散发病例也显示出TP53、CDKN2A、RB1和KRAS等基因的改变，突显了NF-PanNETs的分子复杂性和异质性。虽然基因突变在NF-PanNETs的诊断和治疗中起着关键作用，但越来越多的证据表明表观遗传修饰在塑造其复杂的信号通路中的作用。表观基因组学研究可以为正常细胞功能和病理重编程背后的调控状态转换提供关键见解。特别是，利用转座酶可及染色质测序（ATAC-seq）技术使得全基因组范围内的染色质可及性分析成为可能，为绘制肿瘤进展过程中的顺式和反式调控活动提供了强有力的方法。然而，批量ATAC-seq捕获的是平均的染色质可及性谱，掩盖了样本中肿瘤亚群之间的遗传和调控异质性。单细胞/核ATAC-seq（sc/snATAC-seq）可以在单细胞分辨率下提供调控状态的新见解，但在样本处理过程中会失去空间背景。值得注意的是，最近的空间ATAC-seq技术允许在完整的组织结构内进行染色质可及性分析，恢复了空间背景同时保留了高分辨率的表观基因组数据。这使得可以将表观遗传状态映射到特定的解剖位置和微环境龛穴中，从而提高了我们对空间相互作用如何影响细胞身份、可塑性、命运决定和肿瘤-微环境动态的理解。因此，将研究范围从孤立的单细胞扩展到相互作用的多细胞枢纽，可以进行更精细的表型异质性和调控程序分析，这是仅靠sc/snATAC无法实现的。

在这项研究中，我们使用snATAC-seq和空间ATAC-seq进行了综合染色质可及性分析，以绘制未经治疗的NF-PanNETs在不同肿瘤分级中的表观基因组图谱。我们的发现确定了与肿瘤分级相关的转录调控程序和与细胞异质性相关的空间组织化的肿瘤微环境龛穴，为未来研究分子分化和潜在的治疗脆弱性提供了基础。

结果：
我们对九名未经治疗的NF-PanNET患者的原发肿瘤标本进行了Chromium平台（10x Genomics）的snATAC-seq分析，共获得了38,198个核。此外，还对十二个NF-PanNET标本进行了空间ATAC-seq分析，其中六个标本与snATAC-seq分析的标本相匹配（图1A）。snATAC-seq数据通过≥500个独特核片段/核和转录起始位点（TSS）富集得分≥3的阈值进行过滤，以排除低质量核（图S1A）。在所有snATAC-seq样本中，我们识别出62,896个可及染色质区域（ACRs），其中大多数位于启动子区域（46.3%）、内含子区域（29.69%）和远端基因间区域（18.98%）。

图1 NF-PanNET患者细胞的染色质可及性分析
- （A）NF-PanNET队列的概览，显示每个样本的肿瘤分级以及snATAC-seq和空间ATAC-seq数据的可用性。
- （B）snATAC-seq谱型的UMAP投影，显示主要的注释细胞群。右侧的条形图显示每个细胞群中的核数量。
- （C）细胞类型特异性的DACRs热图。颜色表示跨细胞群的平均标准化可及性（z缩放）。
- （D）代表性谱系标记基因的基因活性得分点图。点大小表示具有可检测基因活性的核百分比，颜色表示平均基因活性得分。
- （E）选定谱系标记基因的基因活性得分UMAP投影。
- （F）基因组浏览器覆盖图，显示跨细胞群的代表性标记基因位点的染色质可及性。
- （G）基于chromVAR偏差得分的细胞群间TF基序活性热图。颜色表示z缩放的平均chromVAR偏差得分。
- （H）RFX5在细胞群中的TF足迹分析。

无监督聚类根据染色质可及性谱型识别出了四种主要的细胞群——肿瘤细胞、免疫细胞、成纤维细胞和外分泌细胞（图1B；图S1B）。在这些细胞群中，我们识别出了5,540个细胞类型特异性的差异ACRs（DACRs）。基于邻近簇特异性顺式调控元件的基因注释揭示了肿瘤微环境（TME）的不同细胞成分（图1C）。例如，簇1在肿瘤相关基因（包括KCNH2和KIF1A）附近的可及性增加。KCNH2已被报道在胰腺癌中上调，并与细胞增殖、迁移、侵袭和上皮-间充质转化（EMT）相关，而KIF1A则与神经内分泌分化有关。簇2在免疫相关基因（包括CXCR4、RUNX3和LAPTM5）附近的ACRs富集。簇3在成纤维细胞相关基因（包括COL1A1和MMP19）附近的可及性增加，支持其成纤维细胞身份。接下来，我们生成了一个基因活性得分矩阵，并根据典型谱系标记的活性进一步注释了细胞群（图1D）。肿瘤细胞显示CHGA、CHGB、CACNA2D2和SYT7的基因活性增加，而成纤维细胞则表现为FN1、COL6A1和COL6A2。免疫细胞显示CD86、CD4和IL2RB等标记物的活性增加，而外分泌细胞则以CFTR和TSPAN8为特征。将代表性标记基因活性得分投影到snATAC-seq均匀流形近似和投影（UMAP）上进一步支持了从顺式调控元件可及性推断出的细胞身份（图1E和1F）。

接下来，我们使用从含有转录因子结合基序的可及区域计算出的chromVAR偏差得分来评估与转录因子（TF）相关的调控活性。肿瘤细胞显示RFX家族基序的偏差得分增加，这些基序的成员之前已被认为与肿瘤相关调控过程有关（图1G；图S1E）。与此观察结果一致，TF足迹分析显示肿瘤细胞中RFX结合基序周围的明显可及性谱型（图1H）。免疫细胞群体显示ETS家族基序的偏差得分增加，这与ETS家族转录因子在多种免疫细胞群（包括T细胞、B细胞、巨噬细胞和树突状细胞）中的已知作用一致（图1G）。最后，对峰值基因组注释和15状态ChromHMM注释的分析表明，免疫细胞和成纤维细胞中的可及区域经常与启动子和增强子状态相关，这与这些肿瘤微环境群体内的活跃顺式调控景观一致（图S1C和S1D）。

为了研究NF-PanNETs中不同肿瘤分级间的染色质可及性动态，我们从snATAC-seq数据中分离出了肿瘤细胞。使用潜在语义索引（LSI）和UMAP进行降维和可视化，我们识别出九个肿瘤亚群（Tumor_1至Tumor_9；图2A；图S2A），每个亚群都具有不同的染色质可及性谱型和基因活性特征（图2B和2C）。值得注意的是，这些亚群在不同肿瘤分级之间的分布有显著差异。1级肿瘤主要由Tumor_1至Tumor_4代表，而2级肿瘤则富含Tumor_5至Tumor_9。先前在分化不良的神经内分泌癌中的研究确定了五种由转录因子定义的分子亚型，其特征是ASCL1、NEUROD1、HNF4A和YAP1的表达。因此，我们检查了在我们这个分化良好的NF-PanNET队列中是否可以检测到相关的调控程序，使用核心转录因子的基因活性、亚型相关模块得分和chromVAR基序活性（图S3）。主要与1级肿瘤相关的Tumor_1–Tumor_3显示出NEUROD1和HNF4A相关的特征，而Tumor_2还显示了ASCL1相关的成分。富含2级肿瘤的Tumor_5显示出最明显的POU2F3相关特征。其他亚群显示出较弱或混合的亚型相关模式，表明NEC衍生的分类并不直接对应于这里识别的NF-PanNET状态。某些2级亚群缺乏典型的亚型标记，而是显示出PCSK1、PTEN、TRIM24和TRIM28等基因的活性增加，表明存在超出经典亚型定义的额外调控多样性。接下来，我们进行了通路富集分析，揭示了跨肿瘤亚群的生物程序的差异富集，包括细胞增殖、EMT和免疫相关通路（图S2B）。1级肿瘤亚群Tumor_4富集了与增殖相关的基因集，包括E2F靶点、G2M检查点和MYC靶点。相反，增殖相关程序富集不明显的Tumor_1和Tumor_5则富集了KRAS信号通路和炎症反应通路，这两种通路都与胰腺炎和胰腺上皮内肿瘤（PanIN）有关。跨肿瘤亚群和分级的DACRs进一步强调了NF-PanNETs的表观遗传异质性（图2D；图S2C）。DACRs的基因组注释和ChromHMM状态分析揭示了它们在基因组分布上的分级相关差异，尤其是在启动子和增强子区域（图S2D和S2E）。

图2 不同肿瘤分级中由染色质可及性动态特征的肿瘤异质性
- （A）基于snATAC-seq谱型的注释肿瘤核的UMAP投影。
- （B）代表性标记基因的基因活性得分点图。点大小表示具有可检测基因活性的核百分比；颜色表示平均活性得分。深红色表示较高的基因活性得分，深灰色表示较低的基因活性得分。
- （C）选定亚群相关标记基因的基因活性得分UMAP投影。
- （D）每个肿瘤亚群的主要DACRs。颜色表示跨肿瘤亚群的z缩放平均标准化可及性。
- （E）跨肿瘤亚群的差异chromVAR偏差得分。颜色表示z缩放的平均chromVAR偏差得分。
- （F）分级相关DACRs中的TF基序富集。

接下来，我们使用从可及区域衍生的偏差得分来检查跨肿瘤亚群的TF活性（图2E）。1级肿瘤中的HNF1A、HNF4A和HNF4G显示出较高的chromVAR基序偏差得分，这与与分化神经内分泌状态相关的调控程序一致。相比之下，2级肿瘤显示出ETS2、E2F2和E2F1基序的较高chromVAR偏差得分，这与增强的增殖和细胞周期相关调控程序一致。足迹分析进一步显示1级肿瘤中NFH1A结合基序周围的染色质可及性更高（图S2F）。为了进一步表征这些可及性差异相关的调控特征，我们对DACRs进行了基序富集分析（图2F；图S2G）。1级肿瘤中KLF和EGR家族成员的基序富集，而低增殖2级肿瘤（Grade2_L，Ki-67 3%–10%）中TFEB和YY1基序占主导。相比之下，高增殖2级肿瘤（Grade2_H，Ki-67 10%–20%）中AP-1基序显著富集。鉴于AP-1家族因子在增殖、分化、凋亡、血管生成和侵袭中的重要作用，AP-1基序的富集与Grade2_H肿瘤中独特的增殖和应激响应调控景观一致。确定了在肿瘤组织学级别间差异的TFs后，我们接下来探讨这些调控状态是否可以沿着连续的染色质可及性轨迹进行组织。因此，我们对所有肿瘤细胞进行了轨迹推断。结果表明，1级肿瘤细胞集中在一个端点，而2级肿瘤细胞集中在另一个端点（图3A）。因此，伪时间轨迹捕捉到了与组织学级别相关的染色质调控状态的连续性。

图3 沿着推断的染色质可及性轨迹的级别相关转录因子调控状态

接下来，我们使用基于chromVAR的TF偏差分数来研究沿推断轨迹的TF基序相关活性的变化。通过对平滑后的分数进行无监督聚类，我们识别出四种主要的活性模式（图3B和3C）。模式1-2包括那些在伪时间上TF偏差分数增加的TFs，如ELK4和E2F2，其相关靶基因在细胞周期相关通路、MAPK/ETS信号通路、TNF信号通路和癌症中的转录调控中富集。相比之下，模式3包含那些在伪时间上TF偏差分数减少的TFs，如FOXA2和HHF1A，其相关靶基因在Notch、Hippo、cGMP-PKG、TGF-β和FoxO通路中富集。模式4显示出一个非单调的轮廓，在中间伪时间状态下TF偏差分数降低。相关靶基因在免疫相关通路中富集，包括PD-L1/PD-1信号通路、Th1/Th2分化和IL-17信号通路。具有不同活性模式的代表性TFs进一步沿伪时间进行了可视化，并投影到snATAC-seq UMAP上（图3C和3D）。HHF1A和FOXA2在轨迹的1级相关端点处显示出更高的TF偏差分数，而ELK4、E2F2和SOX2在2级相关端点处显示出更高的偏差分数。这些发现共同识别了NF-PanNET肿瘤状态中的不同级别相关TF调控程序，但并不暗示级别之间的直接时间或因果转变。

为了表征NF-PanNETs中的免疫微环境，我们将免疫细胞分为髓系和淋巴系亚群（图4A和4B；图S4A）。我们根据ITGAM、CD68和CD163的高基因活性分数对髓系细胞进行分类。在淋巴细胞亚群中，识别出了不同的亚群，包括T细胞（CD8A和CD8B）、B细胞（CD19和CD79A）和自然杀伤（NK）细胞（SPON2、IL2RB和KLRD1）。作为癌症免疫监视的关键参与者，NK细胞通过自然细胞毒性和抗体依赖的细胞毒性发挥抗肿瘤作用。对NK细胞相关DACRs的Kyoto Encyclopedia of Genes and Genomes（KEGG）通路富集分析显示，涉及NK细胞介导的细胞毒性、Fc γ R介导的吞噬作用、抗原处理和呈递以及凋亡的通路富集（图4C）。此外，NK细胞在IL2RB位点的启动子区域显示出增强的染色质可及性（图S4B）。相比之下，B细胞主要通过抗原呈递和抗体分泌贡献于适应性免疫。对B细胞相关DACRs的KEGG通路分析表明，与B细胞受体信号通路、NF-κB信号通路和细胞衰老相关的通路富集。

为了表征NF-PanNETs中的免疫微环境，我们将免疫细胞分为髓系和淋巴系亚群（图4A和4B；图S4A）。我们根据ITGAM、CD68和CD163的高基因活性分数对髓系细胞进行分类。在淋巴细胞亚群中，识别出了不同的亚群，包括T细胞（CD8A和CD8B）、B细胞（CD19和CD79A）和自然杀伤（NK）细胞（SPON2、IL2RB和KLRD1）。NK细胞作为癌症免疫监视的关键参与者，通过自然细胞毒性和抗体依赖的细胞毒性发挥抗肿瘤作用。对NK细胞相关DACRs的Kyoto Encyclopedia of Genes and Genomes（KEGG）通路富集分析显示，涉及NK细胞介导的细胞毒性、Fc γ R介导的吞噬作用、抗原处理和呈递以及凋亡的通路富集。此外，NK细胞在IL2RB位点的启动子区域显示出增强的染色质可及性（图S4B）。相比之下，B细胞主要通过抗原呈递和抗体分泌贡献于适应性免疫。对B细胞相关DACRs的KEGG通路分析表明，与B细胞受体信号通路、NF-κB信号通路和细胞衰老相关的通路富集。

与NK和B细胞相比，髓系细胞在多个免疫和信号通路中显示出更广泛的富集。通路分析显示，髓系细胞在趋化因子信号通路、MAPK信号通路、细胞因子-细胞因子受体相互作用、cAMP信号通路、cGMP-PKG信号通路、Rap1信号通路和Ras信号通路中富集（图4C），突显了髓系细胞在肿瘤微环境（TME）中的多样化调控程序。与这些通路级别的差异一致，对细胞类型特异性DACRs的分析显示，几个髓系相关位点的可及性增加（图4D）。值得注意的例子包括参与白三烯B4降解和炎症消退的酶CYP4F3，以及参与囊泡融合和免疫调节的基因SNAP29。在慢性炎症和晚期癌症背景下产生的髓系来源的抑制性细胞（MDSCs）已知会驱动免疫逃逸和肿瘤进展。MDSC相关的分化已在双信号框架内进行描述：第一个涉及STAT3的激活和包括CEBPA、CEBPB和IRF8在内的转录因子的调控；第二个涉及由炎症介质如TNF-α、各种白细胞介素和HMGB1刺激的NF-κB通路。在我们的数据中，我们观察到髓系细胞中C/EBP家族基序活性的增加以及CEBPB、CEBPE和CEBPG的基因活性分数升高，表明髓系亚群中存在类似MDSC的调控特征（图4E和4H）。对ACRs的基因组注释和ChromHMM分割进一步显示，髓系和T细胞相关的DACRs与其他免疫亚群的启动子和增强子状态的注释频率较低（图4F和4G），表明TME中的免疫细胞群体具有不同的调控结构。

为了进一步剖析髓系的异质性，我们将髓系细胞重新聚类为四个亚群（图5A；图S4C）。然后，我们评估了与M1、M2a、M2c和肿瘤相关巨噬细胞（TAMs）相关的标记物的基因活性（图S4D）。有趣的是，这四个巨噬细胞亚群都没有严格符合经典的M1/M2极化，表明存在更复杂或中间的激活状态。这种复杂性也与之前对PanNETs的单细胞研究一致，这些研究发现原发性和转移性病变中存在异质的免疫群体和动态的髓系程序。在我们的数据中，巨噬细胞亚群1在IRF8和ZBTB46上显示出升高的基因活性（图5B），暗示了一种早期的髓系分化或抗原呈递程序。IRF8在单核吞噬细胞分化和TAM抗原呈递中起成熟作用，尽管其功能后果可能因肿瘤环境而异。巨噬细胞亚群3在HMGB1和CCL18上显示出更高的基因活性（图5B），同时伴有炎症和抑制性通路的富集，这与更活跃的、肿瘤相关的状态一致（图5C）。亚群特定的TF程序进一步支持这一解释。巨噬细胞0中SREBF1/2相关基序活性的增加表明脂质代谢的重编程，这是TAM适应的公认特征，而在巨噬细胞3中观察到的HOX和TEAD家族基序活性表明了经典极化类别所不包含的额外调控异质性（图5D）。总的来说，这些发现通过提供髓系在NF-PanNETs中的染色质调控视图，扩展了之前的转录组观察结果，并支持TME中存在多种中间或特化的髓系状态。

为了区分TME中T细胞的异质性，我们选择了显示经典T细胞标记物CD8A和CD8B基因活性的细胞，识别出三个染色质定义的T细胞簇（图4B和5E）。簇0和2主要在2级肿瘤中富集，而簇1是1级肿瘤中的主要亚群（图5E），表明T细胞簇分布与级别相关。为了进一步表征这些亚群，我们检查了T细胞亚群标记物和免疫调控基因的基因活性（图5F）。值得注意的是，簇0在TOX、CD96、CD44和IKZF2上显示出升高的基因活性分数，这与耗竭相关的调控状态一致。簇0的通路富集分析进一步支持这一解释，Hedgehog信号通路显著富集（图S4E），这条通路先前已被认为与免疫抑制和免疫抑制性肿瘤微环境相关。然而，簇0的耗竭相关特征并未完全符合经典的耗竭特征。具体来说，关键抑制性受体如PDCD1（PD-1）和CTLA4并未显示出显著的基因活性。这一观察结果与先前的发现一致，即PD-1/PD-L1在小肠和PanNETs中的表达通常较低或缺失。这些发现表明NF-PanNET微环境中存在非经典的耗竭或免疫抑制相关调控特征。此外，我们观察到2级肿瘤中几个非经典免疫检查点相关基因的基因活性分数升高（图S4F）。

为了表征T细胞的调控异质性，我们比较了三个T细胞簇间的chromVAR衍生的TF偏差分数（图5G）。这些簇显示出不同的TF程序。簇0在ETS家族基序上显示出更高的chromVAR偏差分数，包括ELK1、ETV6和ELK3。簇1在EOMES和TBR1相关基序上显示出更高的chromVAR偏差分数，而簇2则表现出与众不同的NFIX相关特征。接下来，我们通过比较1级和2级T细胞间的基序富集来研究这些差异背后的染色质调控景观（图5H）。在2级肿瘤中，我们观察到TFEB3和CREM基序的显著富集。TFEB3与TFEB协同作用，调节T细胞中的基因表达，尤其是在应激条件下。CREM已与T细胞耐受性和耗竭相关联，这与2级T细胞中的耗竭相关调控特征一致。相比之下，1级肿瘤中的T细胞显示出POU家族基序的增强可及性，这些基序与细胞身份和干细胞样调控程序相关。总体而言，这些结果表明2级肿瘤中的T细胞与更应激响应和耗竭相关的调控特征相关，而1级T细胞则表现出与静止或干细胞样状态更一致的染色质特征。这些级别相关的T细胞调控状态差异可能反映了NF-PanNET级别间不同的免疫微环境背景。

CAF（癌相关成纤维细胞）是TME的主要组成部分，已被牵涉到肿瘤生长、基质重塑和各种实体恶性肿瘤的疾病进展中。无监督聚类识别出五个不同的CAF亚群（簇0-4），每个亚群都表现出独特的染色质可及性景观（图6A和6B）。计算了典型CAF标记物的基因活性分数以对这些簇进行注释（图6C；图S5A）。簇0和1被鉴定为典型的肌成纤维细胞CAF（myCAFs），其特征是ACTA2、TAGLN和FAP标志物的活性分数升高，这些标志物与收缩性和电子显微镜（ECM）沉积相关，符合肌成纤维细胞的特征。值得注意的是，通路富集分析表明簇1和2中的CAF富集了EMT程序（图6G）。第3组显示出了免疫调节表型，具有高活性的CD74和THY1，这表明抗原呈递的间质细胞（apCAFs）可能参与了免疫反应的调节。相比之下，第2组表现出机械生物学特征，表现为DPT和TPM1的共激活，这与基质组织和收缩相关状态一致。第4组代表了一个基质重塑亚组，其特征是MMP11和转录调节因子HOPX的基因活性升高，这表明这种状态可能与侵袭性特征有关。为了理解间质细胞多样性的调控机制，我们评估了各亚组之间的差异性DNA结合序列（DACRs）和chromVAR衍生的转录因子偏差分数（图6D和6E）。每个亚组都展示了一组独特的DACRs，反映了亚组特异性的调控程序。例如，第0组在代谢基因（如PDK4）附近的染色质可及性增加，表明存在代谢重编程的特征。第3组在DLL4（一种与血管生成调控相关的Notch信号配体）附近的染色质是开放的。第4组在与细胞运动和细胞骨架重塑相关的基因附近的区域显示出染色质可及性增加，进一步支持其与运动和细胞骨架重塑相关的状态。

图6 昵质微环境的染色质景观

图6A CAF snATAC-seq谱型的UMAP投影，显示了标注的CAF亚组和肿瘤分级（右）。

图6B 堆叠条形图显示了每个CAF亚组的相对分级组成。

图6C 点状图显示了CAF亚组中代表性标记基因的基因活性分数。点的大小代表具有可检测基因活性的细胞核的百分比，颜色表示平均基因活性分数。深红色表示较高的基因活性分数，而深灰色表示较低的基因活性分数。

图6D 每个CAF亚组的主要DACRs的热图。颜色代表z缩放的平均标准化可及性。

图6E CAF亚组中chromVAR衍生的转录因子偏差分数的热图。颜色代表z缩放的平均chromVAR偏差分数。

图6F 选定的TF偏差分数在CAF中的UMAP投影。

图6G 与DACRs相关的显著标志通路富集，这些通路在每个CAF亚组中显示出增加的可及性。气泡大小代表基因比例，颜色表示调整后的p值。

图6H 分级相关CAF DACRs内的转录因子基序富集。

为了进一步研究CAF内部与分级相关的调控差异，我们对CAF中的分级相关DACRs进行了基序富集分析。这揭示了不同肿瘤分级之间基序相关调控程序的显著差异（图6H）。在2级肿瘤中，我们观察到与EMT和侵袭性程序相关的转录因子（包括SNAI1、ZEB1和TCF家族成员（如TCF4和TCF12）的基序显著富集，这些都与EMT和侵袭有关。这些发现表明2级CAF与富含EMT和侵袭相关转录因子基序的染色质可及性程序相关。相反，1级CAF则富含与应激反应、分化和中度增殖程序相关的基序（包括FOSL2/JUNB（AP-1复合体）、CREB3L1和BHLHE41）（图6H）。总体而言，这些结果突显了CAF内部显著的表观遗传异质性，并确定了与分级相关的EMT和基质重塑相关调控程序的差异。这些发现表明，不同分级之间的间质染色质状态有所不同，可能反映了肿瘤微环境（TME）中更广泛的组织重塑差异。

空间染色质可及性揭示了NF-PanNETs中的分级相关调控状态

空间ATAC-seq能够映射原生组织结构内的染色质可及性，补充了snATAC-seq的见解。我们对12个原发性NF-PanNETs的组织进行了空间ATAC-seq分析（图S6和S7）。H&E染色的组织切片由病理学家进行了标注（图S8）。整合12个NF-PanNET样本的空间ATAC-seq数据后，发现了六个由染色质可及性谱型定义的不同空间簇，以及93,119个DACRs（图7A；图S9A）。1级和2级肿瘤主要由c0-c4簇代表，而c5簇主要与3级肿瘤相关（图7A，红色圆圈和7B）。使用簇特异性顺式元件及其附近的基因来标注这些簇，揭示了3级肿瘤的独特染色质可及性特征（图7C）。将c5簇相关的空间位置映射回组织切片后，发现它们主要位于病理学家标注的肿瘤区域内（图7D；图S8L）。随后，我们专注于3级肿瘤，并确定了三个具有肿瘤相关、免疫相关和成纤维细胞相关染色质可及性谱型的主要空间簇（图S10J）。基序富集分析显示，在样本PNET15中TCF3和TCF12基序的过度表达（图7E）。这些转录因子已被证明与几种实体瘤（包括胰腺癌、胶质母细胞瘤和乳腺癌）中的EMT、侵袭和干性相关调控程序有关。我们进一步计算了基因活性分数，发现YAP1在3级肿瘤中的基因活性分数显著高于1级和2级肿瘤（图7F和7H）。YAP1是Hippo通路的下游效应因子，在组织修复、再生、细胞外基质重塑和肿瘤发生中起作用。先前的研究已将YAP1活性升高与胰腺导管腺癌的侵袭行为和不良临床结果联系起来。我们的发现与这些观察结果一致，并将其扩展到NF-PanNETs，表明3级肿瘤中YAP1相关的基因活性和染色质可及性增加，提示YAP相关调控程序的激活可能是高级疾病的一个特征。浏览器轨迹可视化进一步显示3级肿瘤中YAP1位置的染色质可及性增加（图7G）。综合来看，这些发现表明3级肿瘤表现出一种独特的染色质可及性景观，富含之前与高级或侵袭性肿瘤状态相关的调控特征。

图7 十二个NF-PanNETs的空间染色质可及性景观

图7A 昵质ATAC-seq谱型的UMAP投影，显示了标注的空间簇。一个红色圆圈突出显示了与3级肿瘤相关的c5簇。

图7B 条形图显示了每个空间簇在各个样本中的相对比例。

图7C 与3级空间簇相关的顶级DACRs的热图。颜色代表z缩放的平均标准化可及性。

图7D 样本PNET5中c5簇的空间定位。

图7E 在样本PNET15中鉴定的富集转录因子基序。

图7F 条形图表示空间位点中的平均YAP1基因活性分数。使用双尾Mann-Whitney U检验评估统计显著性；***p < 0.0001。

图7G 基因组浏览器覆盖图显示了不同肿瘤分级中YAP1位置的染色质可及性。

图7H 代表性样本中YAP1基因活性分数的空间分布。

图7I 和 J 根据空间亚组（I）和肿瘤分级（J）对肿瘤相关空间ATAC-seq谱型进行的UMAP投影。

图7K 分级相关肿瘤染色质状态的伪时间轨迹。

图7L HNF1A和HNF4A基序的chromVAR偏差分数的UMAP投影。

图7M PNET1和PNET8中FOXA2基因活性分数的空间图。

图7N PNET1和PNET8中FOXA2 chromVAR偏差分数的空间图。

接下来，我们整合了1级和2级NF-PanNETs的肿瘤相关空间簇，并确定了10个亚组，命名为Tumor_1至Tumor_10（图7I和7J）。1级肿瘤主要由Tumor_1至Tumor_5代表，而2级肿瘤则富含Tumor_6至Tumor_10。1级肿瘤中的HNF1A、HNF4A和HNF4G基序的chromVAR偏差分数较高，这与snATAC-seq分析中的观察结果一致（图7L；图S9B和S9C）。为了检查1级和2级肿瘤相关的空间肿瘤状态是否可以沿着连续的染色质可及性轨迹进行组织，我们使用空间ATAC-seq谱型进行了伪时间推断。1级相关的肿瘤位点位于推断轨迹的一端，而2级相关的肿瘤位点位于另一端，表明存在分级相关的染色质可及性状态的连续性（图7K）。这种模式与snATAC-seq数据集中观察到的分级相关连续性一致（图3A）。与此观察结果一致的是，FOXA2在snATAC-seq轨迹的1级相关端显示出更高的chromVAR偏差分数，并且PNET1（1级）中的FOXA2 chromVAR偏差分数和基因活性分数均高于PNET8（2级）（图7M和7N）。这些发现进一步支持了snATAC-seq和空间ATAC-seq数据集中分级相关的FOXA2调控模式的一致性。

NF-PanNETs中肿瘤可及性状态与myCAF相关区域之间的空间关联

先前的研究报道，肿瘤间质比例较高的患者预后较差，特别是当间质呈现浸润性生长模式时，间质细胞和肿瘤细胞混合且没有明确的边界。在间质成分中，通常以FAP和ACTA2表达为特征的myCAFs被认为位于肿瘤细胞附近，并与EMT、基质重塑和促进肿瘤的间质状态相关。我们在snATAC-seq数据集中识别了myCAFs，并使用空间ATAC-seq研究了肿瘤相关和成纤维细胞相关可及性状态之间的空间关系。经过降维和聚类后，根据标记基因活性和通路富集分析对空间簇进行了标注（图S10A–S10I）。有趣的是，在样本PNET8中识别出两种不同的肿瘤相关可及性状态，分别命名为Tumor1和Tumor2（图8A–8C）。空间ATAC-seq分析显示，Tumor1和Tumor2相关的空间位点都位于与成纤维细胞相关区域附近，其中Tumor2特别接近ACTA2和FAP的高基因活性评分区域，这与它们与myCAF相关区域的空间接近性一致（图8D和8E）。为了量化这一观察结果，我们计算了myCAF相关区域与两种肿瘤相关可及性状态之间的空间距离，发现Tumor2相关的位点在空间上更接近myCAF相关区域（图8F）。通路富集分析进一步显示，Tumor2在增殖信号通路（如E2F靶标、MYC靶标和mTORC1信号通路）以及氧化磷酸化方面富集。相比之下，Tumor1在KRAS信号通路相关的基因方面富集（图8C）。

图8 通过空间染色质可及性揭示的肿瘤内和肿瘤间异质性

图8A PNET8中标注的空间簇的空间分布。

图8B 点状图显示了PNET8空间亚组中代表性标记基因的基因活性分数。点的大小表示具有基因活性的位点的比例；颜色强度反映了平均基因活性分数。深红色表示较高的基因活性分数，深蓝色表示较低的基因活性分数。

图8C 与DACRs相关的富集标志通路，在每个PNET8空间亚组中显示出增加的可及性。

图8D 在PNET8中无监督聚类和Tumor1、Tumor2及成纤维细胞相关空间簇的空间定位。

图8E PNET8中ACTA2和FAP的空间基因活性图。

图8F PNET8中Tumor1-和Tumor2相关位点到成纤维细胞相关位点的最近空间距离分布。使用双尾Mann-Whitney U检验评估统计显著性（***p < 0.0001）。

图8G PNET8中E2F1和E2F3的染色质可及性谱型。

图8H 样本PNET9中E2F3位置的染色质可及性谱型。

图8I 样本PNET9中标注的空间簇的空间分布。

图8J PNET9中ACTA2和FAP的空间基因活性图。

图8K PNET9中Tumor1-和Tumor2相关空间亚组中富集的转录因子基序。

我们将这种空间分析扩展到其他样本（图S11A–S11C）。PNET1、PNET4和PNET6中的肿瘤相关区域位于ACTA2和FAP基因活性较高的空间区域附近，并且富含增殖通路，类似于PNET8中的Tumor2。相比之下，PNET7和PNET10中的肿瘤相关区域与较低的ACTA2和FAP基因活性相关，并且在KRAS通路方面富集。值得注意的是，PNET9中也观察到了肿瘤内的异质性，其中识别出两种不同的肿瘤相关可及性状态（图8I和8J；图S11D）。与PNET8类似，PNET9中的Tumor2与myCAF相关区域在空间上相关，并且富含增殖通路。总之，与myCAF相关区域空间接近的肿瘤相关区域倾向于表现出更强的增殖相关调控程序的富集，而空间上距离较远的肿瘤相关区域则表现出更强的KRAS相关程序的富集。与肿瘤相关区域相邻的成纤维细胞相关区域也表现出EMT和缺氧相关程序的富集。进一步分析显示，Tumor2中的E2F1和E2F3基因活性分数高于Tumor1（图8G）。浏览器轨迹可视化显示Tumor2中的染色质可及性增加（图8H），这与在Tumor2中识别出的增殖相关程序一致。此外，跨样本的基序富集分析显示，FOX基序在高度增殖的肿瘤中富集，这与它们在调节细胞周期中的作用一致（图S11E）。相比之下，具有更强KRAS通路富集的肿瘤相关区域显示了SNAI1、SNAI2和SNAI3基序的富集（图8K；图S11F和S11G）。KRAS激活可以诱导SNAI1和SNAI2的表达，从而抑制上皮基因表达并促进EMT。这些基序模式与之前描述的KRAS相关信号通路、SNAI家族调控程序、EMT和侵袭性肿瘤表型之间的联系一致。这些发现共同强调了肿瘤间和肿瘤内的染色质可及性异质性，以PNET8和PNET9为例，在不同的基质环境中观察到了不同的肿瘤相关可及性状态。与myCAF相关区域邻近的肿瘤相关区域显示出更强烈的增殖相关调控程序富集，而空间上更为远离的区域则显示出更多的KRAS相关程序富集。这些观察结果揭示了肿瘤调控状态与成纤维细胞相关区域之间的空间关联，并强调了空间ATAC-seq在解析NF-PanNET组织内局部调控异质性方面的实用性。

我们的snATAC-seq和空间ATAC-seq的综合分析提供了NF-PanNETs染色质调控景观的高分辨率视图。通过以单核分辨率分析染色质的可及性，并在天然组织结构内进行，我们识别出肿瘤、免疫和基质组分中的分级相关调控状态，并进一步研究了这些状态在不同空间环境中的组织方式。最近对PanNETs的单细胞转录组学研究揭示了与组织学分级、基因组不稳定性和转移状态相关的显著肿瘤和微环境异质性。我们的发现扩展了这一转录组学框架，表明NF-PanNET的异质性也体现在染色质可及性、顺式调控元件的使用和转录因子（TF）相关调控活性中。我们在同一组织学分级中识别出多种染色质定义的肿瘤状态，包括在低级别肿瘤中更为突出的HNF1A/FOXA2相关程序，以及在高级别肿瘤中更为突出的E2F-、AP-1-、SNAI-和YAP/TEAD相关程序，同时还存在异质的髓系细胞、T细胞和成纤维细胞调控状态。我们的研究与以往转录组学研究在细胞分辨率上的差异可能反映了检测方法、队列组成和采样深度的差异。因此，我们的工作为之前通过单细胞转录组学分析描述的细胞异质性提供了补充性的调控和空间视角，而不是对先前定义的细胞分类的直接复制。

在不同肿瘤细胞中，我们观察到了与分级相关的转录因子相关染色质程序的显著差异。1级肿瘤显示出较高的HNF1A、HNF4A和FOXA2等谱系相关因子的活性，而2级肿瘤则显示出更强的增殖和应激反应调控因子的活性，包括E2F和AP-1家族的成员。3级肿瘤还显示出YAP1基因活性的增加以及YAP1位点的可及性增强，同时伴随着之前与EMT、干性和侵袭性肿瘤表型相关的调控程序的富集。这些发现支持NF-PanNET分级之间存在不同的调控状态。然而，重要的是，我们的样本是横断面收集的，因此伪时间分析代表的是染色质状态的推断连续性，而不是从1级到2级或3级疾病的时间顺序转变的直接证据。空间ATAC-seq分析进一步显示，肿瘤相关的染色质状态在不同的基质环境中是有组织的。在代表性的肿瘤如PNET8和PNET9中，靠近成纤维细胞和myCAF相关区域的肿瘤相关区域显示出更强的增殖相关程序富集，包括E2F、MYC和mTORC1信号通路，而空间上更为远离的肿瘤相关区域则显示出更多的KRAS和SNAI相关调控程序富集。靠近肿瘤相关区域的成纤维细胞相关区域也显示出EMT和缺氧相关通路的富集。这些空间关系揭示了肿瘤和基质调控状态之间的局部耦合，但没有确定交互作用的因果方向。特别是，myCAF相关模式可能反映了肿瘤诱导的反应性基质或更广泛的组织重塑，而不是表明成纤维细胞直接诱导了相关的肿瘤调控程序。此外，由于空间ATAC-seq测量代表的是多细胞的空间位置而不是单个细胞，基因活性和通路评分可能反映了单个位点内的细胞组成和细胞状态的变化。

我们的研究为理解NF-PanNETs中的调控异质性提供了染色质水平和空间框架。涉及HNF/FOXA2、E2F/AP-1、YAP/TEAD和SNAI相关程序的反复出现的分级相关模式为未来的机制研究提供了假设，而肿瘤可及性状态与富含成纤维细胞的区域之间的空间关联强调了在解释肿瘤调控状态时考虑局部组织环境的重要性。

**研究的局限性**
应承认本研究的几个局限性。首先，样本量较小，尤其是3级肿瘤的代表性有限，这限制了分级相关模式的普遍性。其次，空间ATAC-seq测量是在位点水平上进行的，由于数据的稀疏性和相关肿瘤及基质状态的相似性，使得定量细胞类型的稳健估计具有挑战性，因此没有进行正式的空间反卷积。因此，空间分析是基于聚类、基因活性、调控特征和空间邻近性，而不是推断的细胞比例。第三，尽管snATAC-seq和空间ATAC-seq的整合识别出了可重复的分级相关和空间有序的调控程序，但这些发现尚未通过多重免疫荧光、空间转录组学或功能扰动等正交方法独立验证。因此，这里识别的候选调控程序应被视为探索性的和假设生成的。未来将染色质可及性与正交的空间和功能检测、更大的独立队列以及纵向临床数据相结合的研究将对确立这些调控状态的生物学、机制学和临床相关性具有重要意义。

**资源可用性**
如需进一步信息、资源或试剂，请联系联系人Yang Liu（[e-mail]）。

**数据及代码可用性**
• 本研究中生成的snATAC-seq和空间ATAC-seq数据集可通过NCBI基因表达组学数据库（GEO）公开获取，访问号为GSE304813，该数据集也列在关键资源表中。详细的质量控制指标、注释的CRE信息、TF活性数据和支持图表的数据来源在Data S1中提供。
• 用于数据预处理、分析和可视化的所有自定义脚本均可在https://github.com/Xiangjun99/NF-PanNETs-ATAC-project免费下载，并在关键资源表中列出。
• 如需重新分析本文报告的数据，可向联系人索取额外信息。

**致谢**
我们感谢耶鲁大学医学院病理学系和神经学系的支持。同时感谢耶鲁病理组织服务。本研究得到了Maximizing Investigators Research Award (MIRA) for Early Stage Investigators R35 GM150838和NIH R01 HL173271（给Y.L.）的支持。此外，本研究还部分得到了美国癌症协会的#IRG-21-132-60-IRG项目的资助。我们还感谢Leslie H. Warner Postdoctoral Fellowship对D.W.的支持。

**作者贡献**
概念制定：Y. Liu., J.D., 和 D.W.；方法学：D.W., X.D., Y. Liu., F.G., G.L., J.D., 和 Y.J.；研究：D.W., X.D., 和 Y. Liu.；可视化：X.D. 和 Y. Liu.；撰写初稿：X.D. 和 D.W.；审稿和编辑：X.D., D.W., 和 Y. Liu.；资源协调：F.G., G.L., L.L., S.H., D.Z., Y. Liang, J.X., Y.J., J.L., P.K., 和 J.K.；监督：J.D. 和 Y. Liu.；资金争取：J.D. 和 Y. Liu.

**利益声明**
作者声明没有竞争利益。

**STAR★ 方法**
**关键资源表**
**试剂或资源**
**来源标识符**
**生物样本**
人类原发性胰腺神经内分泌肿瘤（NF-PanNET）标本
耶鲁大学医学院
Yale HIC编号：1203009817

**化学品、肽和重组蛋白**
pA-Tn5转座酶 – 未加载
C01070010-20
Diagenode
Tagmentation缓冲液（2x）
C01019043-5000
Diagenode
Digitonin
G9441
Promega
5M NaCl
J60434.AK
Thermo Fisher Scientific
paraformaldehyde溶液（16%）
043368.9M
Thermo Fisher Scientific
Glycine
50046
Sigma-Aldrich
Proteinase K
EO0491
Thermo Fisher Scientific
EDTA（0.5 M）
AM9260G
Thermo Fisher Scientific
Fisher Healthcare Tissue-Plus? O.C.T. Compound
23-730-571
Fisher Scientific
Triton X-100
T8787
Sigma-Aldrich
T4 DNA连接酶
M0202L
New England Biolabs
T4 DNA连接酶反应缓冲液
B0202S
New England Biolabs
NEB
Next High-Fidelity 2X PCR Master Mix
M0541L
New England Biolabs
DPBS
14190144
Thermo Fisher Scientific
DNA Clean & Concentrator-5
D4014
Zymo Research
dNTP混合物
R0192
Thermo Fisher Scientific
EvaGreen染料20X（水溶）
31000-T
Biotium
SPRIselect微球
B23318
Beckman coulter
Poly-L-Lysine Slide
63478-AS
Selection microscopy sciences
D5000 ScreenTape
5067-5588
Agilent
D5000试剂
5067-5589
Agilent
Silicon wafer
C04004
WaferPro
Photoresist SU-8 2010
SU-8 2010
Microchem Laboratory
Tween 20? 10%
113324
65001
Sigma-Aldrich
Tris-HCl，1M溶液，pH 8.0
J22638.K2
Thermo Fisher Scientific
UltraPure? DNase/RNase-Free蒸馏水
109770
15
Thermo Fisher Scientific
DAPI
40043
Biotium
引物、连接接头、DNA条形码
IDT
见表S1
Chromium核分离试剂盒

**本研究生成的snATAC-seq数据**
GSE304813

**本研究生成的空间ATAC-seq数据**
GSE304813

**分析代码**
https://github.com/Xiangjun99/NF-PanNETs-ATAC-project

**软件和算法**
Cell Ranger ATAC
10x Genomics
v2.0
RR for Statistical Computing
v4.3.2
Signac
Stuart等人
v1.1
2.0
Seurat
Hao等人
v5.0.1
Harmony
Korsunsky等人
v1.2.0
chromVAR
Schep等人
v1.2.0
chromVAR
motifs
Bioconductor/authors
v0.2.0
clustree
Zappia和Oshlack
v0.5.1
ChIPpeak
Zhu等人
v3.3
6.0
clusterProfiler
Wu等人
v4.10.0
ArchR
Granja等人
v1.0.1
ggplot2
Wickham等人
v3.5.1
MATLAB
MathWorks
R2020b

**实验模型和研究参与者详情**
人类原发性胰腺神经内分泌肿瘤（NF-PanNET）标本来自耶鲁大学医学院接受手术切除的15名未经治疗的患者。样本通过标准病理学验证，并根据组织病理学肿瘤分级进行分类，包括6个1级、8个2级和1个3级肿瘤。9个样本接受了snATAC-seq分析，12个样本接受了空间ATAC-seq分析，其中6个样本使用了这两种技术。参与者的 demographic和临床特征（包括年龄、性别、Hispanic/Latina种族和肿瘤分级）在表S4中提供。所有人类样本的收集均得到了耶鲁大学医学院机构审查委员会（IRB）的知情同意（Yale HIC批准号：1203009817）。组织收集和OCT包埋是在耶鲁大学IRB的批准和Tissue Resource Oversight Committee的监督下进行的。参与的书面知情同意是根据赫尔辛基宣言获得的。每个样本的处理均遵循HIPAA法规、大学研究政策、病理学部门诊断要求和医院章程。新鲜切除的组织首先在OCT中包埋，然后用液氮快冻，再进行后续处理。

**方法细节**
**组织切片制备和染色**
新鲜解剖的人类样本被浸入OCT中，并用液氮快冻。切片前，将冷冻组织块加热到冷冻切片机（?20 °C）。然后将组织块切成7–10 μm厚，并放置在涂有poly-L-lysine的玻璃载玻片（electron microscopy sciences，CatLog编号63478-AS）的中央。同时收集连续的组织切片以进行DBiT-ATAC-seq和H&E染色。相邻的切片经过标准的H&E染色程序。所有染色均在耶鲁大学医学院病理学系和YPTS进行。这些程序遵循临床实验室改进修正案（CLIA）认证的实验室协议以及YPTS的严格标准协议，确保了组织样本分析的精确性和准确性。对于每个样本，加入100微升的转座混合液，其中包含50微升2×标记缓冲液、33微升1× DPBS、1微升10% Tween 20、1微升1% digitonin、5微升转座体以及10微升无核酸酶的水，然后在37°C下孵育30分钟。标记完成后，立即向组织中加入100微升40 mM EDTA以停止转座过程。最后，去除EDTA，并用DPBS洗涤组织三次，再用清水冲洗并风干以备成像。空间条形码在组织处理前几天就制备好了。条形码（100 μM）和连接子（100 μM）在2X退火缓冲液（20 mM Tris-HCl pH 8.0、100 mM NaCl、2 mM EDTA）中以1:1的比例进行退火，退火程序如下：95°C加热5分钟，然后以-0.1°C/秒的速率缓慢冷却至20°C，接着在12°C下保持3分钟。条形码连接分为三个连续步骤：全图像采集、第一次PDMS装置附着、成像、加载条形码A（A1-A50），以及第二次PDMS装置附着、成像、加载条形码B（B1-B50）。为了进行后续的比对和分析，使用EVOS M7000成像系统以10倍物镜扫描组织三次。第一张图像是在放置第一个PDMS装置之前拍摄的，覆盖整个组织区域；第二和第三张图像分别在放置第一个和第二个PDMS装置后拍摄，只覆盖感兴趣区域（ROI）。微流控条形码技术是通过将PDMS装置通道定位在ROI上来实现的，并使用自制的丙烯酸夹具将PDMS装置和组织载片牢固固定。条形码A和B的连接混合液由69.5微升无核酸酶的水、27微升T4 DNA连接酶缓冲液（10X，New England Biolabs）、11微升T4 DNA连接酶（400 U/μL，New England Biolabs）和5.4微升Triton X-100（5%）组成。将含4微升连接混合液和1微升25 μM DNA条形码A（A1-A50）或条形码B（B1-B50）的5微升溶液加入每个通道的入口，并用真空泵吸入不到3分钟。反应在室温下持续30分钟。反应结束后，移除夹具，小心地将PDMS装置剥离，然后将载片在含有50 mL DPBS的Falcon管中洗涤5分钟，风干后进行下一次条形码处理或裂解。第二次条形码处理后，在ROI上放置一个方形孔PDMS垫片。根据垫片的大小，在垫片形成的储液池中加入50–100 μL新配制的裂解液（0.4 mg/mL蛋白酶K、1 mM EDTA、50 mM Tris-HCl pH 8.0、200 mM NaCl、1% SDS）。将载片在湿盒中于55°C孵育2小时。随后将储液池中的裂解物收集到1.5 mL低结合能力的DNA管中，并在55°C下过夜孵育。使用Zymo DNA Clean & Concentrator-5纯化裂解物中的基因组片段，并将其洗脱到22 μL无核酸酶的水中。ATAC文库的构建是通过直接PCR完成的。PCR混合物包括20 μL DNA溶液、25 μL 2X NEBNext Master Mix、2.5 μL 25 μM索引引物i7、2.5 μL 25 μM P5 PCR引物以及2.5 μL 20X EvaGreen，使用qPCR机器按照以下程序进行：72°C加热5分钟，98°C加热30秒，然后是98°C加热10秒、63°C加热30秒、72°C加热1分钟的循环。当信号开始稳定时停止反应。使用SPRIselect beads按1:1的比例纯化PCR产物，并将其洗脱到20 μL无核酸酶的水中。最后，使用Agilent Bioanalyzer High Sensitivity Chips检查文库的质量，并在Illumina NovaSeq 6000测序仪上进行下一代测序（NGS）。

**snATAC数据处理**
单核ATAC-seq（snATAC-seq）数据使用Cell Ranger ATAC（v2.0，10x Genomics，https://support.10xgenomics.com）进行处理，包括读取对齐、去除PCR重复序列、峰值调用以及为每个样本生成基于峰值的片段计数矩阵。在合并数据之前，首先将所有样本的峰值合并为一个统一的峰值集，然后对每个数据集进行量化。对于后续分析，将filtrated_peak_bc_matrix.h5、singlecell.csv和fragments.tsv.gz文件导入R v4.3.2中，并使用Signac v1.12.0和Seurat v5.0.1进行分析。从Cell Ranger ATAC的每个细胞峰值计数矩阵和相应的片段文件中构建ChromatinAssay。基因注释来自EnsDb.Hsapiens.v98注释包。染色体命名规则转换为UCSC风格，并分配到hg38基因组版本。使用Signac计算质量控制指标，排除了含有少于500个独特片段或转录起始位点富集分数低于3的细胞核，从而排除在后续分析中的使用。为了校正样本间的批次效应，应用Harmony算法并使用默认参数。为了降维，使用RunTFIDF进行基于术语频率-逆文档频率的转换。使用FindTopFeatures和min.cutoff = “q0”选择信息丰富的特征，接着使用RunSVD进行潜在语义索引。排除第一个LSI组分，因为它与测序深度相关。LSI维度2-10用于使用FindNeighbors构建最近邻图，并使用RunUMAP生成UMAP嵌入。通过将PDMS装置通道放置在ROI上，并用自制丙烯酸夹具固定PDMS装置和载片来进行条形码A和B的连接。连接混合液的制备使用69.5微升无核酸酶的水、27微升T4 DNA连接酶缓冲液（10X，New England Biolabs）、11微升T4 DNA连接酶（400 U/μL，New England Biolabs）和5.4微升Triton X-100（5%）。将含4 μL连接混合液和1 μL 25 μM DNA条形码A（A1-A50）或条形码B（B1-B50）的5 μL溶液加入每个通道的入口，并用真空泵吸入不到3分钟。反应在室温下的湿盒中持续30分钟。反应结束后，移除夹具，小心地剥离PDMS装置，然后将载片在含有50 mL DPBS的Falcon管中洗涤5分钟，风干后进行下一次条形码处理或裂解。第二次条形码处理后，在两个条形码都流动过的ROI上放置一个方形孔PDMS垫片。根据垫片的大小，在储液池中加入50–100 μL新配制的裂解液（0.4 mg/mL蛋白酶K、1 mM EDTA、50 mM Tris-HCl pH 8.0、200 mM NaCl、1% SDS）。载片在湿盒中于55°C孵育2小时。然后将储液池中的裂解物收集到1.5 mL DNA低结合能力试管中，并在55°C下过夜孵育。使用Zymo DNA Clean & Concentrator-5纯化裂解物中的基因组片段，并将其洗脱到22 μL无核酸酶的水中。ATAC文库的构建是基于直接PCR完成的。PCR混合物包括20 μL DNA溶液、25 μL 2X NEBNext Master Mix、2.5 μL 25 μM索引引物i7、2.5 μL 25 μM P5 PCR引物以及2.5 μL 20X EvaGreen，在qPCR机器中按照以下程序进行：72°C加热5分钟，98°C加热30秒，接着是98°C加热10秒、63°C加热30秒、72°C加热1分钟的循环。当信号开始稳定时停止反应。使用SPRIselect beads按1:1的比例纯化PCR产物，并将其洗脱到20 μL无核酸酶的水中。最后，使用Agilent Bioanalyzer High Sensitivity Chips检查文库的质量，并在Illumina NovaSeq 6000测序仪上进行下一代测序（NGS，双端测序，150碱基对模式）。

**SNATAC数据处理**
单核ATAC-seq（snATAC-seq）数据使用Cell Ranger ATAC（v2.0，10x Genomics，https://support.10xgenomics.com）进行处理，包括读取对齐、去除PCR重复序列、峰值调用以及为每个样本生成基于峰值的片段计数矩阵。在合并数据之前，通过聚合所有样本的峰值构建一个统一的峰值集，然后对每个数据集进行量化。为了后续分析，将filtered_peak_bc_matrix.h5、singlecell.csv和fragments.tsv.gz文件导入R v4.3.2，并使用Signac v1.12.0和Seurat v5.0.1进行分析。从Cell Ranger ATAC的每个细胞峰值计数矩阵和对应的片段文件中构建ChromatinAssay。基因注释来自EnsDb.Hsapiens.v98注释包。染色体命名规则转换为UCSC风格，并分配到hg38基因组版本。使用Signac计算质量控制指标，排除了含有少于500个独特片段或转录起始位点富集分数低于3的细胞核，从而排除在后续分析中的使用。为了校正样本间的批次效应，应用Harmony算法并使用默认参数。为了降维，使用RunTFIDF进行基于术语频率-逆文档频率的转换。使用FindTopFeatures和min.cutoff = “q0”选择信息丰富的特征，接着使用RunSVD进行潜在语义索引。排除第一个LSI组分，因为它与测序深度相关。LSI维度2-10用于使用FindNeighbors构建最近邻图，并使用RunUMAP生成UMAP嵌入。使用Seurat的FindClusters函数和algorithm = 1实现Louvain社区检测算法来识别细胞簇。对于每个主要细胞区室，评估了从0.1到1.0的不同聚类分辨率，并使用clustree（图S12）进行比较。最终选择的分辨率基于簇结构在相邻分辨率间的持久性、UMAP嵌入中的分离度、样本间的一致性、典型谱系标记的可访问性以及避免生物连贯群体的过度片段化。肿瘤细胞的分辨率为0.6，整体免疫细胞区室的分辨率为0.3，T细胞的分辨率为0.4，髓系细胞的分辨率为0.4，癌相关纤维细胞的分辨率为0.4。除非另有说明，否则所有其他参数都设置为它们的包默认值。

**空间ATAC-seq数据预处理**
空间ATAC-seq数据使用先前发布的流程（https://github.com/dyxmvp/Spatial_ATAC-seq）进行预处理。简要地说，使用Linker 1和Linker 2过滤Read 1。过滤后的序列被重新格式化为Cell Ranger ATAC（10x Genomics）输入格式。在新格式中，基因组序列保留在Read 1中，而条形码A和条形码B被合并到Read 2中。然后将得到的FASTQ文件与人类参考基因组（GRCh38）对齐，接着使用Cell Ranger ATAC v2.0去除重复序列并计数片段。为后续分析生成片段文件，每个条目包含基因组坐标和相应的空间条形码信息（条形码A × 条形码B）。

**空间ATAC-seq数据处理和整合**
使用Signac和Seurat分别处理由Cell Ranger ATAC生成的每个样本的空间ATAC-seq计数矩阵。使用Read10X_h5导入原始的基于条形码的峰值矩阵。使用预定义的条形码查找表将空间条形码转换为相应的二维坐标。从样本特定的空间坐标文件确定的组织区域外的条形码在后续分析前被去除。然后使用CreateChromatinAssay和min.cells = 1以及min.features = 1从过滤后的基于条形码的峰值矩阵和相应的fragments.tsv.gz文件创建ChromatinAssay。排除位于非标准染色体上的峰值。从EnsDb.Hsapiens.v98获取基因注释，转换为UCSC染色体命名规则，分配到hg38基因组版本，并添加到染色质分析中。对于每个样本，使用RunTFIDF进行基于术语频率-逆文档频率的转换来标准化染色质可访问性数据。使用FindTopFeatures（min.cutoff = “q0”）保留特征，接着使用RunSVD进行潜在语义索引。由于第一个LSI组分与测序深度相关，因此排除它，LSI维度2-10用于使用RunUMAP和FindNeighbors构建最近邻图和生成UMAP嵌入。使用FindClusters和algorithm = 1的Louvain社区检测算法识别细胞簇。对于每个主要细胞区室，评估了从0.1到1.0的不同聚类分辨率，并使用clustree进行比较（图S12）。最终分辨率的选择基于簇结构在相邻分辨率间的持久性、UMAP嵌入中的分离度、样本间的一致性、典型谱系标记的可访问性以及避免生物学上连贯群体的过度片段化。选择的分辨率分别为肿瘤细胞0.6、整体免疫细胞区室0.3、T细胞0.4、髓系细胞0.4和癌相关纤维细胞0.4。除非另有说明，否则所有其他参数都设置为它们对应的包默认值。

**空间ATAC-seq数据预处理**
空间ATAC-seq数据使用之前发布的流程（https://github.com/dyxmvp/Spatial_ATAC-seq）进行预处理。简要来说，首先使用Linker 1和Linker 2过滤Read 1。然后将过滤后的序列重新格式化为Cell Ranger ATAC（10x Genomics）输入格式。在新格式中，基因组序列保留在Read 1中，而条形码A和条形码B被合并到Read 2中。将得到的FASTQ文件与人类参考基因组（GRCh38）对齐，接着使用Cell Ranger ATAC v2.0去除重复序列并计数片段。为后续分析生成片段文件，每个条目包含基因组坐标和相应的空间条形码信息（条形码A × 条形码B）。

**SNATAC数据后期处理**
使用Signac和Seurat分别处理每个样本由Cell Ranger ATAC生成的基于条形码的计数矩阵。使用Read10X_h5导入原始的基于条形码的峰值矩阵。使用预定义的条形码查找表将空间条形码转换为对应的二维坐标。从样本特定的空间坐标文件确定位于组织区域外的条形码在后续分析前被去除。然后使用CreateChromatinAssay和min.cells = 1以及min.features = 1从过滤后的基于条形码的峰值矩阵和相应的fragments.tsv.gz文件创建ChromatinAssay。排除位于非标准染色体上的峰值。从EnsDb.Hsapiens.v98获取基因注释，转换为UCSC染色体命名规则，分配到hg38基因组版本，并添加到染色质分析中。对于每个样本，使用RunTFIDF进行基于术语频率-逆文档频率的转换来标准化染色质可访问性数据。使用FindTopFeatures（min.cutoff = “q0”）保留特征，接着使用RunSVD进行潜在语义索引。由于第一个LSI组分与测序深度相关，因此排除它，LSI维度2-10用于UMAP可视化和使用RunUMAP和FindNeighbors构建最近邻图。使用FindClusters和smart local moving算法（algorithm = 3）以及resolution = 0.4进行初始的样本内聚类。通过检查典型肿瘤、外分泌细胞、纤维细胞的基因活性分数来分配簇身份。使用GeneActivity生成基因活性矩阵，并作为RNA分析添加到每个Seurat对象中，然后使用LogNormalize进行标准化。随后使用条形码查找表将簇身份映射回它们的空间坐标。

**跨样本分析**
在分配样本和组织学级别元数据后，合并了12个单独处理的space ATAC-seq对象。使用RunTFIDF、FindTopFeatures（min.cutoff = 20）和RunSVD重新处理合并后的峰值矩阵。使用Harmony在LSI嵌入上校正样本间的批量相关变异，以样本身份作为分组变量。使用Harmony校正后的嵌入使用FindNeighbors和维度2-25以及n.trees = 100构建最近邻图。使用FindClusters和algorithm = 1进行Louvain聚类。根据样本组成、肿瘤级别、典型谱系标记的基因活性以及它们在组织切片中的空间分布来评估得到的簇。

**基因活性和转录因子动态分析**
使用Cicero计算基因活性分数，该工具汇总了基因 promoter及其顺式可访问远端调控元件的可访问性信号。使用Signac中的RunChromVAR函数根据chromVARmotifs R包中的位置权重矩阵（PWMs）评估转录因子的活性。使用FindAllMarkers函数（min.pct = 0.1，logfc.threshold = 0.25）识别差异可访问的染色质区域（DACRs）和转录因子基序（DACRs）。

**细胞簇注释**
通过比较已知标记基因的基因活性分数来注释细胞簇。肿瘤细胞由ARX、IRX2、KRT8、MEN1、GCK、CHGA、CHGB、SYP、CPLX2、IQSEC3、KCNH2、PSD、SPTBN4、SYNGR4、SYT7和ITPK1的表达定义。外分泌细胞由CPB1、CFTR和TSPAN8标记。纤维细胞通过FN1、COL6A1、COL6A2、CD9和EMP1识别。免疫细胞类型分别定义为：B细胞（CD83、CD86、CD19、MS4A1、CD79A、CD79B）、T细胞（CD8A、CD8B、CD4）、NK细胞（IL2RB、SPON2、KLRD1、CD160）、常规树突状细胞（cDCs；FLT3、IRF8）和巨噬细胞（CD68、CD163、LAMP2、ITGAM、LGALS3）。CAF亚型定义为：炎症性基因集富集是使用clusterProfiler包中的enricher功能针对Molecular Signatures Database（MSigDB）中整理的50个标志性基因集进行的。空间ATAC数据可视化：首先使用MATLAB 2020b（https://github.com/edicliuyang/Hip-lex_proteome）从明场图像中确定空间像素位置。为了可视化染色质特征（例如，独特的核片段）的空间分布，使用ArchR（v1.0.1）处理片段文件以生成ArchRProject，然后将其导入Seurat（v5.0.1）中进行组织级映射。使用了pt.size.factor参数来调整像素大小以实现最佳可视化效果。对于基因活性得分或转录因子（TF）偏差得分，从Signac对象中提取元数据并相应地空间投影到组织切片上。

量化与统计分析：统计分析使用R（v4.3.2）进行。两个独立组之间的比较采用非配对双尾Mann–Whitney U检验进行。差异性染色质可及性和 motif富集分析分别使用Seurat和Signac中实现的统计程序完成。除非另有说明，p < 0.05被视为具有统计显著性。当统计显著性用星号表示时，****表示p < 0.0001。具体统计测试和绘制的摘要度量的定义在相应的图例中提供。

补充信息（2）：
- 下载全部内容：PDF（10.50 MB）：文档S1，包括图S1–S12和表S1–S4
- 电子表格：（12.76 MB）：数据S1，包含所示图表背后的质量控制指标和数值源数据