view Intchecks/Script_intensity_check.R @ 0:c2c2e1be904a draft

Uploaded

####################################################################
# SCRIPT INTENSITY CHECK
# V1: Fold and NA
#
# Input: Data Matrix, VariableMetadata, SampleMetadata
# Output: VariableMetadata, Graphics (barplots and boxplots)
#
# Dependencies: RcheckLibrary.R 
#
####################################################################


# Parameters (for dev)
if(FALSE){
  
  rm(list = ls())
  setwd("Y:\\Developpement\\Intensity check\\Pour tests")
  
  DM.name <- "DM_NA.tabular"
  SM.name <- "SM_NA.tabular"
  VM.name <- "vM_NA.tabular"
  type <- "One_class"
  class.col <- "2"
  class1 <- "Blanks"
  VM.output <- "new_VM.txt"
  graphs.output <- "Barplots_and_Boxplots.pdf"
}

intens_check <- function(DM.name, SM.name, VM.name, type, class.col, class1, VM.output, graphs.output){
  # This function allows to check the intensities considering classes with a fold calculation, the number and
  # the proportion of NA
  #
  # Two options: 
  # - one class (selected by the user) against all the remaining samples ("One_class")
  # - tests on each class ("Each_class")
  #
  # Parameters:
  # DM.name, SM.name, VM.name: dataMatrix, sampleMetadata, variableMetadata files access
  # type: "One_class" or "Each_class"
  # class.col: number of the sampleMetadata's column with classes
  # class1: name of the class if type="One_class"
  # VM.output: output file's access (VM with the new columns)
  # graphs.output: pdf with barplots for the proportion of NA and boxplots with the folds values
  
  
  
  
  # Input ---------------------------------------------------------
  
  DM <- read.table(DM.name, header=TRUE, sep="\t", check.names=FALSE)
  SM <- read.table(SM.name, header=TRUE, sep="\t", check.names=FALSE)
  VM <- read.table(VM.name, header=TRUE, sep="\t", check.names=FALSE)
  
  
  
  # Table match check with Rchecklibrary
  table.check <- match3(DM, SM, VM)
  check.err(table.check)
  
  
  rownames(DM) <- DM[,1]
  var_names <- DM[,1]
  DM <- DM[,-1]
  DM <- data.frame(t(DM))
  
  class.col <- colnames(SM)[as.numeric(class.col)]
  
  # check class.col, class1 and the number of classes---------------
  
  if(!(class.col %in% colnames(SM))){
    stop("\n- - - - - - - - -\n", "The column ",class.col, " is not a part of the specify sample Metadata","\n- - - - - - - - -\n")
  }
  
  c_class <- SM[,class.col]
  c_class <- as.factor(c_class)
  nb_class <- nlevels(c_class)
  classnames <- levels(c_class)
  
  
  if((nb_class > (nrow(SM))/3)){
    class.err <- c("\n There are too many classes, think about reducing the number of classes and excluding those
                   with few samples \n")
    cat(class.err)
  }

  if(type == "One_class"){  
    if(!(class1 %in% classnames)){
      list.class1 <- c("\n Classes:",classnames,"\n")
      cat(list.class1)
      err.class1 <- c("The class ",class1, " does not appear in the column ", class.col)
      stop("\n- - - - - - - - -\n", err.class1,"\n- - - - - - - - -\n")
    }
  }
  
  #If type is "one_class", change others classes in "other"
  if(type == "One_class"){
    for(i in 1:length(c_class)){
      if(c_class[i]!=class1){
        c_class <- as.character(c_class)
        c_class[i] <- "Other"
        c_class <- as.factor(c_class)
        nb_class <- nlevels(c_class)
        classnames <- levels(c_class)
        
      }
    }
  }

  DM <- cbind(DM,c_class)
  
  # fold -------------------------------------------------------
  n <- 1
  fold <- data.frame()
  for(j in 1:(nb_class-1)){
    for(k in (j+1):nb_class) {
      for (i in 1:(length(DM)-1)){
        fold[i,n] <- mean(DM[which(DM$c_class==classnames[k]),i], na.rm=TRUE)/
          mean(DM[which(DM$c_class==classnames[j]),i], na.rm=TRUE)}
      names(fold)[n] <- paste("fold",classnames[k],"VS", classnames[j], sep="_")
      n <- n + 1}
  }
  
  # NA ---------------------------------------------------------
  calcul_NA <- data.frame()
  pct_NA <- data.frame()
  for (i in 1:(length(DM)-1)){
    for (j in 1:nb_class){
      n <- 0
      new_DM <- DM[which(DM$c_class==classnames[j]),i]
      for(k in 1:length(new_DM)){
        if (is.na(new_DM[k])){
          n <- n + 1}
        calcul_NA[i,j] <- n
        pct_NA[i,j] <- (calcul_NA[i,j]/length(new_DM))*100}
    }
  }
  names(calcul_NA) <- paste("Nb_NA",classnames, sep="_")
  names(pct_NA) <- paste("Pct_NA", classnames, sep="_")
  
  # Alert message if there is no NA
  
  sumNA <- colSums(calcul_NA)
  sum_total <- sum(sumNA)
  
  alerte <- NULL
  if(sum_total==0){
    alerte <- c(alerte, "Data Matrix contains no NA.\n")
  }
  
  if(length(alerte) != 0){ 
    cat(alerte,"\n")
  }
  
  table_NA <- cbind(calcul_NA, pct_NA)
  
  # check columns names ----------------------------------------
  
  # Fold
  VM.names <- colnames(VM)
  fold.names <- colnames(fold)
  
  for (i in 1:length(VM.names)){
    for (j in 1:length(fold.names)){
      if (VM.names[i]==fold.names[j]){
        fold.names[j] <- paste(fold.names[j],"2", sep="_")
      }
    }
  }
  colnames(fold) <- fold.names
  
  # NA
  NA.names <- colnames(table_NA)
  
  for (i in 1:length(VM.names)){
    for (j in 1:length(NA.names)){
      if (VM.names[i]==NA.names[j]){
        NA.names[j] <- paste(NA.names[j],"2", sep="_")
      }
    }
  }
  colnames(table_NA) <- NA.names
  
  #for NA barplots ---------------------------------------------
  
  data_bp <- data.frame()
  
  for (j in 1:ncol(pct_NA)){
    Nb_NA_0_20 <- 0
    Nb_NA_20_40 <- 0
    Nb_NA_40_60 <- 0
    Nb_NA_60_80 <- 0
    Nb_NA_80_100 <- 0
    for (i in 1:nrow(pct_NA)){
      
      if ((0<=pct_NA[i,j])&(pct_NA[i,j]<20)){
        Nb_NA_0_20=Nb_NA_0_20+1
      }
      
      if ((20<=pct_NA[i,j])&(pct_NA[i,j]<40)){
        Nb_NA_20_40=Nb_NA_20_40+1}
      
      if ((40<=pct_NA[i,j])&(pct_NA[i,j]<60)){
        Nb_NA_40_60=Nb_NA_40_60+1}
      
      if ((60<=pct_NA[i,j])&(pct_NA[i,j]<80)){
        Nb_NA_60_80=Nb_NA_60_80+1}   
      
      if ((80<=pct_NA[i,j])&(pct_NA[i,j]<=100)){
        Nb_NA_80_100=Nb_NA_80_100+1}   
    }
    data_bp[1,j] <- Nb_NA_0_20
    data_bp[2,j] <- Nb_NA_20_40
    data_bp[3,j] <- Nb_NA_40_60
    data_bp[4,j] <- Nb_NA_60_80
    data_bp[5,j] <- Nb_NA_80_100
  }
  rownames(data_bp) <- c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%")
  colnames(data_bp) <- classnames
  data_bp <- as.matrix(data_bp)
  
  
  # Output--------------------------------------------------------
  
  VM <- cbind(VM,fold,table_NA)
  
  write.table(VM, VM.output,sep="\t", quote=FALSE, row.names=FALSE)
  
  #graphics pdf
  
  pdf(graphs.output)
  
  #Boxplots for NA
  
  par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
  
  bp=barplot(data_bp, col=rainbow(nrow(data_bp)), main="Proportion of NA", xlab="Classes", ylab="Variables")
  legend("topright", fill=rainbow(nrow(data_bp)),rownames(data_bp), inset=c(-0.3,0))
  
  
  stock=0
  for (i in 1:nrow(data_bp)){
    text(bp, stock+data_bp[i,]/2, data_bp[i,], col="white")
    stock <- stock+data_bp[i,]
  }
  
  
  #Boxplots for fold test
  for (j in 1:ncol(fold)){
    title=paste(fold.names[j])
    boxplot(fold[j], main=title)
  }
  
  dev.off()

}  


# Function call---------------

#setwd("Y:\\Developpement\\Intensity check\\Pour tests")
#intens_check("DM_NA.tabular", "SM_NA.tabular", "VM_NA.tabular", "One_class", "class", "Blanks", "VM_oneclass_NA.txt", 
#"Barplots_and_Boxplots")