diff Intchecks/Script_intensity_check.R @ 0:51c39ea1fd54 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Intchecks/Script_intensity_check.R	Mon Jan 14 08:37:27 2019 -0500
@@ -0,0 +1,327 @@
+#########################################################################
+# SCRIPT INTENSITY CHECK                                                #
+#                                                                       #
+# 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"
+  class.col <- "2"
+  type <- "One_class"
+  class1 <- "Blanks"
+  fold.frac <- "Top"
+  logarithm <- "log2"
+  VM.output <- "new_VM.txt"
+  graphs.output <- "Barplots_and_Boxplots.pdf"
+}
+
+
+
+
+intens_check <- function(DM.name, SM.name, VM.name, class.col, type, class1, fold.frac, logarithm,
+                         VM.output, graphs.output){
+  
+  
+  # This function allows to check the intensities considering classes with a mean fold change calculation,  
+  # the number and the proportion of missing values (NA) in dataMatrix
+  #
+  # 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
+  # class.col: number of the sampleMetadata's column with classes
+  # type: "One_class" or "Each_class"
+  # class1: name of the class, only if type="One_class"
+  # fold.frac: if type="One class": class1/other ("Top") or other/class1 ("Bottom")
+  # logarithm: "log2", "log10" or "none" for log mean fold change
+  # VM.output: output file's access (VM with new columns)
+  # graphs.output: pdf file's access 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 < 2){
+    err.1class <- c("\n The column",class.col, "contains only one class, fold calculation could not be executed \n")
+    cat(err.1class)  
+  }
+  
+  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 <- c(class1,"Other")
+        
+      }
+    }
+  }
+
+  DM <- cbind(DM,c_class)
+  
+  
+  
+  # fold calculation -------------------------------------------------------------------------------------------
+  
+  if(nb_class >= 2){
+    
+    
+    fold <- data.frame()
+    n <- 1
+    ratio1 <- NULL
+    ratio2 <- NULL
+    
+    if(type=="Each_class"){
+      fold.frac <- "Top"
+    }
+    
+    for(j in 1:(nb_class-1)){
+      for(k in (j+1):nb_class) {
+        
+        if(fold.frac=="Bottom"){
+          ratio1 <- classnames[k]
+          ratio2 <- classnames[j]
+        }else{
+          ratio1 <- classnames[j]
+          ratio2 <- classnames[k]
+        }
+        
+        for (i in 1:(length(DM)-1)){
+          fold[i,n] <- mean(DM[which(DM$c_class==ratio1),i], na.rm=TRUE)/
+            mean(DM[which(DM$c_class==ratio2),i], na.rm=TRUE)
+          if(logarithm=="log2"){
+            fold[i,n] <- log2(fold[i,n])
+          }else if(logarithm=="log10"){
+            fold[i,n] <- log10(fold[i,n])
+          }
+        }
+        names(fold)[n] <- paste("fold",ratio1,"VS", ratio2, sep="_")
+        if(logarithm != "none"){
+          names(fold)[n] <- paste(logarithm,names(fold)[n], sep="_")
+        }
+        n <- n + 1}
+    }
+    
+  }
+  
+  # number and proportion of 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("NA",classnames, sep="_")
+  names(pct_NA) <- paste("Pct_NA", classnames, sep="_")
+  
+  # Alert message if there is no NA in data matrix
+  
+  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 ---------------------------------------------------------------------------------------
+  
+
+  VM.names <- colnames(VM)
+  
+  # Fold
+  
+  if(nb_class >=2){
+    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
+    
+    VM <- cbind(VM,fold)
+  }
+  
+  # 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
+  VM <- cbind(VM,table_NA)
+  
+  
+  #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 ---------------------------------------------------------------------------------------------------
+  
+  
+  write.table(VM, VM.output,sep="\t", quote=FALSE, row.names=FALSE)
+  
+  #graphics pdf
+  
+  pdf(graphs.output)
+  
+  #Barplots 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", cex=0.7)
+    stock <- stock+data_bp[i,]
+  }
+  
+  
+  #Boxplots for fold test
+  
+  if(nb_class >= 2){
+    
+    clean_fold <- fold
+    for(i in 1:nrow(clean_fold)){
+      for(j in 1:ncol(clean_fold)){
+        if(is.infinite(clean_fold[i,j])){
+          clean_fold[i,j] <- NA
+        }
+      }
+    }
+    for (j in 1:ncol(clean_fold)){
+      title <- paste(fold.names[j])
+      boxplot(clean_fold[j], main=title)
+    }
+    }
+  
+  dev.off()
+
+}  
+
+  
\ No newline at end of file