作者:顾祖光(z.gu@dkfz.de）
日期: 2018-03-13
通常基因组信号被表示为数值(如WGBS的甲基化率或ChIP-seq的组蛋白修饰强度)或二进制值(如CpG岛在给定位置的存在)，然而，基因组信号有时也可以表示为分类值或离散值。一个非常典型的例子是染色质状态分割ChromHMM．ChromHMM基本上将染色质状态(例如主动转录状态或抑制转录状态)分配给基因组中的一个给定窗口，并且不同窗口中的状态分配总是相互排斥的(这意味着一个窗口只能有一种状态)。在这个小插图中，我们将演示如何可视化围绕某些基因组靶标的各种染色质状态的富集，以及如何与其他表观基因组信号整合。
在下面的例子中，我们使用路线图的数据集作为示例数据集。首先，我们展示了使用一个样本(样本id: E003，胚胎干细胞，H1细胞系)的基本用法。
负载包:
库(genome . ranges)库(data.table)库(EnrichedHeatmap)库(圈)染色质状态分割的训练和应用五种组蛋白修饰．下面正在使用的文件来自http://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/coreMarks/jointModel/final/E003_15_coreMarks_mnemonics.bed.gz．
states_bed = fread("zcat ~/EnrichedHeatmap_test/ e003_15_cocomments_mnemonics .bed.gz") states = GRanges(seqnames = states_bed[[1]]， ranges = IRanges(states_bed[[2]] + 1, states_bed[[3]])， states = states_bed[[4]]) unique(states_bed[[4]])# #[1]“十五_quies”“9 _het”“14 _reprpcwk”“13 _reprpc”“10 _tssbiv”“7 _enh”# #[7]“1 _tssa”“11 _bivflnk”“2 _tssaflnk”“5 _txwk”“4 _tx”“8 _znf / rpt”# #[13]“6 _enhg”“十二_enhbiv”“3 _txflnk”在这15个州中，有一些州彼此相似，比如1 _tssa(活动TSS)和2 _tssaflnk(侧翼主动TSS)。为了降低分析的复杂性，我们合并了一些相似的状态。
map = c("1_TssA" = "TssActive"， "2_TssAFlnk" = "TssActive"， "3_TxFlnk" = "Transcript"， "4_Tx" = "Transcript"， "5_TxWk" = "Transcript"， "6_EnhG" = "增强剂"，"7_Enh" = "增强剂"，"8_ZNF/Rpts" = "异染色质"，"9_Het" = "异染色质"，"10_TssBiv" = "TssBivalent"， "11_BivFlnk" = "TssBivalent"， "12_EnhBiv" = "增强剂"，"13_ReprPC" = " repression "， "14_ReprPCWk" = " repression "， "15_Quies" = "Quiescent") states$states_simplified = map[states$states]我们还设置了7个合并状态的颜色。
states_col = c("TssActive" = "Red"， "Transcript" = "Green"， "Enhancer" = "Yellow"， "Heterochromatin" = "PaleTurquoise"， "TssBivalent" = "Orange"， " repression " = "Grey"， "Quiescent" = "black")接下来，我们将演示如何将染色质状态正常化为TSS。首先加载转录组并提取TSS区域。转录组注释来自http://egg2.wustl.edu/roadmap/data/byDataType/rna/annotations/gen10.long.gtf.gz我们只使用蛋白质编码基因。转录组数据库文件由GenomicFeatures包(makeTxDbFromGFF ()功能)。
library(GenomicFeatures) txdb = loadDb("~/EnrichedHeatmap_test/gencode19_protein_coding_txdb.sqlite") g =基因(txdb) tss =启动子(g，上游= 0，下游= 1)为了减少运行时间，这里我们只使用1号染色体。分类信号的归一化与数值信号的归一化基本相同。这里我们只需要为分类信号指定列名。
tss_chr1 = tss[seqnames(tss) == "chr1"] mat_states = normalizeToMatrix(states, tss_chr1, value_column = "states_simplified") mat_states##正常化状态到tss_chr1: ##上游5000 bp(50个窗口)##下游5000 bp(50个窗口)##包括目标区域(宽度= 1)## 2065目标区域##信号是分类的(7级)在消息的最后一行，它明确了它是7级分类信号。
归一化分类信号的实现实际上非常简单。在内部,n_states生成规范化矩阵，其中每个矩阵对应一个染色质状态，每个窗口中的值是窗口与状态(与染色质状态)重叠的比例w0或加权意味着模式)。自从ith行和jth所有矩阵中的列都对应着围绕同一个目标的同一个窗口，如果有多个状态与该窗口重叠，在对所有状态进行汇总时，将重叠分数最大的状态赋给该窗口。
为分类信号设计了一些特殊的可视化，其中在顶部注释中分别总结和显示每个状态。在这里states_col必须是命名向量，其中名称应对应于分类信号的级别。
EnrichedHeatmap(mat_states, name = "states"， col = states_col)
你可能会觉得行顺序有点乱。虽然这些信号是分类的，但在内部，它们被编码为整数。就像因素在R中存储，每个分类级别对应一个整数。如果信号列只是一个字符向量，则整数的赋值基于该字符向量的自然顺序。因此，对信号的数字编码州对象是:
Data.frame (states = unique(states$states_simplified)， value = 1:7)##状态值## 1静态1 ## 2异染色质2 ## 3抑制3 ## 4 tss双价4 ## 5增强子5 ## 6 TssActive 6 ## 7转录7如果一个窗口没有任何状态与之重叠，并且它被用白色绘制，则将0赋给该窗口。
数值编码可以通过将信号设置为因子变量来控制，因子的电平顺序控制相应的编码。在下面的代码中，我们进行转换states_simplified作为指定级别顺序的因子(注意这里的顺序states_name也反映了不同州之间的密切关系。)
Data.frame (state = states_name, value = 1:7)##状态值## 1 TssActive 1 ## 2 Transcript 2 ## 3 Enhancer 3 ## 4 Heterochromatin 4 ## 5 TssBivalent 5 ## 6 repression 6 ## 7 quiet 7states$states_simplified = factor(states$states_simplified, levels = states_name) mat_states = normalizeToMatrix(states, tss_chr1, value_column = "states_simplified") EnrichedHeatmap(mat_states, name = "states"， col = states_col)
我们建议应用层次聚类来重新排序行。在这里，数字编码特别重要，因为它影响行之间距离的计算。
EnrichedHeatmap(mat_states, name = "states"， col = states_col, cluster_rows = TRUE)
由于归一化矩阵是数值的，所以k-means聚类也可以应用。
EnrichedHeatmap(mat_states, name = "states"， col = states_col, cluster_rows = TRUE, km = 2)
仔细观察上面的热图，我们发现TSS周围的状态一直比较丰富，而在侧翼区域，状态略有不同。实际上，这表明为了增强TSS相关状态的模式，我们可以只对TSS附近的窗口应用k-means聚类。下面我们只聚类TSS的上下游1kb的矩阵子集(TSS的默认扩展为上下游5kb，归一化矩阵中的列数为100，因此从40th第60列th列为TSS上游和下游1kb的窗口)。
当然我们需要提前计算这个分区。
split = kmeans(mat_states[， 40:60]， centers = 2)$cluster EnrichedHeatmap(mat_states, name = "states"， col = states_col, cluster_rows = TRUE, split = split)
现在很高兴看到活跃TSS的基因都聚集在簇2中，而在簇1中，基因要么没有功能，要么有二价TSS功能。
类似地，我们可以看到染色质状态在基因体中是如何丰富的。由于基因体的宽度不相等，我们添加了一个额外的点图来显示基因的宽度。
g_chr1 = g[seqnames(g) == "chr1"] mat_states_2 = normalizeToMatrix(states, g_chr1, value_column = "states_simplified") EnrichedHeatmap(mat_states_2, name = "states"， col = states_col, cluster_rows = TRUE) + rowAnnotation(gene_len = row_anno_points(log10(width(g_chr1) + 1)， size = unit(1， "mm")， gp = gpar(col = "#00000040")， axis = TRUE)， width = unit(4， "cm"))
接下来我们添加TSS周围的甲基化以及相应基因的表达。表达式数据来自http://egg2.wustl.edu/roadmap/data/byDataType/rna/expression/57epigenomes.RPKM.pc.gz甲基化数据来自http://egg2.wustl.edu/roadmap/data/byDataType/dnamethylation/WGBS/FractionalMethylation.tar.gz
接下来，原始甲基化数据已被平滑bsseq包中。的旧版本生成的rds文件bsseq这就是为什么我们使用一种奇怪的方法来提取CpG位点的位置以及平滑的甲基化数据，可以使用农庄()而且getMeth ()的最新版本应用平滑时，则使用bsseq包中。
expr = read.table("~/EnrichedHeatmap_test/57epigenome . rpkm .pc.gz"， row.names = 1, header = TRUE) obj = readRDS("~/EnrichedHeatmap_test/chr1_roadmap_merged_bsseq.rds") meth = obj@rowData meth_mat = obj@trans(obj@assays$data$coef) mcols(meth) = meth_mat .gz我们取存在于tss_chr1而且expr．
Names (tss_chr1) = gsub("\\。"\\d+$"， ""， names(tss_chr1)) cn = intersect(names(tss_chr1)， rownames(expr)) tss_chr1 = tss_chr1[cn] expr = expr[cn，]我们正常化染色质状态以及甲基化到染色体1上的TSS。
mat_states = normalizeToMatrix(states, tss_chr1, value_column = "states_simplified") mat_meth = normalizeToMatrix(meth, tss_chr1, value_column = "E003"， mean_mode = "absolute"， smooth = TRUE)同样，在TSS上游和下游1kb的染色质状态上应用k-means。
分= kmeans (mat_states(40:6 0),中心= 2)$集群meth_col_fun = colorRamp2 (c(0、0.5、1)、c(“蓝色”,“白色”,“红”))EnrichedHeatmap (mat_states name =“州”,上校= states_col cluster_rows = TRUE, =分裂,分裂top_annotation = HeatmapAnnotation(丰富= anno_enriched (gp = gpar (lty = 1:2)))) + EnrichedHeatmap (mat_meth name =“冰毒”,上校= meth_col_fun top_annotation = HeatmapAnnotation(丰富= anno_enriched (gp = gpar (lty = 1:2)))) +热图(log2 (expr[名称(tss_chr1),"E003"] + 1)， name = "expr"， show_row_names = FALSE, width = unit(5， "mm"))
与更多的表观基因组信号相对应，我们可以看到具有活跃TSS状态的启动子总是与低甲基化和高基因表达相关，而标记为二价状态的启动子通常表达较低。
认为在胚胎干细胞(ESC)中，有大量的基因在启动子上显示出二价状态。它们被称为二价状态，因为H3K4me3是活跃转录的组蛋白标记H3K27me3是抑制转录的组蛋白标记。活性组蛋白标记或抑制性组蛋白标记的丧失/获得对组织发育和分化至关重要，这些活性组蛋白标记或抑制性组蛋白标记转化为成熟细胞。
为了证明这一点，我们分析了一个成熟组织(肺组织，路线图样本id: E096)的这种转变。
states_bed = fread("zcat ~/EnrichedHeatmap_test/ e096_15_coremarks_mnemontics .bed.gz") states_lung = GRanges(seqnames = states_bed[[1]]， ranges = IRanges(states_bed[[2]] + 1, states_bed[[3]])， states = states_bed[[4]]) states_lung$states_simplified = factor(map[states_lung$states]， levels = states_name)ChromHMM的窗口为200bp，因此我们将整个基因组以200bp的窗口进行分割，并为每个窗口分配相应的状态，ESC和肺都是如此。
window = makeWindows(states, w = 200) # makeWindow()来自EnrichedHeatmap包mtch = as。matrix(finoverllaps (window, states)) window$ESC_states = states$states_simplified[mtch[， 2]] mtch = as。matrix(finoverlaps (window, states_lung)) window$lung_states = states_lung$states_simplified[mtch[， 2]]我们只使用带有注释的窗口TssBivalentESC或肺部的状态。
window = window[window$ESC_states == "TssBivalent" | window$lung_states == "TssBivalent"]我们首先制作一个和弦图来显示转换是如何发生的。
transition_mat = table(mcols(window)[， c("ESC_states"， "lung_states")]) * 200 class(transition_mat) = "matrix" transition_mat = transition_mat[states_name, states_name] transition_mat .TssActive转录增强器异染色质TssActive双价抑制性静默## TssActive 000 0 136200 00 ##转录000 00 8200 00 ##增强器000 00 38800 00 ##异染色质000 00 55800 00 ## TssActive双价4485600 665400 3711800 653000 406000 1799200 766200 ##抑制性000 000 00 23800 000转换矩阵中的值表示有多少碱基对从一种染色质状态转变为另一种状态。
行和列来自不同的单元格，我们为行和列名添加不同的前缀。
rownames(transition_mat) = paste0("ESC_"， rownames(transition_mat)) colnames(transition_mat) = paste0("lung_"， colnames(transition_mat)))现在我们用circlize包中。
circos.par(差距。After = c(rep(1, n_states - 1)， 5, rep(1, n_states - 1)， 5))网格。col = c(states_col, states_col) names(grid.col) = c(rownames(transition_mat)， colnames(transition_mat)) chordDiagram(transition_mat, grid.col)Col =网格。col, annotationTrack = c("grid"， "axis")， directional = TRUE) circos.clear() text(0.5， -1， "ESC") text(0.5, 1， "lung") legend("left"， pch = 15, col = states_col, legend = names(states_col))
从Chord图中我们可以看到，ESC中很多具有双价tss状态的区域在肺中已经过渡到活跃或抑制状态。
为了更深入地了解这种转变是怎样的，我们可以制作丰富的热图来观察不同表观基因组信号之间的联系。
自TssBivalent状态基本上是TSS相关区域的状态，我们只保留在TSS上游和下游1kb内必须有窗口的基因TssBivalent带注释的状态。同样，我们只使用1号染色体进行演示。
mat_bivtss = normalizeToMatrix(states[states$states_simplified == "TssBivalent"]， tss_chr1) l = rowsum (mat_bivtss[， 40:60]) > 0 # recall it is 1kb upstream and downstream tss_biv = tss_chr1[l] tss_biv# #农庄与443范围和对象元数据列:# # seqnames范围链| gene_id # # < Rle > < IRanges > < Rle > | <人物> # # ENSG00000000938 chr1(27961788, 27961788)——| ENSG00000000938.8 # # ENSG00000006555 chr1(55266940, 55266940)——| ENSG00000006555.6 # # ENSG00000007968 chr1(23857712, 23857712)——| ENSG00000007968.6 # # ENSG00000008118 chr1(209757062、209757062)+ | ENSG00000008118.5 # # ENSG00000009709 chr1(18957500、18957500)+ | ENSG00000009709.7 ## ... ... ... ... . ...## ENSG00000243749 chr1 [35450954, 35450954] - | ENSG00000243749.1 ## ENSG00000248485 chr1 [161228517,161228517] + | ENSG00000248485.1 ## ENSG00000249087 chr1 [23695711, 23695711] + | ENSG00000249087.3 ## ENSG00000251246 chr1 [155036224, 155036224] + | ENSG00000251246.1 ## ENSG00000253304 chr1 [29450447, 29450447] - | ENSG00000253304.1 ## ------- # seqinfo:来自未指定基因组的25个序列;没有seqlengths我们将染色质状态和甲基化正常化到tss_biv．
mat_states_ESC = normalizeToMatrix(states, tss_biv, value_column = "states_simplified") mat_states_lung = normalizeToMatrix(states_lung, tss_biv, value_column = "E003"， mean_mode = "absolute"， smooth = TRUE) mat_meth_lung = normalizeToMatrix(meth, tss_biv, value_column = "E096"， mean_mode = "absolute"， smooth = TRUE)正常情况下，甲基化变化发生在低甲基化区边界。为了增强甲基化变化的效果，我们直接可视化了甲基化差异。
mat_meth_diff = mat_meth_ESC - mat_meth_lung meth_diff_col_fun = colorRamp2(c(-0.25, 0,0.25)， c("#3794bf"， "#FFFFFF"， "#df8640"))现在我们构建热图列表。
ht_list = EnrichedHeatmap (name = " states_ESC mat_states_ESC,坳= states_col top_annotation = HeatmapAnnotation(丰富= anno_enriched (gp = gpar (lty = 1:2), ylim = c (0,1))), column_title =“ESC”)+ EnrichedHeatmap (name = " states_lung ", mat_states_lung坳= states_col top_annotation = HeatmapAnnotation(丰富= anno_enriched (gp = gpar (lty = 1:2), ylim = c (0,1))), show_heatmap_legend = FALSE, column_title =“州肺”)ht_list = ht_list + EnrichedHeatmap (mat_meth_ESC,name = " meth_ESC ",坳= meth_col_fun top_annotation = HeatmapAnnotation(丰富= anno_enriched (gp = gpar (lty = 1:2), ylim = c (0,1))), column_title =“冰毒ESC”)+ EnrichedHeatmap (name = " meth_lung mat_meth_lung,坳= meth_col_fun top_annotation = HeatmapAnnotation(丰富= anno_enriched (gp = gpar (lty = 1:2), ylim = c (0,1))), show_heatmap_legend = FALSE, column_title =“冰毒肺”)+ EnrichedHeatmap (name = " meth_diff mat_meth_diff,坳= meth_diff_col_fun,top_annotation = HeatmapAnnotation(enrich = anno_enrich (gp = gpar(lty = 1:2, pos_col = "#df8640"， neg_col = "#3794bf"))))， column_title = "Meth ESC - lung")实际上我们可以离散这个数字mat_meth_diff．这里我们分配超甲基化差异大于0.2和的视窗海波到甲基化差小于-0.2的Windows。在这里离散化()通过提供区间列表，将具有连续信号的矩阵转换为分类信号。
mat_meth_diff_离散=离散(mat_meth_diff, rule = list("hypo" = c(-Inf， -0.2)， "hyper" = c(0.2, Inf))##正常化到tss_biv: ##上游5000bp(50个窗口)##下游5000bp(50个窗口)##包括目标区域(宽度= 1)## 443目标区域##信号是分类的(2级)我们将继续添加更多的热图。
ht_list = ht_list + EnrichedHeatmap(mat_meth_diff_discrete, name = "meth_diff_discrete"， col = c("hyper" = "#df8640"， "hypo" = "#3794bf")， top_annotation = HeatmapAnnotation(enrich = anno_enrich (gp = gpar(lty = 1:2)))) . txt))， ht_list = ht_list + EnrichedHeatmap(mat_meth_diff_discrete, name = "meth_diff_discrete"， col = c("hyper" = "#df8640"， "hypo" = "#3794bf")最后是基因表达的热图。
e = log2(expr[names(tss_biv)， c("E003"， "E096")] + 1) ht_list = ht_list + Heatmap(e, name = "expr"， show_row_names = FALSE, width = unit(1， "cm")， cluster_columns = FALSE)所有热图的行排序来自于合并的归一化状态矩阵上的层次聚类，该矩阵对应于ESC和肺中TSS的上游和下游1kb。所有的热图按行分为两组，一组是肺中的高基因表达，另一组是肺中的低基因表达。
row_order = hclust(dist(cbind(mat_states_ESC[， 40:60]， mat_states_lung[， 40:60])))$order split = ifelse(e[， "E096"] > e[， "E003"]， "激活"，"抑制")draw(ht_list, row_order = row_order, split = split)
还有一些基因组信号是分类的例子。
sessionInfo ()## R版本3.3.1(2016-06-21)##平台:x86_64-pc-linux-gnu(64位)##运行在:openSUSE 13.1 (Bottle) (x86_64) ## ## locale: ## [1] LC_CTYPE=en_GB。UTF-8 LC_NUMERIC=C LC_TIME=en_GB。UTF-8 ## [4] LC_COLLATE=en_GB。utf - 8 LC_MONETARY = en_GB。utf - 8 LC_MESSAGES = en_GB。UTF-8 ## [7] LC_PAPER=en_GB。UTF-8 LC_NAME=C LC_ADDRESS= c# # [10] lc_phone =C LC_MEASUREMENT=en_GB。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:##[1]网格并行stats4统计图形grDevices utils数据集方法##[10]基础## ##其他附加包:## [1] genomicfeature_1 .26.4 annotationdbiase_2.34.0 circlize_0.4.3 ## [5] EnrichedHeatmap_1.9.2 ComplexHeatmap_1.17.1 data.table_1.10.4 GenomicRanges_1.26.4 ## [9] GenomeInfoDb_1.10.3 IRanges_2.8.2 S4Vectors_0.12.2 BiocGenerics_0.20.0 ## [13] markdown_0.8 knitr_1.18 colorout_1.1-3 ## ##通过命名空间加载(并且没有附加):## [1] Rcpp_0.12.15 RColorBrewer_1.1-2 XVector_0.14.1 bitops_1.0-6 ## [7] tools_3.3.1 zlibbioc_1.20.0 biomaRt_2.30.0 ## [10] digest_0.6.15 bit_1. 1.1-12 memoise_1.1.0 ## [13] evaluate_0.10.1 RSQLite_2.0 tibble_1.4.2 ## [19] Matrix_1.2-7.1 DBI_0.8 rtracklayer_1.34.2 ## [22] string_1 .2.0 Biostrings_2.42.0 GlobalOptions_0.0.12 ## [25] locfit_1. 1.5-9.1 bit64_0.9-7 GetoptLong_0.1.8 ## [28] BiocParallel_1.8.0 xml_3 . 1.5- 1.5## [34] matrixStats_0.53.1 SummarizedExperiment_1.4.0 shape_1.4.3 ## [37] colorspace_1.3-2 stringi_1.1.2 RCurl_1.95-4.10 ## [40] rjson_0.2.15