R/metrics.R

In bzhanglab/OmicsEV: A tool for large scale omics datasets evaluation

#
import_network_data=function(file=NULL,type="protein",species="human"){
    if(species=="human"){
        if((type=="gene" || type=="protein") && is.null(file)){
            x <- readRDS(system.file("extdata/corum.rds",package = "OmicsEV"))
        }else if(type=="phosphopeptide"){
            x <- readRDS(system.file("extdata/psp.rds",package = "OmicsEV"))
        }
    }else if(species=="drosophila"){
        if((type=="gene" || type=="protein") && is.null(file)){
            x <- readRDS(system.file("extdata/flybase_ppi.rds",package = "OmicsEV"))
        }
    }else{
        cat("species:",species,"\n")
        stop("The species is not supported!")
        x <- NULL
    }
    return(x)
}

import_network_data_from_table=function(file){
    a <- readr::read_tsv(file)
    x <- a$complex
    n_p <- sapply(x,function(y){
        a=unique(unlist(strsplit(split=";",x=y)));
        length(a)})
    # filter complex only having one protein
    x <- x[n_p>=2]
    x <- sapply(x, function(y){
        a=unique(unlist(strsplit(split=";",x=y)));
        a=sort(a);paste(a,sep=";",collapse = ";")}) %>%
        unique
    aa <- sapply(x, function(y){
        t(combn(unlist(strsplit(split = ";",x=y)),2))})
    names(aa) <- NULL
    aa <- do.call(rbind,aa)
    res <- list()
    res$complex <- apply(aa,1,function(xx){
        sort(str_replace_all(xx,pattern=" ",replacement=""))}) %>%
        t %>%
        tibble::as_tibble() %>%
        distinct() %>%
        filter(V1!="",V2!="")
    res$complex <- res$complex %>% filter(V1!=V2)

    ## random complex
    rnd <- c(res$complex$V1,res$complex$V2) %>%
        unique %>%
        sort %>%
        combn(m=2) %>%
        t %>%
        tibble::as_tibble()

    yyy <- paste(res$complex$V1,res$complex$V2,sep="|")
    rnd_p <- paste(rnd$V1,rnd$V2,sep="|")
    rnd <- rnd[ !(rnd_p %in% yyy), ]
    res$rnd_complex <- rnd

    names(res) <- c("network","rnd_network")
    return(res)
}


import_data_from_corum=function(file){
    a <- read.csv(file,stringsAsFactors = FALSE)
    x <- a$subunits.Gene.name.
    n_p <- sapply(x,function(y){
        a=unique(unlist(strsplit(split=";",x=y)));
        length(a)})
    # filter complex only having one protein
    x <- x[n_p>=2]
    x <- sapply(x, function(y){
        a=unique(unlist(strsplit(split=";",x=y)));
        a=sort(a);paste(a,sep=";",collapse = ";")}) %>%
        unique
    aa <- sapply(x, function(y){
        t(combn(unlist(strsplit(split = ";",x=y)),2))})
    names(aa) <- NULL
    aa <- do.call(rbind,aa)
    res <- list()
    res$complex <- apply(aa,1,function(xx){
        sort(str_replace_all(xx,pattern=" ",replacement=""))}) %>%
        t %>%
        tibble::as_tibble() %>%
        distinct() %>%
        filter(V1!="",V2!="")
    res$complex <- res$complex %>% filter(V1!=V2)

    ## random complex
    rnd <- c(res$complex$V1,res$complex$V2) %>%
        unique %>%
        sort %>%
        combn(m=2) %>%
        t %>%
        tibble::as_tibble()

    yyy <- paste(res$complex$V1,res$complex$V2,sep="|")
    rnd_p <- paste(rnd$V1,rnd$V2,sep="|")
    rnd <- rnd[ !(rnd_p %in% yyy), ]
    res$rnd_complex <- rnd

    names(res) <- c("network","rnd_network")

    return(res)
}

# wget http://www.droidb.org/data/DroID_v2018_08/flybase_ppi.txt
# file = "inst/extdata/flybase_ppi.txt"
# import_data_from_flybase_ppi(file)
import_data_from_flybase_ppi=function(file){

    a <- read.delim(file,stringsAsFactors = FALSE)
    x <- a %>% select(SYMBOL1,SYMBOL2)
    names(x) <- c("V1","V2")
    x <- x %>% filter(V1!=V2) %>% distinct() %>%
        apply(1,function(y) {sort(y)}) %>%
        t %>%
        as.data.frame(stringsAsFactors=FALSE) %>%
        distinct()
    all_genes <- sort(unique(c(x$V1,x$V2)))

    res <- list()
    res$network <- x

    ## random complex
    rnd <- sort(all_genes) %>%
        combn(m=2) %>%
        t %>%
        as.data.frame(stringsAsFactors=FALSE)

    yyy <- paste(res$network$V1,res$network$V2,sep="|")
    rnd_p <- paste(rnd$V1,rnd$V2,sep="|")
    rnd <- rnd[ !(rnd_p %in% yyy), ]
    res$rnd_network <- rnd
    return(res)
}


calc_network_corr=function(x,y,sample_class=NULL,missing_value_ratio=0.00001,cpu=0){

    res <- list()
    res$use_sample_class <- sample_class

    if(!is.null(sample_class)){
        cat("Only use samples from class: ",
            paste(sample_class,collapse = ","),"\n")
        x <- lapply(x, function(xx){
                rm_samples <- xx@sampleList %>% filter(!(class %in% sample_class))
                yy <- removeSample(xx,rsamples = rm_samples$sample)
                cat("Use samples:\n")
                out <- yy@peaksData %>% dplyr::select(sample,class) %>%
                    dplyr::distinct() %>%
                    dplyr::group_by(class) %>%
                    dplyr::summarise(n=n())
                print(out)
                return(yy)
        })
    }else{
        cat("Use all samples!\n")
    }

    res$missing_value_ratio <- missing_value_ratio

    if(missing_value_ratio>0){
        cat("Filter missing value with cutoff: ",missing_value_ratio,"\n")
        x <- lapply(x, function(xx){
            cat("dataset:",xx@ID,"\n")
            yy <- xx %>% filterPeaks(ratio=missing_value_ratio)
            return(yy)
        })
    }else{
        cat("No missing value filtering!\n")
    }


    for(i in 1:length(x)){
        if(i == 1){
            common_ids <- x[[1]]@peaksData$ID
        }else{
            common_ids <- intersect(common_ids,x[[i]]@peaksData$ID)
        }
    }


    common_ids <- unique(common_ids)
    res$common_ids <- common_ids
    res$n_used_pairs <- nrow(y$network)

    cat("Common genes:",length(common_ids),"\n")
    cat("Used complexes:",nrow(y$network),"\n")
    if(nrow(y$rnd_network) > 10*nrow(y$network)){
        set.seed(1900)
        n_rnd <- nrow(y$rnd_network)
        n <- nrow(y$network)
        ind <- sample(1:n_rnd,2*n)
        y$rnd_network <- y$rnd_network[ind,]

        cat("Total random complexes:",n_rnd,"\n")
        cat("Used random complexes:",nrow(y$rnd_network),"\n")
    }

    if(cpu==0){
        cpu <- detectCores()
    }
    if(cpu > length(x)){
        cpu <- length(x)
    }
    cat("Used cpus:",cpu,"\n")
    cl <- makeCluster(getOption("cl.cores", cpu))
    clusterExport(cl, c("calc_pair_cor"),envir=environment())
    clusterExport(cl, c("calcCor"),envir=environment())
    #clusterEvalQ(cl,library("dplyr"))

    tmp_res <- parLapply(cl,x,fun = calc_pair_cor,pdat=y,idset=common_ids)
    res$result <- tmp_res

    stopCluster(cl)

    return(res)
}


calc_pair_cor=function(para,pdat,idset=NULL){
    ## filter genes
    valid_ids <- para@peaksData$ID %>% unique
    pdat$network <- pdat$network %>%
        filter(V1 %in% valid_ids, V2 %in% valid_ids)
    pdat$rnd_network <- pdat$rnd_network %>%
        filter(V1 %in% valid_ids, V2 %in% valid_ids)

    ## only consider the genes provided by users
    if(!is.null(idset)){
        valid_ids <- idset
        pdat$network <- pdat$network %>%
            filter(V1 %in% valid_ids, V2 %in% valid_ids)
        pdat$rnd_network <- pdat$rnd_network %>%
            filter(V1 %in% valid_ids, V2 %in% valid_ids)
    }

    pdat$network <- pdat$network %>% mutate(nid=1:nrow(pdat$network))
    dat <- pdat$network %>% gather(key="network",value="ID",-nid)
    m <- inner_join(dat,para@peaksData,by="ID")
    res <- list()
    res$network <- m %>% group_by(nid) %>%
        do(cor=calcCor(.)) %>%
        ungroup() %>%
        as.data.frame() %>%
        mutate(cor=unlist(cor))

    ## random complex
    pdat$rnd_network <- pdat$rnd_network %>% mutate(nid=1:nrow(pdat$rnd_network))
    dat <- pdat$rnd_network %>% gather(key="network",value="ID",-nid)
    m <- inner_join(dat,para@peaksData,by="ID")
    res$rnd_network <- m %>% group_by(nid) %>%
        do(cor=calcCor(.)) %>%
        ungroup() %>%
        as.data.frame() %>%
        mutate(cor=unlist(cor))
    return(res)
}

calcCor=function(x){
    x <- x %>% dplyr::select(ID,sample,value,nid) %>%
        spread(key=ID,value=value) %>%
        as.data.frame()
    res <- cor(x[,3],x[,4],use = "com",method = "spearman")
    return(res)
}


plot_network_cor=function(x, out_dir="./",prefix="test"){

    res <- list()

    dataset_name <- names(x)
    for(i in 1:length(dataset_name)){
        x[[i]]$network$nid <- dataset_name[i]
        x[[i]]$rnd_network$nid <- dataset_name[i]
    }

    ks_res <- x %>% lapply(function(x){
        x1 <- unlist(x$network$cor)
        x2 <- unlist(x$rnd_network$cor)
        # rr <- ks.test(x1,x2)
        # rr$statistic
        roc_input <- dplyr::bind_rows(data.frame(class="network",cor=x1,stringsAsFactors = FALSE),
                                      data.frame(class="rnd_network",cor=x2,stringsAsFactors = FALSE))
        roc.obj <- pROC::roc(roc_input$class,roc_input$cor,percent = FALSE)
        auroc_value <- roc.obj$auc
        return(auroc_value)
    })
    res$ks <- ks_res
    #ks_file <- paste(out_dir,"/",prefix,"-ks.txt",sep="")
    #write(ks_res,file = ks_file)
    #ppi_cor_res$res11$complex$nid <- "RNA"

    network_data_list <- lapply(x, function(y){
        return(y$network)
    })

    network_data_df <- bind_rows(network_data_list) %>% mutate(type="IntraComplex")

    rnd_network_data_list <- lapply(x, function(y){
        return(y$rnd_network)
    })

    rnd_network_data_df <- bind_rows(rnd_network_data_list) %>% mutate(type="InterComplex")


    dat <- rbind(network_data_df,rnd_network_data_df)
    save(dat,file = "dat.rda")
    fig <- paste(out_dir,"/",prefix,"-complex_all_boxplot.png",sep="")
    res$network_boxplot <- fig
    png(fig,width = 1000,height = 500,res=160)
    aaaa <- dat %>% group_by(nid) %>%
        filter(type=="InterComplex") %>%
        summarize(c=median(cor))
        #arrange(c)
    dat -> xxxx
    dat$nid <- factor(dat$nid,levels = aaaa$nid )
    gg <- ggplot(dat,aes(x=nid,y=cor,colour=type))+
        theme_bw()+
        #geom_violin(position="dodge")+
        geom_boxplot(position="dodge")+
        xlab("data table")+
        ylab("Correlation")+
        theme(axis.text.x = element_text(angle = 90, hjust = 1))
    print(gg)
    dev.off()
    save(dat,file="complex_res.rda")

    cor_res <- dat %>% group_by(nid,type) %>%
        summarise(cor=median(cor)) %>%
        ungroup() %>%
        spread(key="type",value="cor") %>%
        mutate(diff=IntraComplex-InterComplex) %>%
        rename(dataSet=nid)

    ks_df <- data.frame(dataSet=names(ks_res),ks=ks_res %>% unlist())
    res$cor <- merge(cor_res,ks_df,sort=FALSE)


    ## save plot data
    fig_data_file <- paste(out_dir,"/",prefix,"-complex_all_boxplot.rds",sep="")
    fig_data <- list()
    fig_data$plot_object <- gg
    fig_data$plot_data <- dat
    saveRDS(fig_data,file = fig_data_file)


    return(res)
}


## protein rna correlation
# "pearson", "kendall", "spearman"
# This function is used to calculate the correlation between protein and mRNA.
calcCorBetweenProteinAndRNA=function(para1,para2,log=TRUE,select_by=1,top_n=1000,
                                     geneset = NULL,
                                     outdir = "./",use_class=NULL,
                                     reg=FALSE,
                                     cor_method="spearman",valueID="value"){

    prefix <- para1@ID

    if(!is.null(use_class)){
        cat("Only use class ",use_class,"\n")
        if(reg == FALSE){
            para1@peaksData <- para1@peaksData %>% dplyr::filter(class==use_class)
            para2@peaksData <- para2@peaksData %>% dplyr::filter(class==use_class)
        }else{
            para1@peaksData <- para1@peaksData %>% dplyr::filter(grepl(pattern = use_class,x = class))
            para2@peaksData <- para2@peaksData %>% dplyr::filter(grepl(pattern = use_class,x = class))
        }
        cat("Used samples:",length(unique(para1@peaksData$sample)),length(unique(para2@peaksData$sample)),"\n")
    }else{
        cat("Use all samples!\n")
    }

    if(!is.null(geneset)){
        cat("Only use specified gene set\n")
        para1@peaksData <- para1@peaksData %>% dplyr::filter(ID %in% geneset)
        para2@peaksData <- para2@peaksData %>% dplyr::filter(ID %in% geneset)
    }

    x1 <- para1@peaksData %>% dplyr::select(ID,sample,!!valueID)
    x2 <- para2@peaksData %>% dplyr::select(ID,sample,!!valueID)
    m <- merge(x1,x2,by=c("ID","sample"))

    ## firstly, plot a scatterplot: x axis is sd, y axis is cor
    names(m)[3:4] <- c("x","y")
    m$x[m$x <= 0] <- NA
    m$y[m$y <= 0] <- NA
    if(log==TRUE){
        m$x <- log2(m$x)
        m$y <- log2(m$y)
    }

    ## filter genes with all values are missing
    na_row <- apply(m[,3:4],1,function(x){any(is.na(x))})
    m <- m[!na_row,]

    ## filter genes occurring less than 3 samples
    n3 <- m %>% group_by(ID) %>% summarise(n=sum(!is.na(x) & !is.na(y))) %>% filter(n>=3)
    m <- m %>% filter(ID %in% n3$ID)

    ## Sample wise correlation
    save(m,file = "m.rda")
    sample_wise_cor <- m %>% group_by(sample) %>%
        dplyr::summarise(cor=cor(x,y,use = "com",method = cor_method)) %>%
        dplyr::mutate(label=prefix)

    ## Gene wise correlation
    res <- m %>% group_by(ID) %>% dplyr::summarise(sd_x=sd(x,na.rm = TRUE),
                                                   sd_y=sd(y,na.rm=TRUE),
                                                   cor=cor(x,y,use = "com",method = cor_method),
                                                   pvalue=get_cor_pvalue(x,y,cor_method = cor_method))
    #save(para1,para2,valueID,na_row,m,res,file="test111.rda")

    png(paste(outdir,"/",prefix,"-cor-sd.png",sep=""),width = 600,height = 1200,res=120)
    cor_name <- paste("Cor (",cor_method,")",sep="")
    par(mfrow=c(4,2),mar=c(3,3,2,1),mgp=c(1.6,0.6,0))
    plot(res$sd_x,res$cor, xlab="SD (1)",ylab=cor_name,col=rgb(0.5,0.4,0.7,0.4),pch=15,cex=0.6)
    plot(res$sd_y,res$cor, xlab="SD (2)",ylab=cor_name,col=rgb(0.5,0.4,0.7,0.4),pch=15,cex=0.6)
    plot(res$sd_x,res$sd_y, xlab="SD (1)",ylab="SD (2)",col=rgb(0.5,0.4,0.7,0.4),pch=15,cex=0.6)
    hist(res$sd_x,nclass=50,xlab="SD (1)",main="")
    hist(res$sd_y,nclass=50,xlab="SD (2)",main="")
    cor_mean <- sprintf("%.4f",mean(res$cor,na.rm=TRUE))
    cor_md <- sprintf("%.4f",median(res$cor,na.rm=TRUE))
    hist(res$cor,nclass=50,xlab=cor_name,main=paste("Data point: ",nrow(res),"\nmean = ",cor_mean,", median = ",cor_md,sep = ""))
    dev.off()

    rr <- data.frame(n=nrow(res),
                     n5=sum(res$cor>=0.5,na.rm = TRUE),
                     n6=sum(res$cor>=0.6,na.rm = TRUE),
                     n7=sum(res$cor>=0.7,na.rm = TRUE),
                     n8=sum(res$cor>=0.8,na.rm = TRUE),
                     median_cor=median(res$cor,na.rm=TRUE),
                     stringsAsFactors = FALSE)

    fres <- list()
    fres$feature_wise_cor_data <- res
    fres$feature_wise <- rr
    fres$sample_wise <- sample_wise_cor

    fres$hist_fig <- plot_mRNA_protein_hist(fres$feature_wise_cor_data,out_dir = outdir,prefix = prefix)

    return(fres)
}

get_cor_pvalue=function(x,y,cor_method="spearman"){
    cor_res <- cor.test(x,y,method = cor_method,use = "com")
    return(cor_res$p.value)
}

plot_mRNA_protein_hist=function(x,out_dir="./",prefix="test"){
    br <- c(seq(min(x$cor,na.rm = TRUE)-0.1,-0.000000001,by=0.01),
            seq(0,max(x$cor,na.rm = TRUE),by=0.01))
    max_count <- max(table(cut(x$cor,breaks = br)))
    positive_cor_ratio <- sum(x$cor>0,na.rm = TRUE)/nrow(x)

    x$pvalue <- p.adjust(x$pvalue,method = "BH")
    sig_positive_cor_ratio <- sum(x$pvalue <= 0.01,na.rm = TRUE)/nrow(x)

    mean_cor <- mean(x$cor,na.rm = TRUE)
    n_pairs <- nrow(x)

    positive_cor_ratio <- sprintf("%.2f%%",100*positive_cor_ratio)
    sig_positive_cor_ratio <- sprintf("%.2f%%",100*sig_positive_cor_ratio) #%>% paste("*\\'%\\'~",sep = "")

    gg <- x %>% mutate(col=ifelse(cor>0,"red","blue")) %>%
        ggplot(aes(x=cor,fill=col)) +
        geom_histogram(breaks=br,color="black",size=0.1)+
        xlab("Spearman’s correlation")+
        ylab("Frequency")+
        ggpubr::theme_pubr()+
        theme(legend.position = "none")+
        scale_fill_manual(breaks = c("red", "blue"),
                          values=c("red", "blue"))+
        geom_vline(xintercept = mean_cor,linetype=2)+
        annotate("text",x=mean_cor,y=1.05*max_count,label=paste("Mean = ",sprintf("%.2f",mean_cor),sep = ""),
                 hjust=-0.05,size=3.5)+
        annotate("text",x=min(br)+0.01,y=1.05*max_count,
                 label=paste0(n_pairs," pairs\n",
                             positive_cor_ratio," positive correlation\n",
                             sig_positive_cor_ratio," significant\npositive correlation\n(adjusted P <= 0.01)"),
                 vjust=1,hjust=0,size=3.5)
    fig <- paste(out_dir,"/",prefix,"-mRNA-protein-hist.png",sep = "")
    png(fig,width = 650,height = 500,res=150)
    print(gg)
    dev.off()
    gg_obj <- paste(out_dir,"/",prefix,"-mRNA-protein-hist_ggplot.rds",sep = "")
    saveRDS(gg,file = gg_obj)
    return(fig)
}


calc_protein_rna_corr=function(x,rna,sample_class=NULL,out_dir="./",cpu=0,
                               missing_value_ratio=0.5){

    x <- lapply(x, function(xx){
        yy <- metaX::filterPeaks(xx,ratio = missing_value_ratio)
        return(yy)
    })

    if(!is.null(sample_class)){
        cat("Only use samples from class: ",
            paste(sample_class,collapse = ","),"\n")
        x <- lapply(x, function(xx){
            rm_samples <- xx@sampleList %>% filter(!(class %in% sample_class))
            yy <- removeSample(xx,rsamples = rm_samples$sample)
            cat("The number of samples used in correlation analysis:",yy@peaksData$sample %>% unique %>% length,"\n")
            return(yy)
        })
    }else{
        cat("Use all samples!\n")
    }

    for(i in 1:length(x)){
        if(i == 1){
            common_ids <- x[[1]]@peaksData$ID
        }else{
            common_ids <- intersect(common_ids,x[[i]]@peaksData$ID)
        }
    }

    common_ids <- intersect(common_ids,rna@peaksData$ID)

    common_ids <- unique(common_ids)
    cat("Common IDs:", length(common_ids),"\n")

    if(cpu==0){
        cpu <- detectCores()
    }
    if(cpu > length(x)){
        cpu <- length(x)
    }
    cat("Use cpus:",cpu,"\n")
    cl <- makeCluster(getOption("cl.cores", cpu))

    res <- parLapply(cl,x,fun = calcCorBetweenProteinAndRNA,
                     para2=rna,
                     log=TRUE,
                     select_by=1,
                     top_n=1000,
                     cor_method="spearman",
                     valueID="value")

    stopCluster(cl)

    dataset_name <- names(res)
    for(i in 1:length(dataset_name)){
        res[[i]]$feature_wise$dataSet <- dataset_name[i]
        res[[i]]$sample_wise$dataSet <- dataset_name[i]
        res[[i]]$feature_wise_cor_data$dataSet <- dataset_name[i]
    }

    ## feature wise
    cor_res <- lapply(res, function(y){y$feature_wise}) %>% bind_rows()
    cor_res <- bind_cols(cor_res %>% select(dataSet),cor_res %>% select(-dataSet))


    ## sample wise
    sample_wise_cor_res <- lapply(res, function(y){y$sample_wise}) %>% bind_rows()
    sample_wise_cor_res_table <- sample_wise_cor_res %>% group_by(dataSet) %>%
        summarise(median_cor=median(cor,na.rm = TRUE))

    fres <- list(feature_wise_cor_table = cor_res,
                 sample_wise_cor_table = sample_wise_cor_res_table,
                 data=res,
                 use_sample_class = sample_class)

    ## plot
    gg <- ggplot(sample_wise_cor_res,aes(x=dataSet,y=cor,color=dataSet)) +
        geom_violin()+
        geom_boxplot(width=0.1)+
        geom_text(data = sample_wise_cor_res_table, aes(x = dataSet, y = median_cor, label = sprintf("%.4f",median_cor)),
                  size = 5, vjust = -0.2, color = "black",)+
        xlab("data table")+
        ylab("Sample wise correlation")+
        theme(legend.position="none")

    sample_wise_cor_fig <- paste(out_dir,"/sample_wise_cor.png",sep="")
    png(sample_wise_cor_fig,width = 800,height = 400,res=120)
    print(gg)
    dev.off()
    fres$sample_wise_cor_fig <- sample_wise_cor_fig

    # generate CDF and boxplot plots for gene-wise correlation analysis
    ## boxplot
    gene_wise_cor_res_all <- lapply(res, function(y){y$feature_wise_cor_data}) %>% bind_rows()
    fres$gene_wise_cor_boxplot_fig <- plot_feature_wise_cor_boxplot(gene_wise_cor_res_all,out_dir = out_dir,prefix = "OmicsEV")
    ## CDF
    fres$gene_wise_cor_cdf_fig <- plot_feature_wise_cor_cdf(gene_wise_cor_res_all,out_dir = out_dir,prefix = "OmicsEV")

    return(fres)

}

plot_feature_wise_cor_cdf=function(x, out_dir="./",prefix="test"){

    fig <- paste(out_dir,"/",prefix,"-feature_wise_cor_cdf.png",sep="")
    png(fig,width = 550,height = 550,res=150)
    gg <- ggplot(x,aes(cor,group=dataSet,color=dataSet))+
        stat_ecdf(geom = "step",size=0.3)+
        theme(legend.position = c(0, 1),
              legend.justification = c(0, 1),
              legend.key = element_blank(),
              legend.background=element_blank())+
        ylab("CDF")+
        xlab("Correlation")
    print(gg)
    dev.off()

    ## save plot data
    fig_data_file <- paste(out_dir,"/",prefix,"-feature_wise_cor_cdf.rds",sep="")
    fig_data <- list()
    fig_data$plot_object <- gg
    fig_data$plot_data <- x
    saveRDS(fig_data,file = fig_data_file)

    return(fig)
}

plot_feature_wise_cor_boxplot=function(x, out_dir="./",prefix="test"){

    median_cor <- x %>% group_by(dataSet) %>% summarise(median_cor=median(cor,na.rm = TRUE)) %>% arrange(desc(median_cor))
    x$dataSet <- factor(x$dataSet,levels = median_cor$dataSet)

    fig <- paste(out_dir,"/",prefix,"-feature_wise_cor_boxplot.png",sep="")
    png(fig,width = 800,height = 400,res=120)
    gg <- ggplot(x,aes(x=dataSet,y=cor))+
        geom_boxplot(width=0.5,outlier.size=0.2)+
        ylab("Correlation")+
        theme(axis.text.x = element_text(angle = 90, hjust = 1))
    print(gg)
    dev.off()
    return(fig)

}


calc_basic_metrics=function(x,class_color=NULL,out_dir="./",cpu=0){

    if(!is.null(class_color)){
        class_col <- read_tsv(class_color) %>% filter(!is.na(col))
    }else{
        class_col <- NULL
    }

    dataset_name <- names(x)
    xx <- list()
    for(i in 1:length(dataset_name)){
        xx[[i]] <- list(name=dataset_name[i],data=x[[i]])
    }
    names(xx) <- dataset_name

    ## output
    fres <- list()

    plist <- list()
    plist$scale <- "pareto"
    plist$t <- 3
    plist$missValueRatioQC <- 0.5
    plist$missValueRatioSample <- 0.5
    plist$classCol <- class_col
    plist$pcaLabel <- NA
    plist$center <- TRUE

    plist$pngWidth <- 6.5
    plist$pngHeight <- 6.5
    plist$pdfWidth <- 6
    plist$pdfHeight <- 6
    plist$legendRowBatch <- 4
    plist$legendRowClass <- 4

    if(length(dataset_name) >=2){
        fig <- plot_upset(x,out_dir = out_dir, prefix = "total_identification_upset")
        fres$datasets_id_overlap <- fig
    }
    density_fig <- plot_density(x,out_dir = out_dir, prefix = "sample_wise_density")


    if(cpu==0){
        cpu <- detectCores()
    }
    ncpu = cpu
    if(cpu > length(dataset_name)){
        cpu <- length(dataset_name)
    }
    cat("Use cpu:",cpu,"\n")

    if(cpu > 1){
        cl <- makeCluster(getOption("cl.cores", cpu))
        clusterExport(cl, c("metaXpipe"),envir=environment())


        res <- parLapply(cl,xx,fun = run_basic_metrics,plist=plist,out_dir=out_dir)
        stopCluster(cl)
    }else{
        res <- lapply(xx,run_basic_metrics,plist=plist,out_dir=out_dir)
    }


    fres$parameters <- plist
    fres$density_plot <- density_fig
    fres$datasets <- res

    ## distribution
    data_dv <- calc_metrics_for_data_distribution_roc(x,cpu=ncpu)
    fres$quant_median_metric <- data_dv

    return(fres)
}

min_max_scale = function(x){
    return((x- min(x)) /(max(x)-min(x)))
}


plot_upset=function(x,out_dir="./",prefix="test"){

    a <- lapply(x, function(y){
        return(y@peaksData$ID %>% unique())
    })

    fig <- paste(out_dir,"./",prefix,".png",sep = "")
    png(fig, width = 9.5, height = 5.5,res=120, units = "in")

    m  <- ComplexHeatmap::make_comb_mat(a)
    ss = ComplexHeatmap::set_size(m)
    cs = ComplexHeatmap::comb_size(m)
    ht <- ComplexHeatmap::UpSet(m,
               set_order = order(ss),
               comb_order = order(ComplexHeatmap::comb_degree(m), -cs),
               top_annotation = HeatmapAnnotation(
                   "Intersection size" = anno_barplot(cs,
                                                        ylim = c(0, max(cs)*1.1),
                                                        border = FALSE,
                                                        gp = grid::gpar(fill = "black"),
                                                        height = grid::unit(4, "cm")
                   ),
                   annotation_name_side = "left",
                   annotation_name_rot = 90))


    ht <- ComplexHeatmap::draw(ht)
    od <- ComplexHeatmap::column_order(ht)
    ComplexHeatmap::decorate_annotation("Intersection size", {
        grid::grid.text(cs[od], x = seq_along(cs), y = grid::unit(cs[od], "native") + grid::unit(3, "pt"),
                  default.units = "native", just = "bottom", gp = grid::gpar(fontsize = 9))
    })

    dev.off()
    return(fig)
}




run_basic_metrics=function(x,plist,out_dir="./"){

    res <- list()

    para <- x$data
    dat_name <- x$name
    c_out_dir <- paste(out_dir,"/",dat_name,sep="")
    para@outdir <- c_out_dir
    para@prefix <- dat_name

    if(is.null(para@sampleList) || is.na(para@sampleList) ||
       nrow(para@sampleList) ==0){
        sampleList  <- read.delim(para@sampleListFile,stringsAsFactors = FALSE)
    }else{
        sampleList  <- para@sampleList
    }

    checkSampleList(sampleList)

    if(is.null(para@ratioPairs)){
        stop("Please set the value of ratioPairs!")
    }

    ##
    plsdaPara <- new("plsDAPara")
    plsdaPara@scale <- plist$scale
    plsdaPara@t <- plist$t
    plsdaPara@do <- FALSE
    missValueRatioQC <- plist$missValueRatioQC
    missValueRatioSample <- plist$missValueRatioSample
    classCol <- plist$classCol
    pcaLabel <- plist$pcaLabel
    center <- plist$center

    ## directory structure, data (all figures and other files), report.html
    raw_outdir <- para@outdir
    ## all figures and other files
    para@outdir <- paste(para@outdir,"/data",sep="")
    makeDirectory(para)

    if(is.null(para@sampleList) || is.na(para@sampleList) ||
       nrow(para@sampleList) ==0){
        expDesign  <- read.delim(para@sampleListFile,stringsAsFactors = FALSE)
    }else{
        expDesign  <- para@sampleList
    }

    expDesign$class <- as.character(expDesign$class)
    expDesign$class[is.na(expDesign$class)] <- "QC"

    expClass <- expDesign %>% group_by(class) %>% summarise(n=length(sample))

    expBatch <- expDesign %>% group_by(batch) %>% summarise(n=length(sample))

    message(date(),"\treset the peaksData...")
    para <- reSetPeaksData(para = para)

    ##
    res$total_features <- para@peaksData$ID %>% unique() %>% length()

    ## 2. pre-processing, filter noise peaks
    if(hasQC(para)){
        message(date(),"\tfilter peaks in QC sample:>=",missValueRatioQC)
        pre_para <- para
        para <- filterQCPeaks(para = para,ratio = missValueRatioQC)
        rpeak <- length(unique(pre_para@peaksData$ID)) - length(unique(para@peaksData$ID))

        paste("Remove peaks in QC samples with missing
            value greater than ",100*missValueRatioQC, " percent:",rpeak,".",sep="")
        rm(pre_para)

    }
    message(date(),"\tfilter peaks in non-QC sample:>=",missValueRatioSample)
    pre_para <- para
    para <- filterPeaks(para = para,ratio = missValueRatioSample)
    rpeak <- length(unique(pre_para@peaksData$ID)) - length(unique(para@peaksData$ID))

    ##
    res$valid_features <- para@peaksData$ID %>% unique() %>% length()

    paste("Remove peaks in non-QC samples with missing
            value greater than ",100*missValueRatioSample," percent:",rpeak,".")
    rm(pre_para)

    message("The CV distribution before normalization: ")
    printCV(para,valueID = "value")

    ## 3. Quality control
    message("plot peak number distribution...")

    fig <- plotPeakNumber(para,legendRowBatch=10)
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")

    ##
    res$features_number_distribution <- fig$fig

    message("plot CV distribution...")
    class_sample <- para@peaksData %>% dplyr::select(class,sample) %>% distinct()
    class_2 <- which(table(class_sample$class)>=2) %>% names()
    cv_para <- para
    cv_para@peaksData <- cv_para@peaksData %>% dplyr::filter(class %in% class_2)
    cv_para@sampleList <- cv_para@sampleList %>% dplyr::filter(class %in% class_2)

    fig <- metaX::plotCV(cv_para)
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")
    ##
    res$cv_distribution <- fig$fig
    res$cv_stat <- get_cv(cv_para@peaksData) %>% mutate(dataSet=cv_para@ID)

    message("plot missing value distribution...")
    fig <- plotMissValue(para,height = 3.7,width = 9.15)
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")

    ##
    res$missing_value_distribution <- fig$fig
    res$non_missing_value_ratio <- data.frame(dataSet=para@ID,
                                              ratio=get_missing_value_ratio(para@peaksData),
                                              stringsAsFactors = FALSE)

    message("plot peak intensity distribution...")
    fig <- plotIntDistr(para,width = 9.15)
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")

    ##
    res$features_quant_distribution <- fig$fig

    message("plot TIC distribution...")
    fig <- plotPeakSumDist(para,valueID = "value",height = 3.7,width = 5.5)
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")

    ticTable <- para@peaksData %>% group_by(sample,batch,class) %>%
        dplyr::summarise(tic=sum(value,na.rm = TRUE)) %>%
        group_by(batch) %>%
        dplyr::summarise(n00=quantile(tic)[1],
                         n25=quantile(tic)[2],
                         n50=quantile(tic)[3],
                         n75=quantile(tic)[4],
                         n100=quantile(tic)[5])
    print(ticTable)

    message("plot average intensity distribution...")
    fig <- plotPeakMeanDist(para,valueID = "value")
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")

    meanTable <- para@peaksData %>% group_by(sample,batch,class) %>%
        dplyr::summarise(meanIntensity=mean(value,na.rm = TRUE)) %>%
        group_by(batch) %>%
        dplyr::summarise(n00=quantile(meanIntensity)[1],
                         n25=quantile(meanIntensity)[2],
                         n50=quantile(meanIntensity)[3],
                         n75=quantile(meanIntensity)[4],
                         n100=quantile(meanIntensity)[5])
    print(meanTable)

    message(date(),"\tmissing value inputation...")
    para <- missingValueImpute(para,cpu=1)
    paste("The missing value were imputed by ",para@missValueImputeMethod,".",sep="")


    save(para,file = "para.rds")
    fig <- plotCorHeatmap(para = para,valueID = "value",anno = TRUE,
                          samples = NULL,
                          height = 6,width = 6,sortBy="batch",cluster = FALSE,
                          classCol=classCol)

    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")
    ##
    res$sample_wise_cor_heatmap <- fig$fig

    if(hasQC(para)){
        message("plot correlation heatmap...")
        #saveRDS(para,file = "para.rda")
        prefix_w <- para@prefix
        para@prefix <- paste(prefix_w,"-qc",sep = "")
        fig <- plotCorHeatmap(para = para,valueID = "value",anno = TRUE,
                              samples = NA, ## QC sample
                              height = 6,width = 6,cluster = FALSE,
                              classCol=classCol)
        #fig1 <- paste("data/",basename(fig$fig),sep="")
        #fig2 <- paste("data/",basename(fig$highfig),sep="")
        para@prefix <- prefix_w
        res$qc_sample_wise_cor_heatmap <- fig$fig

    }

    #save(para,t,file="test.rda")
    save(para,plist,scale,center,pcaLabel,classCol,file=paste(para@outdir,"/","pca.rda",sep=""))
    ppca <- transformation(para,method = plist$t,valueID = "value")
    ppca <- metaX::preProcess(ppca,scale = plist$scale,center = plist$center,
                              valueID = "value")
    fig <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = TRUE,
                          rmQC = FALSE,label = pcaLabel,classColor = classCol,
                          legendRowClass = plist$legendRowClass,
                          legendRowBatch = plist$legendRowBatch,
                          pngWidth = plist$pngWidth,
                          pngHeight = plist$pngHeight,
                          pdfWidth = plist$pdfWidth,
                          pdfHeight = plist$pdfHeight)
    prefix_bak <- ppca@prefix
    ppca@prefix <- paste(ppca@prefix,"-nobatch",sep="")
    fig_nobatch <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = FALSE,
                                  rmQC = FALSE,label = pcaLabel,classColor = classCol,
                                  legendRowClass = plist$legendRowClass,
                                  legendRowBatch = plist$legendRowBatch,
                                  pngWidth = plist$pngWidth,
                                  pngHeight = plist$pngHeight,
                                  pdfWidth = plist$pdfWidth,
                                  pdfHeight = plist$pdfHeight)

    ##
    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")
    res$pca_with_batch <- fig$fig
    res$pca_plot_data <- fig$plotdata


    ## plot heatmap
    save(ppca,classCol,file = "heatmap.rda")
    fig <- plotHeatMap(ppca,valueID="value",log=FALSE,rmQC=FALSE,
                        #scale="row",
                        #clustering_distance_rows="euclidean",
                        #clustering_distance_cols="euclidean",
                        #clustering_method="ward.D2",
                        seriation = FALSE,
                        classCol=classCol,
                        show_colnames=FALSE)

    ppca@prefix <- prefix_bak


    #fig1 <- paste("data/",basename(fig$fig),sep="")
    #fig2 <- paste("data/",basename(fig$highfig),sep="")
    res$cluster_heatmap <- fig$fig


    ## remove QC
    ppca@prefix <- paste(ppca@prefix,"-noqc",sep="")
    fig <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = TRUE,
                          rmQC = TRUE,label = pcaLabel,classColor = classCol,
                          legendRowClass = plist$legendRowClass,
                          legendRowBatch = plist$legendRowBatch,
                          pngWidth = plist$pngWidth,
                          pngHeight = plist$pngHeight,
                          pdfWidth = plist$pdfWidth,
                          pdfHeight = plist$pdfHeight)
    ## no QC, no batch
    ppca@prefix <- paste(ppca@prefix,"-nobatch",sep="")
    fig <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = FALSE,
                          rmQC = TRUE,label = pcaLabel,classColor = classCol,
                          legendRowClass = plist$legendRowClass,
                          legendRowBatch = plist$legendRowBatch,
                          pngWidth = plist$pngWidth,
                          pngHeight = plist$pngHeight,
                          pdfWidth = plist$pdfWidth,
                          pdfHeight = plist$pdfHeight)
    #plotPLSDA(para)

    if(!is.null(para@ratioPairs) && para@ratioPairs != ""){
        ## only perform this analysis when a user provides comparison group
        ## information.
        para@peaksData$valueNorm <- para@peaksData$value
        #para <- peakStat(para = para,plsdaPara = plsdaPara,doROC=FALSE,
        #                 pcaLabel=pcaLabel,classColor = pcaColor)

        qf <- paste("data/",para@prefix,"-quant.txt",sep="")
    }

    saveRDS(para,file = paste(para@outdir,"/",para@prefix,"-result.rds",
                                  sep=""))

    message("Print information about the current R session:")
    writeLines(capture.output(sessionInfo()),
               paste(para@outdir,"/",para@prefix,"-sessionInfo.txt",sep=""))

    return(res)
}


get_sample_data=function(x){
    a <- lapply(x, function(y){
        yy <- y@peaksData %>% dplyr::select(sample,class,batch) %>%
            dplyr::distinct() %>%
            dplyr::mutate(dataSet=y@ID)
        yy$sample <- as.character(yy$sample)
        yy$class <- as.character(yy$class)
        return(yy)
    })
    xx <- dplyr::bind_rows(a)
    class_samples <- xx %>% group_by(class,dataSet) %>%
        dplyr::summarise(samples=dplyr::n()) %>%
        ungroup() %>%
        tidyr::spread(key=dataSet,value=samples)

    return(class_samples)
}

get_identification_summary_table=function(x,format=TRUE){
    x1 <- get_metrics(x$basic_metrics$datasets,metric="total_features")
    x2 <- get_metrics(x$basic_metrics$datasets,metric="valid_features")
    a <- data.frame(dataSet=names(x1),x1=x1,x2=x2,stringsAsFactors = FALSE)
    if(format==TRUE){
        a <- a %>% mutate(x1=cell_spec(x1, bold = ifelse(x1 >= max(x1), TRUE, FALSE), color = ifelse(x1 >= max(x1), "red", "black"))) %>%
            mutate(x2=cell_spec(x2, bold = ifelse(x2 >= max(x2), TRUE, FALSE), color = ifelse(x2 >= max(x2), "red", "black")))
            # dplyr::mutate(x1 = formattable::color_bar("lightgreen")(x1),
            #              x2 = formattable::color_bar("lightgreen")(x2)) #%>%
            #dplyr::mutate(x11 = ifelse(x1 >= max(x1),
            #                           cell_spec(x1, color = "red", bold = T),
            #                           cell_spec(x1, color = "green", italic = T)),
            #              x22 = ifelse(x2 >= max(x2),
            #                           cell_spec(x2, color = "red", bold = T),
            #                           cell_spec(x2, color = "green", italic = T))) %>%
            #rename(x1=NULL,x2=NULL)


    }
    names(a) <- c("dataSet","#identified features","#quantifiable features")
    rownames(a) <- NULL
    return(a)
}

get_full_path=function(x){
    for(i in 1:length(x)){
        x[i] <- normalizePath(x[i])
    }
    return(x)
}


run_kbet=function(x,out_dir="./",prefix="test"){
    res <- lapply(x, function(y){
        xx <- metaX:::getPeaksTable(y,value="value")
        png(paste(y@ID,".png",sep=""))

        if(min(table(xx$batch)) <= 10){
            dat <- kBET(df=xx[,-c(1:4)],batch = xx$batch,plot = TRUE,
                        heuristic = FALSE,k0 = 3)
        }else{
            dat <- kBET(df=xx[,-c(1:4)],batch = xx$batch,plot = TRUE)
        }
        ff <- list(kbet=dat,name=y@ID)
        dev.off()
        return(ff)
    })

    res_table <- lapply(res, function(a){
        y <- a$kbet$summary[1,] %>% mutate(dataSet=a$name) %>% select(dataSet,everything())
    })
    res_table <- bind_rows(res_table)

    dd <- lapply(res, function(a){
        y <- a$kbet
        dat <- data.frame(value=c(y$stats$kBET.expected,y$stats$kBET.observed),
                          Test=c(rep("Expected(random)",length(y$stats$kBET.expected)),
                                 rep("Observed(kBET)",length(y$stats$kBET.observed))),
                          dataSet=a$name,
                          stringsAsFactors = FALSE)
    })

    gdat <- bind_rows(dd)
    fig <- paste(out_dir,"/",prefix,"-kBET-boxplot.png",sep="")
    png(fig,width = 800,height = 500,res=130)
    gg <- ggplot(gdat,aes(x=Test,y=value,color=dataSet))+
        geom_boxplot()+
        ylab("Rejection rate")
        #facet_grid(.~dataSet)
    print(gg)
    dev.off()

    fres <- list(kBET_boxplot=fig,table=res_table,kbet=res)
    return(fres)
}


calc_silhouette_width=function(x,transform_method="3",scale_method="pareto",
                               center=TRUE,missing_value_ratio=0.5){
    res <- lapply(x, function(y){
        y <- metaX::filterPeaks(y,ratio = missing_value_ratio)
        y@missValueImputeMethod <- "knn"
        y <- metaX::missingValueImpute(y,valueID = "value",method="knn")
        ppca <- metaX::transformation(y,method = transform_method,valueID = "value")
        ppca <- metaX::preProcess(ppca,scale = scale_method,center = center,
                                  valueID = "value")
        xx <- metaX:::getPeaksTable(y,value="value")
        #data: a matrix (rows: samples, columns: features (genes))
        #batch: vector or factor with batch label of each cell
        pca.data <- prcomp(xx[,-c(1:4)], center=FALSE) #compute PCA representation of the data
        batch.silhouette <- batch_sil(pca.data, xx$batch)
        return(batch.silhouette)
    })
    res <- res %>% unlist()
    return(res)
}


calc_pca_batch_regression=function(x,transform_method=3,scale_method="pareto",
                               center=TRUE,missing_value_ratio=0.5,format=FALSE,
                               top_pc=10){
    res <- lapply(x, function(y){
        y <- metaX::filterPeaks(y,ratio = missing_value_ratio)
        y@missValueImputeMethod <- "knn"
        y <- metaX::missingValueImpute(y,valueID = "value",method="knn")
        ppca <- metaX::transformation(y,method = transform_method,valueID = "value")
        ppca <- metaX::preProcess(ppca,scale = scale_method,center = center,
                                  valueID = "value")
        xx <- metaX:::getPeaksTable(ppca,value="value")
        #data: a matrix (rows: samples, columns: features (genes))
        #batch: vector or factor with batch label of each cell
        pca.data <- prcomp(xx[,-c(1:4)], center=FALSE) #compute PCA representation of the data
        yy <- pcRegression(pca.data,batch = xx$batch,n_top = top_pc)
        save(pca.data,yy,xx,top_pc,file = "pca_batch_reg.rda")
        return(list(pcr=yy,name=y@ID))
        return(yy)
    })

    save(res,file="res.rda")

    fres <- lapply(res, function(a){
        r2 <- as.data.frame(a$pcr$r2)
        r2 <- r2[1:min(top_pc,nrow(r2)),]
        #r2$R.squared <- cell_spec(r2$R.squared, bold = T,
        #                          color = ifelse(r2$p.value.lm <= 0.05, "red", "black"),
        #                          font_size = spec_font_size(r2$R.squared))
        #r2$R.squared <- cell_spec(r2$R.squared, )
        r2$PC <- row.names(r2)
        r2$PC <- as.integer(str_replace_all(r2$PC,pattern = "PC",replacement = ""))
        r2$dataSet <- a$name
        r2$ExplainedVar <- a$pcr$ExplainedVar[1:nrow(r2)]
        return(r2)
    })

    fres <- bind_rows(fres)
    fres$R.squared <- sprintf("%.3f",fres$R.square) %>% as.numeric()
    fres$R.squared <- cell_spec(fres$R.squared, bold = T,
                              color = ifelse(fres$p.value.lm <= 0.05, "red", "black"),
                              font_size = spec_font_size(fres$R.squared))
    fres2 <- fres %>% select(PC,R.squared,dataSet) %>% spread(key=dataSet,value=R.squared)

    ## explained_var
    explained_var <- fres %>% select(PC,ExplainedVar,dataSet) %>% spread(key=dataSet,value=ExplainedVar)

    ## pcRegscale
    pcRegscale <- sapply(res, function(a){
        return(a$pcr$pcRegscale)
    })
    pcRegscale <- data.frame(dataSet=names(pcRegscale),
                             pcRegscale=pcRegscale,
                             stringsAsFactors = FALSE)
    row.names(pcRegscale) <- NULL

    return(list(pcr=res,
                table=fres2,
                explained_var=explained_var,
                pcRegscale=pcRegscale,
                top_pc=top_pc))
}

get_pcr_table=function(x,top_pc=10){
    fres <- lapply(x, function(a){
        r2 <- as.data.frame(a$pcr$r2)
        r2 <- r2[1:min(top_pc,nrow(r2)),]
        #r2$R.squared <- cell_spec(r2$R.squared, bold = T,
        #                          color = ifelse(r2$p.value.lm <= 0.05, "red", "black"),
        #                          font_size = spec_font_size(r2$R.squared))
        #r2$R.squared <- cell_spec(r2$R.squared, )
        r2$PC <- row.names(r2)
        r2$dataSet <- a$name
        return(r2)
    })

    fres <- bind_rows(fres)
    return(fres)
}


calc_batch_effect_metrics=function(x,out_dir="./",prefix="test",missing_value_cutoff=0.5){

    cat("Missing value imputation ...\n")
    x <- lapply(x, function(y){
        y <- metaX::filterPeaks(y,ratio = missing_value_cutoff)
        y <- missingValueImpute(y)
        return(y)
    })
    #kbet_res <- run_kbet(x)
    sil_res <- calc_silhouette_width(x)
    pcr_res <- calc_pca_batch_regression(x)
    #res <- list(kbet=kbet_res,sil=sil_res,pcr=pcr_res)
    res <- list(kbet=NULL,sil=sil_res,pcr=pcr_res)
    return(res)
}


plot_pca=function(x,out_dir="./",prefix="test",pointSize=0.8,labelSize=4,
                  legendRowBatch = 10,
                  legendRowClass = NULL,show_class=FALSE,
                  pc="12",ncol_sub_fig=3,fig_res=120, pngWidth=6.5, legend_pos=NULL){
    # pc = "12" or "13"

    plotData <- lapply(x,function(y){y$pca_plot_data}) %>% bind_rows()

    if(pc == "12"){
        ggobj <- ggplot(data = plotData,aes(x=x,y=y,color=batch))
        ggobj <- ggobj + geom_hline(yintercept=0,colour="gray")+
            geom_vline(xintercept=0,colour="gray")+
            #geom_point()+
            #xlab(paste("PC1"," (",sprintf("%.2f%%",100*pca.res@R2[1]),") ",sep=""))+
            #ylab(paste("PC2"," (",sprintf("%.2f%%",100*pca.res@R2[2]),") ",sep=""))+
            xlab("PC1")+
            ylab("PC2")+
            theme_bw()+
            theme(#legend.position = legendPosition,
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank())
    }else{
        ggobj <- ggplot(data = plotData,aes(x=x,y=z,color=batch))
        ggobj <- ggobj + geom_hline(yintercept=0,colour="gray")+
            geom_vline(xintercept=0,colour="gray")+
            #geom_point()+
            #xlab(paste("PC1"," (",sprintf("%.2f%%",100*pca.res@R2[1]),") ",sep=""))+
            #ylab(paste("PC2"," (",sprintf("%.2f%%",100*pca.res@R2[2]),") ",sep=""))+
            xlab("PC1")+
            ylab("PC3")+
            theme_bw()+
            theme(#legend.position = legendPosition,
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank())
    }
    ggobj <- ggobj + stat_ellipse(aes(),geom = "path")
    if(show_class==FALSE){
        ggobj <- ggobj + geom_point(size=pointSize)
    }else{
        ggobj <- ggobj + geom_point(aes(shape=class),size=pointSize)+
            scale_shape_manual(values=1:n_distinct(plotData$class))
    }

    ggobj <- ggobj + facet_wrap(.~dataSet,ncol = ncol_sub_fig,scales="free")

    if(!is.null(legendRowBatch)){
        ggobj <- ggobj + guides(color = guide_legend(nrow = legendRowBatch))

    }

    if(!is.null(legend_pos)){
        ggobj <- ggobj + theme(legend.position=legend_pos)
    }

    #pdf(file = highfig,width = pdfWidth,height = pdfHeight)
    #pngWidth <- 6.5
    ht <- ceiling(length(unique(plotData$dataSet))/ncol_sub_fig) * pngWidth/ncol_sub_fig
    fig <- paste(out_dir,"/",prefix,"-pca_batch.png",sep="")
    png(fig,width = pngWidth,height = ht,res = fig_res,units = "in")

    print(ggobj)
    dev.off()

    return(list(fig=fig,data=plotData))

}


plot_density=function(x,out_dir="./",prefix="test",filter_by_quantile=0){
    a <- lapply(x, function(y){
        dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
            as_tibble()
        dat$value[dat$value<=0] <- NA
        return(dat)
    })

    if(filter_by_quantile>0){
        a <- lapply(a, function(y){
            b <- quantile(y$value,probs = c(filter_by_quantile,1-filter_by_quantile),na.rm = TRUE)
            y <- filter(y,value >= b[1], value <= b[2])
            return(y)
        })
    }

    dat <- bind_rows(a)
    #dat$value[dat$value<=0] <- NA

    gg <- dat %>% ggplot(aes(x=log2(value),color=sample)) +
        #stat_density(geom="line")+
        geom_line(stat="density",alpha=0.4)+
        #geom_density()+
        facet_wrap(.~dataSet,ncol = 3,scales="free")+
        guides(color=FALSE)+
        theme_bw()+
        theme(#legend.position = legendPosition,
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank())

    pngWidth <- 6.5
    ht <- ceiling(length(unique(dat$dataSet))/3) * pngWidth/3
    fig <- paste(out_dir,"/",prefix,"-density.png",sep="")
    png(fig,width = pngWidth,height = ht,res = 120,units = "in")
    print(gg)
    dev.off()
    return(fig)
}


calc_ml_metrics=function(x,sample_list=NULL,sample_class=NULL,use_all=TRUE,
                         missing_value_cutoff=0.5, valueID="value",cpu=0,
                         n_repeats=20){

    if(missing_value_cutoff > 0){
        x <- lapply(x, filterPeaks, ratio=missing_value_cutoff)
    }

    if(use_all==FALSE){
        ## use common genes/proteins
        for(i in 1:length(x)){
            if(i == 1){
                common_ids <- x[[1]]@peaksData$ID
            }else{
                common_ids <- intersect(common_ids,x[[i]]@peaksData$ID)
            }
        }
        common_ids <- unique(common_ids)
        cat("Common IDs:", length(common_ids),"\n")

        x <- lapply(x, function(y){
            y@peaksData <- y@peaksData %>% filter(ID %in% common_ids)
            return(y)
        })

    }

    if(!is.null(sample_list)){
        sam <- read.delim(sample_list,stringsAsFactors = FALSE) %>%
            dplyr::select(sample,class)
        x <- lapply(x, function(y){
            y@sampleList <- read.delim(sample_list,stringsAsFactors = FALSE)
            y@peaksData$class<- NULL
            y@peaksData <- merge(y@peaksData,sam,by="sample")
            return(y)
        })

        class_group <- sam$class %>% unique()
    }else{
        if(!is.null(sample_class)){
            class_group <- sample_class
        }else{
            class_group <- x[[1]]@peaksData$class %>% unique()
        }
    }

    cat("Use class:",paste(class_group,collapse = ","),"\n")

    res <- lapply(x, function(y){

        if(requireNamespace("doMC",quietly = TRUE)){
            if(cpu==0){
                library(doMC)
                cpu = detectCores()
                registerDoMC(cores = cpu)
            }else if(cpu>=2){
                library(doMC)
                registerDoMC(cores = cpu)
            }
        }
        if(any(is.na(y@peaksData$value))){
            y <- metaX::missingValueImpute(y)
        }else{
            cat("No missing value!\n")
        }


        dd <- lapply(1:n_repeats,function(i){
            cat("Repeat:", i,"\n")
            dat <- featureSelection(y,group=class_group,method = "rf",
                                    valueID = valueID, fold = 5,
                                    resampling_method = "LOOCV",
                                    repeats = 10,
                                    plot_roc = FALSE,
                                    ratio = 2/3, k = 100,
                                    metric = "ROC",
                                    sizes = length(unique(y@peaksData$ID)),
                                    plotCICurve = FALSE,verbose=TRUE)
            dat$results$Repeat_ID = i
            return(dat)

        })

        eval_res <- lapply(dd,function(x){ x$results %>% select(everything())}) %>%
            bind_rows()

        return(eval_res)
        #return(dat)
    })

    name_datasets <- names(res)

    for(i in 1:length(res)){
        res[[i]]$dataSet <- name_datasets[i]
    }

    res_table <- bind_rows(res)

    ftable <- res_table %>% group_by(dataSet) %>% dplyr::summarise(mean_ROC=mean(ROC),median_ROC=median(ROC),sd_ROC=sd(ROC))

    fres <- list(data=res_table,table=ftable,class_group=class_group)

    return(fres)
}

plot_ml_boxplot=function(x, out_dir="./",prefix="test"){

    mean_roc <- x %>% group_by(dataSet) %>% summarise(mean_ROC=mean(ROC)) %>% arrange(desc(mean_ROC))
    x$dataSet <- factor(x$dataSet,levels = mean_roc$dataSet)



    fig <- paste(out_dir,"/",prefix,"-ml_boxplot.png",sep="")
    png(fig,width = 800,height = 400,res=120)
    gg <- ggplot(x,aes(x=dataSet,y=ROC))+
        geom_boxplot(width=0.5,outlier.size=0.2)+
        theme(axis.text.x = element_text(angle = 90, hjust = 1))+
        stat_summary(fun.y=mean, colour="darkred", geom="point",
                     shape=18, size=3, show.legend = FALSE) +
        geom_text(data = mean_roc %>% mutate(mean_roc_label=sprintf("%.2f",mean_ROC)),
                  aes(label = mean_roc_label, y = mean_ROC + 0.1*(max(x$ROC)-min(x$ROC))),colour="darkred")
    print(gg)
    dev.off()

    ## save plot data
    fig_data_file <- paste(out_dir,"/",prefix,"-ml_boxplot.rds",sep="")
    fig_data <- list()
    fig_data$plot_object <- gg
    fig_data$plot_data <- x
    saveRDS(fig_data,file = fig_data_file)

    return(fig)

}


calc_metrics_for_data_distribution=function(x,cpu=0){
    # x is list of datasets
    if(cpu==0){
        cpu = detectCores()
    }
    a <- lapply(x, function(y){
        # cat("Dataset:",y@ID,"\n")
        dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
            as_tibble()
        dat$value[dat$value<=0] <- NA

        # log2 scale
        dat$value <- log2(dat$value)
        all_samples <- dat$sample %>% unique()
        ss <- combn(all_samples,m=2)

        cl <- makeCluster(getOption("cl.cores", cpu))
        clusterExport(cl, c("get_ks_statistic"),envir=environment())
        ksv <- parLapply(cl,1:ncol(ss),fun = get_ks_statistic,ss=ss,dat=dat)
        stopCluster(cl)
        ksv <- unlist(ksv)
        ks <- median(ksv,na.rm = TRUE)
        res <- data.frame(dataSet=y@ID,quant_median_ks=ks,n=length(ksv),stringsAsFactors = FALSE) %>%
            as_tibble()
        return(res)
    }) %>% bind_rows()
    a$raw_quant_median_ks <- a$quant_median_ks
    #a$quant_median_ks <- 1 - min_max_scale(a$quant_median_ks)
    a$quant_median_ks <- 1 - a$quant_median_ks/(a$quant_median_ks+1)
    return(a)
}

calc_metrics_for_data_distribution_roc=function(x,cpu=0){
    # x is list of datasets
    if(cpu==0){
        cpu = detectCores()
    }
    a <- lapply(x, function(y){
        # cat("Dataset:",y@ID,"\n")
        dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
            as_tibble()
        dat$value[dat$value<=0] <- NA

        # log2 scale
        dat$value <- log2(dat$value)
        all_samples <- dat$sample %>% unique()
        ss <- combn(all_samples,m=2)

        cl <- makeCluster(getOption("cl.cores", cpu))
        clusterExport(cl, c("get_two_sample_auroc"),envir=environment())
        ksv <- parLapply(cl,1:ncol(ss),fun = get_two_sample_auroc,ss=ss,dat=dat)
        stopCluster(cl)
        ksv <- unlist(ksv)
        ks <- median(ksv,na.rm = TRUE)
        res <- data.frame(dataSet=y@ID,quant_median_metric=ks,n=length(ksv),stringsAsFactors = FALSE) %>%
            as_tibble()
        return(res)
    }) %>% bind_rows()
    #a$raw_quant_median_ks <- a$quant_median_ks
    #a$quant_median_ks <- 1 - min_max_scale(a$quant_median_ks)
    #a$quant_median_ks <- 1 - a$quant_median_ks/(a$quant_median_ks+1)
    return(a)
}



get_ks_statistic = function(s,ss,dat){
    ii <- as.character(ss[,s])
    x1 <- dat$value[dat$sample==ii[1]]
    x2 <- dat$value[dat$sample==ii[2]]
    b <- ks.test(x1,x2)
    return(abs(b$statistic))
}

get_two_sample_auroc = function(s,ss,dat){
    ii <- as.character(ss[,s])
    x <- dat %>% dplyr::filter(sample %in% ii[1:2]) %>%
        dplyr::filter(!is.na(value))
    x$sample <- factor(x$sample,levels = ii[1:2])
    roc.obj <- pROC::roc(x$sample,x$value,percent = FALSE)
    auroc_value <- roc.obj$auc
    metric <- 1-2*abs(auroc_value-0.5)
    return(metric)
}

# generate overview table and overview radar figure
generate_overview_table=function(x,highlight_top_n=3,min_auc=0.8){
    ## generate an overview table which contains different metrics
    ## input: x =>

    ## two data tables, one for radar plot real value, one for overview table
    total_features <- x$input_parameters$total_features
    dat <- get_identification_summary_table(x,format = FALSE)
    show_dat <- dat
    if("#identified features" %in% names(dat) && "#quantifiable features" %in% names(dat)){
        dat <- dat %>% dplyr::mutate(`#identified features`=`#identified features`/total_features,
                          `#quantifiable features`=`#quantifiable features`/total_features)
        show_dat <- show_dat %>% dplyr::mutate(`#identified features` = paste(`#identified features`,"\n(",format_number(dat$`#identified features`),")",sep=""),
                                   `#quantifiable features`=paste(`#quantifiable features`,"\n(",format_number(dat$`#quantifiable features`),")",sep=""))
    }


    ## add missing value metric
    missing_value_table <- lapply(x$basic_metrics$datasets,function(y)y[["non_missing_value_ratio"]]) %>%
        dplyr::bind_rows() %>%
        dplyr::rename(non_missing_value_ratio=ratio)
    dat <- merge(dat,missing_value_table %>% dplyr::select(dataSet,non_missing_value_ratio),by="dataSet")
    show_dat <- merge(show_dat,missing_value_table %>%
                          dplyr::select(dataSet,non_missing_value_ratio),by="dataSet") %>%
        dplyr::mutate(non_missing_value_ratio=format_number(non_missing_value_ratio))

    ## data distribution metric
    if(!("quant_median_metric" %in% names(x$basic_metrics)) && length(x$input_parameters$datasets) >= 2){
        ncpu <- ifelse("cpu" %in% names(x$input_parameters),x$input_parameters$cpu,0)
        data_dv <- calc_metrics_for_data_distribution_roc(x$input_parameters$datasets,cpu=ncpu)
        dat <- merge(dat,data_dv %>%
                         dplyr::select(dataSet,quant_median_metric) %>%
                         dplyr::rename(data_dist_similarity=quant_median_metric),by="dataSet")
        show_dat <- merge(show_dat,data_dv %>%
                         dplyr::select(dataSet,quant_median_metric) %>%
                         dplyr::mutate(data_dist_similarity=format_number(quant_median_metric)) %>%
                         dplyr::select(-quant_median_metric),by="dataSet")
    }else if("quant_median_metric" %in% names(x$basic_metrics)){
        data_dv <- x$basic_metrics$quant_median_metric
        dat <- merge(dat,data_dv %>%
                         dplyr::select(dataSet,quant_median_metric) %>%
                         dplyr::rename(data_dist_similarity=quant_median_metric),by="dataSet")
        show_dat <- merge(show_dat,data_dv %>%
                              dplyr::select(dataSet,quant_median_metric) %>%
                              dplyr::mutate(data_dist_similarity=format_number(quant_median_metric)) %>%
                              dplyr::select(-quant_median_metric),by="dataSet")
    }

    if(!is.null(x$batch_effect_metrics)){

        # kBET
        if(!is.null(x$batch_effect_metrics$kbet)){
            dat <- merge(dat,x$batch_effect_metrics$kbet$table %>%
                             select(dataSet,kBET.observed))

            dat$kBET <- 1 - dat$kBET.observed
            dat$kBET.observed <- NULL

            show_dat <- merge(show_dat,x$batch_effect_metrics$kbet$table %>%
                             select(dataSet,kBET.observed))
            show_dat$kBET <- paste(format_number(show_dat$kBET.observed),"\n(",format_number(1-show_dat$kBET.observed),")",sep = "")
            show_dat$kBET.observed <- NULL
        }
        #dat$kBET.observed <- cell_spec(dat$kBET.observed,
        #                               color = ifelse(y >= max(y), "red", "black"))


        # Silhouette width
        sil <- data.frame(dataSet=names(x$batch_effect_metrics$sil),
                          silhouette_width=x$batch_effect_metrics$sil)
        dat <- merge(dat,sil)
        dat$silhouette_width <- 1 - abs(dat$silhouette_width)

        show_dat <- merge(show_dat,sil)
        show_dat$silhouette_width <- paste(format_number(show_dat$silhouette_width),"\n(",format_number(1 - abs(show_dat$silhouette_width)),")",sep = "")


        # PCR
        # pcr <- get_pcr_table(x$batch_effect_metrics$pcr$pcr)
        # sig_pc <- pcr %>% dplyr::filter(p.value.lm<=0.05) %>%
        #     dplyr::select(PC) %>%
        #     dplyr::distinct()
        # pcr <- pcr %>% dplyr::filter(PC %in% sig_pc$PC) %>%
        #     dplyr::select(dataSet,PC,R.squared) %>%
        #     dplyr::mutate(PC=paste("batch_effect_",PC,sep="")) %>%
        #     tidyr::spread(key = PC,value = R.squared)
        pcRegscale <- x$batch_effect_metrics$pcr$pcRegscale
        #pcRegscale$pcRegscale <- 1 - pcRegscale$pcRegscale

        dat <- merge(dat,pcRegscale)
        dat$pcRegscale <- 1 - dat$pcRegscale

        show_dat <- merge(show_dat,pcRegscale)
        show_dat$pcRegscale <- paste(format_number(show_dat$pcRegscale),"\n(",format_number(1 - show_dat$pcRegscale),")",sep = "")
    }

    ## complex
    if(!is.null(x$network_table)){
        dat <- merge(dat,x$network_table$cor %>% dplyr::select(dataSet,ks) %>%
                         dplyr::rename(complex_auc=ks))
        show_dat <- merge(show_dat,x$network_table$cor %>% dplyr::select(dataSet,ks) %>%
                         dplyr::rename(complex_auc=ks)) %>%
            dplyr::mutate(complex_auc=format_number(complex_auc))

    }

    ## function prediction
    if(!is.null(x$fun_pred)){
        f_res <- get_func_pred_meanAUC(x$fun_pred$data,min_auc = min_auc)
        dat <- merge(dat,f_res)
        show_dat <- merge(show_dat,f_res) %>%
            dplyr::mutate(func_auc=format_number(func_auc))
    }

    ## class prediction
    if(!is.null(x$ml)){
        dat <- merge(dat,x$ml$table %>%
                         dplyr::select(dataSet,mean_ROC) %>%
                         dplyr::rename(class_auc=mean_ROC))

        show_dat <- merge(show_dat,x$ml$table %>%
                         dplyr::select(dataSet,mean_ROC) %>%
                         dplyr::rename(class_auc=mean_ROC)) %>%
            dplyr::mutate(class_auc=format_number(class_auc))
    }

    ## median CV: QC samples
    cv_table <- get_cv_table(x$basic_metrics$datasets,"cv_stat")
    if("QC" %in% cv_table$class){
        cv_table <- cv_table %>%
        dplyr::filter(class=="QC") %>%
        dplyr::select(dataSet,cv30) %>%
        dplyr::rename(CV30=cv30) # %>%
        #dplyr::mutate(scaled_median_CV=1-median_CV/max(median_CV))

        dat <- merge(dat,cv_table,by="dataSet") # %>%
            #dplyr::select(-median_CV) %>%
            #dplyr::rename(median_CV=scaled_median_CV)
        show_dat <- merge(show_dat,cv_table,by="dataSet") %>%
            dplyr::mutate(CV30=format_number(CV30)) #%>%
            #dplyr::select(-scaled_median_CV)

    }

    ## mRNA-protein correlation
    if(!is.null(x$input_parameters$x2)){
        dat <- merge(dat,x$protein_rna$feature_wise_cor_table %>%
                         dplyr::select(dataSet,median_cor) %>%
                         dplyr::rename(gene_wise_cor=median_cor))

        show_dat <- merge(show_dat,x$protein_rna$feature_wise_cor_table %>%
                         dplyr::select(dataSet,median_cor) %>%
                         dplyr::rename(gene_wise_cor=median_cor)) %>%
            dplyr::mutate(gene_wise_cor=format_number(gene_wise_cor))

        dat <- merge(dat,x$protein_rna$sample_wise_cor_table %>%
                         dplyr::select(dataSet,median_cor) %>%
                         dplyr::rename(sample_wise_cor=median_cor))

        show_dat <- merge(show_dat,x$protein_rna$sample_wise_cor_table %>%
                         dplyr::select(dataSet,median_cor) %>%
                         dplyr::rename(sample_wise_cor=median_cor)) %>%
            dplyr::mutate(sample_wise_cor=format_number(sample_wise_cor))
    }


    return(list(dat=dat,show_dat=show_dat))

}



plot_radar=function(x){

    x <- x[,c(1:2,ncol(x):3)]
    radar_dat <- x

    p <- plot_ly(
        type = 'scatterpolargl',
        #fill = 'toself',
        #alpha = 0.1,
        #line = list(smoothing = 1, shape = "spline"),
        mode = 'lines+markers'#,
        #marker = list(color = "white",
        #              size = 4,
        #              line = list(width = 2)),

    )
    #
    for(i in 1:nrow(radar_dat)){
        p <- p %>% add_closed_trace(
            r = radar_dat[i,-1] %>% as.numeric(),
            theta = names(radar_dat)[-1],
            #line_close = TRUE,
            name = x$dataSet[i]
            #line = list(color = "#d3d3d3", dash = "3px")
            #mode = "lines"


        )
    }
    p <- p %>%
        layout(
            polar = list(
                radialaxis = list(
                    visible = T,
                    range = c(0,1)
                )
            )
        )

    return(p)
}

add_closed_trace <- function(p, r, theta, ...){
    plotly::add_trace(p, r = c(r, r[1]), theta = c(theta, theta[1]), ...)
}

## Without missing value imputation
noise_signal_analysis=function(x, qc_sample=NULL,bio_sample=NULL, out_dir="./"){

    a <- lapply(x, function(y){
        dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
            as_tibble()
        dat$value[dat$value<=0] <- NA

        dat$value <- log2(dat$value)

        dat <- dat %>% filter(class %in% c(qc_sample, bio_sample))

        dat$class <- ifelse(dat$class %in% qc_sample,"QC","Bio")

        b <- dat %>% dplyr::group_by(ID,class) %>%
            dplyr::summarise(sd=sd(value,na.rm = TRUE),.groups = 'drop') %>%
            tidyr::spread(key = class, value = sd)

        b <- b %>% dplyr::mutate(nsr=QC/Bio, dataSet=y@ID)

        return(b)
    })

    snr_res <- dplyr::bind_rows(a)

    res <- list()
    res$table <- snr_res %>% group_by(dataSet) %>%
        dplyr::summarise(nsr=median(nsr,na.rm=TRUE),n=n()) %>%
        dplyr::arrange(nsr)
    res$raw_data <- snr_res

    snr_res$dataSet <- factor(snr_res$dataSet,levels = res$table$dataSet)

    fig <- paste(out_dir,"/snr.png",sep="")
    png(fig,height = 500,width = 700,res=130)
    gg <- ggplot(snr_res,aes(x=dataSet,y=nsr,color=dataSet)) +
        geom_boxplot(outlier.size = 0.5) +
        ylab("Noise to signal ratio")+
        theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust=0.5),legend.position = "none")
    print(gg)
    dev.off()
    res$fig <- fig


    ## MAD

    mad_figs <- lapply(x, function(y){
        dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
            as_tibble()
        dat$value[dat$value<=0] <- NA

        dat$value <- log2(dat$value)

        b <- dat %>% group_by(sample,class) %>%
            dplyr::summarise(mad=mad(value,na.rm = TRUE)) %>%
            ungroup() %>%
            dplyr::arrange(mad)
        b$i <- 1:nrow(b)

        fig <- paste(out_dir,"/mad",y@ID,".png",sep = "")
        png(fig,width = 1000,height = 400,res=150)
        gg <- ggplot(b,aes(x=i,y=mad,color=class)) +
            geom_point() +
            geom_hline(yintercept = median(b$mad)) +
            xlab("Rank") +
            ylab("MAD (sample)")
        print(gg)
        dev.off()
        return(fig)
    })

    names(mad_figs) <- names(x)

    res$mad_figs <- mad_figs

    return(res)

}


get_cv=function(peaksData){
    cvstat <- peaksData %>%
        dplyr::group_by(class,ID) %>%
        dplyr::summarize(cv=sd(value,na.rm=TRUE)/mean(value,na.rm=TRUE)) %>%
        dplyr::ungroup() %>%
        dplyr::group_by(class) %>%
        dplyr::summarize(median_cv=median(cv,na.rm=TRUE),
                         mean_cv=mean(cv,na.rm=TRUE),
                         n=n(),
                         cv30=sum(cv<=0.3,na.rm = TRUE)/n()) %>%
        as.data.frame()
    return(cvstat)
}

get_cv_table=function(x, metric = "total_features"){
    res <- lapply(x,function(y)y[[metric]]) %>% dplyr::bind_rows()
    return(res)
}

get_missing_value_ratio=function(x){
    n_samples <- length(unique(x$sample))
    n_features <- length(unique(x$ID))
    n_v <- n_samples * n_features
    n_valid_v <- sum(!is.na(x$value))
    ratio <- n_valid_v/n_v
    return(ratio)
}

get_cv_table=function(x, metric = "total_features"){
    res <- lapply(x,function(y)y[[metric]]) %>% dplyr::bind_rows()
    return(res)
}

format_overview_table=function(ov_table){
    a <- ov_table$dat
    show_a <- ov_table$show_dat
    nm <- names(ov_table)
    for(i in 2:ncol(a)){
        y <- a[,i]
        y_show <- show_a[,i]
        y_show <- cell_spec(y_show, bold = ifelse(y >= max(y), TRUE, FALSE),color = ifelse(y >= max(y), "red", "black"))
        show_a[,i] <- y_show
    }
    show_ov_table <- as.data.frame(t(show_a[,-1]))
    names(show_ov_table) <- show_a$dataSet
    show_ov_table$metric <- rownames(show_ov_table)
    rownames(show_ov_table) <- NULL
    show_ov_table <- show_ov_table %>% dplyr::select(metric,everything())
    return(show_ov_table)
}

format_number=function(x,digit=4){
    format_str <- paste("%.",digit,"f",sep = "")
    y <- sprintf(format_str,x)
    return(y)
}

get_formated_id_table=function(final_res){
    id_table <- get_identification_summary_table(final_res,format = FALSE) %>% dplyr::rename(`data table`=dataSet)
    total_features <- final_res$input_parameters$total_features
    show_id_table <- id_table
    if("#identified features" %in% names(id_table) && "#quantifiable features" %in% names(id_table)){
        id_table <- id_table %>% dplyr::mutate(`#identified features`=`#identified features`/total_features,
                                               `#quantifiable features`=`#quantifiable features`/total_features)
        show_id_table <- show_id_table %>% dplyr::mutate(`#identified features` = paste(`#identified features`,"\n(",sprintf("%.2f%%",100*id_table$`#identified features`),")",sep=""),
                                                         `#quantifiable features` = paste(`#quantifiable features`,"\n(",sprintf("%.2f%%",100*id_table$`#quantifiable features`),")",sep=""))
    }
    y <- id_table$`#identified features`
    show_id_table$`#identified features` <- cell_spec(show_id_table$`#identified features`, bold = ifelse(y >= max(y), TRUE, FALSE),color = ifelse(y >= max(y), "red", "black"))
    y <- id_table$`#quantifiable features`
    show_id_table$`#quantifiable features` <- cell_spec(show_id_table$`#quantifiable features`, bold = ifelse(y >= max(y), TRUE, FALSE),color = ifelse(y >= max(y), "red", "black"))

    return(list(id_table=id_table,show_id_table=show_id_table))
}