R/enrichTest.R

In takahirosuzuki1980/InfiniumDiffMetMotR: Motif Overrepresentation analysis for Illumina Methylation Arrays

#' 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)
}