#' Significant Calling and summary output
#'
#' Judge whether PWM is significantly overrepresented at the DMP regions.
#' The Judgement is made by Poisson p-value at DMP position and enrichment score at DMP position which is more than Q3+IQR*3 of enrichment score of analyzed range.
#'
#' @param significant_ranges generated by SigRange
#' @param motif_counts_matrix generated by countMatrix
#' @param seq_range Analyzed range from DMP
#' @param enrichment_scores enrichScoreDist
#' @param windowSize windowsize used for p-value coumutation etc... (should be same as windowSize parameter of countMatrix)
#'
#' @importFrom stats IQR quantile
#'
#' @return summary of analysis and signicifant call
#' @keywords summary, significant call
#' @export
#'
enrichTest <- function(significant_ranges, motif_counts_matrix, seq_range, enrichment_scores, windowSize){
if(!is.null(significant_ranges)){ #if significant range exist
centSigRangeInd <- which(significant_ranges[,1] < 0 & significant_ranges[,2] > 0) #extract index (row numbers) of significanrly enriched rgion which is overlapped with CpG (center)
centSigRange <- significant_ranges[centSigRangeInd,] #center peak range
}else{
centSigRangeInd <- NULL
}
if(length(centSigRangeInd) != 0){
cat (" ----Poison Test: Significant----\n\n")
##Max value (FC) in the center peak
rowNoMaxClass <- which(motif_counts_matrix[,3] == max(motif_counts_matrix[which(as.numeric(rownames(motif_counts_matrix)) >= centSigRange[1] & as.numeric(rownames(motif_counts_matrix)) <= centSigRange[2]),3])) #row index of max FC
if( length(rowNoMaxClass) > 1){ #if the min Pvalue is not single, adopt a class which is the closest to 0
rowNoMaxClass <- which(rownames(motif_counts_matrix) == min(rownames(motif_counts_matrix)[rowNoMaxClass]))
}
maxClass <- rownames(motif_counts_matrix)[rowNoMaxClass]
maxFC <- motif_counts_matrix[rowNoMaxClass, 3]
poisPvals <- min(motif_counts_matrix[which(as.numeric(rownames(motif_counts_matrix)) > significant_ranges[centSigRangeInd,1] & as.numeric(rownames(motif_counts_matrix)) < significant_ranges[centSigRangeInd,2]),4])
poispeakPosiEvals <- "Overlapped"
##TEST: Mixture distribution based test
ranks <- seq(seq_range[1], seq_range[2], length=1001)
centPeakProb <- enrichment_scores[which(ranks >centSigRange[1] & ranks < centSigRange[2])] #"enrichments" (values of mixture distribution) of center peak
max_value <- max(centPeakProb) #max of "enrichment" in the center peak
max_ratio <- (max_value - quantile(enrichment_scores)[4])/IQR(enrichment_scores)
outlier <- as.numeric(quantile(enrichment_scores)[4]) + (IQR(enrichment_scores)*3) #windowsize used for p-value coumutation (should be same as windowSize parameter of countMatrix)# set of outlier cut-off (Q3+IQR*3)
enrichment_Q1<- quantile(enrichment_scores)[2]
enrichment_Med <-quantile(enrichment_scores)[3]
enrichment_Q3 <- quantile(enrichment_scores)[4]
enrichment_IQR <- IQR(enrichment_scores)
centPeakRatio <- max_ratio
centPeakOutL <- outlier
centPeakMax <- max_value
##main body of judgement of center peak test
if (max_value > outlier){
peakMaxEval <- "significant"
cat (" ----local enrichment: Significant----\n\n")
} else {
peakMaxEval <- "NS"
cat (" ----local enrichment: Insignificant----\n\n")
}
}else{
##case of no center peak
cat (" ----Poison Test: Insignificant----\n\n")
poisPvals <- min(motif_counts_matrix[which(as.numeric(rownames(motif_counts_matrix)) > (0-(windowSize/2)) & as.numeric(rownames(motif_counts_matrix)) < (0+(windowSize/2))),4]) #put min P of cnter window
poispeakPosiEvals <- "Not-Overlapped" #Judgement of center peak
maxClass <- "NA"
maxFC <- "NA"
centPeakMax <- "NA"
enrichment_Q1<- "NA"
enrichment_Med <-"NA"
enrichment_Q3 <- "NA"
enrichment_IQR <- "NA"
centPeakRatio <- "NA"
centPeakOutL <- "NA"
peakMaxEval <- "NA"
}
parameters <- c(poispeakPosiEvals, maxClass, maxFC, poisPvals, enrichment_Q1, enrichment_Med, enrichment_Q3, enrichment_IQR, centPeakOutL, centPeakMax, centPeakRatio, peakMaxEval)
names(parameters) <- c("Posison.Test", "max.class", "mac.FC", "Poison.Test.P.value", "Q1", "Median", "Q3", "IQR", "Outlier","max.peak", "peak.IQR.ratio", "peak.Test")
return (parameters)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.