R/norm.quantile.R

In xuz1/ENmix: Quality control and analysis tools for Illumina DNA methylation BeadChip

normalize.quantiles2 <-function(mat){
    if(any(is.na(mat))){
	    xm=is.na(mat)
	    mat[is.na(mat)]=mean(mat,na.rm=TRUE)
            ref=apply(apply(mat, 2, sort),1,mean)
            mat=apply(mat,2,function(x)ref[rank(x)])
            mat[xm]=NA
            mat}else{
            ref=apply(apply(mat, 2, sort),1,mean)
            apply(mat,2,function(x)ref[rank(x)])
    }
}

normalize.quantiles.use.target2 <-function(mat,target){
    if(any(is.na(mat))){
            xm=is.na(mat)
            mat[is.na(mat)]=mean(mat,na.rm=TRUE)
            ref=sort(target)
            mat=apply(mat,2,function(x)ref[rank(x)])
            mat[xm]=NA
            mat}else{
            ref=sort(target)
            apply(mat,2,function(x)ref[rank(x)])
    }
}


normalize.q <- function(x)
{
    if(sum(is.na(x))>0)
    {stop("Multi-array quantile normalization does not allow missing value")}
    rname <- rownames(x)
    cname <- colnames(x)
    x <- normalize.quantiles2(x)
    rownames(x) <- rname
    colnames(x) <- cname
    x
}

norm.quantile <- function(mdat,method="quantile1")
{
#   if(!require("preprocessCore")){stop("Can not load preprocessCore package")}
    if(!is(mdat, "methDataSet") & !is(mdat, "MethylSet")){
      stop("object needs to be of class 'methDataSet' or 'MethylSet'")}
    cat("Analysis is running, please wait...!","\n")

    if(is(mdat,"MethylSet")){
      anno <- getAnnotation(mdat)
    }else if(is(mdat,"methDataSet")){
      anno=rowData(mdat)
      names(anno)[which(names(anno)=="Infinium_Design_Type")]="Type"
      anno$Type[anno$Type %in% "snpI"]="I"
      anno$Type[anno$Type %in% "snpII"]="II"
    }

    if (method == "quantile1")
    {
    mdat_I <- mdat[anno$Type == "I",]
    mdat_II <- mdat[anno$Type == "II",]
    assays(mdat_I)$Meth<-normalize.q(assays(mdat_I)$Meth)
    assays(mdat_I)$Unmeth<-normalize.q(assays(mdat_I)$Unmeth)
    assays(mdat_II)$Meth<-normalize.q(assays(mdat_II)$Meth)
    assays(mdat_II)$Unmeth<-normalize.q(assays(mdat_II)$Unmeth)
    methData <- rbind(assays(mdat_I)$Meth,assays(mdat_II)$Meth)
    unmethData <- rbind(assays(mdat_I)$Unmeth,assays(mdat_II)$Unmeth)
    CpGID <- rownames(mdat)
    SampleID <- colnames(mdat)
    methData <- methData[CpGID,]
    unmethData <- unmethData[CpGID,]
    assays(mdat)$Meth<-methData
    assays(mdat)$Unmeth<-unmethData
    }
    else if(method == "quantile2")
    {
    mdat_I <- mdat[anno$Type == "I",]
    mdat_II <- mdat[anno$Type == "II",]
    mat <- rbind(assays(mdat_I)$Meth,assays(mdat_I)$Unmeth)
    mat  <-  normalize.q(mat)
    assays(mdat_I)$Meth <- mat[1:(nrow(mat)/2),]
    assays(mdat_I)$Unmeth <- mat[((nrow(mat)/2)+1):nrow(mat),]
    mat <- rbind(assays(mdat_II)$Meth,assays(mdat_II)$Unmeth)
    mat  <-  normalize.q(mat)
    assays(mdat_II)$Meth <- mat[1:(nrow(mat)/2),]
    assays(mdat_II)$Unmeth <- mat[((nrow(mat)/2)+1):nrow(mat),]
    methData <- rbind(assays(mdat_I)$Meth,assays(mdat_II)$Meth)
    unmethData <- rbind(assays(mdat_I)$Unmeth,assays(mdat_II)$Unmeth)
    CpGID <- rownames(mdat)
    SampleID <- colnames(mdat)
    methData <- methData[CpGID,]
    unmethData <- unmethData[CpGID,]
    assays(mdat)$Meth<-methData
    assays(mdat)$Unmeth<-unmethData

    }
    else if(method == "quantile3")
    {
    ##quantile normalization of combined intensity values
    ##similar with lumi lumimethyN
    mat <- rbind(assays(mdat)$Meth,assays(mdat)$Unmeth)
    mat  <-  normalize.q(mat)
    methData <- mat[1:(nrow(mat)/2),]
    unmethData <- mat[((nrow(mat)/2)+1):nrow(mat),]
    assays(mdat)$Meth<-methData
    assays(mdat)$Unmeth<-unmethData
    }
    if(is(mdat, "methDataSet")){
      annotation=metadata(mdat)$preprocessMethod
      metadata(mdat)$preprocessMethod=c(annotation,paste("Normalization: ", method,sep=""))
    }else if(is(mdat, "MethylSet")){
      annotation=mdat@preprocessMethod
      mdat@preprocessMethod <- c(annotation,paste("Normalization: ", method,sep=""))
    }
    mdat
}