In wzthu/wzzw: An Easy-to-use Systematic pipeline for ATACseq data analysis

knitr::opts_chunk$set(echo = TRUE)

load("Report2.Rdata")

Case Control Peak Compare

The following is the peak comparation result for case and control.

library(VennDiagram)
library(grid)
PeakNumInCC <- getReportVal(comp_result$peakCom, "venn.data")
venn.plot <- draw.pairwise.venn(
    area1 = PeakNumInCC[1] + PeakNumInCC[3],
    area2 = PeakNumInCC[2] + PeakNumInCC[3],
    cross.area = PeakNumInCC[3],
    category = c("Case", "Control"),
    fill = c("skyblue", "mediumorchid")
)

The total peaks in case are r PeakNumInCC[4], specific peaks in case are r PeakNumInCC[1].

The total peaks in control are r PeakNumInCC[5], specific peaks in control are r PeakNumInCC[2].

The overlap peaks are merged to r PeakNumInCC[3] peaks.

GO Analysis

Gene ontology analysis for all genes around specific peak regions of case and control.

Significant GO terms of genes associated with case specific peaks

go_path <- getReportVal(comp_result$goAna.case, "goOutput")
go_data <- read.table(file = go_path, header = TRUE, sep = "\t")
go_data <- subset(go_data, select = c("ID", "Description", "GeneRatio", "pvalue", "qvalue"))
go_data$pvalue <- signif(go_data$pvalue, digits = 3)
go_data$pvalue <- as.character(go_data$pvalue)
go_data$qvalue <- signif(go_data$qvalue, digits = 3)
go_data$qvalue <- as.character(go_data$qvalue)
if(nrow(go_data)==0){
    print("No GO terms found: empty table")
}else if(nrow(go_data) < 15){
    knitr::kable(go_data, align = "l")
}else{
    knitr::kable(go_data[1:15, ], align = "l")
}

Click to Visit Case Differential Go Analysis file

Significant GO terms of genes associated with control specific peaks

go_path <- getReportVal(comp_result$goAna.ctrl, "goOutput")
go_data <- read.table(file = go_path, header = TRUE, sep = "\t")
go_data <- subset(go_data, select = c("ID", "Description", "GeneRatio", "pvalue", "qvalue"))
go_data$pvalue <- signif(go_data$pvalue, digits = 3)
go_data$pvalue <- as.character(go_data$pvalue)
go_data$qvalue <- signif(go_data$qvalue, digits = 3)
go_data$qvalue <- as.character(go_data$qvalue)
if(nrow(go_data)==0){
    print("No GO terms found: empty table")
}else if(nrow(go_data) < 15){
    knitr::kable(go_data, align = "l")
}else{
    knitr::kable(go_data[1:15, ], align = "l")
}

Click to Visit Ctrl Differential Go Analysis file

Motif Enrichment Analysis

This function takes case and control specific peaks as foreground respectively, overlap peaks between case and control as background.

Significant motifs of case specific peaks

motif_enrich.case <- getReportVal(comp_result$mout, "rdsOutput.peak1")
motif_enrich.case <- motif_enrich.case[, c(1, 3, 4)]
colnames(motif_enrich.case) <- c("motif", "motif length", "p_value")
motif_enrich.case <- motif_enrich.case[order(motif_enrich.case$p_value), ]
rownames(motif_enrich.case) <- seq(nrow(motif_enrich.case))
motif_enrich.case$p_value <- signif(motif_enrich.case$p_value, digits = 3)
motif_enrich.case$p_value[motif_enrich.case$p_value < 1e-300] <- 0
motif_enrich.case$p_value <- as.character(motif_enrich.case$p_value)
if(nrow(motif_enrich.case) < 15){
    knitr::kable(motif_enrich.case, align = "l")
}else{
    knitr::kable(motif_enrich.case[1:15, ], align = "l")
}

Significant motifs of control specific peaks

motif_enrich.ctrl <- getReportVal(comp_result$mout, "rdsOutput.peak2")
motif_enrich.ctrl <- motif_enrich.ctrl[, c(1, 3, 4)]
colnames(motif_enrich.ctrl) <- c("motif", "motif length", "p_value")
motif_enrich.ctrl <- motif_enrich.ctrl[order(motif_enrich.ctrl$p_value), ]
rownames(motif_enrich.ctrl) <- seq(nrow(motif_enrich.ctrl))
motif_enrich.ctrl$p_value <- signif(motif_enrich.ctrl$p_value, digits = 3)
motif_enrich.ctrl$p_value[motif_enrich.ctrl$p_value < 1e-300] <- 0
motif_enrich.ctrl$p_value <- as.character(motif_enrich.ctrl$p_value)
if(nrow(motif_enrich.ctrl)==0){
    print("No motif found: empty table")
}else if(nrow(motif_enrich.ctrl) < 15){
    knitr::kable(motif_enrich.ctrl, align = "l")
}else{
    knitr::kable(motif_enrich.ctrl[1:15, ], align = "l")
}

Genomic Footprint

The following is the footprint for motif occurance of case&control peaks.

par(mfrow=c(1,2))

footprint_data <- getReportVal(caselist$atacProcs$footprint, "footprint.data")
if("CTCF" %in% names(footprint_data)){
    footprint_figure.name <- "CTCF"
    footprint_figure.data <- as.vector(footprint_data$CTCF)
}else{
    footprint_figure.name <- names(footprint_data[1])
    footprint_figure.data <- as.vector(footprint_data[[1]])
}
footprint_figure.length <- length(footprint_figure.data) - 200
footprint_text <- paste(footprint_figure.name, "(Case)", sep = "")
plot(footprint_figure.data, type = "l", col = "blue", lwd = 2,
     main = footprint_text,
    xlab = "Relative Distance From Motif (bp)",
    ylab = "Cut Site Count", xaxt = "n", yaxt = "n")
axis(1, at = seq(1, 100, len = 3),
    labels = -(100 + 1 - seq(1, 100 + 1, len = 3)),
    padj = -1.0, tck = -0.01)
axis(1, at = 100 + footprint_figure.length + seq(1, 100, len = 3),
    labels = seq(0, 100, len = 3),
    padj = -1.0, tck = -0.01)
axis(2, padj = 1.0,tck = -0.02)
abline(v = c(100, 100 + footprint_figure.length + 1), lty = 2)

footprint_data <- getReportVal(ctrllist$atacProcs$footprint, "footprint.data")
if("CTCF" %in% names(footprint_data)){
    footprint_figure.name <- "CTCF"
    footprint_figure.data <- as.vector(footprint_data$CTCF)
}else{
    footprint_figure.name <- names(footprint_data[1])
    footprint_figure.data <- as.vector(footprint_data[[1]])
}
footprint_figure.length <- length(footprint_figure.data) - 200
footprint_text <- paste(footprint_figure.name, "(Control)", sep = "")
plot(footprint_figure.data, type = "l", col = "blue", lwd = 2,
     main = footprint_text,
    xlab = "Relative Distance From Motif (bp)",
    ylab = "Cut Site Count", xaxt = "n", yaxt = "n")
axis(1, at = seq(1, 100, len = 3),
    labels = -(100 + 1 - seq(1, 100 + 1, len = 3)),
    padj = -1.0, tck = -0.01)
axis(1, at = 100 + footprint_figure.length + seq(1, 100, len = 3),
    labels = seq(0, 100, len = 3),
    padj = -1.0, tck = -0.01)
axis(2, padj = 1.0,tck = -0.02)
abline(v = c(100, 100 + footprint_figure.length + 1), lty = 2)

case.dir <- getReportVal(caselist$atacProcs$footprint, "pdf.dir")
ctrl.dir <- getReportVal(ctrllist$atacProcs$footprint, "pdf.dir")

For all footprint of case sample, The absolute path is r R.utils::getAbsolutePath(case.dir).

For all footprint of control sample, The absolute path is r R.utils::getAbsolutePath(ctrl.dir).

The following is the footprint for motif occurance of case&control specific peaks.

par(mfrow=c(1,2))

footprint_data <- getReportVal(comp_result$footprint.case, "footprint.data")
if("CTCF" %in% names(footprint_data)){
    footprint_figure.name <- "CTCF"
    footprint_figure.data <- as.vector(footprint_data$CTCF)
}else{
    footprint_figure.name <- names(footprint_data[1])
    footprint_figure.data <- as.vector(footprint_data[[1]])
}
footprint_figure.length <- length(footprint_figure.data) - 200
footprint_text <- paste(footprint_figure.name, "(Case specific)", sep = "")
plot(footprint_figure.data, type = "l", col = "blue", lwd = 2,
     main = footprint_text,
    xlab = "Relative Distance From Motif (bp)",
    ylab = "Cut Site Count", xaxt = "n", yaxt = "n")
axis(1, at = seq(1, 100, len = 3),
    labels = -(100 + 1 - seq(1, 100 + 1, len = 3)),
    padj = -1.0, tck = -0.01)
axis(1, at = 100 + footprint_figure.length + seq(1, 100, len = 3),
    labels = seq(0, 100, len = 3),
    padj = -1.0, tck = -0.01)
axis(2, padj = 1.0,tck = -0.02)
abline(v = c(100, 100 + footprint_figure.length + 1), lty = 2)

footprint_data <- getReportVal(comp_result$footprint.ctrl, "footprint.data")
if("CTCF" %in% names(footprint_data)){
    footprint_figure.name <- "CTCF"
    footprint_figure.data <- as.vector(footprint_data$CTCF)
}else{
    footprint_figure.name <- names(footprint_data[1])
    footprint_figure.data <- as.vector(footprint_data[[1]])
}
footprint_figure.length <- length(footprint_figure.data) - 200
footprint_text <- paste(footprint_figure.name, "(Control specific)", sep = "")
plot(footprint_figure.data, type = "l", col = "blue", lwd = 2,
     main = footprint_text,
    xlab = "Relative Distance From Motif (bp)",
    ylab = "Cut Site Count", xaxt = "n", yaxt = "n")
axis(1, at = seq(1, 100, len = 3),
    labels = -(100 + 1 - seq(1, 100 + 1, len = 3)),
    padj = -1.0, tck = -0.01)
axis(1, at = 100 + footprint_figure.length + seq(1, 100, len = 3),
    labels = seq(0, 100, len = 3),
    padj = -1.0, tck = -0.01)
axis(2, padj = 1.0,tck = -0.02)
abline(v = c(100, 100 + footprint_figure.length + 1), lty = 2)

casesp.dir <- getReportVal(comp_result$footprint.case, "pdf.dir")
ctrlsp.dir <- getReportVal(comp_result$footprint.ctrl, "pdf.dir")