看到可視化的包，真的停不下來，就想先翻譯一邊先。

參考鏈接

APLS 全稱 AnaLysis routines for ePigenomicS data, 顧名思義，是一個(gè)分析表觀數(shù)據(jù)的工具。產(chǎn)生可用發(fā)表的可視化文件，主要針對全基因組表觀基因組學(xué)數(shù)據(jù)，例如 ChIP-seq，ATAC-seq等。

功能簡介：

• 計(jì)算基因組區(qū)域的富集
• 計(jì)算樣品相關(guān)性
• Plot enrichments across groups
• 繪制 UCSC 基因組瀏覽器圖
• 注釋基因組區(qū)間
• 基因組注釋相關(guān)可視化

我覺得這個(gè)包值得學(xué)習(xí)的是它底層的代碼，有機(jī)會可以看看。

能得到什么圖？

Bigwig 文件

Bigwig 文件通常用于存儲全基因組范圍的定量數(shù)據(jù)，如 ChIP-seq，ATAC-seq，WGBS，GRO-seq等。下圖說明了 Bigwig 文件的示例。

怎么產(chǎn)生 Bigwig 文件

• 使用 UCSC 中的生成 bigwig 方法 的 deeptools 軟件的 bamCoverage 功能將 BAM → bigwig。一旦生成了標(biāo)準(zhǔn)化的 bigwig 文件并從 BAM 文件中識別出 Peak，就很少在后續(xù)流程中再次使用 BAM 文件（ 啊哈哈，也就是說是為 bw 文件而生 ）。所有下游過程的要求都可以通過標(biāo)準(zhǔn)化的 bigwig 文件來滿足，例如，量化峰值或啟動子的標(biāo)準(zhǔn)化 reads 數(shù)，通過基因組瀏覽器或 IGV 中的進(jìn)行展示。

• 鑒定出 Peak 后，隨后的步驟將是對鑒定出的 Peak 處的標(biāo)準(zhǔn)化 reads 數(shù)進(jìn)行定量，以進(jìn)行探索性數(shù)據(jù)分析（ explorative data analysis: EDA ），無監(jiān)督聚類分析 PCA，以識別所考慮樣品中的模型并產(chǎn)生新穎的生物學(xué)見解。

ALPS 流程

• ALPS 包的設(shè)計(jì)理念是從最少的輸入開始，并從數(shù)據(jù)中獲取豐富的見解。大多數(shù)時(shí)候，ALPS 中的大多數(shù)功能都只需要一個(gè)數(shù)據(jù)表格，該表具有指向 bigwig 文件的路徑和相關(guān)的示例元信息。各種功能將利用此數(shù)據(jù)表格并生成下游輸出。該包產(chǎn)生可發(fā)表的可視化效果，其中大部分可以使用 ggplot2 在 R 中進(jìn)行自定義。
• 以下是 ALPS 流程和可用功能的概述：
  

實(shí)戰(zhàn)演練

• 安裝包，注意：R3.6.1

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ALPS")

• 計(jì)算基因組區(qū)域的富集

• 表觀基因組學(xué)中的大多數(shù)探索性分析都始于對一組基因組區(qū)域的富集或甲基化進(jìn)行量化，比如: 啟動子區(qū)域 或 Peak 區(qū)域。量化后的文件將用作下游分析（ 例如 PCA，聚類 ）的輸入。multiBigwig_summary() 函數(shù)使用帶有 bigwig 路徑 和 bed 文件路徑 的樣本文件表格計(jì)算基因組區(qū)域的富集。此功能是 rtracklayer 包對 bigwig 實(shí)用程序的包裝。在計(jì)算富集之前，該功能會同時(shí)根據(jù)輸入數(shù)據(jù)表中存在的所有 bed 文件生成共同的 Peak。
• 從 ALPS 包中讀取數(shù)據(jù)表格
  可以看到這一步就是得到上面所說的表格 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 （我覺得這里用英文表述會好點(diǎn)，就不翻譯了）
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

然后運(yùn)行multiBigwig_summary() 函數(shù)，通過同時(shí)從 bed_pa??th 列中的所有 bed 文件準(zhǔn)備共同的 Peak 來計(jì)算所有 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

幾乎不費(fèi)吹灰之力，multiBigwig_summary() 函數(shù)的輸出可以很容易地與其他 R/bioconductor 包裝集成在一起，以進(jìn)行探索性分析，PCA 或 聚類 。

• get_variable_regions() 獲取 multiBigwig_summary 或??類似格式的輸出，并返回 n 個(gè)可縮放的可變區(qū)域，可以直接與 ComplexHeatmap （ 顧大佬的包真的強(qiáng) ）等工具一起使用。
• 以下是有關(guān)如何通過 get_variable_regions() 函數(shù)將 multiBigwig_summary() 函數(shù)輸出結(jié)果集成到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))

• 計(jì)算樣品相關(guān)性

樣品間相關(guān)性如何是我們大家都所關(guān)心的。plot_correlation() 函數(shù)就是用來做這個(gè)事情的。該函數(shù)與multiBigwig_summary() 函數(shù)的輸出以及與其他格式相似的工具兼容。

plot_correlation(enrichments_df = enrichments, 
                 log_transform = TRUE, 
                 plot_type = "replicate_level", 
                 sample_metadata = chr21_data_table)

除了在組內(nèi)和組之間進(jìn)行重復(fù)相關(guān)之外，還可以在對組內(nèi)的所有樣本取平均后進(jìn)行組與組之間相關(guān)性計(jì)算。plot_correlation() 函數(shù)中的參數(shù) 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 的來進(jìn)一步修改 copy_level 或group_level 圖的外觀。

• Plot enrichments across groups

• 一旦使用各種差異富集包確定了特定于組的基因組區(qū)域（ 或 Peak ），例如 DESeq2，diffBind，QSEA ，可能希望可視化所有組所有樣本中的富集程度，以顯示富集差異的大小。plot_enrichments() 函數(shù)采用來自 multiBigwig_summary() 或??類似格式的數(shù)據(jù)富集的 data.frame，并以箱形圖和小提琴圖的組合形式繪制富集圖。該功能受論文（ Allen M 2018 ）和 ggplot2 擴(kuò)展包 gghalves ( Frederik 2019 )的所啟示。有兩種方法可以繪制富集差異，一種方法是在對每個(gè)區(qū)域的一組中的所有樣本取平均后直接繪制組水平富集，另一種方法是繪制每組的 paired conditions，例如對于 處理和未經(jīng)處理的轉(zhuǎn)錄因子 富集。在這兩種情況下，函數(shù)都需要一個(gè)sample_metadata 表以及 enrichments data.frame。
• 以下示例說明了在不同設(shè)置中對 plot_enrichments() 函數(shù)的用法。如果 plot_type = "separate"，函數(shù)將繪制組與組之間的富集水平

## plot_type == "separate"

plot_enrichments(enrichments_df = enrichments, 
                 log_transform = TRUE, 
                 plot_type = "separate", 
                 sample_metadata = chr21_data_table)

• 如果 plot_type = "overlap"，函數(shù)會繪制箱形圖和重疊小提琴，以顯示成對條件下的分布。
  這些圖的 sample_metadata 還需要另外一列來描述樣本狀態(tài)。請參閱以下示例

## 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"))

有單獨(dú)和重疊的其他參數(shù)可用于修改外觀（請檢查 ?plot_enrichments ），此外，該函數(shù)返回 ggplot2對象，使用戶能夠更改圖的其他組件。

• 繪制 UCSC 基因組瀏覽器圖

• 在任何全基因組的表觀基因組學(xué)分析中，通常有趣的是檢查某些基因組基因座上的富集，例如基因組區(qū)域的各種組蛋白修飾，它們定義了染色質(zhì)狀態(tài)或啟動子轉(zhuǎn)錄因子的共結(jié)合或者 增強(qiáng)子元件。實(shí)現(xiàn)此任務(wù)的經(jīng)典方法是將所有 bigwig 文件加載到 IGV 或在 UCSC 中創(chuàng)建數(shù)據(jù)中心，然后導(dǎo)航到感興趣的區(qū)域。這并不總是可行的，需要大量的人工工作，另外，還需要一臺 UCSC 基因組瀏覽器服務(wù)器才能使用未發(fā)布的數(shù)據(jù)完成此任務(wù)。
• 為了解決這個(gè)問題，設(shè)計(jì)了幾種 R/bioconductor 包（例如 Gviz，karyoploteR ）。即使在 R 環(huán)境中，也需要投入大量精力來創(chuàng)建 UCSC 基因組瀏覽器（如圖）。ALPS 包中的 plot_broswer_tracks() 函數(shù)需要最少的數(shù)據(jù)表和基因組區(qū)域輸入，并產(chǎn)生像 plot 這樣的發(fā)布質(zhì)量瀏覽器。該函數(shù)使用 Gviz 包中的實(shí)用程序生成可視化效果（ 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")

• 注釋基因組區(qū)間

• 全基因組表觀基因組分析的常見任務(wù)之一是確定 Peak/binding 的基因組位置。概述了特定轉(zhuǎn)錄因子經(jīng)常結(jié)合的位置或觀察到特定類型的組蛋白修飾的位置。函數(shù) get_genomic_annotations() 利用上述數(shù)據(jù)表并返回在每個(gè)基因組特征（如啟動子，UTR，基因間區(qū)域等）中發(fā)現(xiàn)的 Peak/binding 位點(diǎn)的百分比。此函數(shù)是 ChIPseeker 的 annotatePeak() 函數(shù)（ Yu, Wang, He 2015 ）。

• 函數(shù)還提供了帶有 merge_level 的選項(xiàng)，用于合并來自不同級別的不同樣本的 overlaping Peak。
  • 全部通過合并 data_table 中存在的所有樣本的 overlaping Peak 來創(chuàng)建一個(gè)共同的 Peak 集
  • group_level 創(chuàng)建組與組之間共同識別的 Peak 集。表示每組所有樣本的overlaping Peak 將被合并
  • none 不會創(chuàng)建任何共同識別的 Peak 集。將返回每個(gè)樣本的基因組注釋。

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

• 基因組注釋相關(guān)可視化

從 get_genomic_annotations() 返回的結(jié)果可以直接傳遞到函數(shù) plot_genomic_annotations() 以可視化每個(gè)特征中的 Peak 百分??比。該功能可以產(chǎn)生可視化的條形圖或熱圖

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, 這里是在數(shù)據(jù)中沒有的，所以出不來的。
plot_motif_logo(motif_path = myc_transfac, 
                database = "transfac", 
                plot_type = "logo")