NutriTools.R

From Organic Design wiki
Revision as of 03:29, 9 May 2006 by Sven (talk | contribs) (newer versions of Wikitable.R functions)

createPlots <- function(RG, plotDir=getwd(), 

                       bc.correct     = c("none","subtract","minimum", "movingmin"),
                       method.within  = c("none", "loess", "printtiploess"),
                       method.between = c("none","scale", "quantile")) {

 if(exists("RG")) {
  1. 1) Create directories
   dir.create(file.path(plotDir, "Densityplots"), recursive=TRUE)
   for(i in c("ImagePlots","MAplots")) {
     dir.create(file.path(plotDir, i), recursive=TRUE)
     if(i == "ImagePlots") {
       for(x in bc.correct) {
         for(y in method.within) {
           for(z in method.between) {
             for( j in c("R","G", "M","A")) {
               dir.create(file.path(plotDir, i, paste(x,"-",y,"-",z, sep=""),j), recursive=TRUE)
               if(x=="none" & y =="none" & z =="none") {
                 for(k in c("Rb", "Gb")) {
                   dir.create(file.path(plotDir, i, paste(x,"-",y,"-",z, sep=""),k), recursive=TRUE)
                 }
               }
             }
           }
         }
       }
     }
     if(i == "MAplots") {
       for(x in bc.correct) {
         for(y in method.within) {
           for(z in method.between) {
             dir.create(file.path(plotDir, i, paste(x,"-",y,"-",z, sep="")), recursive=TRUE)
           }
         }
       }
     }
   }
  1. 2) Create plots
   for(i in bc.correct) {
     for( j in method.within) {
       MA <- normalizeWithinArrays(RG, method=j, bc.method=i)
       for (k in method.between) {
         MA.Between <- normalizeBetweenArrays(MA, method=k)
         print(paste("[ Graphing:",i,j,k,"]", sep=" "))
         # plotDensities
         png(file.path(plotDir, "Densityplots", paste(i,"-",j,"-",k,".png", sep="")))
         plotDensities(MA.Between)
         dev.off()
         # plotMA3by2
         plotMA3by2(MA.Between, path=file.path(plotDir, "MAplots", paste(i,"-",j,"-",k, sep="")))
         # Imageplots
         imageplot3by2(MA.Between, z="M", path = file.path(plotDir,"Imageplots",  paste(i,"-",j,"-",k, sep=""), "M"))
         imageplot3by2(MA.Between, z="A", path = file.path(plotDir,"Imageplots",  paste(i,"-",j,"-",k, sep=""), "A"))
         imageplot3by2(RG.MA(MA.Between), z="G", path = file.path(plotDir,"Imageplots", paste(i,"-",j,"-",k, sep=""),"G"),
                       low="black", high="green")
         imageplot3by2(RG.MA(MA.Between), z="R", path = file.path(plotDir,"Imageplots", paste(i,"-",j,"-",k, sep=""),"R"),
                       low="black", high="red")
         if(i=="none" & j=="none" & k=="none") {
           imageplot3by2(RG, z="Gb", path = file.path(plotDir,"Imageplots", paste(i,"-",j,"-",k, sep=""),"Gb"),
                         low="black", high="green")
           imageplot3by2(RG, z="Rb", path = file.path(plotDir,"Imageplots", paste(i,"-",j,"-",k, sep=""),"Rb"),
                         low="black", high="red")
         }
       }
     }
   }
 }
 return("done")

}

dir.exists <- function(x, curdir = ".", rootdir) {

  1. -------------------------------------------------------------------------------------- #
  2. Author : Marcus Davy
  3. This is a recursive program to take in a string such as:
  4. foo/fodda/fi, and create subdirectories foo, fodda, and fi
  5. from the current directory location (if they don't already exist)
  6. -------------------------------------------------------------------------------------- #
 if( missing(rootdir)) {
   rootdir <- getwd()
 }
 setwd(curdir)

 x <- unlist(strsplit(x,"/"))
 if( !file.exists(x[1]) ){
   dir.create(x[1])
   print(paste("Creating", file.path(getwd(), x[1])))
 }
 curdir <- x[1]

 if(length(x)==1) {
   setwd(rootdir)
 } else {
   return(dir.exists(x[-1], curdir = curdir,  rootdir = rootdir))
 }

}

fixLabel <- function(x, string, sep=".") {

  1. ----------------------------------------------------------------------------- #
  2. Function for fixing labels for limma package functions:
  3. modelMatrix and make Contrasts e.g. 14DAFB converts to DAFB.14
  4. Author: Marcus Davy
  5. ----------------------------------------------------------------------------- #
 beginning <- regexpr(string, x)
 ending    <- nchar(x)

 labelstart   <- substr(x, beginning, ending)
 labelend <- substr(x, 0, beginning - 1)
 newlabel <- paste(labelstart, labelend, sep=sep)
 
 return(newlabel)

}

"gprDiagnostics" <- 

 function(fnames, path = ".", annotation = c("Block","Column","Row", "Name", "ID")) {
   skipcol <- c()
   dimensions <- data.frame(matrix(NA,nc=2, nr=length(fnames)))
   dimnames(dimensions) <- list(fnames, c("nrows", "ncols"))
   diagnostics <- list()
   fullfnames <- file.path(path, fnames)
   for(i in 1:length(fullfnames)){
     y <- readLines(fullfnames[i], n = 100)
     skipcol[fnames[i]] <- grep(annotation[1], y)[1] - 1
     h <- scan(fullfnamesi, quiet=TRUE, what=character(1), sep="\t", skip = skipcol[fnames[i]], quote="\"", nlines=1)    
     names(h) <- h
     h[] <- rep("NULL", length(h))
     h[annotation] <- NA
     dimensions[fnamesi,] <- dim(read.table(fullfnames[i], sep="\t", quote="\"", colClasses=h, skip=skipcol[fnames[i]], header=T))
   }
   list(skipcol, dimensions)
 }

NApadding <- function(OBJ) { 

 if (is(OBJ, "RGList")) {
   OBJ$R[OBJ$weights==0]  <- NA
   OBJ$G[OBJ$weights==0]  <- NA
   OBJ$Gb[OBJ$weights==0] <- NA
   OBJ$Rb[OBJ$weights==0] <- NA
 } else {
   OBJ$M[OBJ$weights==0]  <- NA
   OBJ$A[OBJ$weights==0]  <- NA
 }
 return(OBJ)

}

padGPR <- function(RG) { 

 mapGPRind <- function(RG){
   N <- apply(RG$genes[,c("Block","Column","Row")],2, max)
   ind <- N["Column"] * (RG$genes[,"Row"]-1) + RG$genes[,"Column"] + (RG$genes[,"Block"] - 1) * (N["Row"] * N["Column"])
   return(ind)
 }

 OUT         <- new("RGList")
 OUT$R       <- matrix(NA, nr=prod(unlist(RG$printer)), nc=ncol(RG))
 OUT$Rb      <- matrix(NA, nr=prod(unlist(RG$printer)), nc=ncol(RG))
 OUT$G       <- matrix(NA, nr=prod(unlist(RG$printer)), nc=ncol(RG))
 OUT$Gb      <- matrix(NA, nr=prod(unlist(RG$printer)), nc=ncol(RG))
 OUT$weights <- matrix(NA, nr=prod(unlist(RG$printer)), nc=ncol(RG))
 OUT$genes   <- data.frame(matrix(NA, nr=prod(unlist(RG$printer)), nc=ncol(RG$genes)))  # data.frame coersion critical
 
 OUT$targets       <- RG$targets
 OUT$source        <- RG$source
 OUT$printer       <- RG$printer
 ind               <- mapGPRind(RG)
 OUT$R[ind,]       <- RG$R
 OUT$Rb[ind,]      <- RG$Rb
 OUT$G[ind,]       <- RG$G
 OUT$Gb[ind,]      <- RG$Gb
 OUT$weights[ind,] <- RG$weights
 OUT$genes[ind,]   <- RG$genes
 return(OUT)

}

readAgilent <- function(dataDir="/Volumes/HD2/Data/Nutrigenomics/MultipleScans/Agilent"){ 

 require(limma)
  1. 1) Read in data using limma
 if(exists("targets", inherits=FALSE)) {
   RGsub <- read.maimages(targets$Filename, path = dataDir, source="agilent", name=targets$Filetype)
 } else {
   RGsub <- read.maimages(dir(dataDir, pattern=".txt"), path = dataDir, source="agilent")
 }
  1. 2) Find Agilent indices from x ,y coords

getInd <- function(x,y, nx, byrow=T) { 

 return( (x-1)* nx + y)

} 

 ind <- getInd(RGsub$genes$Row, RGsub$genes$Col, nx=215)
  1. 3) Construct complete RGList
 RGfull <- new("RGList")
 RGfull$R <-RGfull$G <-RGfull$Rb <-RGfull$Gb <- matrix(NA, nc=ncol(RGsub), nr=max(ind))
 RGfull$genes <- data.frame(matrix(NA, nc=7, nr=max(ind), dimnames = list(NULL, names(RGsub$genes))))

 matrixNames <- dimnames(RGsub$R)
  1. 4) Populate complete list
 RGfull$R[ind,]     <- RGsub$R
 RGfull$Rb[ind,]    <- RGsub$Rb
 RGfull$G[ind,]     <- RGsub$G
 RGfull$Gb[ind,]    <- RGsub$Gb
 RGfull$genes[ind,] <- RGsub$genes
 RGfull$targets <- RGsub$targets

 dimnames(RGfull$R) <- dimnames(RGfull$Rb) <- dimnames(RGfull$G) <- dimnames(RGfull$Gb) <- matrixNames
  1. Checking the rows that needed padding
 rowstopad <- as.numeric(rownames(RGsub$genes) [ !(rownames(RGsub$genes) %in% ind) ])
 rowstopad
  1. Indexing looks ok
 RGfull[rowstopad,]

 RGfull$printer <- structure(list(ngrid.r=1, ngrid.c=1, nspot.r=105, nspot.c=215), class = "PrintLayout")
 rm(RGsub)
 return(RGfull)

}


"readGPR" <- 

 function (sourceType = "median", dataDir = "/Volumes/HD2/Data/Nutrigenomics/MultipleScans/GPR", 
           ...)

{ 

 require(limma)
 print(sessionInfo())
 source.type <- switch(sourceType, mean = "genepix", median = "genepix.median")
 print(paste("source =", source.type))
 if (exists("targets")) {
   RG <- read.maimages(targets$Filename, path = dataDir, 
                       source = source.type, wt.fun = wtflags(0), names = targets$Filename, ...)
 }
 else {
   RG <- read.maimages(dir(dataDir, pattern = ".gpr"), path = dataDir, 
                       source = source.type, wt.fun = wtflags(0), ...)
 }
 RG$printer <- getLayout(RG$genes)
 return(RG)

}

"wikitable" <- 

 function (genelist, filename, title, othernames, table.head, 
           table.center = TRUE, repository = "ll")

{ 

 nl <- "\n"
 outfile <- file(filename, "w")
 if (!(missing(title) | is.null(title))) {
   cat("=", title, "=", nl, file = outfile, sep = "")
 }
 cat("{| border=\"0\"", file = outfile, sep = "\n")
 if (!(missing(table.head) | is.null(table.head))) {
   headout <- paste("| ", "", table.head, "", sep="")
   cat("|-  align=\"center\" style=\"background:#006699; color:white\"", headout, file = outfile, sep = "\n")
   }
 if (is.list(genelist)) {
   if (all.equal(rep(length(genelist1), length(genelist)), 
                 unlist(lapply(genelist, length), use.names = FALSE))) {
     nrows <- length(genelist1)
   } else {
     stop("The lists of genes to annotate must all be of equal length.")
   }
 } else {
   if(!is.null(genelist)) { # nrows from either genelist or othernames
     nrows <- length(genelist)
   } else {
     nrows <- length(othernames1)
   }
 }
 if(!is.null(repository)) {
   if (is.list(repository)) {
     rows <- ""
     for (i in seq(along = repository)) {
       rows <- paste(rows, getWikiRows(genelisti, repositoryi), sep="\n")
     }
   } else {
     rows <- getWikiRows(genelist, repository, sep="\n")
   }
 } else {
   rows <- NULL
 }
   if (!missing(othernames)) {
     if (is.list(othernames)) {
       others <- ""
       for (nm in othernames) {
         if (is.matrix(nm)) {
           for (i in 1:dim(nm)[2]) {
             others <- paste(others, "\n| ", nm[, i], sep = "")
           }
       }
       if (is.list(nm)) {
         out <- vector()
         for (j in seq(along = nm)) {

out[j] <- paste("

", nmj, "

", sep = "", 

                           collapse = "")
         }
         out <- paste("\n| ", out, sep = "")
         others <- paste(others, out, sep = "")
       }
       if ((is.vector(nm) || is.factor(nm)) && !is.list(nm)) 
         others <- paste(others, "\n| ", nm, sep = "")
     }
   } else {
     others <- NULL
   }
   if( !(is.null(rows) | is.null(others))) {
     if (length(rows) != length(others)) 
       stop(paste("There are", length(rows), "rows in your genelist, but", 
                  length(others), "rows in othernames.\n This will not give", 
                  "good results!\n"))
  }
     rows <- paste(rows, others)
     
 }
 for (i in 1:nrows) cat("|- bgcolor=\"#f2f2f2\"", rows[i], nl, file = outfile, 
                        sep = "")
 cat("|}", file = outfile, sep="")
 close(outfile)

}


"getWikiRows" <- 

 function (ids, repository = "ug")

{ 

 out <- paste("|", getWikiCells(ids, repository), sep = " ")
 return(out)

}

"getWikiCells" <- 

 function (ids, repository = "ug")

{ 

 if (is.list(ids)) {
   out <- vector()
   temp <- lapply(ids, getQueryLink, repository = repository)
   for (i in seq(along = ids)) {
     if (temp[i] != " ") {

out[i] <- paste("

[", tempi, " ", idsi, "]

", sep = "", collapse = "") 

     }
     else out[i] <- temp[i]
   }
 }
 else {
   temp <- getQueryLink(ids, repository)
   blanks <- temp == " "
   out <- paste("[", temp, " ", ids, "]", sep = "")
   out[blanks] <- " "
 }
 return(out)

}

formatTopTable <- function(OBJ, ...) { 

 if(!is(fit2, "MArrayLM")) {
   return("OBJ is not of MArrayLM class")
 }
 genelist <- topTable(OBJ, ...)
 # Remove the B statistic column
 genelist <- genelist[,-grep("^B$", colnames(genelist))]
 rawt <- fit2$coef / fit2$stdev.unscaled / fit2$sigma
 colnames(genelist)[grep("^t$", colnames(genelist))] <- "Mod.t"
 return(cbind(genelist[,1:(grep("^Mod.t$", colnames(genelist))-1)],
              "Raw.t" = rawt[as.numeric(rownames(genelist))],
              genelist[,grep("^Mod.t$", colnames(genelist)):length(colnames(genelist))])
        )

}

Retrieved from "https://organicdesign.nz/wiki/index.php?title=NutriTools.R&oldid=57073"