看到可视化的包，真的停不下来，就想先翻译一边先。
APLS 全称 AnaLysis routines for ePigenomicS data, 顾名思义，是一个分析表观数据的工具。产生可用发表的可视化文件，主要针对全基因组表观基因组学数据，例如 ChIP-seq，ATAC-seq等。
功能简介：
- 计算基因组区域的富集
- 计算样品相关性
- Plot enrichments across groups
- 绘制 UCSC 基因组浏览器图
- 注释基因组区间
- 基因组注释相关可视化
我觉得这个包值得学习的是它底层的代码，有机会可以看看。
能得到什么图？
Bigwig 文件
Bigwig文件通常用于存储全基因组范围的定量数据，如ChIP-seq，ATAC-seq，WGBS，GRO-seq等。下图说明了Bigwig文件的示例。
怎么产生 Bigwig 文件
使用 UCSC 中的生成 bigwig 方法 的
deeptools软件的bamCoverage功能将BAM→bigwig。一旦生成了标准化的bigwig文件并从BAM文件中识别出Peak，就很少在后续流程中再次使用BAM文件（ 啊哈哈，也就是说是为 bw 文件而生 ）。所有下游过程的要求都可以通过标准化的bigwig文件来满足，例如，量化峰值或启动子的标准化reads数，通过基因组浏览器或IGV中的进行展示。鉴定出
Peak后，随后的步骤将是对鉴定出的Peak处的标准化reads数进行定量，以进行探索性数据分析（explorative data analysis: EDA），无监督聚类分析PCA，以识别所考虑样品中的模型并产生新颖的生物学见解。
ALPS 流程
-
ALPS包的设计理念是从最少的输入开始，并从数据中获取丰富的见解。大多数时候，ALPS中的大多数功能都只需要一个数据表格，该表具有指向bigwig文件的路径和相关的示例元信息。各种功能将利用此数据表格并生成下游输出。该包产生可发表的可视化效果，其中大部分可以使用ggplot2在R中进行自定义。 - 以下是
ALPS流程和可用功能的概述：
实战演练
- 安装包，
注意：R3.6.1
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ALPS")
- 计算基因组区域的富集
- 表观基因组学中的大多数探索性分析都始于对一组基因组区域的富集或甲基化进行量化，比如:
启动子区域或Peak 区域。量化后的文件将用作下游分析（ 例如PCA，聚类）的输入。multiBigwig_summary()函数使用带有bigwig 路径和bed 文件路径的样本文件表格计算基因组区域的富集。此功能是rtracklayer包对bigwig实用程序的包装。在计算富集之前，该功能会同时根据输入数据表中存在的所有bed文件生成共同的Peak。- 从
ALPS包中读取数据表格
可以看到这一步就是得到上面所说的表格chr21_data_table
chr21_data_table <- system.file("extdata/bw", "ALPS_example_datatable.txt", package = "ALPS", mustWork = TRUE)
## attach path to bw_path and bed_path （我觉得这里用英文表述会好点，就不翻译了）
d_path <- dirname(chr21_data_table)
# dirname 返回文件 chr21_data_table 的路径部分
chr21_data_table <- read.delim(chr21_data_table, header = TRUE)
chr21_data_table$bw_path <- paste0(d_path, "/", chr21_data_table$bw_path)
chr21_data_table$bed_path <- paste0(d_path, "/", chr21_data_table$bed_path)
chr21_data_table %>% head
sample_id group color_code
1 ACCx_1 ACCx #8DD3C7
2 ACCx_2 ACCx #8DD3C7
3 BRCA_1 BRCA #FFED6F
4 BRCA_2 BRCA #FFED6F
5 GBMx_1 GBMx #B3DE69
6 GBMx_2 GBMx #B3DE69
bw_path
1 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/ACCx_025FE5F8_T1.bw
2 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/ACCx_025FE5F8_T2.bw
3 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/BRCA_000CFD9F_T2.bw
4 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/BRCA_01112370_T1.bw
5 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/GBMx_09C0DCE7_T1.bw
6 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/GBMx_09C0DCE7_T2.bw
bed_path
1 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/ACCx_1.bed
2 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/ACCx_2.bed
3 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/BRCA_1.bed
4 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/BRCA_2.bed
5 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/GBMx_1.bed
6 C:/Users/ql/Documents/R/win-library/3.6/ALPS/extdata/bw/GBMx_2.bed
然后运行
multiBigwig_summary()函数，通过同时从bed_path列中的所有bed文件准备共同的Peak来计算所有bigwig文件的富集
enrichments <- multiBigwig_summary(data_table = chr21_data_table,
summary_type = "mean",
parallel = FALSE)
enrichments %>% head
chr start end ACCx_1 ACCx_2 BRCA_1 BRCA_2 GBMx_1 GBMx_2 LGGx_1 LGGx_2
1 chr21 45639550 45640549 3.4293169 4.30173375 6.229088 6.053632 6.2633310 8.964782 8.0253673 6.3168158
2 chr21 45640994 45641493 5.1058709 4.39355741 2.412814 3.262252 6.2341950 5.569840 2.2929616 4.2814719
3 chr21 45642859 45643914 118.9869371 132.23839670 63.987536 82.130318 23.8171352 29.146589 57.2475263 50.2498071
4 chr21 45668163 45668662 0.3810360 0.02440864 4.359150 1.732686 0.6438128 1.015109 0.8307302 0.9188698
5 chr21 45689906 45690405 0.4229491 0.81362319 2.638012 0.923278 1.1011821 0.723354 2.4658742 1.0604877
6 chr21 45698386 45698910 46.5824657 54.70859597 1.141300 1.746901 24.9284784 21.264368 2.8639639 2.1893534
几乎不费吹灰之力，
multiBigwig_summary()函数的输出可以很容易地与其他R/bioconductor包装集成在一起，以进行探索性分析，PCA或聚类。
get_variable_regions()获取multiBigwig_summary或类似格式的输出，并返回n个可缩放的可变区域，可以直接与ComplexHeatmap（ 顾大佬的包真的强 ）等工具一起使用。- 以下是有关如何通过
get_variable_regions()函数将multiBigwig_summary()函数输出结果集成到ComplexHeatmap的示例
enrichments_matrix <- get_variable_regions(enrichments_df = enrichments,
log_transform = TRUE,
scale = TRUE,
num_regions = 100)
suppressPackageStartupMessages(require(ComplexHeatmap))
suppressPackageStartupMessages(require(circlize))
Heatmap(enrichments_matrix, name = "enrichments",
col = colorRamp2(c(-1, 0, 1), c("green", "white", "red")),
show_row_names = FALSE,
show_column_names = TRUE,
show_row_dend = FALSE,
column_names_gp = gpar(fontsize = 8))
- 计算样品相关性
样品间相关性如何是我们大家都所关心的。
plot_correlation()函数就是用来做这个事情的。该函数与multiBigwig_summary()函数的输出以及与其他格式相似的工具兼容。
plot_correlation(enrichments_df = enrichments,
log_transform = TRUE,
plot_type = "replicate_level",
sample_metadata = chr21_data_table)
除了在组内和组之间进行重复相关之外，还可以在对组内的所有样本取平均后进行组与组之间相关性计算。
plot_correlation()函数中的参数plot_type = "group_level"正是这样做的。
plot_correlation(enrichments_df = enrichments,
log_transform = TRUE,
plot_type = "group_level",
sample_metadata = chr21_data_table)
可以分别使用传递给
corrplot :: corrplot或GGally :: ggpairs的来进一步修改copy_level或group_level图的外观。
- Plot enrichments across groups
- 一旦使用各种差异富集包确定了特定于组的基因组区域（ 或 Peak ），例如
DESeq2，diffBind，QSEA，可能希望可视化所有组所有样本中的富集程度，以显示富集差异的大小。plot_enrichments()函数采用来自multiBigwig_summary()或类似格式的数据富集的data.frame，并以箱形图和小提琴图的组合形式绘制富集图。该功能受论文（ Allen M 2018 ）和ggplot2扩展包gghalves( Frederik 2019 )的所启示。有两种方法可以绘制富集差异，一种方法是在对每个区域的一组中的所有样本取平均后直接绘制组水平富集，另一种方法是绘制每组的 paired conditions，例如对于处理和未经处理的转录因子富集。在这两种情况下，函数都需要一个sample_metadata表以及enrichments data.frame。- 以下示例说明了在不同设置中对
plot_enrichments()函数的用法。如果plot_type = "separate"，函数将绘制组与组之间的富集水平
## plot_type == "separate"
plot_enrichments(enrichments_df = enrichments,
log_transform = TRUE,
plot_type = "separate",
sample_metadata = chr21_data_table)
- 如果
plot_type = "overlap"，函数会绘制箱形图和重叠小提琴，以显示成对条件下的分布。
这些图的sample_metadata还需要另外一列来描述样本状态。请参阅以下示例
## plot_type == "overlap"
enrichemnts_4_overlapviolins <- system.file("extdata/overlap_violins",
"enrichemnts_4_overlapviolins.txt",
package = "ALPS", mustWork = TRUE)
enrichemnts_4_overlapviolins <- read.delim(enrichemnts_4_overlapviolins, header = TRUE)
## metadata associated with above enrichments
data_table_4_overlapviolins <- system.file("extdata/overlap_violins",
"data_table_4_overlapviolins.txt",
package = "ALPS", mustWork = TRUE)
data_table_4_overlapviolins <- read.delim(data_table_4_overlapviolins, header = TRUE)
## enrichments table
enrichemnts_4_overlapviolins %>% head
chr start end s1 s2 s3 s4 s5 s6 s7 s8
1 chr21 5101659 5102227 0.25526578 0.4611994 0.21438043 0.1573575 -0.05081961 0.5747147 0.56832485 0.3590694
2 chr21 5128223 5128741 0.82057469 0.9359073 0.66987324 0.8718381 0.31443426 1.1137688 -0.20168400 0.8110701
3 chr21 5154593 5155112 0.80378039 0.8452937 0.61775645 0.5620449 0.29282731 0.8309985 -0.04550986 0.9127605
4 chr21 5220488 5221005 0.04583463 0.7161867 0.04999978 0.5506898 0.41634277 0.8684639 0.37864897 0.3802162
5 chr21 5221136 5221912 0.87553172 0.1238063 0.94939878 0.5923653 0.24742289 0.3206298 -0.13936118 0.7587155
6 chr21 5223327 5223826 0.12800909 0.5948275 0.58255203 0.7278129 -0.14432328 0.3536091 0.48731594 0.8615313
## metadata table
data_table_4_overlapviolins %>% head
bw_path sample_id sample_status group color_code
1 path.bw s1 untreated tf1 gray
2 path.bw s2 untreated tf2 gray
3 path.bw s3 untreated tf1 gray
4 path.bw s4 untreated tf2 gray
5 path.bw s5 treated tf1 red
6 path.bw s6 treated tf2 red
plot_enrichments(enrichments_df = enrichemnts_4_overlapviolins,
log_transform = FALSE,
plot_type = "overlap",
sample_metadata = data_table_4_overlapviolins,
overlap_order = c("untreated", "treated"))
有单独和重叠的其他参数可用于修改外观（请检查
?plot_enrichments），此外，该函数返回ggplot2对象，使用户能够更改图的其他组件。
- 绘制 UCSC 基因组浏览器图
- 在任何全基因组的表观基因组学分析中，通常有趣的是检查某些基因组基因座上的富集，例如基因组区域的各种组蛋白修饰，它们定义了
染色质状态或启动子转录因子的共结合或者增强子元件。实现此任务的经典方法是将所有bigwig文件加载到IGV或在UCSC中创建数据中心，然后导航到感兴趣的区域。这并不总是可行的，需要大量的人工工作，另外，还需要一台UCSC基因组浏览器服务器才能使用未发布的数据完成此任务。- 为了解决这个问题，设计了几种
R/bioconductor包（例如Gviz，karyoploteR）。即使在R环境中，也需要投入大量精力来创建UCSC基因组浏览器（如图）。ALPS包中的plot_broswer_tracks()函数需要最少的数据表和基因组区域输入，并产生像plot这样的发布质量浏览器。该函数使用Gviz包中的实用程序生成可视化效果（ Hahne, Ivanek 2016 ）。- 以下代码段说明了如何使用此功能:
## gene_range
gene_range = "chr21:45643725-45942454"
plot_browser_tracks(data_table = chr21_data_table,
gene_range = gene_range,
ref_gen = "hg38")
- 注释基因组区间
- 全基因组表观基因组分析的常见任务之一是确定
Peak/binding的基因组位置。概述了特定转录因子经常结合的位置或观察到特定类型的组蛋白修饰的位置。函数get_genomic_annotations()利用上述数据表并返回在每个基因组特征（如启动子，UTR，基因间区域等）中发现的Peak/binding位点的百分比。此函数是ChIPseeker的annotatePeak()函数（ Yu, Wang, He 2015 ）。
- 函数还提供了带有
merge_level的选项，用于合并来自不同级别的不同样本的overlaping Peak。
- 全部通过合并
data_table中存在的所有样本的overlaping Peak来创建一个共同的Peak集group_level创建组与组之间共同识别的Peak集。表示每组所有样本的overlaping Peak将被合并none不会创建任何共同识别的Peak集。将返回每个样本的基因组注释。
g_annotations <- get_genomic_annotations(data_table = chr21_data_table,
ref_gen = "hg38",
tss_region = c(-1000, 1000),
merge_level = "group_level")
Warning message:
In data.table::dcast(., Feature ~ .id, value.var = "Frequency") :
The dcast generic in data.table has been passed a data.frame and will attempt to redirect to the reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now deprecated as well. Please do this redirection yourself like reshape2::dcast(.). In the next version, this warning will become an error.
g_annotations %>% head
Feature ACCx BRCA GBMx LGGx
1 Promoter 45.762712 31.2977099 41.666667 41.666667
2 5' UTR 1.694915 0.7633588 NA NA
3 3' UTR 6.779661 3.8167939 6.250000 5.000000
4 1st Exon 1.694915 1.5267176 4.166667 5.000000
5 Other Exon 8.474576 9.1603053 6.250000 6.666667
6 1st Intron 5.084746 3.0534351 NA 3.333333
- 基因组注释相关可视化
从
get_genomic_annotations()返回的结果可以直接传递到函数plot_genomic_annotations()以可视化每个特征中的Peak百分比。该功能可以产生可视化的条形图或热图
plot_genomic_annotations(annotations_df = g_annotations, plot_type = "heatmap")
plot_genomic_annotations(annotations_df = g_annotations, plot_type = "bar")
Warning message:
Removed 4 rows containing missing values (position_stack).
## logo plot
## emm, 这里是在数据中没有的，所以出不来的。
plot_motif_logo(motif_path = myc_transfac,
database = "transfac",
plot_type = "logo")
