基于scRNA-seq的GRN分析三阴性乳腺癌的肿瘤异质性

三阴性乳腺癌即TNBC是一种肿瘤异质性高的乳腺癌亚型。最近的研究表明，TNBC患者可能包含具有不同分子亚型的细胞。此外，基于scRNA-seq数据构建的GRN已经证明了对关键调控因子研究的重要性。作者使用scRNA-seq对TNBC患者的GRN进行了全面分析。从scRNA-seq数据鉴定出545个恶性细胞。根据PAM50模型对恶性细胞进行亚型划分。细胞间通讯分析显示巨噬细胞在肿瘤微环境中起主导作用。接下来，通过整合基因共表达和富集结合基序构建各个亚型的GRN。然后，根据基因的中心性度量来确定关键基因。重要的是，关键基因ETV6在所有亚型中都普遍上调，但它通过调节不同的靶基因在每个亚型中发挥不同的作用。综上所述，基于scRNA-seq数据构建GRN有助于剖析TNBC的肿瘤内异质性，并鉴定TNBC的关键基因。

来自：Single-cell RNA-seq dissects the intratumoral heterogeneity of triple-negative breast cancer based on gene regulatory networks

数据收集和预处理

从GEO（https:// www.ncbi.nlm.nih.gov/geo/）数据库（GEO: GSE118390）下载了6例TNBC患者（PT039、PT058、PT081、PT084、PT089和PT126）的scRNA-seq数据，该数据库包含1189个高质量细胞，并注释了细胞类型。既往研究表明，乳腺癌的恶性细胞起源于上皮细胞，因此，仅保留868个上皮细胞用于恶性细胞的鉴定。此外，下载了240个正常乳腺上皮细胞的基因表达谱，其中39个作为参考（正常）细胞用于拷贝数估算和正常乳腺上皮细胞GRN构建。此外，从TCGA数据库中获取105例正常和115例TNBC样本的基因表达谱，进行差异基因表达分析。METABRIC数据库收集了186例TNBC患者的临床信息，用于生存分析。

恶性上皮细胞和乳腺癌亚型的鉴定

如先前的研究所述，恶性细胞总是伴随着高变CNV（拷贝数变异）。从基因表达谱中推断出868个上皮细胞的CNV。将240个正常上皮细胞作为CNV推断的参考集。根据基因组位置对基因进行排序，以上游和下游基因的平均基因表达量作为推定CNV。在本研究中，使用了intercnv R包进行CNV估计。在从基因表达谱推断出CNV后，进行层次聚类以识别恶性细胞簇。在亚型鉴定中，CNV高变量的细胞簇被认为是恶性细胞。

根据肿瘤细胞的基因表达谱，采用PAM50模型将肿瘤细胞分为5个亚型（normal-like, basal-like, Her2+， LumA，和LumB），该模型基于50个标记基因（PAM50基因特征）广泛用于乳腺癌亚型鉴定。然后，对基因表达谱进行了log2转化。使用genefu R包的SubPred_pam50函数分配每个乳腺癌细胞的亚型。

细胞-细胞通信和GSVA分析

为了研究不同细胞类型（包括5种分子亚型）的细胞间通讯，首先，从之前的研究中获得了细胞类型注释，包括B细胞、T细胞、巨噬细胞、基质细胞和内皮细胞然后使用python CellPhoneDB（https://www.cellphonedb.org/）包估计两种细胞类型间，配体-受体对在不同细胞类型中的意义，p < 0.05的配体受体对被认为是存在显著相互作用的。

为了测试5种分子亚型的差异激活基因集，首先从分子特征数据库（https://www.gsea-msigdb.org/gsea/ msigdb）下载了50个标记基因集。然后，使用GSVA R包对5个分子亚型进行GSVA分析（基因集变异分析）。

GRN的构建于关键节点识别

利用基因共表达和转录结合基序信息构建了乳腺癌各亚型和正常上皮细胞的GRN。首先，这些TF是从先前的研究中获得的，该研究包含1797个特定的TF，接下来，进行共表达分析，以估计TF和靶基因之间的潜在调节强度，然后，在每个亚型和正常上皮细胞中保留富含TF结合基序的TF-靶基因对。使用python pySCENIC软件包进行共表达和TF结合基序富集分析。

在获得正常上皮细胞和5个分子亚型的GRN后，采用节点中心性指标构建网络中节点的重要性。在本研究中，使用度degree、间度betweenness、特征值eigenvalue、PageRank和接近度closeness来评估节点的中心性。这些指标使用R的igraph包计算。

• 度是对应节点的相邻节点的个数。高度节点通常被认为是基本的枢纽。
• 中间度是根据经过相应节点的最短路径数来计算的。高中间度的节点可能成为GRN的瓶颈。
• 特征值通过考虑邻居的重要性来度量节点的重要性。
• PageRank是对应节点随机游走的概率。PageRank类似于特征值，而PageRank引入了阻尼因子，默认为0.85。
• 接近度是对应节点到所有节点的平均距离。接近度高的节点表示该节点位于网络的中心位置。

在计算GRN的每个节点的五个中心性指标后，使用Q统计对每个GRN的五个节点的中心性指标进行整合：$$Q(r_1,r_2,...,r_N)=N!V_N$$$$V_k=\sum _{j=1}^k(-1{)}^{j-1}\frac{V_{k-j}}{j!}{r}_{N-k+1}^j$$其中，$V_0=1,N=5$，$r_i$为第$i$个中心度量指标，$r^j$代表$r$的$j$次方。最后，得到了相应GRN中每个节点的Q统计量。前1%的Q统计量被认为是每个GRN中的关键节点。

功能富集分析和生存分析

为了获得TFs在5种分子亚型特异性GRN中的功能，使用DAVID（Database for Annotation, Visualization and Integrated Discovery）在线工具进行GO注释。此外，从以前的研究中下载了与癌症标志物相关的GO术语。为了可视化，只保留了与癌症标志相关的术语。在本研究中，还研究了ETV6在五种分子亚型中的表达以及TNBC患者的临床预后。采用R生存和survminer包评估基因表达和预后意义。该研究流程图如图所示。

• 首先，使用scRNA-seq数据推断CNV。根据推断的CNV鉴定出恶性细胞。然后利用PAM50模型对每个细胞进行亚型划分，并对每个亚型构建GRN。其次，计算五个中心性指标来衡量每个GRN中节点的重要性。然后使用Q统计量来整合这些中心性指标。最后，对常见关键基因的不同作用进行了评估。对关键基因进行差异表达分析和生存分析。

• 彩色节点表示关键基因。共6类亚型对应的GRN，正常上皮normal epithelial A，normal-like B，basal-like C，Her2+ D，LumA E 和 LumB F。

• A：ETV6在每个亚型中调控不同的基因。线条颜色表示不同亚型的ETV6调控。
• B：各亚型中，被ETV6调控基因的GO注释。彩色点阵表示ETV6调控基因的GO注释。
• C：ETV6在正常上皮细胞和5种分子亚型中的表达。
• D：TCGA正常样本和TNBC样本中ETV6的表达。
• E：基于ETV6表达的METABRIC数据集中TNBC患者的生存分析。低表达的ETV6表现出较好的临床效果。