2.1转录网络推断(TNI)

的交通噪音指数管道从泛型函数开始tni.constructor ()并创建一个TNI-class对象，提供了从高通量基因表达数据进行TRN推断的方法。的tni.constructor采用归一化基因表达的矩阵和相应的基因和样本注释，以及待评估的调控因子列表。在这里，基因表达矩阵和注释可在tniData数据集，提取，预处理和缩小弗莱彻等人(2013)而且柯蒂斯等人(2012)．此数据集包含一个包含3个对象的列表:一个命名的基因表达矩阵(tniData expData美元)，一个带有基因注释的数据帧(tniData rowAnnotation美元)，以及带有示例注释的数据帧(tniData colAnnotation美元)．我们将使用这个数据集来演示一个只有5个规则的小型TRN的构造。但总的来说，我们建议为给定物种注释的所有tf构建规则;请参阅案例研究部分为其他建议。

的tni.constructor方法将检查所有给定参数的一致性。的交通噪音指数然后分三步执行pipeline: (i)计算调节器和所有潜在目标之间的MI，通过排列分析去除不显著的关联;(ii)通过自举去除不稳定的相互作用;(iii)应用ARACNe算法。示例中还提供了其他注释。

库(研制数据(tniData)

#输入1:“expData”，一个命名的基因表达矩阵(基因在行，样本在cols);#输入3:“rowAnnotation”，一个可选的数据帧，带有基因注释#输入4:“colAnnotation”，一个可选的数据帧，带有示例注释TFs <- c(“FOXM1”，“E2F2”，“E2F3”，“RUNX2”，“PTTG1”)rtni <- tni。构造函数(expData = tniData$expData, regulatoryElements = tfs, rowAnnotation = tniData$rowAnnotation, colAnnotation = tniData$colAnnotation) # p.s.或者，'expData'可以是一个' summarizeexperiment '对象

然后tni.permutation ()函数接受预处理TNI-class对象，并返回由排列分析推断出的TRN(对多个假设检验进行更正)。

#请设置nPermutations >= 1000 rtni <- tni。排列(rtni, nPermutations = 100)

不稳定的相互作用，随后通过bootstrap分析使用tni.bootstrap ()函数，它创建了一个共识引导网络，这里称为refnet(引用网络)。

Rtni <- tni.bootstrap(Rtni)

在此阶段，TRN中的每个目标可能与多个tf相关联。由于调节可以通过直接(TF-target)和间接交互(TF-TF-target)发生，因此管道接下来应用ARACNe算法(Margolin, Nemenman, et al. 2006)去除由两个tf和一个共同靶基因组成的任何三元组中最弱的相互作用，保留优势tf -靶对(Meyer, Lafitte, and Bontempi 2008)．ARACNe算法使用数据处理不等式(DPI)定理来丰富具有直接TF-target相互作用的规则，创建DPI过滤的TRN，这里称为tnet(转录网络)。详情请参阅马戈林，王，等人(2006)而且弗莱彻等人(2013)．简单地，考虑三个随机变量，X，Y而且Z形成一个网络三元组X相互作用Z只有通过Y（即。，交互网络为X - Y > > Z)，两者之间不存在其他路径X而且Z)．DPI定理说明了信息之间的传递Y而且Z总是大于传递之间的信息X而且Z．基于此假设，ARACNe算法扫描由两个调节器和一个目标组成的所有三元组，并去除每个三元组中MI值最小的边，认为这是一种冗余关联。

Rtni <- tni.dpi.filter(Rtni)

对于结果监管网络的摘要，我们建议使用tni.regulon.summary ()函数。从下面的总结中，我们可以看到规则的数量，推断目标的数量和规则大小分布。

tni. regulation .summary(rtni) ##由5条规则组成的监管网络。##—dpi过滤网络:## regulatoryElements目标边缘## 5 1217 1380 ##最小第一曲。中值平均第三曲。最大## 76 236 348 276 351 369 ##——参考网络:## regulatoryElements目标边缘## 5 1217 2544 ##最小第一曲。中位数平均第三曲。最大## 76 475 616 509 670 707 ##——

的tni.regulon.summary ()函数还可以让我们获得关于特定规则的详细信息，包括正目标和负目标的数量。请使用此信息作为评估规则集的初始指南。通常小的规则(<15个目标)对于大多数下游方法是没有用的，高度不平衡的规则(如。那些只有积极目标)提供不稳定的活动读数。

tni.regulon。## FOXM1规则有236个目标，这是一个大而平衡的规则。##——DPI过滤网络目标:##总正负## 236 135 101 ##——参考网络目标:##总正负## 616 365 251 ##——相互信息的调节器:## E2F2, E2F3, PTTG1 ## ##——

所有结果可在TNI-class对象可以使用tni.get ()辅助功能。例如，设置What = ' rules .and.mode'将返回一个包含规则的列表，包括为每个交互分配的权重。

规则<- tni。Get (rtni, what = " rules .and。mode"， idkey = "SYMBOL") head(规则$FOXM1) ## C1orf106 HBB MAPT ZNF385B FOSB CBX2 ## 0.1321503 -0.1378999 -0.2063462 -0.1092636 -0.1271561 0.1762198

权重的绝对值代表MI值，而符号(+/-)表示基于调节器与其目标之间的Pearson相关性的预测作用模式。我们还可以将推断的规则检索到igraph-class对象(Csardi 2019)使用tni.graph ()函数，该函数为可视化设置一些基本图形属性红色的R包(Castro et al. 2012)．

G <- tni。图(rtni, regulatoryElements = c(“FOXM1”，“E2F2”，“PTTG1”))

下一部分将展示如何绘制igraph-class对象使用红色的（图1)．

library(reader) rdp <- RedPort() call (rdp) addGraph(rdp, g, layout=NULL) addLegend. dll (rdp, g, layout=NULL)color(rdp, g, type="edge")shape(rdp, g) relax(rdp, ps = TRUE)

图1所示。的三种规则的图表示研制包中。

2．2转录网络分析(TNA)

在前一节中，我们概述了交通噪音指数管道，它计算TRN并提供有关其规则的信息。在这里，我们描述TNA管道，它提供了在规则列表上进行丰富分析的方法。的TNA管道从函数开始tni2tna.preprocess ()，将TNI-class成一个TNA-class对象。然后，规则将测试与给定的基因表达签名的关联，在这里提供tnaData数据集。请注意，本插图中的数据集仅用于演示目的。它由一个包含3个对象的列表组成:一个来自差异基因表达分析的具有log2倍变化的命名数字向量，这里称为“显型”(美元tnaData表型)，一个列出差异表达基因的特征向量(tnaData美元打击)，以及映射到表型的基因注释数据帧(tnaData phenoIDs美元)．

#输入1:“对象”，一个有规则的TNI对象#输入2:“表型”，一个命名的数字向量，通常是log2差异表达水平#输入3:“hits”，一个字符向量，通常是一组差异表达基因#输入4:“phenoIDs”，一个可选的数据帧，基因注释映射到表型数据(tnaData) rtna <- tni2tna。预处理(object = rtni, phenotype = tnaData$phenotype, hits = tnaData$hits, phenoIDs = tnaData$phenoIDs)

的tna.mra ()函数取TNA-class对象，并运行主调节器分析(MRA)(Carro et al. 2010)在一个规则列表(修正多个假设检验)。该MRA评估的重叠之间的每个规则和基因列在支安打论点。

运行MRA方法rtna <- tna.mra(rtna)

所有结果可在TNA-class对象可以使用tna.get ()辅助函数;设置What = 'mra'将返回一个数据帧，列出TRN的基因总数(宇宙。大小)，表示每个规则的目标数目(调节子。大小)中所列的基因数目支安打参数(总计支安打)，即某一规定与“命中次数”(预期。支安打)，即所观察到的给定规则与“命中次数”(观察到。支安打)，计算所观测到的重叠的统计意义phyper ()函数(Pvalue)，调整后的p值(调整。Pvalue)．

#获得MRA结果;# . .设置'ntop = -1'将返回所有结果，无论阈值mra <- tna。get(rtna, what="mra"， ntop = -1) head(mra) ##规则宇宙。大小调节子。总大小。支安打预期。支安打观察到。支安打## 1870 E2F2 5304 348 660 43.30 85 ## 2305 FOXM1 5304 236 660 29.37 62 ## 9232 PTTG1 5304 369 660 45.92 75 ## 1871 E2F3 5304 351 660 43.68 60 ## 860 RUNX2 5304 76 660 9.46 7 ## Pvalue Adjusted.Pvalue ## 1870 1.4e-10 7.0e-10 ## 2305 2.7e-09 6.7e-09 ## 9232 5.8e-06 9.7e-06 ## 1871 5.3e-03 6.6e-03 ## 860 8.5e-01 8.5e-01

作为一种补充的方法，tna.gsea1 ()函数运行单尾基因集富集分析(GSEA-1T)，以寻找与特定“响应”相关的规则，该响应由由差异基因表达签名生成的排序基因列表表示(即。包括在tni - tna预处理步骤中的“表型”)。在这里，调节的目标基因被认为是一个基因集，它是根据表型来评估的。GSEA-1T使用基于排名的评分指标，以测试基因集和表型之间的关联(Subramanian et al. 2005)．

#执行GSEA方法#请设置nPermutations >= 1000 rtna <- tna。gsea1 (rtna nPermutations = 100)

设置What = 'gsea1'在tna.get ()函数将检索列出GSEA统计信息的数据帧，并且可以使用tna.plot.gsea1 ()函数(图2)．

#获取GSEA结果gsea1 <- tna。get(rtna, what="gsea1"， ntop = -1) head(gsea1) ##规则规则。观察到的大小。得分Pvalue已调整。Pvalue## 2305 FOXM1 236 0.68 0.009901 0.012376 ## 1870 E2F2 348 0.67 0.009901 0.012376 ## 9232 PTTG1 368 0.64 0.009901 0.012376 ## 1871 E2F3 351 0.61 0.009901 0.012376 ## 860 RUNX2 76 0.46 0.118810 0.118810

#绘制GSEA结果。gsea1(rtna, labPheno="abs(log2 fold changes)"， ntop = -1)

图2。GSEA分析显示每个调控基因(竖条)按差异基因表达分析(表型)排序。类的输出TNA由tna.gsea1 ()函数。

然而，GSEA-1T并没有表明这种联系的方向。接下来,TNA管道使用双尾GSEA (GSEA- 2t)方法测试该调控与表型是否呈正相关或负相关。的tna.gsea2 ()函数将该调控分为正靶和负靶(基于TF与其靶之间的Pearson相关性)，然后评估靶在排序基因列表中的分布。

#执行GSEA-2T方法#请设置nPermutations >= 1000 rtna <- tna。gsea2(rtna, nPermutations = 100)

设置What = 'gsea2'在tna.get ()函数将检索列出GSEA- 2t统计信息的数据帧，并且可以使用tna.plot.gsea2 ()函数。

#获取GSEA-2T结果gsea2 <- tna。get(rtna, what = "gsea2"， ntop = -1) head(gsea2$differential) ##观察到的大小。得分Pvalue已调整。Pvalue## 9232 PTTG1 368 1.11 0.009901 0.049505 ## 1870 E2F2 348 0.38 0.217820 0.277230 ## 2305 FOXM1 236 0.41 0.257430 0.277230 ## 1871 E2F3 351 0.11 0.267330 0.277230 ## 860 RUNX2 76 0.67 0.277230 0.277230

在GSEA-2T中(图3)，规定的正面目标和负面目标分别被视为独立的pos而且负的基因集，然后根据表型进行评估。对于每个基因集(pos而且负的)沿着排名列表逐级进行。当在基因集中找到一个基因时，它在排序列表中的位置就会被标记出来。连续求和，显示为粉红色和蓝色(pos而且负的基因集，分别)系，增加时，该位置的基因属于基因集，减少时，它不属于。每个运行和到x轴的最大距离表示富集分数。GSEA-2T产生两个每表型富集评分(ES)，其差异(dES = ESpos- - - - - - ES负的)表示规则活动。目标是评估目标基因是否在阳性或阴性差异表达的基因中过度表达。大的阳性dES表明一个诱导(激活)的规则，而大的阴性dES表明一个抑制的规则。详情请参阅坎贝尔等人(2016)而且坎贝尔等人(2018)说明使用这种方法的案例;将GSEA-2T扩展到单个样品卡斯特罗等人(2016)而且Groeneveld等人(2019)．

GSEA-2T结果图。gsea2(rtna, labPheno="log2 fold changes"， tfs="PTTG1")

图3。双尾GSEA分析显示了通过差异基因表达分析(表型)对调控基因的阳性或阴性靶标(红/蓝竖条)进行排序。类的输出TNA由tna.gsea2 ()函数(坎贝尔等人(2016)而且坎贝尔等人(2018)提供如何解释该方法的结果的示例)。

3．1上下文

方法的输入数据准备研制使用公开的mRNA-seq数据，以及TCGA-BRCA队列的临床/分子数据。我们展示了如何从基因组数据共享(GDC)下载统一的GRCh38/hg38数据TCGAbiolinks包(Colaprico et al. 2016)．预处理将生成一个SummarizedExperiment对象，包含基因表达数据，然后我们将使用它来计算brca特定的规则。

3.2包装和数据要求

在继续之前，请确保已安装所有库。

library(RTN) library(TCGAbiolinks) library(summarizeexperiment) library(TxDb.Hsapiens.UCSC.hg38.knownGene) library(snow)

3．3数据预处理

我们将使用Bioconductor包TCGAbiolinks查询及下载。我们需要TCGA-BRCA队列的GRCh38/hg38协调、标准化RNA-seq数据。TCGAbiolinks将在你的工作目录中创建一个名为GDCdata的目录，并将从GDC下载的文件保存在其中。每个患者的文件将在单独的文件中下载，我们将使用500个病例的子集仅用于演示目的。由于需要下载大量的rna -seq文本文件，总计>140 MB，下载需要一段时间。然后,GDCprepare ()函数将文件编译为类的R对象RangedSummarizedExperiment．的RangedSummarizedExperiment有6个插槽。最重要的是rowRanges(基因注释),colData(示例注释)，和化验，其中包含基因表达矩阵。

#为TCGA-BRCA队列设置GDCquery #基因表达数据将根据hg38查询进行比对<- GDCquery(project = "TCGA-BRCA"， data. query)category =“转录组分析”，数据。type = "基因表达量化"，工作流。type = "HTSeq - FPKM-UQ"， sample。type = c("原发实体瘤"))#获取一个子集用于演示(n = 500个病例)病例<- getResults(query, cols = "cases")病例<- sample(cases, size = 500) query <- GDCquery(project = "TCGA-BRCA"， data. txt)category =“转录组分析”，数据。type = "基因表达量化"，工作流。type = "HTSeq - FPKM-UQ"， sample。type = c(“原发实体瘤”)，barcode =病例)GDCdownload(查询)tcgaBRCA_mRNA_data <- GDCprepare(查询)

从GDC下载的对象包含基因级表达数据，包括编码和非编码基因(如。lncRNAs)。我们将过滤这些，只保留UCSC hg38基因列表中注释的基因(~30,000个基因)。

# subsetByOverlaps(tcgaBRCA_mRNA_data, geneRanges)

nrow(rowData(tcgaBRCA_mRNA_data)) ## [1] ~ 30000

最后，我们将简化队列注释中的名称，以便在后续RTN函数中更好地总结。运行此步骤后，tcgaBRCA_mRNA_data对象已为交通噪音指数管道。请确保保存tcgaBRCA_mRNA_data对象在适当的文件夹中。

#修改基因注释中的列名，以便更好地总结colnames(rowData(tcgaBRCA_mRNA_data)) <- c("ENSEMBL"， "SYMBOL"， "OG_ENSEMBL")

#保存预处理数据，以便后续分析Save (tcgaBRCA_mRNA_data, file = "tcgaBRCA_mRNA_data_preprocessed.RData")

3.4推断转录调节网络

接下来，我们将为TCGA-BRCA队列生成规则，遵循中描述的步骤快速启动部分。对于规则重构，我们将使用可在tfsData数据集，由1612个tf编译兰伯特等人(2018)．该管道应用于包含至少100个转录组谱的数据集。这代表了ARACNe算法的经验下界(Margolin, Wang, et al. 2006)．

#加载TF注释数据("tfsData")

#检查TF注释:#从Lambert et al.(2018)与基因注释#从TCGA-BRCA队列geneannot <- rowData(tcgaBRCA_mRNA_data) regulatoryElements <- Intersect (tfsData$Lambert2018$SYMBOL, geneannot$SYMBOL)

#运行TNI构造函数rtni_tcgaBRCA <- TNI。构造函数(expData = tcgaBRCA_mRNA_data, regulatoryElements = regulatoryElements)

我们可以提供一些一般性的实际指导。为了计算一个大型的管理网络，我们建议使用多线程模式雪包中。作为最少的计算资源，我们建议处理器具有>= 4核和RAM >= 8 GB每核(您应该调整特定的例程以适应您的可用资源)。的makeCluster ()函数将设置要在本地机器上创建的节点数量，从而生成一个集群对象的TNI-class方法。例如，在2.9 GHz Core i9-8950H工作站和32GB DDR4 RAM上运行时，从约30,000个基因和500个样本的基因表达矩阵中重建1600个规则需要大约3小时。请注意，并行运行大型数据集可能会消耗所有系统内存。我们建议监控并行化，以避免工作太接近内存上限;的parChunks参数中可用的tni.permutation ()而且tni.bootstrap ()函数可用于调整发送用于并行化的作业大小。最后，对于RNA-seq数据，我们建议使用互信息的非参数估计器(缺省选项)tni.permutation ()功能)。

通过排列和自举分析计算参考监管网络。#请根据您可用的硬件选项设置“spec”(cluster=snow::makeCluster(spec=4，“SOCK”))rtni_tcgaBRCA <- tni。permutation(rtni_tcgaBRCA, pValueCutoff = 1e-7) rtni_tcgaBRCA <- tni.bootstrap(rtni_tcgaBRCA) stopCluster(getOption("cluster"))

接下来我们运行ARACNe算法Eps = 0，它设置了每个三元组中重量最小的边的移除公差。例如，XY在边缘XYZ如果重量低于三联，则移除三联YZ - eps而且XZ - eps(见2.1节和aracne ()函数的文档)。根据经验，应该使用0(无公差)到0.15之间的值(Margolin, Wang, et al. 2006)．作为一种不那么武断的方法，我们建议设置eps = NA，它将从排列和自举步骤中计算的空分布估计阈值。

#计算经dpi过滤的监管网络rtni_tcgaBRCA <- tni.dpi. dpi。filter(rtni_tcgaBRCA, eps = 0)

#保存TNI对象，以便后续分析Save (rtni_tcgaBRCA, file="rtni_tcgaBRCA. rdata ")

的选择pValueCutoff阈值将取决于假阳性和假阴性之间的理想权衡。一个明智的pValueCutoff可以选择阈值先天的通过将期望的假阳性数量除以测试tf -目标相互作用的数量(Margolin, Wang, et al. 2006)．例如，在本例中，我们测试了约5e+7 TF-target相互作用(即。约1,600个tf vs.约30,000个基因);因此pValueCutoff1e-7的阈值应该会导致大约5个假阳性。

请注意，某些程度的缺失注释是预期的，因为不是所有的基因符号列在regulatoryElements可能存在于TCGA-BRCA预处理数据中。此外，与计算不一致的数据可以在tni.constructor进行预处理。例如，如果一个基因的表达在队列中没有变化，则不可能将该基因的表达与队列中其他基因的表达联系起来。作为一个极端的例子，没有变异的基因(如。并不是所有样本中都有表达)排除在分析之外。对于结果监管网络的摘要，我们建议使用tni.regulon.summary ()函数(请参阅快速启动部分)。

还要注意，MI指标是基于基因在队列中的表达变化。大量的肿瘤样本通常包含多个分子亚型，并且通常为构建调控提供良好的表达可变性。相比之下，更均匀的样本集可能更具有挑战性的探索规则，这可能是正常的，非癌症样本集的情况。我们不建议对低表达可变性的样本集计算规则。

3．5从不同队列中推断规律时，如何选择显著性水平

比较不同队列的法规信息的一个重要挑战是在每个队列中在相似的统计条件下生成法规。例如，队列一个有300个样本的基因表达数据吗B100个样本。在每个队列中，我们会推断规则研制通过连续运行四个步骤:tni.constructor ()检查输入数据;tni.permutation ()，在这里我们设置了pValueCutoff为基因对表达谱之间的互信息(MI);tni.bootstrap ()，用于评估MI的稳定性;最后tni.dpi.filter ()，这增强了直接调节-目标相互作用。

如果我们设置一个排列pValueCutoff = 1e-5推断…的规则一个，我们应该选择什么pValueCutoff为B？这里的问题是p值是样本容量的逆函数;也就是说，对于较小的队列，我们应该使用较大(不太严格)的p值。理想情况下，应该确定第二个队列的样本量先天的，但这很少是一种选择，因为规则通常是从大型回顾性队列研究中重建的。

对于RTN调控，“效应量”是调控因子和潜在靶基因之间的MI值，我们的零假设是调控因子和潜在靶基因的表达谱在统计上无显著关系;因此，对于这个调控因子，转录网络不应该保留这个靶基因。一般而言，在RTN的MI排列计算中，p值阈值控制了I型错误率(即。当我们拒绝一个真实的零假设，产生一个假阳性的结果)，但它也会影响第二类错误的发生率(当我们未能拒绝一个虚假的零假设，产生一个假阴性的结果)米勒(2019)．类型I和类型II错误率分别用α和β表示，统计检验的幂值定义为1 - β。α减少β增加，反之亦然(Mudge et al. 2012)．因为α决定了检测给定效应大小的能力，a的选择pValueCutoff对β做出隐性决策;也就是真正的监管互动研制可能会认为不重要，所以会拒绝。当涉及到评估从两个不同队列重建的规则之间的监管共性时(如。一个而且B)，但只在一个应该被认为是严重的错误，错误地发现一个不存在的影响，只在B．为了将此类问题的风险降至最低，我们建议通过观察I型和II型错误之间的权衡来调整α。Mudge et al. (2012)开发R代码来计算相关性的最佳α水平;我们通过实现tni.alpha.adjust ()功能，以协助选择pValueCutoff当分析不同样本数量的数据集时。

例如，在给定'nA'， 'alphaA'和'betaA'的情况下，估计'nB'中的'alphaB'， alphaB <- tni.alpha。调整(nB = 100, nA = 300，阿尔法a = 1e-5，贝塔aa = 0.2)阿尔法b # [1] 0.029

对于第二组B，我们建议使用接近的值alphaB为pValueCutoff在tni.permutation ()功能(和队列一个，alphaA)，以便研制将返回两个队列的I型和II型错误之间的相似权衡所生成的规则结果。

4．1上下文

弗莱彻等人(2013)利用METABRIC乳腺癌队列的微阵列转录组数据重建了809个转录因子(tf)的调控(Curtis et al. 2012)．卡斯特罗等人(2016)发现其中36个TF调控因子与乳腺癌的遗传风险相关。风险tf分为两个不同的集群。“簇1”风险tf与雌激素受体阳性(ER+)乳腺癌风险相关，包括ESR1、FOXA1和GATA3，而“簇2”风险tf与雌激素受体阴性(ER-)、基底样乳腺癌相关。我们在这里的目标是演示如何生成单个肿瘤样本的规则活动概况，使用由弗莱彻等人(2013)36个风险tf。

4.2包装和数据要求

在继续之前，请确保已安装所有库。的Fletcher2013b数据包可从R/Bioconductor存储库获得。安装并加载此包将提供本案例研究所需的所有数据。的rtni1st数据集是预处理的TNI-class对象重构的规则弗莱彻等人(2013)METABRIC队列1 (n=997个原发肿瘤)的基因表达矩阵。

library(RTN) library(Fletcher2013b) library(pheatmap)

#加载“rtni1st”数据对象，其中包括规则和表达式配置文件数据(“rtni1st”)

#感兴趣的转录因子列表(这里有36个风险相关的转录因子)。tfs < - c(“AFF3”、“基于“增大化现实”技术”,“ARNT2”、“BRD8”、“CBFB”、“CEBPB”、“E2F2”、“E2F3”、“三”、“ESR1”、“FOSL1”、“FOXA1”、“GATA3”、“GATAD2A”、“LZTFL1”、“MTA2”、“MYB”、“MZF1”、“独立”、“PPARD”、“RARA”、“RB1”、“RUNX3”、“SNAPC2”、“SOX10”、“SPDEF”、“TBX19”、“TCEAL1”、“TRIM29”、“XBP1”、“YBX1”、“YPEL3”、“ZNF24”、“ZNF434”、“ZNF552”、“ZNF587”)

4.３规则活动概要

调控活性谱(rap)描述了一个队列中样本之间调控程序的相似性和差异。为了评估大量的样本，我们实现了一个函数来计算整个队列的双尾GSEA(更多细节在Groeneveld等人(2019))．简单地说，对于每个规则，tni.gsea2 ()函数中每个样本的规则活动分数TNI-class对象。对于样本中的每个基因，从其在样本中的表达相对于其在队列中的平均表达计算差异基因表达;然后将基因排序为表示该样本中差异基因表达签名的排序列表，该列表用于运行GSEA-2T，如中所述tna.gsea2方法(参见2.2节)．

#计算单个样本rtni1st <- tni的规则活动。gsea2(rtni1st, regulatoryElements = risk.tfs) metabric_regact <- tni. tfs。get(rtni1st, what = "regulonActivity")

的tni.gsea2返回一个包含计算富集分数的列表:ESpos,西文负的dES表示单个样品的调控活性。接下来,pheatmap包用于生成热图，显示rap和一些示例属性(图4)．

#从'rtni1st'数据集metabric_annot <- tni获取样本属性。get(rtni1st， "colAnnotation") #为pheatmap属性获取ER+/-和PAM50属性<- c("LumA"，"LumB"，"Basal"，"Her2"，"Normal"，"ER+"，"ER-") metabric_annot <- metabric_annot[， atbs]

# Plot regular activity profiles pheatmap(t(metabric_regact$dif)， main=" metabric_cohort 1 (n=977 samples)"， annotation_col = metabric_annot, show_colnames = FALSE, annotation_legend = FALSE, clustering_method =" ward. "D2"， fontsize_row = 6, clustering_distance_rows = "correlation"， clustering_distance_cols = "correlation")

图4。来自METABRIC队列1的单个肿瘤样本的调控活性谱(对于使用该方法的示例，请参见图4A和图5A罗伯逊等人(2017))．