3.1.1获得ctss

CAGEfightR从简单的ctss开始分析，这是在基因组中映射到每个碱基对(bp)的CAGE标签5 '末端的数量。ctss通常存储为bed -file、bedgraph -file或bigwig -file。由于CTSSs是稀疏的(只有一小部分bps是CTSSs)，这些文件相对较小，因此易于共享，许多研究将通过在线存储库提供CTSSs，例如地理．

当从新的CAGE库准备ctss时，标签通常首先进行条形码分割、修剪和过滤，然后使用标准命令行工具映射到参考基因组。确切的步骤将取决于所使用的给定CAGE协议。CTSSs随后可以从生成的bam文件中每次提取一个库，例如使用genomecov从bedtools与5设置。TSRchitect而且篮球选手也包括函数(inputToTSS ()，getCTSSs ()，分别)用于从R内的bam文件中获取ctss篮球选手在映射时可以选择校正g偏置。

CAGEfightR可以分析许多类型的5 '端数据，只要它们可以用类似于CTSSs的格式表示。

3.1.2进口ctss

CAGEfightR使用方便BigWigFile而且BigWigFileList用于处理存储在bigwig文件中的ctss的容器(每条链一个文件)，因为这些容器允许对存储在文件中的基因组信息进行快速随机访问、检查和总结(例如通过进口()，seqinfo ()而且总结())．首先，我们需要分辨CAGEfightR在哪里可以找到硬盘上包含CTSSs的bigwig文件。通常，我们会提供每个文件的路径(例如。/ CAGEdata BigWigFiles / Sample1_plus.bw)，但这里我们将使用纳米管数据包:

#在硬盘bw_plus <- system. exe上设置文件路径。bw_minus <- system. file("extdata"， nanotubes$BigWigPlus, package = "nanotubes"， mustWork = TRUE)file("extdata"， nanotubes$BigWigMinus, package = "nanotubes"， mustWork = TRUE) #保存为命名的BigWigFileList bw_plus <- BigWigFileList(bw_plus) bw_minus <- BigWigFileList(bw_minus) names(bw_plus) <- names(bw_minus) <- nanotubes$Name

第一步是量化所有样本的CTSS使用情况。这通常是最耗时的步骤之一CAGEfightR但它可以通过使用多核来加速(如果可用，请参阅材料和方法)。我们还需要指定基因组，我们可以从BSgenome.Mmusculus.UCSC.mm9基因组计划:

CTSSs <- quantifyCTSSs(plusStrand = bw_plus, minusStrand = bw_minus, genome = seqinfo(bsg)， design = nanotubes) #>正在检查设计…#>检查提供的基因组兼容性…#>在11个样本中使用3个worker迭代1个基因组瓦…从正链导入CTSSs…从负链导入CTSSs…#>合并股…格式化输出…#> ### CTSS摘要### #>样本数量:11 #> CTSS数量:933.39万#>稀疏性:81.68% #>类型的行wranges: StitchedGPos #>最终对象大小:280,883,032 B

大约有900万ctss被存储为RangedSummarizedExperiment，这是Bioconductor高通量实验的标准容器。我们可以检查范围和CTSS计数:

#>类:rangedsummarize实验#> dim: 9338802 11 #>元数据(0):#> assays(1):计数#> rownames: NULL #> rowData names(0): #> colnames(11): C547 C548…N17 N18 #> colData names(4): Class Name BigWigPlus BigWigMinus # Extract CTSS positions rowRanges(CTSSs) #> StitchedGPos object with 9338802 positions and 0 metadata columns: #> seqnames pos strand #>    # b[1] [1] chr1 3024556 + #> [2] chr1 3025704 + #> [3] chr1 3025705 + #> [4] chr1 3028283 + #> [5] chr1 3146133 + #> ... ... ... ...#> [9338798] chrUn_random 5810899 - #> [9338799] chrUn_random 5813784 - #> [9338800] chrUn_random 5880838 - #> [9338801] chrUn_random 5893536 - #> [9338802] chrUn_random 5894263 - #> ------- #> seqinfo: 35个序列(1个循环)从mm9基因组#提取CTSS计数分析(CTSSs，“计数”)%>%头#> 6 x 11稀疏矩阵类“dgCMatrix”#>[[压缩11列名称'C547'， 'C548'， 'C549'…#> #> [1，] . .1 . . . . . . . .#>[2，]…1 . . . . . . .#> [3，] . . . .1 . . . . . .#> [4，] . . . . 1 . . . . . . #> [5,] . . . . . . 1 . . . . #> [6,] . 1 . . . . . . . . .

3.1.3用于寻找TSS和增强器候选的单向和双向聚类:

CAGEfightR我们首先将每个样本中的CTSS计数归一化为每百万标签(TPM)值，然后将所有样本中的TPM值相加:

CTSSs <- CTSSs %>% calcTPM() %>% calcPooled() #>计算库大小…计算TPM…

这将添加一些新的信息ctss:每个库中的标签总数，这是一种新的分析方法TPM，以及每个CTSS的合并信号。

我们可以用单向聚类来定位单向集群，通常简称为标签集群(tc)，它们是tss的候选。的quickTSSs将在一个函数调用中定位和量化tc:

TCs <- quicktss (CTSSs) #>使用现有的分数列!#> #> - Running clusterunidirectional: #> split by strand…#>切片缩减找到集群…#>统计…准备输出…#>标签集群摘要:# b>宽度计数百分比# b> Total 3602099 1e+ 02% # b> >=1 2983433 8e+01 % #> >=10 577786 2e+01 % #> >=100 40842 1e+00 % #> >=1000 38 1e- 03% #> # b> - Running quantifyClusters: #> Finding overlaps…#>在集群内聚合…

注意:quickTSSs运行CAGEfightR使用默认设置。如果你有更大或更嘈杂的数据集，你很可能想用不同的设置进行更稳健的分析。看到CAGEfightR有关更多信息的小插图。

许多鉴定出的tc表达非常低。为了获得可能与生物学相关的tss，我们只保留至少5个样本中表达在1个TPM以上的TCs(5个样本是最小的实验组的大小):

tss <- TCs %>% calcTPM() %>% subsetBySupport(inputAssay = "TPM"， unexpression = 1, minSamples = 4) #>计算支持…计算库大小…在calcTotalTags中#>警告(object = object, inputAssay = inputAssay, outputColumn #> = outputColumn): object已经在colData中有一个名为totalTags的列:它将#>被覆盖!计算TPM…# >构造子集……从3602099个区域中移除3573214 (99.2%)

这除去了大量低表达的tc，留给我们近3万个tc进行分析。为简单起见，我们将这些tc称为TSS候选人，因为每个TC都可以被视为一个转录本或基因的TSS的位置和活性的测量。请注意，这是一种简化，因为从技术上讲，TC将几个bp精确的ctss组合在一起。

然后我们转向双向聚类用于识别双向集群(bc)。类似地，我们可以用quickEnhancers来定位和量化bc (bc通过聚类的两条链上的求和标签进行量化):

BCs <- quickEnhancers(CTSSs) #>使用现有的分数列!#> #> - Running clusterbidirectional: #>预过滤双向候选区域…#>保留分析:68% #>分流…计算正链上的窗口覆盖率…计算负链上的窗口覆盖率…#>计算平衡得分…#>切片减少寻找双向集群…#>统计…准备输出…#> #双向集群汇总:#>双向集群数量:106779 #>最大平衡评分:1 #>最小平衡评分:0.950001090872574 #>最大宽度:1866 #>最小宽度:401 #> #> - Running subsetByBidirectionality: #>计算双向… #> Subsetting... #> Removed 73250 out of 106779 regions (68.6%) #> #> - Running quantifyClusters: #> Finding overlaps... #> Aggregating within clusters...

注意:quickEnhancers运行CAGEfightR使用默认设置。如果你有更大或更嘈杂的数据集，你很可能想用不同的设置进行更稳健的分析。看到CAGEfightR有关更多信息的小插图。

同样，我们通常对低表达的bc不感兴趣。由于与tc相比，bc总体上表达较低，我们将简单地过滤掉5个样本中不需要至少1个计数的bc:

BCs <- subsetBySupport(BCs, inputAssay = "counts"， unexpression = 0, minSamples = 4) #>计算支持度…# >构造子集……从33529个区域中移除20017 (59.7%)

3.1.4用文本模型注释集群

在定位了单向和双向簇后，我们可以根据已知的转录本和基因模型对它们进行注释。这些类型的注释是通过存储的TxDb- Bioconductor中的对象。这里我们将简单地使用UCSC成绩单所包含的TxDb.Mmusculus.UCSC.mm9.knownGene包装，但是CAGEfightR小插图包括如何获得TxDb对象从其他来源(GFF/GTF文件，AnnotationHub等)。

从TSS候选对象开始，我们不仅可以注释带有重叠文本的TSS，还可以注释在什么地方部分TSS是通过使用层次注释方案实现的。由于一些候选TSS可能非常宽，我们只使用TSS峰值用于注释:

#注释TSSs文本id <- assignTxID(TSSs, txModels = txdb, swap = "thick") #>提取文本…#>发现层次重叠…#> ###重叠摘要:### #>特征重叠转录本:87.65% #>唯一转录本数量:31898 #注释转录本上下文TSSs <- assignTxType(TSSs, txModels = txdb, swap = "thick") #>发现交换范围的层次重叠…#> ##重叠总结:### #> txType计数百分比#> 1启动子13395 46.4 #> 2近端2246 7.8 #> 3 fiveUTR 2112 7.3 #> 4 threeUTR 1200 4.2 #> 5 CDS 3356 11.6 #> 6外显子161 0.6 #> 7内含子2810 9.7 #> 8反义1294 4.5 #> 9基因间2311 8.0

几乎一半的tss是在注释启动子中发现的，而另一半与UCSC已知的转录本相比是新的。

转录本注释对于增强子分析特别有用，因为双向聚类也可以检测双向启动子。因此，常用的筛选方法是只考虑基因间区或内含子区的bc作为增强子候选:

#用文本上下文BCs注释<- assignTxType(BCs, txModels = txdb, swap = "thick") #>查找交换范围的层次重叠…重叠摘要:### #> txType计数百分比#> 1启动子766 5.7 #> 2近端1649 12.2 #> 3 fiveUTR 67 0.5 #> 4 threeUTR 596 4.4 #> 5 CDS 420 3.1 #> 6外显子71 0.5 #> 7内含子6815 50.4 #> 8反义0 0.0 #> 9基因间3128 23.1 #只保留非外显子BCs作为增强子候选增强子<-子集(BCs, txType %in% c(“基因间”，“内含子”))

剩下近10000个bc需要分析。同样，为了简单起见，我们将这些非外显子bc称为增强剂的候选人对于工作流程的剩余部分。

3.1.5合并到单个数据集

对于许多下游分析，特别是规范化和差分表达式，将TSS和候选增强子组合到一个数据集中是有用的。这确保了集群不会重叠，因此每个CAGE标记只被计数一次。

我们必须首先确保增强程序和TSS候选程序具有相同的附加信息，因为CAGEfightR只允许合并具有相同样本和集群信息的集群:

# Clean colData TSS $totalTags <- NULL Enhancers$totalTags <- NULL # Clean rowData rowData(TSS)$balance <- NA rowData(TSS)$bidirectionality <- NA rowData(Enhancers)$txID <- NA #添加标签，使以后的检索更容易rowData(TSS)$clusterType <- "TSS" rowData(Enhancers)$clusterType <- "Enhancer"

然后集群可以合并:由于增强子在技术上可以被检测为两个独立的TSS，我们只在TSS和增强子候选重叠时保留增强子:

RSE <- combineClusters(object1 = tss, object2 = Enhancers, removeifoverloverl重叠= "object1") #>从object1: 374中删除重叠特征#>保持分析:计数#>保持列:score, thick, support, txID, txType, balance, bidirecality, clusterType #>合并元数据…堆叠和重排序…

我们最终计算出最终合并数据集的标签总数和tpm缩放计数:

RSE <- calcTPM(RSE) #>计算库大小…计算TPM…

3.2.1之上tss和增强子的基因组浏览图

首先，我们可以简单地在基因组浏览器风格的图中绘制tss和增强子的一些例子Gviz包(Hahne and Ivanek 2016)．它需要一些代码来设置，但产生的轨道可以在后面的示例中重用:

#基因组轨迹axis_track <- GenomeAxisTrack() #注释轨迹tx_track <- GeneRegionTrack(txdb, name = "Gene Models"， col = NA, fill = "bisque4"， shape = "arrow"， showId = TRUE)

快速生成基因组浏览器图的一个良好的一般策略是首先定义一个感兴趣的区域，然后仅使用该区域内的数据subsetByOverlaps．下面的代码使用第一个候选TSS演示了这一点:

#提取第一个TSS plot_region周围的100 bp <- RSE %>% rowRanges() %>%子集(clusterType == "TSS") %>% .[1] %>% add(100) %>% unstrand() # CTSS track ctss_track <- CTSSs %>% rowRanges() %>% subsetByOverlaps(plot_region) %>% trackCTSS(name = "CTSSs") #>分割pooled信号的链…#>准备轨道…#集群跟踪cluster_track <- RSE %>% subsetByOverlaps(plot_region) %>% trackClusters(name = "Clusters"， col = NA, showId = TRUE) #>设置thick和thin特性…#>合并排序…#>准备轨道…# Plot tracks together (list(axis_track, ctss_track, cluster_track, tx_track)， from = start(plot_region)， to = end(plot_region)，染色体= as.character(seqnames(plot_region)))

图1:TSS候选的基因组浏览器示例

顶部轨道显示了合并的CTSS信号，中间轨道显示了识别的TSS候选信号。粗条表示TSS候选峰值(TSS候选中使用最多的总体ctss)。底部的轨迹显示了该基因组位置的已知转录模型。在这种情况下，cage定义的TSS候选与注释很好地对应。

我们还可以绘制第一个增强子候选:

#在第一个增强器plot_region周围提取100 bp <- RSE %>% rowRanges() %>%子集(clusterType == " enhancer ") %>% .[1] %>% add(100) %>% unstrand() # CTSS track ctss_track <- CTSSs %>% rowRanges() %>% subsetByOverlaps(plot_region) %>% trackCTSS(name = "CTSSs") #>通过链分割pooled信号…#>准备轨道…#集群跟踪cluster_track <- RSE %>% rowRanges %>% subsetByOverlaps(plot_region) %>% trackClusters(name = "Clusters"， col = NA, showId = TRUE) #>设置thick和thin特性…#>合并排序…#>准备轨道…# Plot tracks together (list(axis_track, ctss_track, cluster_track, tx_track)， from = start(plot_region)， to = end(plot_region)，染色体= as.character(seqnames(plot_region)))

图2:增强子候选的基因组浏览器示例

这里我们看到了主动增强子的双向模式特征。增强子候选在中间的轨道中。厚的中点标志着增强子候选中的最大平衡点。

3.2.2tss和增强子的定位和表达

除了查看TSS和增强子候选的单个示例之外，我们还希望获得与转录注释相关的集群的数量和表达的概述。中提取所有必要的数据RangedSummarizedExperiment变成一个普通的data.frame：

cluster_info <- RSE %>% rowData() %>% as.data.frame()

然后我们使用ggplot2绘制每个注释类别中的聚类的数量和表达级别:

#集群数量ggplot(cluster_info, aes(x = txType, fill = clusterType)) + geom_bar(alpha = 0.75, position = "dodge"， color = "black") + scale_fill_colorblind("Cluster type") + labs(x = "Cluster annotation"， y = "Frequency") + theme(axis.text. text)。X = element_text(角度= 90,hjust = 1))

图3:每个注释类别中的集群数量

#集群表达式ggplot(cluster_info, aes(x = txType, y = log2(score / ncol(RSE))， fill = clusterType)) + geom_violin(alpha = 0.75, draw_quantiles = c(0.25, 0.50, 0.75)) + scale_fill_colorblind("Cluster type") + labs(x = "Cluster annotation"， y = "log2(TPM)") + theme(axis.text. txt)。x = element_text(angle = 90, hjust = 1)) #>正则化警告。值(x, y, ties, missing(ties)， na。Rm = na.rm): #>崩溃为唯一的'x'值

图4:每个注释类别中的集群的表达式

我们发现TSS候选基因在注释启动子上普遍高表达。大多数新型候选TSS表达水平较低，除5 ' -UTRs中部分候选TSS表达水平较低。增强子候选表达水平远低于TSS候选表达水平。

3.2.3分析TSS的形状和序列

CAGE数据的经典分析是将tss划分为锋利的而且广泛的类，显示不同的核心启动子区域和不同的组织表达模式(Carninci et al. 2006)．

CAGEfightR能算几个形状统计总结了TSS的形状。分位数范围(IQR)可以用来发现尖锐和广泛的tss。IQR测量bp-distance保持，例如TSS候选中10-90%的池化表达，这抑制了可能的离散标记的影响，这些标记可以极大地扩展TSS候选的宽度，而对表达没有多大贡献。由于低表达的候选TSS由于其低宽度和低数量的标签而不能表现出很大的形状变化，我们在这里关注高表达的候选TSS(平均TPM >= 10):

#选择高表达的TSS hightss <-子集(RSE, clusterType == 'TSS' & score / ncol(RSE) >= 10) #计算IQR为10%-90%的区间hightss <- calcShape(hightss, pooled = CTSSs, shapeFunction = shapeIQR, lower = 0.10, upper = 0.90) #>拆分#>应用功能到每个集群…准备输出输出…

然后我们可以绘制出iqr的双峰分布。我们使用一个放大面板来突出这两个类之间的区别:

highTSSs %>% rowData() %>% as.data.frame() %>% ggplot(aes(x = IQR)) + geom_histogram(binwidth = 1, fill = "hotpink"， alpha = 0.75) + geom_vline(x截取= 10,linetype = "虚线"，alpha = 0.75, color = "黑色")+ facet_zoom(xlim = c(0,100)) +实验室(x = "10-90% IQR"， y = "频率")

图5:高表达tss四分位范围(IQRs)的双峰分布

我们发现大多数候选TSS的IQR都低于或高于10 bp(虚线)，所以我们使用这个截断线将候选TSS分为Sharp和Broad:

#划分组rowData(hightss)$shape <- ifelse(rowData(hightss)$IQR < 10， "Sharp"， "Broad") #计算组大小表(rowData(hightss)$shape) #> #> Broad Sharp #> 9555 812

我们现在可以研究两类候选TSS的核心启动子序列。我们首先需要提取每个候选TSS的周围启动子序列:我们将其定义为TSS候选峰值-40/+10 bp，并使用BSgenome.Mmusculus.UCSC.mm9基因组计划:

promoter_seqs <- hightss %>% rowRanges() %>% swapRanges() %>% promoters(upstream = 40, downstream = 10) %>% getSeq(bsg， .)

返回一个DNAStringSet-对象，我们可以绘制作为序列标志(Schneider and Stephens 1990)通过ggseqlogo包(Wagih 2017)：

Promoter_seqs %>% as。character %>% split(rowData(hightss)$shape) %>% ggseqlogo(data = .， ncol = 2, nrow = 1) + theme_logo() + theme(轴.title.)x = element_blank()， axis.text.x = element_blank()， axis.ticks.x = element_blank()) #>警告:' guides( = FALSE) '已弃用。请使用“guides( = #> "none")”代替。

图6:Sharp和Broad类tss核心启动子区域序列标识

正如预期的那样，我们观察到Sharp TSS候选者与Broad TSS候选者相比，在TSS峰值的上游往往有一个更强的TATA-box。

3.2.4寻找相互作用tss和增强子的候选

除了简单地识别增强子之外，尝试并推断它们可能调控的基因也很有趣。CAGE数据本身不能提供增强器与TSS物理相互作用的直接证据。这将需要专门的染色质确认捕获测定，如HiC, 4C, 5C等。然而，以往的研究表明，相互接近且表达高度相关的tss和增强子更容易相互作用。因此，我们可以使用tss和增强子之间的距离和表达相关性来识别tss -增强子链接作为物理相互作用的候选(Andersson et al. 2014)．

为了做到这一点CAGEfightR，我们首先需要将这两种类型的聚类表示为一个具有两个层次的因子:

rowData(RSE)$clusterType <- RSE %>% rowData() %>% use_series("clusterType") %>% as_factor() %>% fct_relevel("TSS")

然后，我们可以计算tss和增强子之间在50 kbp距离内的所有成对相关性。在这里，我们使用非参数Kendall 's tau作为相关性的度量，但可以提供其他用于计算相关性的函数(例如，可以在对数转换TPM值上计算Pearson 's r，仅捕获线性关系):

#查找链接并计算相关性all_links <- RSE %>% swapRanges() %>% findLinks(maxDist = 5e4L, directional = "clusterType"， inputAssay = "TPM"， method = "kendall") #>查找从TSS到Enhancer的方向链接…计算41311对相关性…准备输出…#> #链路汇总:#>链路数量:41311 #>成对距离汇总:#>最小第1段，中位数，平均第3段，最大值。#> 205 8832 21307 22341 35060 50000 #显示链接all_links #> GInteractions对象与41311个交互和4个元数据列:#> seqnames1 ranges1 seqnames2 ranges2 | orientation #>     |  # b> [1] chr1 6204746—chr1 6226837 | downstream #> [2] chr1 7079251—chr1 7083527 | downstream #> [3] chr1 9535519—chr1 9554735 | downstream #> [4] chr1 9538162—chr1 9554735 | downstream #> [5] chr1 20941781—chr1 20990601 | downstream #> ... ... ... ... ... ... . ...#> [41307] chr9_random 193165——chr9_random 217926 | upstream #> [41309] chr9_random 223641——chr9_random 217926 | downstream #> [41310] chr9_random 223641——chr9_random 242951 | upstream #> [41311] chrUn_random 3714359——chrUn_random 3718258 | upstream #> distance estimate .value #>    #> [1] 22090 -0.0603023 0.805434 #> [2] 4275 0.3659942 0.128613 #> [3] 19215 -0.2132007 0.392330 #> [4] 1657248819 0.1407053 - 0.565461 0.3411211 - 0.171112 # > [5]  #> ... ... ... ...#>[41307] 24760 0.4770843 0.0423302 #>[41308] 49785 0.1809068 0.4599290 #>[41309] 5714 -0.0366988 0.8758961 #>[41310] 19309 -0.2613098 0.2857948 #>[41311] 3898 -0.1705606 0.4937737 #> ------- #>区域:38454范围和8个元数据列#> seqinfo:来自mm9基因组的35个序列(1个循环)

输出为GInteractions-object fromInteractionSet包(Lun, Perry, and Ing-Simmons 2016):对于每个TSS增强器链接，除了相关估计和p值外，还计算距离和方向(相对于TSS的上游/下游)。如果要提取一组高度相关的链接进行进一步分析，则可以适当地对p值进行多次测试，例如使用p.adjust ()．现在，我们只对顶部的正相关感兴趣，所以我们对链接进行子集和排序:

#子集仅为正相关cor_links <-子集(all_links，估计> 0)#基于相关性cor_links <- cor_links[order(cor_links$估计，递减= TRUE)]

然后，我们可以可视化整个基因组区域的相关模式，这里使用最相关的tss增强子链接:

#提取感兴趣的链接周围的区域plot_region <- cor_links[1] %>% boundingBox() %>% linearize(1:2) %>% add(1000L) #集群轨道cluster_track <- RSE %>% subsetByOverlaps(plot_region) %>% trackClusters(name = "Clusters"， col = NA, showId = TRUE) #>设置粗细特征…#>合并排序…#>准备轨道…#链接轨道link_track <- cor_links %>% subsetByOverlaps(plot_region) %>% trackLinks(name = "Links"，交互。Measure = "p.value"， interactive .dimension.transform = "log"， col.outside = "grey"， plot。anchors = FALSE, col.interactions = "black") # Plot tracks together (list(axis_track, link_track, cluster_track, tx_track)， from = start(plot_region)， to = end(plot_region)，染色体= as.character(seqnames(plot_region)))

图7:tss增强子的基因组浏览器示例链接

最上面的轨迹显示了3个TSS候选基因围绕Atp1b1基因的相关性。最显著的相关性见于上游TSS候选基因和最远端增强子候选基因之间。

对于以这种方式计算TSS和增强子之间的相关性，需要注意的是:当我们计算来自两种条件(Ctrl和Nano)的生物复制之间的表达相关性时，高相关性也可能出现在TSS和增强子候选对处理的响应方向相同时。这意味着在不同条件下组合所有样本时观察到的相关性可能不同于在每个条件下计算的相关性(这种不直观的现象被称为辛普森悖论)。为了避免包含这样的情况，可以分别分析每个条件，以在每个状态中找到tss增强器链接。作为这种方法的扩展，人们还可以查看tss增强子链接，在不同状态下显示不同的相关性强度。这种类型的分析被称为微分共表达式分析．

3.2.5寻找增强剂的拉伸

一些研究发现，与单例增强子相比，紧密间隔的增强子组或延伸倾向于显示不同的染色质特征和功能。这类增强子通常被称为“超级增强子”或“拉伸增强子”。(Pott and Lieb 2015)．

CAGEfightR可以检测到增强器延伸基于CAGE数据。CAGEfightR将附近的增强子分组，并计算它们之间的平均成对相关性，如下所示(再次使用Kendall 's tau):

#子集到只增强子增强子<-子集(RSE, clusterType == "Enhancer") #查找12.5 kbp范围内的拉伸<- find拉伸(enhancers, inputAssay = "TPM"， mergeDist = 12500L, minSize = 5L, method = "kendall") #>查找拉伸…计算相关性…#> #拉伸汇总:#>拉伸数量:95 #>拉伸内部的集群总数:587 / 9943 #>最小集群:5 #>最大集群:15 #>最小宽度:7147 #>最大宽度:92531 #>平均对相关性的总结:#>最小第一曲。中位数。平均第三曲。最大。#> -0.10038 0.01351 0.08107 0.09097 0.16171 0.37105

类似于tss和增强子，我们也可以根据它们与已知转录本的关系来注释拉伸:

#注释转录模型拉伸<- assignTxType(拉伸，txModels = txdb) #>查找层次重叠…# > # # #重叠简介:# # # # > txType数百分比# > 1启动子50 52.6 # 0.0 # > > 2近端0 3 fiveUTR 6 # 6.3 # > 4 threeUTR 5 5.3 > 3.2 5 cd 3 # > 6外显子2 2.1 # > 7基因内区15 # # > > 8反义0 0.0 15.8 14.7 9基因间14 #按相关性排序延伸< -延伸(顺序(延伸aveCor美元,减少= TRUE)] #显示结果延伸# >农庄与95范围和4元数据对象列:#> seqnames ranges | #>    b| #> chr11 98628005-98647506 * | #> chr7:139979437-140003112 chr7 139979437-140003112 * | #> chr15:31261340-31279984 chr15 31261340-31279984 * | #> chr11 117733009-117752208 * | #> chr7:97167988-97188451 * | #> ... ... ... ... .#> chr16:91373912-91399202 chr16 91373912-91399202 * | #> chr7:132619265-132644381 * | #> chr15:79181690-79208915 chr15 79181690-79208915 * | #> chr10:94708643-94729408 chr10 94708643-94729408 * | #> revmap nClusters aveCor txType #>     #> chr11:98628005-98647506 6600,6601,6602，…6 0.371053启动子#> chr7:139979437-140003112 4220,4221,4222，…5 0.328631 promoter #> chr15:31261340-31279984 7962,7963,7964，…5 0.301604内含子#> chr11:117733009-117752208 6785,6786,6787，…6 0.284399启动子#> chr7:97167988-97188451 4022,4023,4024，…6 0.262200启动器#> ... ... ... ... ...#> chr15:101076561-101093429 8320,8321,8322，… 5 -0.0549688 intergenic #> chr16:91373912-91399202 8643,8644,8645,... 7 -0.0598361 fiveUTR #> chr7:132619265-132644381 4160,4161,4162,... 5 -0.0626249 promoter #> chr15:79181690-79208915 8144,8145,8146,... 5 -0.0981772 promoter #> chr10:94708643-94729408 5823,5824,5825,... 5 -0.1003807 intron #> ------- #> seqinfo: 35 sequences (1 circular) from mm9 genome

返回的农庄包含每个拉伸的位置、增强器的数量和平均相关性。拉伸存在于各种环境中，一些是基因间的，另一些横跨基因的各个部分。让我们画出一个顶级的基因间延伸:

#提取一段增强器周围的区域plot_region <- extends ["chr17:26666593-26675486"] + 1000 #集群轨道cluster_track <- RSE %>% subsetByOverlaps(plot_region) %>% trackClusters(name = "Clusters"， col = NA, showId = TRUE) #>设置粗瘦特性…#>合并排序…#>准备轨道…# CTSS track ctss_track <- CTSSs %>% subsetByOverlaps(plot_region) %>% trackCTSS(name = "CTSSs") #>按链分割池信号…#>准备轨道…#拉伸增强器跟踪stretch_track <-拉伸%>% subsetByOverlaps(plot_region) %>% AnnotationTrack(name = "拉伸"，fill = "hotpink"， col = NULL) # Plot tracks together (list(axis_track, stretch_track, cluster_track, ctss_track)， from = start(plot_region)， to = end(plot_region)，染色体= as.character(seqnames(plot_region)))

图8:基因组浏览器增强子拉伸的例子

该片段由至少5个候选增强子组成，每个增强子都显示双向转录。

3.3.1表达规范化和EDA

在对差分表达式(DE)的各种度量进行统计测试之前，重要的是首先进行彻底的探索性数据分析(EDA)，以确定我们需要在最终的DE模型中包括哪些因素。

这里我们将使用DESeq2(Love, Huber, and Anders 2014)用于规范化和EDA，因为它提供了易于使用的函数来执行基本分析。其他流行的工具如刨边机(罗宾逊，麦卡锡，和史密斯2010)而且limma(Ritchie et al. 2015)提供类似的功能，以及更专门的EDA包，如EDASeq．

DESeq2以方差稳定转换的形式提供复杂的计数数据规范化和转换:这在日志转换之前向规范化表达式值添加动态伪计数，以抑制计数数据固有的均值-方差关系。这对于CAGE数据特别有用，因为CAGE甚至可以检测非常低表达的tss和增强子。

注意:由于增强子候选标签整体表达量较低，因此仅占标签总数的一小部分。由于适当的归一化因子估计假定了大量的非DE特征，TSS和增强子候选程序通常都应该包括在DE分析中。

首先，我们拟合方差稳定变换的“盲”版本，因为我们还不知道什么设计适合这个特定的研究:

#用空白设计创建DESeq2对象dds_blind <- DESeqDataSet(RSE, design = ~ 1) #规范化和日志转换vst_blind <- vst(dds_blind, blind = TRUE)

一个非常有用的第一个表示是主成分分析(PCA)图，将整个实验的方差汇总为主成分(PCs):

plotPCA (vst_blind,“阶级”)

图9:方差稳定表达的pca图

我们观察到PC2按照实验组(Nano vs Ctrl)对样品进行分离。然而，PC1也将样本分为两组。这暗示了对表达的一种不必要的但系统的影响，通常被称为a批处理效果．批次效应产生的原因有很多，例如不同的实验室或人员制备的文库或使用稍有不同的试剂池。通常，批处理效应与准备的日期库和Bornholdt同时发生等暗示这是原始研究中批量效应的原因。

当我们测试DE时，我们不想将这种不需要的变异误认为是生物变异。为了防止这种情况，我们可以将批处理效应作为最终DE模型的一个因素。首先，我们定义批处理变量:

#提取PCA结果pca_res <- plotPCA(vst_blind， "Class"， returnData = TRUE) #使用PC1 batch_var <- ifelse(pca_res$PC1 > 0， " a "， "B") #将批处理变量作为因子附加到实验RSE$ batch <- factor(batch_var) #显示新设计RSE %>% colData() %>%子集(select = c(Class, batch)) %>% knitr::kable(标题= "添加新的批处理共变量后的设计矩阵。")

作为手动定义批处理变量的替代方法，工具如股东价值分析而且RUVSeq可用于从数据中估计未知批效应。

3.3.2集群级差异表达式

按照我们上面的简短EDA，我们已经准备好指定实验的最终设计:我们希望同时考虑样本的类别和批次:

#指定设计dds <- DESeqDataSet(RSE, design = ~ Batch + Class) #拟合DESeq2模型dds <- DESeq(dds) #>估计大小因子#>估计分散度#>基因分散估计#>均值-分散关系#>最终分散估计#>拟合模型和测试

现在，我们可以为nano和ctrl的比较提取估计的影响(对数折叠变化)和统计显著性(p值)，隐式地校正批量效应:

#提取结果res <- results(dds, contrast = c("Class"， "Nano"， "Ctrl")， alpha = 0.05, independentFiltering = TRUE, tidy = TRUE) %>% bind_cols(as.data.frame(rowData(RSE))) %>% as_tibble() #显示top hits res %>% top_n(n = -10, wt = padj) %>% dplyr::select(cluster = row, clusterType, txType, baseMean, log2FoldChange, padj) %>% knitr::kable(标题= "顶级差异表达TSS和增强子候选")

检查几块诊断图，以确定是否正确总是一个好主意DESeq2分析很成功。一个这样的例子是ma图(另一个有用的图是p值直方图):

ggplot(res, aes(x = log2(baseMean)， y = log2FoldChange, color = padj < 0.05)) + geom_point(alpha = 0.25) + geom_hline(yintercept = 0, linetype = "虚线"，alpha = 0.75) + facet_grid(clusterType ~ .))

图10:差异表达分析的诊断MA图

我们可以看到，与增强子候选基因相比，我们总体上发现了更多的DE TSS，这是意料之中的，因为TSS候选基因的表达水平更高。许多增强考生在决赛中被淘汰DESeq2分析(“独立过滤”步骤)，因为它们的表达水平太低，无法检测任何DE:这增加了剩余TSS和增强子候选在较高表达水平检测DE的能力。

我们可以将DE TSS和增强子候选的总数制成表格:

table(clusterType = rowRanges(RSE)$clusterType, DE = res$padj < 0.05) #> DE #> clusterType FALSE TRUE #> TSS 22071 6385 #>增强器3034 199

3.3.3校正批处理效果的表达式估计

除了查看每个聚类的估计值和显著性外，我们还可能希望查看一些顶级命中值的单个表达式值。然而，我们还需要纠正表达式估计本身的批处理效果，就像我们对对数折叠变化和p值所做的那样(当然使用相同的模型)。

这里我们使用ComBat(Johnson, Li, and Rabinovic 2007)从股东价值分析包，适用于从小型实验中去除简单的批量效应。对于更高级的设置，removeBatchEffect从limma可以删除任意复杂的批处理效果。的RUVSeq包和fsva从股东价值分析可用于删除未知批处理效果。

我们再次使用方差稳定变换来准备数据战斗(这使得计数数据类似于从微阵列获得的表达式估计，因为ComBat最初是为微阵列开发的):

#引导/非盲方差稳定转换vst_guided <- variancestablizingtransformation (dds, blind = FALSE)

运行战斗我们还需要两个额外的信息:i)描述生物效应或预期效应的设计矩阵;ii)已知但不需要的批量效应。我们首先指定设计矩阵，然后运行战斗：

#设计所需效果的矩阵bio_effects <-模型。matrix(~ Class, data = colData(RSE)) #运行ComBat assay(RSE， "ComBat") <- ComBat(dat = assay(vst_guided)， batch = RSE$ batch, mod = bio_effects) #>在单个批次中发现了253个表达一致的基因(全为零);这些将不会针对批次进行调整。#>发现2批次#>调整1个协变量或协变量水平#>跨基因标准化数据#>拟合L/ s模型并寻找先验#>寻找参数调整#>调整数据

我们可以重做pca图，看看批量效应修正的全局效应:

#覆盖实验分析(vst_guided) <- assay(RSE， "ComBat") # Plot as as before plotPCA(vst_guided， "Class")

图11:批量校正表达式的pca图

现在Nano和Ctrl沿着PC1分开(修正前是PC2)。由于PC1捕获了最多的方差，这表明批量效应已被消除，实验组现在是表达方差的主要贡献者。

然后，我们使用以下方法提取前10个DE增强程序候选tidyverse代码:

#查找前10位DE增强器top10 <- res %>% filter(clusterType == "Enhancer"， padj < 0.05) %>% group_by(log2FoldChange >= 0) %>% top_n(n = 5, wt = abs(log2FoldChange)) %>% pull(row) #以tidy格式提取表达式值tidyEnhancers <- assay(RSE， "ComBat")[top10，] %>% t() %>% as_tibble(rownames = "Sample") %>% mutate(Class = RSE$Class) %>% gather(key = "Enhancer"， value = " expression "， -Sample， -Class, factor_key = TRUE)

最后，我们可以绘制出每个顶级增强子候选的批量修正表达式:

ggplot(tidyEnhancers, aes(x = Class, y = Expression, fill = Class)) + geom_dotplot(stackdir = "center"， binaxis = "y"， dotsize = 3) + facet_wrap(~ Enhancer, ncol = 2, scales = "free_y") #> Bin width默认为数据范围的1/30。用' binwidth '选择更好的值。

图12:前10个差异表达增强子候选表达谱

3.3.4dna结合基序的富集

在鉴别差异表达的TSS和增强子候选之后，一个典型的问题是转录因子(tf)可能参与了它们的调控。为了阐明这个问题，我们可以用JASPAR数据库中的dna结合基元注释tss和增强子(Mathelier et al. 2016)．

首先我们提取tss和增强子周围的序列。这里我们简单定义为TSS候选峰值或增强子候选中点附近+/- 500 bp:

cluster_seqs <- RSE %>% rowRanges() %>% swapRanges() %>% unstrand() %>% add(500) %>% getSeq(bsg, names = .)

其次，我们使用TFBSTools(Tan and Lenhard 2016)包来获取作为位置频率矩阵(PFMs)的基元JASPAR2016数据库:

#提取motif作为PFMs motif_pfms <- getMatrixSet(JASPAR2016, opts = list(species = "10090")) #查看前几个motifs的id和名称:head(name(motif_pfms)) #> MA0004.1 MA0006.1 MA0029.1 MA0063.1 MA0067.1 MA0078.1 #> "Arnt" "Ahr::Arnt" "Mecom" "Nkx2-5" "Pax2" "Sox17"

第三，我们使用motifmatchr包(Schep 2018)要在启动子序列中找到匹配项:

#查找匹配motif_hits <- matchMotifs(motif_pfms, subject = cluster_seqs) #匹配以稀疏矩阵形式返回:motifMatches(motif_hits)[1:5, 1:5] #> 5 x 5稀疏矩阵类“lgCMatrix”#> MA0004.1 MA0006.1 MA0029.1 MA0063.1 MA0067.1 #> [1，] . . . .| #> [2，] . . . . .#>[3，] |。|。#> [4，] . . . . .#>[5，]。|。|。

最后，我们可以做一个简单的Fisher’s Exact测试，看看一个motif与DE TSS和增强子候选的共现次数是否比我们预期的要多。在这里，我们将看看FOS::JUN motif (MA0099.2):

# 2x2 table for Fisher table(FOSJUN = motifMatches(motif_hits)[，"MA0099.2"]， DE = res$padj < 0.05) %>% print() %>% Fisher . Test () #> DE #> FOSJUN FALSE TRUE #> FALSE 22144 5596 #> TRUE 2961 988 #> #> #> Data:。#> p-value = 5.839e-12 #>替代假设:真实优势比不等于1 #> 95%置信区间:#> 1.220330 1.427821 #>样本估计:#>优势比#> 1.320361

高于1的显著优势比表明FOS::JUN是调控纳米管响应的候选TF(或者，更准确地说，是候选TF二聚体)。这并不奇怪，因为FOS::JUN是TNF-alpha炎症通路的一部分(见下文)。Bornholdt也发现了类似的主题等在最初的研究中。

当然，这只是一个非常快速和简单的分析主题丰富。人们可以很容易地使用TSS和增强子候选体周围的不同区域和/或在TSS和增强子之间拆分富集分析。其他Bioconductor包PWMEnrich，rGADEM而且motifcounter实现更先进的统计方法，计算丰富的已知主题。rGADEM，BCRANK而且motifRG也可以用来计算新母题的丰富性，有时被称为主题发现．

3.3.5基因水平差异表达

虽然CAGE数据自然是在聚类水平(TSS和增强子候选)进行分析，但在许多情况下，也有兴趣研究基因水平的表达估计。一个主要的例子是基因本体论(GO)和京都基因和基因组百科全书(KEGG)术语的丰富(Ashburner et al. 2000;基因本体联盟2019;Kanehisa and Goto 2000)它们只在基因水平上被定义。CAGEfightR包括用基因模型注释聚类和总结基因级表达的功能。

我们可以用与转录本id相同的方式注释基因id聚类:

RSE <- assignGeneID(RSE, geneModels = txdb) #>提取基因…#> 84个基因被删除，因为它们在同一参考序列的两条链#>上或在多个参考序列上都有外显子，#>因此不能由单个基因组范围表示。使用'single.strand.genes。only=FALSE'获取一个#> GRangesList对象中的所有基因，或使用suppressMessages()来抑制该消息。考虑到最近TSS的距离时重叠…#>发现层次重叠…#>重叠总结:### #>特征重叠基因:81.34% #>唯一基因数量:13760

然后使用CAGEfightR计算基因内TSS候选基因的总数:

GSE <- RSE %>%子集(clusterType == "TSS") %>% quantifyGenes(genes = "geneID"， inputAssay = "counts") #>量化基因:13567

结果是RangedSummarizedExperiment取值范围是aGRangesList持有每个基因内的TSS候选基因:

#> seqnames ranges strand | score #>    |  #> chr7:80884953-80885056;+ chr7 80884953-80885056 + | 11.0585 #> chr7:80885120-80885677;+ chr7 80885120-80885677 + | 1162.3447 #> thick support txID txType balance #>      #> chr7:80884953-80885056;+ 80885000 5 uc009hrf。2近端NA #> chr7:80885120-80885677;+ 80885256 11 uc009hrf。2启动子N一个#> bidirectionality clusterType geneID #>    #> chr7:80884953-80885056;+ NA TSS 100038347 #> chr7:80885120-80885677;+ NA TSS 100038347 #> ------- #> seqinfo: 35 sequences (1 circular) from mm9 genome

本例中的基因id是Entrez id (Bioconductor软件包广泛使用)。可以将这些系统id转换为更便于人类阅读的符号org.Mm.eg.db注释包:

#翻译符号rowData(GSE)$symbol <- mapIds(odb, keys = rownames(GSE)， column = " symbol "， keytype = "ENTREZID") #> 'select()'返回键和列之间的1:1映射

获得了基因级计数矩阵后，我们现在可以进行基因级DE分析。这里我们使用limma- boom，因为limma便于执行后续的富集分析。其他工具，如DESeq2(上图)或刨边机(见下文)也可以使用。

请注意：limma是对基于计数的数据进行DE分析的强大工具。然而，由于它依赖于日志转换计数，它并不总是适合分析具有非常低计数的特征的数据集。对于基因水平的分析来说，这通常不是问题，但对于增强子来说，这可能是一个问题，增强子通常表达非常低。

类似于DESeq2在上面的分析中，我们首先构建必要的对象，然后规范化表达式值:

#创建DGElist对象dge <- DGElist (counts = assay(GSE， "counts")， genes = as.data.frame(rowData(GSE))) #计算归一化因子dge <- calcNormFactors(dge)

然后我们应用vom -transform来模拟均值-方差趋势，为此我们还需要指定设计矩阵(在这种情况下，设计必须同时包含想要的和不想要的效果!)然后使用相同的设计矩阵拟合基因模型:

#设计矩阵mod <- model。矩阵(~ Batch + Class, data = colData(GSE)) #使用voom v <- voom(dge, design = mod) #拟合和收缩DE模型拟合<- lmFit(v, design = mod) eb <- eBayes(Fit, robust = TRUE) #总结结果dt <- decideTests(eb)

然后，我们可以报告差异基因表达的总体总结，并查看前几个热门:

#全局总结dt %>%总结()%>% knitr::kable(标题=“差异表达基因的全局总结。”)

# Inspect top htis topTable(eb, coef = "ClassNano") %>% dplyr::select(symbol, nClusters, AveExpr, logFC, adj.P.Val) %>% knitr::kable(标题= "顶级差异表达基因。")

3.3.6丰富GO-和kegg -术语

除了查看单个顶级基因，我们还可以查看DE基因与已知基因功能数据库的关系，以深入了解实验中可能影响的生物过程。

limma使得在DE分析之后很容易进行这样的富集分析。由于我们有被Entrez id索引的基因，我们可以直接使用goana找到丰富的go术语:goana使用偏置瓮模型来估计GO-terms的富集，同时考虑DE基因的表达水平:

#查找丰富的GO术语GO <- goana(eb, coef = "ClassNano"， species = "Mm"， trend = TRUE) #显示top hits topGO(GO，本体= "BP"， number = 10) %>% knitr::kable(标题= "顶级丰富或耗尽的GO术语。")

同理，我们可以用KEGG术语kegga：

#找到丰富的KEGG-术语KEGG <- kegga(eb, coef = "ClassNano"，物种= "Mm"，趋势= TRUE) #显示顶部点击topKEGG(KEGG, number = 10) %>% knitr::kable(标题= "顶级丰富或枯竭的KEGG术语。")

这两项分析都表明，与炎症反应和防御反应相关的基因在纳米管暴露后上调(Bornholdt也发现了类似的富集等在最初的研究中)。这支持了纳米管引起类似石棉反应的假设。

kegg术语表示定义良好的路径。我们可以用pathview包(Luo and Brouwer 2013)为了更详细地研究特定富集途径中的基因。例如，我们可以观察TNF信号通路中的基因调控:

# Format DE genes for pathview DE_genes <- dt[， "ClassNano"] %>% as.integer() %>% set_names(rownames(dt)) %>% Filter(f=function(x) x != 0, x=.) #运行pathview;这将保存一个png文件到临时目录路径视图(DE_genes, species = "mmu"，路径。Id = "mmu04668"， kegg。dir = tempdir()) #显示png文件grid.newpage() grid. rster (png::readPNG("mmu04668.pathview.png"))

图13:KEGG TNF信号通路中差异表达基因的详细视图

3.3.7TSS使用差异

在上面的两个分析中，我们研究了单个候选TSS或单个基因是否在不同的实验条件下改变了表达。然而，我们可能还想看看一个基因是否表现出差异TSS使用(DTU):一个基因是否在不同的条件下使用不同的TSS候选。这个问题类似于RNA-Seq中的差异剪接，但研究的是tss而不是转录本/异构体(Soneson, Love, and Robinson 2015)．这里我们将使用刨边机diffSpliceDGE方法来测试DTU，尽管可以使用许多其他包，例如diffSplice从limma，DEXSeq，DRIMSeq等. .

直观地说,diffSpliceDGE测试一个给定的候选TSS是否在同一基因中表现出与其他候选TSS相同的变化，这表明候选TSS在整个基因中受到不同的调控。然而，这并没有考虑给定TSS候选基因的相对组成，例如，TSS候选基因的贡献是占总基因表达的1%-2%还是25%-50%。因此，一个有用的预处理步骤是在分析之前过滤掉对总基因表达贡献很小的候选TSS。

我们使用CAGEfightR去除在5个以上样本中表达不超过总基因表达10%的候选TSS(我们首先删除未分配到基因的候选TSS):

#过滤掉低表达的RSE_filtered <- RSE %>%子集(clusterType == "TSS" & !is.na(geneID)) %>% subsetByComposition(inputAssay = "counts"， genes = "geneID"， unexpression = 0.1, minSamples = 5L) #>计算组成…# >构造子集……24500个区域中删除了8001个(32.7%)

我们只能在有多个tss的基因中检测DTU。因此，有用的第一个可视化方法是看看有多少基因使用了一个以上的候选TSS:

RSE_filtered %>% rowData() %>% as.data.frame() %>% as_tibble() %>% dplyr::count(geneID) %>% ggplot(aes(x = n, fill = n >= 2)) + geom_bar(alpha = 0.75) + scale_fill_colorblind(" multitss ") +实验室(x = "每个基因的tss数量"，y = "频率")

图14:基因内替代TSS使用概述

虽然大多数基因只使用单一的候选TSS，但许多基因使用两个或更多的候选TSS。

同样，我们从构建运行所需的r对象开始刨边机：

#像以前一样用符号注释:rowData(RSE_filtered)$symbol <- mapIds(odb, keys = rowData(RSE_filtered)$geneID, column = " symbol "， keytype = "ENTREZID") #> 'select()'返回键和列之间的1:1映射#提取基因信息TSS_info <- RSE_filtered %>% rowData() %>%子集(select = c(score, txType, geneID, symbol)) %>% as.data.frame() #构建DGEList dge <- DGEList(counts = assay(RSE_filtered， "counts")，基因= TSS_info)

然后，我们使用准似然方法归一化和拟合模型，包括diffSpliceDGE步骤:

#估计标准化因素dge < - calcNormFactors (dge) #估计色散和适合glm disp <——estimateDisp (dge、设计=国防部tagwise = FALSE) QLfit <——glmQLFit (disp、设计= mod,健壮的= TRUE) #应用diffSpliceDGE ds < - diffSpliceDGE (QLfit,系数=“ClassNano geneid =“geneid”)# >外显子总数:16499 # >基因总数:13563 # > 1的基因外显子:11098 # >意味着一个基因的外显子数量:1 # > Max基因的外显子数量:5

我们可以从两个层面观察DTU:单个TSS候选基因是否表现出DTU (TSS水平)或单个基因是否以任何方式表现出DTU(基因水平)。首先，让我们看看个别的TSS候选人(TSS级DTU):

dtu_TSSs <- topSpliceDGE(ds, test = "exon") dplyr::select(dtu_TSSs, txType, geneID, symbol, logFC, FDR) %>% knitr::kable(标题= "顶级差异使用TSSs")

这里对对数折叠变化的解释与以前略有不同:这些对数折叠变化相对于该基因中所有TSS候选基因的总体对数折叠变化。

然后我们可以查看每个基因的结果(基因级DTU):

dtu_genes <- topSpliceDGE(ds, test = "Simes") dplyr::select(dtu_genes, geneID, symbol, NExons, FDR) %>% knitr::kable(row.names = FALSE，标题= "显示任何差异TSS使用的顶级基因。")

我们看到这两个列表是一致的，这并不奇怪，因为基因水平的结果是通过聚合跨基因的tss水平p值得到的。

我们可以通过热图来查看每个候选TSS在Fblim1基因中的批量校正表达(见上文):

RSE_filtered %>%子集(geneID == "74202") %>%化验("ComBat") %>% t() %>% pheatmap(color =岩浆(100)，cluster_cols = FALSE, main="Fblim1")

图15:热图显示Fblim1中tss的表达

Fblim1有3个TSS候选:2个用于Ctrl样本，而Nano样本也使用chr4:141154044-141154185;- TSS(也可以在上面的TSS级别表中看到)。虽然热图对于查看表达变化很有用，但基因组浏览器视图更适合检查每个TSS候选的基因组上下文:

#提取Fblim1 plot_region区域<- RSE_filtered %>% rowRanges() %>%子集(geneID == "74202") %>% GenomicRanges::reduce(min。gapwidth = 1e6L) %>% unstrand() %>% add(5e3L) #集群轨道cluster_track <- RSE_filtered %>% subsetByOverlaps(plot_region) %>% trackClusters(name = "Clusters"， col = NA, showId = TRUE) #>设置厚薄特性…#>合并排序…#>准备轨道…# CTSS跟踪每个组ctrl_track <- CTSSs %>%子集(select = Class == "Ctrl") %>% calcPooled() %>% subsetByOverlaps(plot_region) %>% trackCTSS(name = "Ctrl") #>警告在calcPooled(.):对象已经在rowData中有一个名为score的列:它#>将被覆盖!#>拆分池中的信号…#>准备轨道…nano_track <- CTSSs %>%子集(select = Class == "Nano") %>% calcPooled() %>% subsetByOverlaps(plot_region) %>% trackCTSS(name = "Nano") #>警告在calcPooled(.):对象已经在rowData中有一个名为score的列:它#>将被覆盖!#>拆分池中的信号…#>准备轨道… # Plot tracks together plotTracks(list(axis_track, tx_track, cluster_track, Ctrl=ctrl_track, nano_track), from = start(plot_region), to = end(plot_region), chromosome = as.character(seqnames(plot_region)))

图16:Fblim1中差异TSS使用的基因组浏览器示例

Fblim1基因使用了两个注释的TSS候选基因，但Nano样本也使用了一个新的内含子TSS候选基因。