R/Multiclass-NONE.R

In YingZhang0811/test11: R Package for Assessing Normalization for Metabolomic Data

#' @title Multi-class (N>1) Metabolomic Study with dataset without QCSs and ISs.
#' @description This function handles (1) normalizing the multi-class metabolomic
#' data without QCSs and ISs using 168 methods/strategies,
#' (2) evaluating the normalization performances from multiple perspectives,
#' and (3) enabling the systematic comparison among all methods/strategies
#' based on a comprehensive performance ranking.
#' @param fileName Please find the detail information of the file format from those six sample
#' files in the working directory (in github) “idrblab/METNOR/data”
#' @param impt Input the name of imputation methods.
#' If set 1, method of column mean imputation.
#' If set 2, method of column median imputation.
#' If set 3, method of half of the minimum positive value.
#' If set 4, method of KNN imputation.
#' @param trsf Input the name of transformation methods.
#' If set 1, method of cube root transformation.
#' If set 2, method of log transformation.
#' If set 3, none transformation method.
#' @import DiffCorr affy vsn DT ropls
#' @import e1071 AUC impute MetNorm
#' @import ggsci timecourse multiROC dummies
#' @import ggplot2 ggord ggfortify usethis
#' @import ggrepel ggpubr sampling crmn
#' @rawNamespace import(limma, except=.__C__LargeDataObject)
#' @importFrom grDevices dev.off png rainbow rgb
#' @importFrom systemPipeR overLapper
#' @importFrom graphics abline close.screen legend mtext par points screen split.screen symbols text title
#' @importFrom stats anova as.formula cor dnorm kmeans lm loess loess.control mad median model.matrix na.omit pf pnorm qnorm qt quantile rnorm runif sd var
#' @importFrom utils combn read.csv write.csv write.table
#' @examples allranks_non <- normulticlassnoall(fileName = multi_non_data, impt = "1", trsf = "1")
#' @export normulticlassnoall

normulticlassnoall <- function(fileName, impt = NULL, trsf = NULL){

  cat("METNOR is Running ...","\n")
  cat("\n")

  consistency_M <- function(fold = 3, top = 20) {

    folds <- fold

    DEG <- list()
    for (i in 1:folds) {

      com.x <- test_data[test.fold[[i]], ]

      com.x[sapply(com.x, simplify = 'matrix', is.infinite)] <- 0

      X_matrix <- as.data.frame(com.x[,-1])
      y_label <- as.factor(com.x[,1])

      pos_filter <- OPLSDA_C(X_matrix, y_label)

      DEG[[i]] <- pos_filter
    }
    names(DEG) <- LETTERS[1:folds]

    top.n <-top # Extracting the top n genes.
    DEG.list <- DEG
    for (g in 1:length(DEG.list)) {
      DEG.list[[g]] <- DEG.list[[g]][1:top.n]
    }

    # Calculating consistency score:

    setlist <- DEG.list


    return(setlist) # consistense score

  }

  imput<-function(filter_data2,n){
    matrix<-switch(
      n,
      "1" = t(ImputMean(filter_data2)),#
      "2" = t(ImputMedian(filter_data2)),#
      "3" = t(back(filter_data2)),#
      "4" = t(impute.knn(as.matrix(t(filter_data2)), k = 10,rng.seed = 1024)$data)#
    )
    return(matrix)
  }

  im_nam<-c(
    "Mean",
    "Median",
    "Half Minimum",
    "KNN"
  )

  #Transformation---------------------------------------------------------------------------------
  trans<-function(data,n){
    matrix<-switch(
      n,
      "1" = Cube_root(data),
      "2" = log2(data),
      "3" = data
    )
    return(matrix)
  }

  t_nam<-c(
    "Cube_root",
    "Log",
    "None"
  )

  #normalization---------------------------------------------------------------------------------
  norm_sam<-function(train_data_t,n){
    matrix<-switch(
      n,
      train_data_t,#1
      PQN(train_data_t),#2
      fastlo(train_data_t),#3
      CONTRAST(train_data_t),#4
      QUANTILE(train_data_t),#5
      LINEAR(train_data_t),#6
      LIWONG(train_data_t),#7
      CUBIC(train_data_t),#8
      AUTO(train_data_t),#9
      RANGE(train_data_t),#10
      PARETO(train_data_t),#11
      VAST(train_data_t),#12
      LEVEL(train_data_t),#13
      VSN(train_data_t),#14
      POWER(train_data_t),#15
      MSTUS(train_data_t),#16
      SUM(train_data_t),#17
      MEDIAN(train_data_t),#18
      MEAN(train_data_t),#19
      EIGENMS(train_data_t, sampleLabel)#20
    )
    return(matrix)
  }

  n_nam<-c(
    "NON",
    "PQN",
    "LOE",
    "CON",
    "QUA",
    "LIN",
    "LIW",
    "CUB",
    "AUT",
    "RAN",
    "PAR",
    "VAS",
    "LEV",
    "VSN",
    "POW",
    "MST",
    "SUM",
    "MED",
    "MEA",
    "EIG"
  )

  #-----------------------------------------------------------------
  ###################################################Step-2 调用数据
  data_q <- fileName
  data2 <- data_q[, -c(1:2)]
  data4 <- data_q[, 1:2]
  data2[data2 == 0] <- NA # the zero value has been replaced by NA.

  #filtering-----------------------------------------------------------------------------
  col_f <- apply(data2, 2, function(x) length(which(is.na(x)))/length(x))
  if (length(which(col_f >0.2))==0){
    data2_f <- data2
  }else {
    data2_f <- data2[, -which(col_f > 0.2)]
  }
  filter_data2 <- data2_f

  #---------------------------------------------------------------------------------

  aftetable<-list()
  normal_data <- list()
  train_data_metabolite3 <- NULL

  for (i in as.numeric(impt)){

    afterimpute.table <- NULL
    imput_m <- try(imput(filter_data2,i))
    if(class(imput_m)=="try-error")
    { next }

    imputed_data <- cbind(data4, imput_m)
    afterimpute.table <- imputed_data
    getLabel <-as.character(data_q[, 2])
    data1 <- afterimpute.table
    train_data_t <- t(data1[, -(1:2)])

    for (j in as.numeric(trsf)){

      train_data_Transformation3<-try(trans(train_data_t,j))
      if(class(train_data_Transformation3)=="try-error")
      { next }
      train_data_Transformation3[is.infinite(data.matrix(train_data_Transformation3))]<-NA
      sampleLabel <- as.character(afterimpute.table[, 2])

      imputed_data <- cbind(data4, t(train_data_Transformation3))
      after.table <- imputed_data

      ####
      #sink(file=paste("OUTPUT-METNOR-Record",".txt",sep=""))

      for ( k in c(2:8,16:20)){
        train_data_Preprocess3 <-try(norm_sam(train_data_Transformation3,k))
        if(class(train_data_Preprocess3)=="try-error")
        { next }

        for(h in c(1,9:13,15)){

          train_data_metabolite3<-try(norm_sam(train_data_Preprocess3,h))
          if(class(train_data_metabolite3)=="try-error")
          { next }

          normalized_data3 <- try(t(train_data_metabolite3))
          if(class(normalized_data3)=="try-error")
          { next }

          eva.data3 <- cbind(afterimpute.table[, 1:2], normalized_data3)
          eva.data3 <- eva.data3[, -1]
          eva.data3 <- as.data.frame(eva.data3)
          rownames(eva.data3) <- afterimpute.table[, 1]
          colnames(eva.data3)[1] <- "Group"
          eva_data3<-eva.data3


          normal_data[[paste(n_nam[k],n_nam[h],sep="+")]] <- eva_data3

          eva_data3 <- NULL

          aftetable[[paste(n_nam[k],n_nam[h],sep="+")]] <- after.table


        }
      }
      for(k in c(9:13,15)){

        train_data_Preprocess3 <-try(norm_sam(train_data_Transformation3,k))
        if(class(train_data_Preprocess3)=="try-error")
        { next }
        for(h in c(1,2:8,16:20)){

          train_data_metabolite3<-try(norm_sam(train_data_Preprocess3,h))
          if(class(train_data_metabolite3)=="try-error")
          { next }

          normalized_data3 <- try(t(train_data_metabolite3))
          if(class(normalized_data3)=="try-error")
          { next }

          eva.data3 <- cbind(afterimpute.table[, 1:2], normalized_data3)
          eva.data3 <- eva.data3[, -1]
          eva.data3 <- as.data.frame(eva.data3)
          rownames(eva.data3) <- afterimpute.table[, 1]
          colnames(eva.data3)[1] <- "Group"
          eva_data3<-eva.data3

          normal_data[[paste(n_nam[k],n_nam[h],sep="+")]] <- eva_data3

          eva_data3 <- NULL

          aftetable[[paste(n_nam[k],n_nam[h],sep="+")]] <- after.table


        }
      }

      for(k in c(1,14)){

        train_data_metabolite3<-try(norm_sam(train_data_Transformation3,k))
        if(class(train_data_metabolite3)=="try-error")
        { next }

        normalized_data3 <- try(t(train_data_metabolite3))
        if(class(normalized_data3)=="try-error")
        { next }

        eva.data3 <- cbind(afterimpute.table[, 1:2], normalized_data3)
        eva.data3 <- eva.data3[, -1]
        eva.data3 <- as.data.frame(eva.data3)
        rownames(eva.data3) <- afterimpute.table[, 1]
        colnames(eva.data3)[1] <- "Group"
        eva_data3<-eva.data3

        normal_data[[paste(n_nam[k],n_nam[h],sep="+")]] <- eva_data3

        eva_data3 <- NULL

        aftetable[[paste(n_nam[k],n_nam[h],sep="+")]] <- after.table

      }

    }
  }

  #sink()
  nanmes_right<-names(normal_data)

  Fpmad<-list()
  Fpurity<-list()
  Fscore<-list()
  Fauc<-list()

  for (mmm in 1:length(normal_data)){

    name <- nanmes_right

    ###Fpmad-------------------------------------------------------------------------
    n_data <- as.data.frame(normal_data[mmm],col.names=NULL)
    n_data <- as.matrix(n_data)
    if(sum(is.na(n_data))<length(n_data)/3){
      eva_data3<-as.data.frame(normal_data[mmm],col.names=NULL)
    }else{next}

    pmad3N.log <- eva_data3
    pmad3N <- try(PMAD(pmad3N.log))
    if(class(pmad3N)=="try-error")
    { next }

    Fpmad[names(normal_data[mmm])] <- mean(pmad3N)

    names(after.table)[1] <- "Group"

    pmad3R.log <- after.table[,-1]
    pmad3R <- try(PMAD(pmad3R.log))
    if(class(pmad3R)=="try-error")
    { next }

    C3 <- cbind(pmad3R, pmad3N); colnames(C3) <- c("Before", "After")

    cat(paste("Assessing Method" , paste(mmm,":",sep=""), name[mmm]),"\n")

    cat("   Criterion Ca (reduction of intragroup variation) ...","\n")

    # 2. Fpurity-------------------------------------------------------------------------

    data_kmeans <- eva_data3
    data_kmeans[sapply(data_kmeans, simplify = 'matrix', is.na)] <- 0

    eva_data3_f<-data_kmeans
    del_col <- NULL
    for (m in 1:dim(eva_data3_f)[2]){
      if (sum(eva_data3_f[,m]==0)==nrow(eva_data3_f)) {
        del_col <- c(del_col,m)
      }
    }

    if (is.null(del_col)){
      eva_data3_f <- eva_data3_f
    }else{
      eva_data3_f <- eva_data3_f[,-del_col]
    }

    X_matrix<-eva_data3_f[,-1]
    Group <- as.factor(eva_data3_f$Group)

    sink(file=paste("OUTPUT-METNOR-Record",".txt",sep=""))

    pos_filter <- try(OPLSDA_test(X_matrix, Group, cutoff = 0.8))

    sink()

    if(class(pos_filter)=="try-error")
    { next }

    DEG <- cbind(Group,X_matrix[, pos_filter])
    data_kmeans<-DEG

    clusters <- length(unique(data_kmeans[, 1]))
    obj_kmeans <- try(kmeans(data_kmeans[,-1], centers = clusters, nstart = 1))
    if(class(obj_kmeans)=="try-error")
    { next }

    groups <- factor(data_kmeans[, 1], levels = unique(data_kmeans[,1]))
    unique.groups <- levels(groups)
    cols <- 1:length(unique(data_kmeans[, 1]))
    box_cols <- NULL

    for (ii in 1:length(data_kmeans[, 1])) {

      box_cols[ii] <- cols[match(data_kmeans[, 1][ii],unique.groups)]
    }

    true_label<-box_cols

    pre<-obj_kmeans$cluster
    tru<-true_label
    tmatrix<-table(pre,tru)

    label<-tru
    result<-pre

    accuracy<-try(purity(result,label))
    if(class(accuracy)=="try-error")
    { next }

    Fpurity[names(normal_data[mmm])]<-accuracy

    unique.groups <- levels(as.factor(data_kmeans[,1]))
    T_number <- length(unique.groups)
    cols <- rainbow(length(unique.groups))

    box_cols <- c(rep(NA, length(rownames(data_kmeans))))

    for (ii in 1:length(data_kmeans[, 1])) {
      box_cols[ii] <- cols[match(data_kmeans[, 1][ii],unique.groups)]
    }

    data_kmeans <- as.data.frame(data_kmeans)
    data_kmeans$Color<-box_cols

    data_kmeans$T_label <- data_kmeans[, 1]

    cat("   Criterion Cb (differential metabolic analysis) ...","\n")

    sink(file=paste("OUTPUT-METNOR-Record",".txt",sep=""))
    # 3. Consistency -------------------------------------------------------- #

    test_data <- eva_data3
    test_data[sapply(test_data, simplify = 'matrix',is.nan)] <- 0

    test_data <- test_data[order(test_data$Group),]

    number_labels <- test_data$Group
    folds <- 3
    test.fold <- list()

    test.fold1 <- strata(c("Group"),size=(as.numeric(table(test_data$Group))/3),method="srswor",data=test_data)[,2]
    test.fold[[1]] <- test.fold1

    data.2 <- test_data[-test.fold1,]

    test.fold2 <- strata(c("Group"),size=(as.numeric(table(test_data$Group))/3),method="srswor",data=data.2)[,2]
    test.fold2 <- match(row.names(data.2)[test.fold2],row.names(test_data))

    test.fold[[2]] <- test.fold2

    test.fold[[3]] <- (1:nrow(test_data))[-c(test.fold1,test.fold2)]

    DEG.list <- try(consistency_M(3, 10))

    if(class(DEG.list)=="try-error")
    { next }

    CW_value <- try(CWvalue(DEG.list,Y=(ncol(eva_data3)-1),n=length(DEG.list[[1]])))
    if(class(CW_value)=="try-error")
    { next }

    Fscore[names(normal_data[mmm])]<-CW_value

    setlist3 <- DEG.list
    OLlist3 <- try(overLapper(setlist = setlist3, sep="_", type = "vennsets"))
    if(class(OLlist3)=="try-error")
    { next }

    counts <- list(sapply(OLlist3$Venn_List, length), sapply(OLlist3$Venn_List, length))

    sink()

    cat("   Criterion Cc (consistency in marker discovery) ...","\n")
    # -- 4. AUC value --------------------------------------------------------------------------------- #

    set.seed(3)

    X <- X_matrix[, pos_filter]
    y <- as.factor(Group)

    # cross-validated SVM-probability plot
    folds <- 5
    test.fold <- split(sample(1:length(y)), 1:folds) #ignore warning

    all.pred.tables <-  lapply(1:folds, function(i) {
      test <- test.fold[[i]]
      Xtrain <- X[-test, ]
      ytrain <- as.factor(y[-test])

      sm <- try(svm(Xtrain, ytrain, cost = as.numeric(1), prob = TRUE)) # some tuning may be needed

      prob.benign <- try(attr(predict(sm, X[test,], prob = TRUE), "probabilities")[, 2])

      data.frame(ytest = y[test], ypred = prob.benign) # returning this
    })

    full.pred.table <- try(do.call(rbind, all.pred.tables))
    if(class(full.pred.table)=="try-error")
    { next }
    svm_para3 <- c(1, 5, as.numeric(1))
    roc_data3 <- full.pred.table
    auc.value <- try(auc(roc(full.pred.table[, 2], full.pred.table[, 1])))
    if(class(auc.value)=="try-error")
    { next }

    Fauc[names(normal_data[mmm])] <- auc.value

    cat("   Criterion Cd (classification accuracy) ...","\n")
    cat("\n")

  }

  result<-dplyr::bind_rows("Precision"=unlist(Fpmad),"Cluster_accuracy"=unlist(Fpurity),"Reproducibility"=unlist(Fscore),"Classification"=unlist(Fauc),.id = "id")
  result1<-t(result)
  colnames(result1)<-result1["id",]
  result2<-result1[-1,]
  result2<-data.frame(result2,check.names=FALSE)

  for(i in 1:dim(result2)[2]){result2[,i]=as.numeric(as.character(result2[,i]))}

  Rank<-apply(result2, 2, function(x){rank(-x,ties.method="min",na.last = "keep")})

  if(length(grep("Precision",colnames(Rank)))==1){
    Rank[,"Precision"]<-rank(as.numeric(as.character(result2[,"Precision"])),ties.method="min",na.last = "keep")
  }else{
    Rank<-Rank
  }

  Rank_revision<-apply(Rank, 2, function(x){x[is.na(x)]<-nrow(Rank);return(x)})
  Ranksum0<-apply(Rank_revision, 1, sum)
  Rankres0<-cbind("OverallRank"=Ranksum0,Rank_revision)

  zuihou0<-cbind("Rank"=Rankres0,"Value"=result2)

  zuihou1<-zuihou0[order(Rankres0[,"OverallRank"],decreasing = FALSE),]
  zuihou1[,1]<-rank(zuihou1[,1],ties.method="min")
  zuihou2<-zuihou1
  zuihou3 <- zuihou2
  zuihou3 <- round(zuihou3,4)
  colnames(zuihou3) <- c("Overall-Rank","Criteria.Ca-Rank","Criteria.Cb-Rank","Criteria.Cc-Rank","Criteria.Cd-Rank","Criteria.Ca-Value","Criteria.Cb-Value","Criteria.Cc-Value","Criteria.Cd-Value")

  data_color<-as.data.frame(zuihou2[,-c(1:5)])

  data_color["Value.Precision"][data_color["Value.Precision"]>=0.7]<-1
  data_color["Value.Precision"][data_color["Value.Precision"]<0.7&data_color["Value.Precision"]>=0.3]<-8
  data_color["Value.Precision"][data_color["Value.Precision"]<0.3]<-10

  data_color["Value.Cluster_accuracy"][data_color["Value.Cluster_accuracy"]>=0.8]<-10
  data_color["Value.Cluster_accuracy"][data_color["Value.Cluster_accuracy"]<0.8&data_color["Value.Cluster_accuracy"]>=0.5]<-8
  data_color["Value.Cluster_accuracy"][data_color["Value.Cluster_accuracy"]<0.5]<-1

  data_color["Value.Reproducibility"][data_color["Value.Reproducibility"]>=0.3]<-10
  data_color["Value.Reproducibility"][data_color["Value.Reproducibility"]<0.3&data_color["Value.Reproducibility"]>=0.15]<-8
  data_color["Value.Reproducibility"][data_color["Value.Reproducibility"]<0.15]<-1

  data_color["Value.Classification"][data_color["Value.Classification"]>=0.9]<-10
  data_color["Value.Classification"][data_color["Value.Classification"]<0.9&data_color["Value.Classification"]>=0.7]<-8
  data_color["Value.Classification"][data_color["Value.Classification"]<0.7]<-1

  data_color_m<-as.data.frame(data_color)
  Ranksum_color<-apply(data_color_m, 1, sum)
  data_color_m01<-cbind( "rank_color"=Ranksum_color,data_color_m)
  data_color_m02<-data_color_m01[order(data_color_m01[,"rank_color"],decreasing =T),]
  data_color_m<-data_color_m02

  row<-rownames(data_color_m)
  nfina<-nchar(row[1])
  nstart<-nchar(row[1])-6
  result <- substring(row, nstart,nfina)

  data_heat<-data_color_m[,-1]

  rank_result <- zuihou3[match(row.names(data_heat),row.names(zuihou3)),]
  rank_result[,1] <- 1:nrow(rank_result)

  cat("\n")
  cat("*************************************************************************","\n")
  cat("Congratulations! Assessment Successfully Completed!","\n")
  cat("Thanks for Using METNOR. Wish to See You Soon ...","\n")
  cat("*************************************************************************","\n")
  cat("\n")
  return(rank_result)
}