看到可视化的包，真的停不下来，就想先翻译一边先。

参考链接

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_pa​​th 列中的所有 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")