knitr::opts_chunk$set(echo = TRUE)
load("Report.Rdata")

Concordance

Shared peaks

library(VennDiagram)
library(grid)
if(length(fastqInput1)>5){
    warning("'VennDiagram' will not be generated for more than 5 replicates")
}else{
    gridobj<-venn.diagram(x = peaknumset,filename = NULL)
    grid.draw(gridobj)
}

Correlation

library(corrplot)
coltable <- colorRampPalette(c("red", "white", "blue")) 
corrplot(correlation, method = "color", type = "upper",addCoef.col = "grey", cl.lim = c(0, 1),col = coltable(100))

Reads Alignment Statistics

Fragment size distribution

library(ggplot2)
load("Report.Rdata")
readsCounts<-getReportVal(atacProcs$fragLenDistr,"readsCounts")
ggplot(readsCounts[1:1000,], aes(length,counts))+geom_path(color="Red")+xlab("Fragment length (bp)")+ylab("Read counts") + theme_bw() + theme(panel.grid =element_blank()) + labs(title = "Fragment Length Distribution") + theme(plot.title = element_text(hjust = 0.5))

library(stats)
strength<-Mod(fft(readsCounts$counts))/length(readsCounts$counts)
periodx<-length(readsCounts$counts)/(1:(length(readsCounts$counts)-1))
strength<-strength[2:length(strength)]

rs1<-as.data.frame(cbind(periodx[periodx<20&periodx>2],strength[periodx<20&periodx>2],0))
rs2<-as.data.frame(cbind(periodx[periodx<400&periodx>2],strength[periodx<400&periodx>2],1))
rs<-rbind(rs1,rs2)
colnames(rs)<-c("period","strength","check")

g1<-ggplot(rs[rs["check"]==0,]) + geom_vline(xintercept = 10.4, linetype=2)+ geom_line(aes(x=period,y=strength),color="Red")+ theme_bw() + theme(panel.grid =element_blank()) + annotate("text", x = 10.4, y = max(rs[rs["check"]==0,2]), label = "10.4bp") +xlab("period") + ylab("strength")+ labs(title = "the Pitch of the DNA Helix") + theme(plot.title = element_text(hjust = 0.5))

g2<-ggplot(rs[rs["check"]==1,]) + geom_vline(xintercept = 186, linetype=2)+ geom_line(aes(x=period,y=strength),color="Red")+ theme_bw() + theme(panel.grid =element_blank()) + annotate("text", x = 186, y = max(rs[rs["check"]==1,2]), label = "186bp") +xlab("period") + ylab("strength")+ labs(title = "Nucleosome") + theme(plot.title = element_text(hjust = 0.5))
library(gridExtra)
grid.arrange(g1, g2, ncol=2)

TSS enrichment

The nucleosome free reads (<100bp) and monnucleosome span reads (180~247bp) enrichment around transcription starting site (TSS) are shown below.

library(ggplot2)
load("Report.Rdata")
df<-getReportVal(atacProcs$tssqc100,"tss")
g1<-ggplot(df,aes(pos,counts))+geom_line()+ geom_vline(xintercept = 0, linetype=2)+xlab("upstream<-TSS->downstream")+ylab("reads count")+theme_bw() + theme(panel.grid =element_blank()) + labs(title = "Nucleosome Free Reads") + theme(plot.title = element_text(hjust = 0.5))
df<-getReportVal(atacProcs$tssqc180_247,"tss")
g2<-ggplot(df,aes(pos,counts))+geom_line()+ geom_vline(xintercept = 0, linetype=2)+xlab("upstream<-TSS->downstream")+ylab("reads count")+theme_bw() + theme(panel.grid =element_blank()) + labs(title = "Monnucleosome Span Reads") + theme(plot.title = element_text(hjust = 0.5))
grid.arrange(g1, g2, ncol=2)

Peak Statistics

Blacklist ratio

ifelse(is.null(atacProcs$blacklistQC),"NA",
knitr::kable(x = getReportVal(atacProcs$blacklistQC,"report")))

DHS ratio

ifelse(is.null(atacProcs$DHSQC),"NA",
knitr::kable(x = getReportVal(atacProcs$DHSQC,"report")))

Fraction of reads in peaks (FRiP)

knitr::kable(x = getReportVal(atacProcs$fripQC,"report"))

Peak Annotation

library(ChIPseeker)
peakanno <- getReportVal(atacProcs$Peakanno,"annoOutput.rds")
plotAnnoPie(x = peakanno)

Gene Ontology Analysis

Gene ontology analysis for all genes around peak regions.

go_path <- getReportVal(atacProcs$goAna, "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){
    message("No GO terms found: empty table")
}else if(nrow(go_data) < 15){
    knitr::kable(go_data)
}else{
    knitr::kable(go_data[1:15, ])
}

Click to Visit Go Analysis file

Footprint Analysis

footprint_data <- getReportVal(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, "(Footprinting)", sep = "")

The following figure is r footprint_figure.name footprint.

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)

pdf.dir <- getReportVal(atacProcs$footprint,"pdf.dir")

All motif footprint figures are saved as pdf files. The absolute path is r R.utils::getAbsolutePath(pdf.dir).

Annotation of items in table

For single end sequencing data, esATAC will counts reads number.

For paired end sequencing data, esATAC will counts read pairs or fragment number.

Session Info

sessionInfo()


wzthu/ATACFlow documentation built on Aug. 9, 2022, 2:24 a.m.