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

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.

• Total peaks is the number of peak called by using nucleosome free reads (<100bp)

• Peaks overlaped with union DHS (ratio) is the percentage of called peak overlaped with blacklist. The larger value shows the better quality.

• Peaks overlaped with blacklist (ratio) is the percentage of called peak overlaped with blacklist. The smaller value shows the better quality.

• Fraction of reads in peaks (FRiP) is the fraction of nucleosome free reads (<100bp) in peak. The larger value shows the better quality.

Session Info

sessionInfo()

wzthu/ATACpipe documentation built on Aug. 12, 2022, 7:38 a.m.