Nothing
R/Stage_3_PeakFinderFunctions.R
In MACPET: Model based analysis for paired-end data
Defines functions
FUN
FindRegions_fun
TakeFitAnalysisData_fun
Stage_3_Main_fun
# Helping functions for PeakFinder.R
#' @importFrom plyr dlply . ldply llply ddply
#' @importFrom intervals Intervals clusters interval_included size
#' @importFrom stats ppois p.adjust
#' @importFrom S4Vectors metadata
#' @importFrom GenomeInfoDb seqlengths seqinfo
#' @importFrom BiocParallel bplapply
#' @importClassesFrom S4Vectors Annotated
#' @importClassesFrom intervals Intervals
############################################## Main function for stage 3:
#-------------
#-------------
Stage_3_Main_fun = function(SA_prefix, S3_method, S3_AnalysisDir, S3_Selfobject,
S3_image) {
# global variables for Rcheck:
#--------------------------------------------
#---------------Take and reorder Input:
#--------------------------------------------
# Take time:
Analysis.time.start = Sys.time()
# create folder to save data:
if (!dir.exists(S3_AnalysisDir))
dir.create(S3_AnalysisDir)
#--------------------------------------------
#---------------Keep data you need only:
#--------------------------------------------
FitDataInfo = TakeFitAnalysisData_fun(S3_Selfobject = S3_Selfobject)
PETsData = FitDataInfo$PETsData #all the PETs data, used in inference too.
ChromInf = FitDataInfo$ChromInf #used in inference and fitting.
#--------------------------------------------
#---------------Segment to regions:
#--------------------------------------------
SegmRes = FindRegions_fun(PETsData = PETsData, ChromInf = ChromInf)
SegPETS = SegmRes$SegPETS #segmented PETs data, used in fitting.
RegionCounts = SegmRes$RegionCounts #take the counts
#------------
#-----save the region counts on main data:
#------------
Ordermatch = match(RegionCounts$Chrom, S4Vectors::metadata(S3_Selfobject)$Self_info$Chrom)
S4Vectors::metadata(S3_Selfobject)$Self_info$Region.counts[Ordermatch] = RegionCounts$Counts
#--------------------------------------------
#---------------Fit:
#--------------------------------------------
GlobalPeakCallRes = FitCallGlobal_fun(SegPETS = SegPETS, ChromInf = ChromInf)
#--------------------------------------------
#---------------update the results
#--------------------------------------------
#------classification information
S4Vectors::metadata(S3_Selfobject)$Classification.Info = GlobalPeakCallRes$Classification.Info
#------Peak counts
S4Vectors::metadata(S3_Selfobject)$Self_info$Peak.counts = 0
Match = match(GlobalPeakCallRes$Peak.counts$Chrom, S4Vectors::metadata(S3_Selfobject)$Self_info$Chrom)
S4Vectors::metadata(S3_Selfobject)$Self_info$Peak.counts[Match] = GlobalPeakCallRes$Peak.counts$V1
#--------------------------------------------
#---------------Run inference:
#--------------------------------------------
Peaks.Info = GlobalPeakCallRes$Peaks.Info #peaks found by fit
# run inference
InfRes = SignificanceCall_fun(Peaks.Info = Peaks.Info, PETsData = PETsData, ChromInf = ChromInf,
S3_method = S3_method)
#---add Peak information for peaks found in data with their FDR etc
S4Vectors::metadata(S3_Selfobject)$Peaks.Info = InfRes
#------------------------
#-----update class and save:
#------------------------
class(S3_Selfobject) = "PSFit"
NamepsfitData = paste(SA_prefix, "_psfitData", sep = "")
assign(NamepsfitData, S3_Selfobject) #assign value.
save(list = NamepsfitData, file = file.path(S3_AnalysisDir, NamepsfitData))
#------------------------
#-----plot images
#------------------------
if (S3_image) {
#-------------
# create data:
#-------------
PeakInfo = S4Vectors::metadata(S3_Selfobject)
PeakInfo = PeakInfo$Peaks.Info
PeakInfo$FDR0.05 = "FDR >= 0.05"
PeakInfo$FDR0.05[which(PeakInfo$FDR < 0.05)] = "FDR < 0.05"
PeakInfo$ID = seq_len(nrow(PeakInfo))
# image 1:
Get_image_S3_P1_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
# image 2:
Get_image_S3_P2_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
# image 3:
Get_image_S3_P3_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
# image 4:
Get_image_S3_P4_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
}
#------------------------
# save log:
#------------------------
futile.logger::flog.info("=====================================", name = "SA_LogFile",
capture = FALSE)
futile.logger::flog.info("Stage 3 is done!", name = "SA_LogFile", capture = FALSE)
futile.logger::flog.info(paste("Analysis results for stage 3 are in:\n", S3_AnalysisDir),
name = "SA_LogFile", capture = FALSE)
# save time:
Analysis.time.end = Sys.time()
Total.Time = Analysis.time.end - Analysis.time.start
LogFile = paste("Total stage 3 time:", Total.Time, " ", units(Total.Time))
futile.logger::flog.info(LogFile, name = "SA_LogFile", capture = FALSE)
}
# done
#-------------
#-------------
############## Functions for inputs and data
#-------------
#-------------
# function for converting the data in the analysis format:
TakeFitAnalysisData_fun = function(S3_Selfobject) {
cat("Converting data for analysis...")
#------------
# take Chromosome information(some used in inference):
#------------
ChromInf = GenomeInfoDb::seqlengths(GenomeInfoDb::seqinfo(S3_Selfobject))
ChromInf = data.frame(Chrom = names(ChromInf), size = as.numeric(ChromInf), stringsAsFactors = FALSE)
Self_info = S4Vectors::metadata(S3_Selfobject)$Self_info
MatchChrom = match(Self_info$Chrom, ChromInf$Chrom)
ChromInf$PET.counts[MatchChrom] = Self_info$PET.counts
#------------
# convert to data frame and Swap Anchors:
#------------
PETsData = as.data.frame(S3_Selfobject)
# use Tag mids:
PETsData$UTag = (PETsData$start1 + PETsData$end1)/2
PETsData$DTag = (PETsData$start2 + PETsData$end2)/2
Swap = which(PETsData$UTag > PETsData$DTag)
if (length(Swap) != 0) {
UTagswap = PETsData$UTag[Swap]
DTagswap = PETsData$DTag[Swap]
PETsData$UTag[Swap] = DTagswap
PETsData$DTag[Swap] = UTagswap
}
#------------
# add mainindex
#------------
PETsData$MainIndex = seq_len(nrow(PETsData))
#------------
# keep left part of downstream and right of upstream
#------------
PETsData$Chrom = as.character(PETsData$seqnames1)
PETsData = PETsData[, c("Chrom", "UTag", "DTag", "MainIndex")]
#------------
# return:
#------------
cat("Done\n")
return(list(PETsData = PETsData, ChromInf = ChromInf))
}
# done
#-------------
#-------------
############## Functions for region segmentation:
#-------------
#-------------
# main function for region segmentation:
FindRegions_fun = function(PETsData, ChromInf) {
# rcheck:
Chrom = Region = NULL
#------------
#------------
# split by chromosome to list to pass to bplapply
bppass = plyr::dlply(PETsData, plyr::.(Chrom), function(x) x)
#------------
#------------
#--Segment to regions:
#------------
cat("Segmenting into regions...")
SegmRes = BiocParallel::bplapply(X = bppass, FUN = function(y, ChromInf) {
ChromCur = as.character(unique(y$Chrom))
ChromInfCur = subset(ChromInf, Chrom == ChromCur)
#------------
# Give the region ids:
#------------
y = RegionIds_fun(PETinf = y, ChromCur = ChromCur, ChromInfCur = ChromInfCur)
return(y)
}, ChromInf = ChromInf)
cat("Done\n")
#------------
#---break and rbind:
#------------
# take the region info:
PETinf = lapply(SegmRes, "[[", 1)
PETinf = plyr::ldply(PETinf, function(y) y)
PETinf = subset(PETinf, !is.na(Region)) #keep not NA
if (nrow(PETinf) != 0) {
# update data with the regions:
PETsData$Region = NA
PETsData$Region[PETinf$MainIndex] = PETinf$Region
PETsData = subset(PETsData, !is.na(Region))
} else {
stop("No regions found in data. The data is probably too sparse.", call. = FALSE)
}
#------------
#--take total regions info:
#------------
RegionCounts = lapply(SegmRes, "[[", 2)
RegionCounts = plyr::ldply(RegionCounts, function(y) y)
RegionCounts = RegionCounts[, c("Chrom", "Counts")]
futile.logger::flog.info(paste("Total Regions found:", sum(RegionCounts$Counts)),
name = "SA_LogFile", capture = FALSE)
return(list(SegPETS = PETsData, RegionCounts = RegionCounts))
}
# done
#-------------
#-------------
# Fuction for creating region ids by segmenting into regions:
RegionIds_fun = function(PETinf, ChromCur, ChromInfCur) {
# global variables for Rcheck:
#--------------------------
#------------
# initiate
#------------
PETinf$Region = NA
TotRegions = data.frame(Chrom = ChromCur, Counts = 0)
#------------
# take Pet Pmid to create global intervals:
#------------
PETInt = intervals::Intervals(PETinf[, c("UTag", "DTag")], closed = c(TRUE, TRUE))
#------------
# Take clusters:
#------------
PETclst = intervals::clusters(PETInt, which = TRUE, w = 0)
TotRegions$Counts = length(PETclst) #total regions
if (TotRegions$Counts != 0)
{
#------------
# region ids:
#------------
Region = rep(seq_len(length(PETclst)), lengths(PETclst))
#------------
# pets in the regions:
#------------
TagsIncluded = unlist(PETclst)
#------------
# add them according to the rownames of the PETinf
#------------
PETinf$Region[TagsIncluded] = Region
} #else no region found
PETinf = PETinf[, c("Chrom", "MainIndex", "Region")]
return(list(PETinf = PETinf, TotRegions = TotRegions))
}
# done
#-------------
#-------------
############## Functions for fitting the main model:
#-------------
#-------------
# main function for fitting: breaking each region
FitCallGlobal_fun = function(SegPETS, ChromInf) {
#---R-check
Region = Chrom = NULL
#----
bppass = plyr::dlply(SegPETS, plyr::.(Chrom, Region), function(y) y)
#---------------------------------------
#--------call the function in parallel for fitting regions:
#---------------------------------------
cat("Running peak calling process...")
# the following runs in parallel:
GlobalFitRes = BiocParallel::bplapply(X = bppass, FUN = FitCallLocal_fun_Rcpp,
ChromInf = ChromInf)
cat("Done\n")
#---------------------------
#--------save Classification:
#---------------------------
Peak.Id.Inf = plyr::llply(GlobalFitRes, function(x) x[[1]])
Peak.Id.Inf = do.call(rbind, Peak.Id.Inf)
if (is.null(Peak.Id.Inf)) {
stop("No Peaks found by the algorithm!", call. = FALSE)
}
#---------------------------
#--match classes with SegPETS
#---------------------------
SegPETS$Peak.ID = NA
Match = match(Peak.Id.Inf[, 2], SegPETS$MainIndex)
SegPETS$Peak.ID[Match] = Peak.Id.Inf[, 1]
Classification.Info = SegPETS[, c("MainIndex", "Region", "Peak.ID")]
#---------------------------
#--------save Peaks:
#---------------------------
Peaks.Info = plyr::ldply(GlobalFitRes, function(x) x[[2]])
if (is.null(Peaks.Info)) {
stop("No Peaks found by the algorithm!", call. = FALSE)
}
Peaks.Info = Peaks.Info[, which(colnames(Peaks.Info) != ".id")]
futile.logger::flog.info(paste("Total", nrow(Peaks.Info), "candidate peaks found in data"),
name = "SA_LogFile", capture = FALSE)
#---------------------------
#--------Find counts:
#---------------------------
Peak.counts = plyr::ddply(Peaks.Info, plyr::.(Chrom), nrow)
return(list(Peaks.Info = Peaks.Info, Classification.Info = Classification.Info,
Peak.counts = Peak.counts))
}
# done
#-------------
#-------------
############## Functions for the Inference:
#-------------
#-------------
# Main significance call function:
SignificanceCall_fun = function(Peaks.Info, PETsData, ChromInf, S3_method) {
# global variables for Rcheck:
Chrom = NULL
#--------------------------
#--------------------------------------------
#---------------Take information for the analysis:
#--------------------------------------------
Peaks.Info$Chrom = as.character(Peaks.Info$Chrom)
PETsData$MainIndex = NULL #remove index not needed
#--------------------------------------------
#---------------run inference:
#--------------------------------------------
cat("Splitting data by chromosome for inference...\n")
# split by chromosome all the three inputs: then merge them by chromosome
bppass = Get_SignificanceCall_Data_fun(Peaks.Info = Peaks.Info, PETsData = PETsData,
ChromInf = ChromInf)
#------------
# Run inference(call c++):
#------------
futile.logger::flog.info("Computing p-values...", name = "SA_LogFile", capture = FALSE)
InfRes = BiocParallel::bplapply(X = bppass, FUN = PoissonLocalInference_fun,
windows = c(10, 15))
#------------
# merge and return:
#------------
InfRes = do.call(rbind.data.frame, InfRes)
rownames(InfRes) = NULL
#------
# FDR at 0.05:
#------
futile.logger::flog.info("FDR adjusting p-values...", name = "SA_LogFile", capture = FALSE)
InfRes$FDRUp = stats::p.adjust(p = InfRes$p.valueUp, method = S3_method)
InfRes$FDRDown = stats::p.adjust(p = InfRes$p.valueDown, method = S3_method)
InfRes$FDR = stats::p.adjust(p = InfRes$p.value, method = S3_method)
#------
#---return:
#------
return(InfRes)
}
# done
#-------------
#-------------
# Function for breaking the inputs for significance calling:
Get_SignificanceCall_Data_fun = function(Peaks.Info, PETsData, ChromInf) {
# Rcheck:
Chrom = NULL
#-------------
# break Peaks.Info by chrom:
#-------------
Peaks.Info_list = plyr::dlply(Peaks.Info, plyr::.(Chrom), function(x) x)
#-------------
# break PETsData by chrom:
#-------------
PETsData$MainIndex = NULL #dont need index
PETsData_list = plyr::dlply(PETsData, plyr::.(Chrom), function(x) x)
#-------------
# break ChromInf by chrom:
#-------------
ChromInf_list = plyr::dlply(ChromInf, plyr::.(Chrom), function(x) x)
#-------------
# create the bppass input:
#-------------
bppass = list()
for (chri in seq_along(Peaks.Info_list)) {
# take Peaks_Info_x:
Peaks_Info_x = Peaks.Info_list[[chri]]
Chr_x = Peaks_Info_x$Chrom[1]
# take PETsData_x:
PETsData_x = PETsData_list[[which(names(PETsData_list) == Chr_x)]]
PETsData_x$Chrom = NULL #dont need that anymore
PETsData_x = as.matrix(PETsData_x)
# take ChromInf_x:
ChromInf_x = ChromInf_list[[which(names(ChromInf_list) == Chr_x)]]
ChromInf_x$Chrom = NULL #dont need that
bppass[[chri]] = list(Peaks_Info_x = Peaks_Info_x, PETsData_x = PETsData_x,
ChromInf_x = ChromInf_x)
}
return(bppass)
}
# done
#-------------
#-------------
# Function for Inference:
PoissonLocalInference_fun = function(bppass_x, windows) {
#-------------
# break data:
#-------------
Peaks_Info_x = bppass_x$Peaks_Info_x
PETsData_x = bppass_x$PETsData_x
# take UTag
UTag_x = intervals::Intervals(PETsData_x[, c(1, 1)], closed = c(TRUE, TRUE))
# take DTag
DTag_x = intervals::Intervals(PETsData_x[, c(2, 2)], closed = c(TRUE, TRUE))
# take chromosome info:
ChromSize = bppass_x$ChromInf_x$size
ChromPETs = bppass_x$ChromInf_x$PET.counts
#-------------
# ------ Up window 1:
#-------------
UpW1 = cbind(Peaks_Info_x$CIQ.Up.start - Peaks_Info_x$CIQ.Up.size * windows[1]/2,
Peaks_Info_x$CIQ.Up.start + Peaks_Info_x$CIQ.Up.size * windows[1]/2)
UpW1 = intervals::Intervals(UpW1, closed = c(TRUE, TRUE))
# observed pets in intervals
UpW1_obs = lengths(intervals::interval_included(UpW1, UTag_x))
lambdaUpW1 = UpW1_obs * Peaks_Info_x$CIQ.Up.size/(intervals::size(UpW1) + 1)
#-------------
# ------ Up window 2:
#-------------
UpW2 = cbind(Peaks_Info_x$CIQ.Up.start - Peaks_Info_x$CIQ.Up.size * windows[2]/2,
Peaks_Info_x$CIQ.Up.start + Peaks_Info_x$CIQ.Up.size * windows[2]/2)
UpW2 = intervals::Intervals(UpW2, closed = c(TRUE, TRUE))
# observed pets in intervals
UpW2_obs = lengths(intervals::interval_included(UpW2, UTag_x))
lambdaUpW2 = UpW2_obs * Peaks_Info_x$CIQ.Up.size/(intervals::size(UpW2) + 1)
#-------------
# ------ Down window 1:
#-------------
DownW1 = cbind(Peaks_Info_x$CIQ.Down.start - Peaks_Info_x$CIQ.Down.size * windows[1]/2,
Peaks_Info_x$CIQ.Down.start + Peaks_Info_x$CIQ.Down.size * windows[1]/2)
DownW1 = intervals::Intervals(DownW1, closed = c(TRUE, TRUE))
# observed pets in intervals
DownW1_obs = lengths(intervals::interval_included(DownW1, DTag_x))
lambdaDownW1 = DownW1_obs * Peaks_Info_x$CIQ.Down.size/(intervals::size(DownW1) +
1)
#-------------
# ------ Down window 2:
#-------------
DownW2 = cbind(Peaks_Info_x$CIQ.Down.start - Peaks_Info_x$CIQ.Down.size * windows[2]/2,
Peaks_Info_x$CIQ.Down.start + Peaks_Info_x$CIQ.Down.size * windows[2]/2)
DownW2 = intervals::Intervals(DownW2, closed = c(TRUE, TRUE))
# observed pets in intervals
DownW2_obs = lengths(intervals::interval_included(DownW2, DTag_x))
lambdaDownW2 = DownW2_obs * Peaks_Info_x$CIQ.Down.size/(intervals::size(DownW2) +
1)
#-------------
# ------ genome backgrounds:
#-------------
lambdaBGC_Up = ChromPETs * Peaks_Info_x$CIQ.Up.size/ChromSize
lambdaBGC_Down = ChromPETs * Peaks_Info_x$CIQ.Down.size/ChromSize
#-------------
# ------ find maximum:
#-------------
OptimalLambda = plyr::llply(seq_len(nrow(Peaks_Info_x)), function(i, lambdaUpW1,
lambdaUpW2, lambdaBGC_Up, lambdaDownW1, lambdaDownW2, lambdaBGC_Down) {
UpMax = max(c(2, lambdaUpW1[i], lambdaUpW2[i], lambdaBGC_Up[i]))
DownMax = max(c(2, lambdaDownW1[i], lambdaDownW2[i], lambdaBGC_Down[i]))
return(c(UpMax, DownMax))
}, lambdaUpW1 = lambdaUpW1, lambdaUpW2 = lambdaUpW2, lambdaBGC_Up = lambdaBGC_Up,
lambdaDownW1 = lambdaDownW1, lambdaDownW2 = lambdaDownW2, lambdaBGC_Down = lambdaBGC_Down)
OptimalLambda = do.call(rbind, OptimalLambda)
#-------------
# ------ Find folds and p.values:
#-------------
# Upstream
FoldEnrichUp = Peaks_Info_x$Pets/OptimalLambda[, 1]
pvalueUp = stats::ppois(q = Peaks_Info_x$Pets, lambda = OptimalLambda[, 1], lower.tail = FALSE)
# Downstream
FoldEnrichDown = Peaks_Info_x$Pets/OptimalLambda[, 2]
pvalueDown = stats::ppois(q = Peaks_Info_x$Pets, lambda = OptimalLambda[, 2],
lower.tail = FALSE)
# merged p-value:
pvalue = pvalueUp * pvalueDown
#-------------
# ------Return
#-------------
Peaks_Info_x$lambdaUp = OptimalLambda[, 1]
Peaks_Info_x$FoldEnrichUp = FoldEnrichUp
Peaks_Info_x$p.valueUp = pvalueUp
Peaks_Info_x$lambdaDown = OptimalLambda[, 2]
Peaks_Info_x$FoldEnrichDown = FoldEnrichDown
Peaks_Info_x$p.valueDown = pvalueDown
Peaks_Info_x$p.value = pvalue
return(Peaks_Info_x)
}
#-------------
#-------------
# done Functions for plotting:
#-------------
#-------------
Get_image_S3_P1_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
CIQ.Up.size = NULL
CIQ.Down.size = NULL
FDR0.05 = NULL
# image dir:
S3_P1_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p1_image.jpg",
sep = ""))
# plot
S3_image_p1 = ggplot2::ggplot(PeakInfo, ggplot2::aes(CIQ.Up.size, CIQ.Down.size,
color = factor(FDR0.05))) + ggplot2::geom_point(size = 1) + ggplot2::theme_bw() +
ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"), panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(), panel.background = ggplot2::element_blank()) +
ggplot2::theme(axis.text = ggplot2::element_text(size = 15, color = "black"),
axis.title = ggplot2::element_text(size = 18, color = "black"), plot.title = ggplot2::element_text(size = 18,
color = "black"), legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_text(size = 18)) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_color_discrete(guide = ggplot2::guide_legend(title = "Binding Site FDR")) +
ggplot2::scale_y_continuous(expand = c(0, 0)) + ggplot2::scale_x_continuous(expand = c(0,
0)) + ggplot2::xlab("Upstream Peak Size") + ggplot2::ylab("Downstream Peak Size") +
ggplot2::geom_abline(intercept = 0, slope = 1) + ggplot2::ggtitle("Upstream/Downstream peak sizes")
# save:
ggplot2::ggsave(plot = S3_image_p1, file = S3_P1_image_dir, scale = 2)
}
#-------------
#-------------
Get_image_S3_P2_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
ID = NULL
FDR = NULL
# image dir:
S3_P2_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p2_image.jpg",
sep = ""))
S3_image_p2 = ggplot2::ggplot(PeakInfo, ggplot2::aes(ID, sort(FDR, decreasing = FALSE))) +
ggplot2::geom_line(size = 1) + ggplot2::geom_hline(yintercept = 0.05, color = "red",
size = 1) + ggplot2::theme_bw() + ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank()) + ggplot2::theme(axis.text = ggplot2::element_text(size = 15,
color = "black"), axis.title = ggplot2::element_text(size = 18, color = "black"),
plot.title = ggplot2::element_text(size = 18, color = "black")) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_y_continuous(expand = c(0, 0)) + ggplot2::scale_x_continuous(expand = c(0,
0)) + ggplot2::xlab("Sorted Binding Site IDs") + ggplot2::ylab("FDR") + ggplot2::ggtitle("FDR values for the binding sites")
# save:
ggplot2::ggsave(plot = S3_image_p2, file = S3_P2_image_dir, scale = 2)
}
#-------------
#-------------
Get_image_S3_P3_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
CIQ.Peak.size = NULL
FDR0.05 = NULL
# image dir:
S3_P3_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p3_image.jpg",
sep = ""))
# plot
S3_image_p3 = ggplot2::ggplot(PeakInfo, ggplot2::aes(CIQ.Peak.size, fill = factor(FDR0.05))) +
ggplot2::geom_density(size = 1, alpha = 0.5) + ggplot2::theme_bw() + ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank()) + ggplot2::theme(axis.text = ggplot2::element_text(size = 15,
color = "black"), axis.title = ggplot2::element_text(size = 18, color = "black"),
plot.title = ggplot2::element_text(size = 18, color = "black"), legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_blank()) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_y_continuous(expand = c(0, 0), labels = function(y) format(y,
scientific = FALSE)) + ggplot2::scale_x_continuous(expand = c(0, 0)) +
ggplot2::xlab("Binding Site Sizes") + ggplot2::ylab("Density") + ggplot2::geom_abline(intercept = 0,
slope = 1) + ggplot2::ggtitle("Densities for the binding site sizes")
# save:
ggplot2::ggsave(plot = S3_image_p3, file = S3_P3_image_dir, scale = 2)
}
#-------------
#-------------
Get_image_S3_P4_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
FDRUp = NULL
FDRDown = NULL
FDR0.05 = NULL
# image dir:
S3_P4_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p4_image.jpg",
sep = ""))
# plot
S3_image_p4 = ggplot2::ggplot(PeakInfo, ggplot2::aes(FDRUp, FDRDown, color = factor(FDR0.05))) +
ggplot2::geom_point(size = 1) + ggplot2::theme_bw() + ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank()) + ggplot2::theme(axis.text = ggplot2::element_text(size = 15,
color = "black"), axis.title = ggplot2::element_text(size = 18, color = "black"),
plot.title = ggplot2::element_text(size = 18, color = "black"), legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_text(size = 18)) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_color_discrete(guide = ggplot2::guide_legend(title = "Binding Site FDR")) +
ggplot2::scale_y_continuous(expand = c(0, 0)) + ggplot2::scale_x_continuous(expand = c(0,
0)) + ggplot2::xlab("FDR for Upstream Binding Site Peak") + ggplot2::ylab("FDR for Downstream Binding Site Peak") +
ggplot2::geom_abline(intercept = 0, slope = 1) + ggplot2::ggtitle("Upstream/Downstream peak FDR")
# save:
ggplot2::ggsave(plot = S3_image_p4, file = S3_P4_image_dir, scale = 2)
}
Try the MACPET package in your browser
Any scripts or data that you put into this service are public.
MACPET documentation built on Nov. 8, 2020, 5:47 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions FUN FindRegions_fun TakeFitAnalysisData_fun Stage_3_Main_fun
# Helping functions for PeakFinder.R
#' @importFrom plyr dlply . ldply llply ddply
#' @importFrom intervals Intervals clusters interval_included size
#' @importFrom stats ppois p.adjust
#' @importFrom S4Vectors metadata
#' @importFrom GenomeInfoDb seqlengths seqinfo
#' @importFrom BiocParallel bplapply
#' @importClassesFrom S4Vectors Annotated
#' @importClassesFrom intervals Intervals
############################################## Main function for stage 3:
#-------------
#-------------
Stage_3_Main_fun = function(SA_prefix, S3_method, S3_AnalysisDir, S3_Selfobject,
S3_image) {
# global variables for Rcheck:
#--------------------------------------------
#---------------Take and reorder Input:
#--------------------------------------------
# Take time:
Analysis.time.start = Sys.time()
# create folder to save data:
if (!dir.exists(S3_AnalysisDir))
dir.create(S3_AnalysisDir)
#--------------------------------------------
#---------------Keep data you need only:
#--------------------------------------------
FitDataInfo = TakeFitAnalysisData_fun(S3_Selfobject = S3_Selfobject)
PETsData = FitDataInfo$PETsData #all the PETs data, used in inference too.
ChromInf = FitDataInfo$ChromInf #used in inference and fitting.
#--------------------------------------------
#---------------Segment to regions:
#--------------------------------------------
SegmRes = FindRegions_fun(PETsData = PETsData, ChromInf = ChromInf)
SegPETS = SegmRes$SegPETS #segmented PETs data, used in fitting.
RegionCounts = SegmRes$RegionCounts #take the counts
#------------
#-----save the region counts on main data:
#------------
Ordermatch = match(RegionCounts$Chrom, S4Vectors::metadata(S3_Selfobject)$Self_info$Chrom)
S4Vectors::metadata(S3_Selfobject)$Self_info$Region.counts[Ordermatch] = RegionCounts$Counts
#--------------------------------------------
#---------------Fit:
#--------------------------------------------
GlobalPeakCallRes = FitCallGlobal_fun(SegPETS = SegPETS, ChromInf = ChromInf)
#--------------------------------------------
#---------------update the results
#--------------------------------------------
#------classification information
S4Vectors::metadata(S3_Selfobject)$Classification.Info = GlobalPeakCallRes$Classification.Info
#------Peak counts
S4Vectors::metadata(S3_Selfobject)$Self_info$Peak.counts = 0
Match = match(GlobalPeakCallRes$Peak.counts$Chrom, S4Vectors::metadata(S3_Selfobject)$Self_info$Chrom)
S4Vectors::metadata(S3_Selfobject)$Self_info$Peak.counts[Match] = GlobalPeakCallRes$Peak.counts$V1
#--------------------------------------------
#---------------Run inference:
#--------------------------------------------
Peaks.Info = GlobalPeakCallRes$Peaks.Info #peaks found by fit
# run inference
InfRes = SignificanceCall_fun(Peaks.Info = Peaks.Info, PETsData = PETsData, ChromInf = ChromInf,
S3_method = S3_method)
#---add Peak information for peaks found in data with their FDR etc
S4Vectors::metadata(S3_Selfobject)$Peaks.Info = InfRes
#------------------------
#-----update class and save:
#------------------------
class(S3_Selfobject) = "PSFit"
NamepsfitData = paste(SA_prefix, "_psfitData", sep = "")
assign(NamepsfitData, S3_Selfobject) #assign value.
save(list = NamepsfitData, file = file.path(S3_AnalysisDir, NamepsfitData))
#------------------------
#-----plot images
#------------------------
if (S3_image) {
#-------------
# create data:
#-------------
PeakInfo = S4Vectors::metadata(S3_Selfobject)
PeakInfo = PeakInfo$Peaks.Info
PeakInfo$FDR0.05 = "FDR >= 0.05"
PeakInfo$FDR0.05[which(PeakInfo$FDR < 0.05)] = "FDR < 0.05"
PeakInfo$ID = seq_len(nrow(PeakInfo))
# image 1:
Get_image_S3_P1_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
# image 2:
Get_image_S3_P2_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
# image 3:
Get_image_S3_P3_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
# image 4:
Get_image_S3_P4_fun(S3_AnalysisDir = S3_AnalysisDir, SA_prefix = SA_prefix,
PeakInfo = PeakInfo)
}
#------------------------
# save log:
#------------------------
futile.logger::flog.info("=====================================", name = "SA_LogFile",
capture = FALSE)
futile.logger::flog.info("Stage 3 is done!", name = "SA_LogFile", capture = FALSE)
futile.logger::flog.info(paste("Analysis results for stage 3 are in:\n", S3_AnalysisDir),
name = "SA_LogFile", capture = FALSE)
# save time:
Analysis.time.end = Sys.time()
Total.Time = Analysis.time.end - Analysis.time.start
LogFile = paste("Total stage 3 time:", Total.Time, " ", units(Total.Time))
futile.logger::flog.info(LogFile, name = "SA_LogFile", capture = FALSE)
}
# done
#-------------
#-------------
############## Functions for inputs and data
#-------------
#-------------
# function for converting the data in the analysis format:
TakeFitAnalysisData_fun = function(S3_Selfobject) {
cat("Converting data for analysis...")
#------------
# take Chromosome information(some used in inference):
#------------
ChromInf = GenomeInfoDb::seqlengths(GenomeInfoDb::seqinfo(S3_Selfobject))
ChromInf = data.frame(Chrom = names(ChromInf), size = as.numeric(ChromInf), stringsAsFactors = FALSE)
Self_info = S4Vectors::metadata(S3_Selfobject)$Self_info
MatchChrom = match(Self_info$Chrom, ChromInf$Chrom)
ChromInf$PET.counts[MatchChrom] = Self_info$PET.counts
#------------
# convert to data frame and Swap Anchors:
#------------
PETsData = as.data.frame(S3_Selfobject)
# use Tag mids:
PETsData$UTag = (PETsData$start1 + PETsData$end1)/2
PETsData$DTag = (PETsData$start2 + PETsData$end2)/2
Swap = which(PETsData$UTag > PETsData$DTag)
if (length(Swap) != 0) {
UTagswap = PETsData$UTag[Swap]
DTagswap = PETsData$DTag[Swap]
PETsData$UTag[Swap] = DTagswap
PETsData$DTag[Swap] = UTagswap
}
#------------
# add mainindex
#------------
PETsData$MainIndex = seq_len(nrow(PETsData))
#------------
# keep left part of downstream and right of upstream
#------------
PETsData$Chrom = as.character(PETsData$seqnames1)
PETsData = PETsData[, c("Chrom", "UTag", "DTag", "MainIndex")]
#------------
# return:
#------------
cat("Done\n")
return(list(PETsData = PETsData, ChromInf = ChromInf))
}
# done
#-------------
#-------------
############## Functions for region segmentation:
#-------------
#-------------
# main function for region segmentation:
FindRegions_fun = function(PETsData, ChromInf) {
# rcheck:
Chrom = Region = NULL
#------------
#------------
# split by chromosome to list to pass to bplapply
bppass = plyr::dlply(PETsData, plyr::.(Chrom), function(x) x)
#------------
#------------
#--Segment to regions:
#------------
cat("Segmenting into regions...")
SegmRes = BiocParallel::bplapply(X = bppass, FUN = function(y, ChromInf) {
ChromCur = as.character(unique(y$Chrom))
ChromInfCur = subset(ChromInf, Chrom == ChromCur)
#------------
# Give the region ids:
#------------
y = RegionIds_fun(PETinf = y, ChromCur = ChromCur, ChromInfCur = ChromInfCur)
return(y)
}, ChromInf = ChromInf)
cat("Done\n")
#------------
#---break and rbind:
#------------
# take the region info:
PETinf = lapply(SegmRes, "[[", 1)
PETinf = plyr::ldply(PETinf, function(y) y)
PETinf = subset(PETinf, !is.na(Region)) #keep not NA
if (nrow(PETinf) != 0) {
# update data with the regions:
PETsData$Region = NA
PETsData$Region[PETinf$MainIndex] = PETinf$Region
PETsData = subset(PETsData, !is.na(Region))
} else {
stop("No regions found in data. The data is probably too sparse.", call. = FALSE)
}
#------------
#--take total regions info:
#------------
RegionCounts = lapply(SegmRes, "[[", 2)
RegionCounts = plyr::ldply(RegionCounts, function(y) y)
RegionCounts = RegionCounts[, c("Chrom", "Counts")]
futile.logger::flog.info(paste("Total Regions found:", sum(RegionCounts$Counts)),
name = "SA_LogFile", capture = FALSE)
return(list(SegPETS = PETsData, RegionCounts = RegionCounts))
}
# done
#-------------
#-------------
# Fuction for creating region ids by segmenting into regions:
RegionIds_fun = function(PETinf, ChromCur, ChromInfCur) {
# global variables for Rcheck:
#--------------------------
#------------
# initiate
#------------
PETinf$Region = NA
TotRegions = data.frame(Chrom = ChromCur, Counts = 0)
#------------
# take Pet Pmid to create global intervals:
#------------
PETInt = intervals::Intervals(PETinf[, c("UTag", "DTag")], closed = c(TRUE, TRUE))
#------------
# Take clusters:
#------------
PETclst = intervals::clusters(PETInt, which = TRUE, w = 0)
TotRegions$Counts = length(PETclst) #total regions
if (TotRegions$Counts != 0)
{
#------------
# region ids:
#------------
Region = rep(seq_len(length(PETclst)), lengths(PETclst))
#------------
# pets in the regions:
#------------
TagsIncluded = unlist(PETclst)
#------------
# add them according to the rownames of the PETinf
#------------
PETinf$Region[TagsIncluded] = Region
} #else no region found
PETinf = PETinf[, c("Chrom", "MainIndex", "Region")]
return(list(PETinf = PETinf, TotRegions = TotRegions))
}
# done
#-------------
#-------------
############## Functions for fitting the main model:
#-------------
#-------------
# main function for fitting: breaking each region
FitCallGlobal_fun = function(SegPETS, ChromInf) {
#---R-check
Region = Chrom = NULL
#----
bppass = plyr::dlply(SegPETS, plyr::.(Chrom, Region), function(y) y)
#---------------------------------------
#--------call the function in parallel for fitting regions:
#---------------------------------------
cat("Running peak calling process...")
# the following runs in parallel:
GlobalFitRes = BiocParallel::bplapply(X = bppass, FUN = FitCallLocal_fun_Rcpp,
ChromInf = ChromInf)
cat("Done\n")
#---------------------------
#--------save Classification:
#---------------------------
Peak.Id.Inf = plyr::llply(GlobalFitRes, function(x) x[[1]])
Peak.Id.Inf = do.call(rbind, Peak.Id.Inf)
if (is.null(Peak.Id.Inf)) {
stop("No Peaks found by the algorithm!", call. = FALSE)
}
#---------------------------
#--match classes with SegPETS
#---------------------------
SegPETS$Peak.ID = NA
Match = match(Peak.Id.Inf[, 2], SegPETS$MainIndex)
SegPETS$Peak.ID[Match] = Peak.Id.Inf[, 1]
Classification.Info = SegPETS[, c("MainIndex", "Region", "Peak.ID")]
#---------------------------
#--------save Peaks:
#---------------------------
Peaks.Info = plyr::ldply(GlobalFitRes, function(x) x[[2]])
if (is.null(Peaks.Info)) {
stop("No Peaks found by the algorithm!", call. = FALSE)
}
Peaks.Info = Peaks.Info[, which(colnames(Peaks.Info) != ".id")]
futile.logger::flog.info(paste("Total", nrow(Peaks.Info), "candidate peaks found in data"),
name = "SA_LogFile", capture = FALSE)
#---------------------------
#--------Find counts:
#---------------------------
Peak.counts = plyr::ddply(Peaks.Info, plyr::.(Chrom), nrow)
return(list(Peaks.Info = Peaks.Info, Classification.Info = Classification.Info,
Peak.counts = Peak.counts))
}
# done
#-------------
#-------------
############## Functions for the Inference:
#-------------
#-------------
# Main significance call function:
SignificanceCall_fun = function(Peaks.Info, PETsData, ChromInf, S3_method) {
# global variables for Rcheck:
Chrom = NULL
#--------------------------
#--------------------------------------------
#---------------Take information for the analysis:
#--------------------------------------------
Peaks.Info$Chrom = as.character(Peaks.Info$Chrom)
PETsData$MainIndex = NULL #remove index not needed
#--------------------------------------------
#---------------run inference:
#--------------------------------------------
cat("Splitting data by chromosome for inference...\n")
# split by chromosome all the three inputs: then merge them by chromosome
bppass = Get_SignificanceCall_Data_fun(Peaks.Info = Peaks.Info, PETsData = PETsData,
ChromInf = ChromInf)
#------------
# Run inference(call c++):
#------------
futile.logger::flog.info("Computing p-values...", name = "SA_LogFile", capture = FALSE)
InfRes = BiocParallel::bplapply(X = bppass, FUN = PoissonLocalInference_fun,
windows = c(10, 15))
#------------
# merge and return:
#------------
InfRes = do.call(rbind.data.frame, InfRes)
rownames(InfRes) = NULL
#------
# FDR at 0.05:
#------
futile.logger::flog.info("FDR adjusting p-values...", name = "SA_LogFile", capture = FALSE)
InfRes$FDRUp = stats::p.adjust(p = InfRes$p.valueUp, method = S3_method)
InfRes$FDRDown = stats::p.adjust(p = InfRes$p.valueDown, method = S3_method)
InfRes$FDR = stats::p.adjust(p = InfRes$p.value, method = S3_method)
#------
#---return:
#------
return(InfRes)
}
# done
#-------------
#-------------
# Function for breaking the inputs for significance calling:
Get_SignificanceCall_Data_fun = function(Peaks.Info, PETsData, ChromInf) {
# Rcheck:
Chrom = NULL
#-------------
# break Peaks.Info by chrom:
#-------------
Peaks.Info_list = plyr::dlply(Peaks.Info, plyr::.(Chrom), function(x) x)
#-------------
# break PETsData by chrom:
#-------------
PETsData$MainIndex = NULL #dont need index
PETsData_list = plyr::dlply(PETsData, plyr::.(Chrom), function(x) x)
#-------------
# break ChromInf by chrom:
#-------------
ChromInf_list = plyr::dlply(ChromInf, plyr::.(Chrom), function(x) x)
#-------------
# create the bppass input:
#-------------
bppass = list()
for (chri in seq_along(Peaks.Info_list)) {
# take Peaks_Info_x:
Peaks_Info_x = Peaks.Info_list[[chri]]
Chr_x = Peaks_Info_x$Chrom[1]
# take PETsData_x:
PETsData_x = PETsData_list[[which(names(PETsData_list) == Chr_x)]]
PETsData_x$Chrom = NULL #dont need that anymore
PETsData_x = as.matrix(PETsData_x)
# take ChromInf_x:
ChromInf_x = ChromInf_list[[which(names(ChromInf_list) == Chr_x)]]
ChromInf_x$Chrom = NULL #dont need that
bppass[[chri]] = list(Peaks_Info_x = Peaks_Info_x, PETsData_x = PETsData_x,
ChromInf_x = ChromInf_x)
}
return(bppass)
}
# done
#-------------
#-------------
# Function for Inference:
PoissonLocalInference_fun = function(bppass_x, windows) {
#-------------
# break data:
#-------------
Peaks_Info_x = bppass_x$Peaks_Info_x
PETsData_x = bppass_x$PETsData_x
# take UTag
UTag_x = intervals::Intervals(PETsData_x[, c(1, 1)], closed = c(TRUE, TRUE))
# take DTag
DTag_x = intervals::Intervals(PETsData_x[, c(2, 2)], closed = c(TRUE, TRUE))
# take chromosome info:
ChromSize = bppass_x$ChromInf_x$size
ChromPETs = bppass_x$ChromInf_x$PET.counts
#-------------
# ------ Up window 1:
#-------------
UpW1 = cbind(Peaks_Info_x$CIQ.Up.start - Peaks_Info_x$CIQ.Up.size * windows[1]/2,
Peaks_Info_x$CIQ.Up.start + Peaks_Info_x$CIQ.Up.size * windows[1]/2)
UpW1 = intervals::Intervals(UpW1, closed = c(TRUE, TRUE))
# observed pets in intervals
UpW1_obs = lengths(intervals::interval_included(UpW1, UTag_x))
lambdaUpW1 = UpW1_obs * Peaks_Info_x$CIQ.Up.size/(intervals::size(UpW1) + 1)
#-------------
# ------ Up window 2:
#-------------
UpW2 = cbind(Peaks_Info_x$CIQ.Up.start - Peaks_Info_x$CIQ.Up.size * windows[2]/2,
Peaks_Info_x$CIQ.Up.start + Peaks_Info_x$CIQ.Up.size * windows[2]/2)
UpW2 = intervals::Intervals(UpW2, closed = c(TRUE, TRUE))
# observed pets in intervals
UpW2_obs = lengths(intervals::interval_included(UpW2, UTag_x))
lambdaUpW2 = UpW2_obs * Peaks_Info_x$CIQ.Up.size/(intervals::size(UpW2) + 1)
#-------------
# ------ Down window 1:
#-------------
DownW1 = cbind(Peaks_Info_x$CIQ.Down.start - Peaks_Info_x$CIQ.Down.size * windows[1]/2,
Peaks_Info_x$CIQ.Down.start + Peaks_Info_x$CIQ.Down.size * windows[1]/2)
DownW1 = intervals::Intervals(DownW1, closed = c(TRUE, TRUE))
# observed pets in intervals
DownW1_obs = lengths(intervals::interval_included(DownW1, DTag_x))
lambdaDownW1 = DownW1_obs * Peaks_Info_x$CIQ.Down.size/(intervals::size(DownW1) +
1)
#-------------
# ------ Down window 2:
#-------------
DownW2 = cbind(Peaks_Info_x$CIQ.Down.start - Peaks_Info_x$CIQ.Down.size * windows[2]/2,
Peaks_Info_x$CIQ.Down.start + Peaks_Info_x$CIQ.Down.size * windows[2]/2)
DownW2 = intervals::Intervals(DownW2, closed = c(TRUE, TRUE))
# observed pets in intervals
DownW2_obs = lengths(intervals::interval_included(DownW2, DTag_x))
lambdaDownW2 = DownW2_obs * Peaks_Info_x$CIQ.Down.size/(intervals::size(DownW2) +
1)
#-------------
# ------ genome backgrounds:
#-------------
lambdaBGC_Up = ChromPETs * Peaks_Info_x$CIQ.Up.size/ChromSize
lambdaBGC_Down = ChromPETs * Peaks_Info_x$CIQ.Down.size/ChromSize
#-------------
# ------ find maximum:
#-------------
OptimalLambda = plyr::llply(seq_len(nrow(Peaks_Info_x)), function(i, lambdaUpW1,
lambdaUpW2, lambdaBGC_Up, lambdaDownW1, lambdaDownW2, lambdaBGC_Down) {
UpMax = max(c(2, lambdaUpW1[i], lambdaUpW2[i], lambdaBGC_Up[i]))
DownMax = max(c(2, lambdaDownW1[i], lambdaDownW2[i], lambdaBGC_Down[i]))
return(c(UpMax, DownMax))
}, lambdaUpW1 = lambdaUpW1, lambdaUpW2 = lambdaUpW2, lambdaBGC_Up = lambdaBGC_Up,
lambdaDownW1 = lambdaDownW1, lambdaDownW2 = lambdaDownW2, lambdaBGC_Down = lambdaBGC_Down)
OptimalLambda = do.call(rbind, OptimalLambda)
#-------------
# ------ Find folds and p.values:
#-------------
# Upstream
FoldEnrichUp = Peaks_Info_x$Pets/OptimalLambda[, 1]
pvalueUp = stats::ppois(q = Peaks_Info_x$Pets, lambda = OptimalLambda[, 1], lower.tail = FALSE)
# Downstream
FoldEnrichDown = Peaks_Info_x$Pets/OptimalLambda[, 2]
pvalueDown = stats::ppois(q = Peaks_Info_x$Pets, lambda = OptimalLambda[, 2],
lower.tail = FALSE)
# merged p-value:
pvalue = pvalueUp * pvalueDown
#-------------
# ------Return
#-------------
Peaks_Info_x$lambdaUp = OptimalLambda[, 1]
Peaks_Info_x$FoldEnrichUp = FoldEnrichUp
Peaks_Info_x$p.valueUp = pvalueUp
Peaks_Info_x$lambdaDown = OptimalLambda[, 2]
Peaks_Info_x$FoldEnrichDown = FoldEnrichDown
Peaks_Info_x$p.valueDown = pvalueDown
Peaks_Info_x$p.value = pvalue
return(Peaks_Info_x)
}
#-------------
#-------------
# done Functions for plotting:
#-------------
#-------------
Get_image_S3_P1_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
CIQ.Up.size = NULL
CIQ.Down.size = NULL
FDR0.05 = NULL
# image dir:
S3_P1_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p1_image.jpg",
sep = ""))
# plot
S3_image_p1 = ggplot2::ggplot(PeakInfo, ggplot2::aes(CIQ.Up.size, CIQ.Down.size,
color = factor(FDR0.05))) + ggplot2::geom_point(size = 1) + ggplot2::theme_bw() +
ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"), panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(), panel.background = ggplot2::element_blank()) +
ggplot2::theme(axis.text = ggplot2::element_text(size = 15, color = "black"),
axis.title = ggplot2::element_text(size = 18, color = "black"), plot.title = ggplot2::element_text(size = 18,
color = "black"), legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_text(size = 18)) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_color_discrete(guide = ggplot2::guide_legend(title = "Binding Site FDR")) +
ggplot2::scale_y_continuous(expand = c(0, 0)) + ggplot2::scale_x_continuous(expand = c(0,
0)) + ggplot2::xlab("Upstream Peak Size") + ggplot2::ylab("Downstream Peak Size") +
ggplot2::geom_abline(intercept = 0, slope = 1) + ggplot2::ggtitle("Upstream/Downstream peak sizes")
# save:
ggplot2::ggsave(plot = S3_image_p1, file = S3_P1_image_dir, scale = 2)
}
#-------------
#-------------
Get_image_S3_P2_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
ID = NULL
FDR = NULL
# image dir:
S3_P2_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p2_image.jpg",
sep = ""))
S3_image_p2 = ggplot2::ggplot(PeakInfo, ggplot2::aes(ID, sort(FDR, decreasing = FALSE))) +
ggplot2::geom_line(size = 1) + ggplot2::geom_hline(yintercept = 0.05, color = "red",
size = 1) + ggplot2::theme_bw() + ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank()) + ggplot2::theme(axis.text = ggplot2::element_text(size = 15,
color = "black"), axis.title = ggplot2::element_text(size = 18, color = "black"),
plot.title = ggplot2::element_text(size = 18, color = "black")) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_y_continuous(expand = c(0, 0)) + ggplot2::scale_x_continuous(expand = c(0,
0)) + ggplot2::xlab("Sorted Binding Site IDs") + ggplot2::ylab("FDR") + ggplot2::ggtitle("FDR values for the binding sites")
# save:
ggplot2::ggsave(plot = S3_image_p2, file = S3_P2_image_dir, scale = 2)
}
#-------------
#-------------
Get_image_S3_P3_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
CIQ.Peak.size = NULL
FDR0.05 = NULL
# image dir:
S3_P3_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p3_image.jpg",
sep = ""))
# plot
S3_image_p3 = ggplot2::ggplot(PeakInfo, ggplot2::aes(CIQ.Peak.size, fill = factor(FDR0.05))) +
ggplot2::geom_density(size = 1, alpha = 0.5) + ggplot2::theme_bw() + ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank()) + ggplot2::theme(axis.text = ggplot2::element_text(size = 15,
color = "black"), axis.title = ggplot2::element_text(size = 18, color = "black"),
plot.title = ggplot2::element_text(size = 18, color = "black"), legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_blank()) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_y_continuous(expand = c(0, 0), labels = function(y) format(y,
scientific = FALSE)) + ggplot2::scale_x_continuous(expand = c(0, 0)) +
ggplot2::xlab("Binding Site Sizes") + ggplot2::ylab("Density") + ggplot2::geom_abline(intercept = 0,
slope = 1) + ggplot2::ggtitle("Densities for the binding site sizes")
# save:
ggplot2::ggsave(plot = S3_image_p3, file = S3_P3_image_dir, scale = 2)
}
#-------------
#-------------
Get_image_S3_P4_fun = function(S3_AnalysisDir, SA_prefix, PeakInfo) {
# Rcheck:
FDRUp = NULL
FDRDown = NULL
FDR0.05 = NULL
# image dir:
S3_P4_image_dir = file.path(S3_AnalysisDir, paste(SA_prefix, "_stage_3_p4_image.jpg",
sep = ""))
# plot
S3_image_p4 = ggplot2::ggplot(PeakInfo, ggplot2::aes(FDRUp, FDRDown, color = factor(FDR0.05))) +
ggplot2::geom_point(size = 1) + ggplot2::theme_bw() + ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank()) + ggplot2::theme(axis.text = ggplot2::element_text(size = 15,
color = "black"), axis.title = ggplot2::element_text(size = 18, color = "black"),
plot.title = ggplot2::element_text(size = 18, color = "black"), legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_text(size = 18)) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
ggplot2::scale_color_discrete(guide = ggplot2::guide_legend(title = "Binding Site FDR")) +
ggplot2::scale_y_continuous(expand = c(0, 0)) + ggplot2::scale_x_continuous(expand = c(0,
0)) + ggplot2::xlab("FDR for Upstream Binding Site Peak") + ggplot2::ylab("FDR for Downstream Binding Site Peak") +
ggplot2::geom_abline(intercept = 0, slope = 1) + ggplot2::ggtitle("Upstream/Downstream peak FDR")
# save:
ggplot2::ggsave(plot = S3_image_p4, file = S3_P4_image_dir, scale = 2)
}
Try the MACPET package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.