[BioC] bluefuse export and limma read.maimages

Steve Taylor stephen.taylor at molecular-sciences.ox.ac.uk
Thu Jul 7 16:45:46 CEST 2005


Hi Greg,

That's great! Are there any plans to add this facility to limmaGUI?

Kind Regards,

Steve

> 
> here are some hacks added to read.maimages() function in LIMMA package, 
> in order to read BlueFuse export files.
> It works fine for me.
> 
> Let me know.
> Greg
> 
> 
> #############################
> # read bluefuse export files headers #
> #############################
> my.readBlueFuseHeader <- function(file){
>     con <- file(file, "r")
>     on.exit( close( con ) )
>         error <- 0
>         i<-0
>    
>         # First line contains BlueFuse version    
>     i <- i+1
>     if( length( grep( "BlueFuse", readLines(con, n = 1) ) ) == 0 ){
>         cat( "Humm...I didnt find any BlueFuse expression in the fisrt 
> line !\n" )
>         error <- error + 1
>     }
>        
>     # Second line is a blank line
>     i <- i+1
>     if( ( nbsp <- readLines(con, n = 1) ) != "" ){
>         cat( "file format error founded in line 2 !\n" )
>         error <- error + 1
>     }
>    
>     if( error == 2 )
>         stop( "There is almost likely a problem with your file ; I am 
> stopping now !" )
>        
>    
>         # From the Third line begin headers until a new blank line from 
> which begin data
>         out <- list()
>         while( ( line <- readLines(con, n = 1) ) != "" ){
>             i<-i+1
>         txt <- strsplit( sub( "\"$", "", sub( "^\"", "", line ) ), split 
> = ": ")[[1]]
>           out[[i]] <- txt[2]
>           names(out)[i] <- txt[1]
>     }
>         out$NHeaderRecords <- i + 1
>         out
> }
> 
> 
> ################
> # bluefuse enabled #
> ################
> my.read.maimages <- function (files, source = "spot", path = NULL, ext = 
> NULL, names = NULL,
>     columns = NULL, other.columns = NULL, annotation = NULL,
>     wt.fun = NULL, verbose = TRUE, sep = "\t", quote = "\"",
>     ...)
> {
>     if (missing(files)) {
>         if (missing(ext))
>             stop("Must specify input files")
>         else {
>             extregex <- paste("\\.", ext, "$", sep = "")
>             files <- dir(path = ifelse(is.null(path), ".", path),
>                 pattern = extregex)
>             files <- sub(extregex, "", files)
>         }
>     }
>     ## modif ##
>     source <- match.arg(source, c( "bluefuse", "agilent", "arrayvision", 
> "genepix",
>         "imagene", "quantarray", "smd.old", "smd", "spot", 
> "spot.close.open"))
>     if (source == "imagene")
>         return(read.imagene(files = files, path = path, ext = ext,
>             names = names, columns = columns, wt.fun = wt.fun,
>             verbose = verbose, sep = sep, quote = quote, ...))
>     slides <- as.vector(as.character(files))
>     if (!is.null(ext))
>         slides <- paste(slides, ext, sep = ".")
>     nslides <- length(slides)
>     if (is.null(names))
>         names <- removeExt(files)
>     if (is.null(columns))
>         columns <- switch(source,
>           agilent = list(Gf = "gMeanSignal", Gb = "gBGMedianSignal", Rf 
> = "rMeanSignal", Rb = "rBGMedianSignal"),
>            
>         smd.old = list(Gf = "CH1I_MEAN", Gb = "CH1B_MEDIAN", Rf = 
> "CH2I_MEAN", Rb = "CH2B_MEDIAN"),
>        
>         smd = list(Gf = "Ch1 Intensity (Mean)", Gb = "Ch1 Background 
> (Median)", Rf = "Ch2 Intensity (Mean)",
>                     Rb = "Ch2 Background (Median)"),
>         spot = list(Rf = "Rmean", Gf = "Gmean", Rb = "morphR", Gb = 
> "morphG"),
>            
>         spot.close.open = list(Rf = "Rmean", Gf = "Gmean", Rb = 
> "morphR.close.open", Gb = "morphG.close.open"),
>            
>         genepix = list(Rf = "F635 Mean", Gf = "F532 Mean", Rb = "B635 
> Median", Gb = "B532 Median"),
>        
>         quantarray = list(Rf = "ch2 Intensity", Gf = "ch1 Intensity", Rb 
> = "ch2 Background", Gb = "ch1 Background"),
> 
>     # modif #   
>         bluefuse = list( Rf="AMPCH2", Gf="AMPCH1", Rb="AMPCH2", 
> Gb="AMPCH1" )
>             )
>            
>     fullname <- slides[1]
>     if (!is.null(path))
>         fullname <- file.path(path, fullname)
>     switch(source, quantarray = {
>         firstfield <- scan(fullname, what = "", sep = "\t", flush = TRUE,
>             quiet = TRUE, blank.lines.skip = FALSE, multi.line = FALSE)
>         skip <- grep("Begin Data", firstfield)
>         if (length(skip) == 0)
>             stop("Cannot find \"Begin Data\" in image output file")
>         nspots <- grep("End Data", firstfield) - skip - 2
>         obj <- read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", nrows = nspots, ...)
>     }, arrayvision = {
>         skip <- 1
>         cn <- scan(fullname, what = "", sep = sep, quote = quote,
>             skip = 1, nlines = 1, quiet = TRUE)
>         fg <- grep(" Dens - ", cn)
>         if (length(fg) != 2)
>             stop(paste("Cannot find foreground columns in", fullname))
>         bg <- grep("Bkgd", cn)
>         if (length(fg) != 2)
>             stop(paste("Cannot find background columns in", fullname))
>         columns <- list(Rf = fg[1], Rb = bg[1], Gf = fg[2], Gb = bg[2])
>         obj <- read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", ...)
>         nspots <- nrow(obj)
>     }, genepix = {
>         skip <- readGPRHeader(fullname)$NHeaderRecords
>         obj <- read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", fill = TRUE, ...)
>         nspots <- nrow(obj)
>     }, smd.old = {
>         skip <- readSMDHeader(fullname)$NHeaderRecords
>         obj <- read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", fill = TRUE, ...)
>         nspots <- nrow(obj)
>     }, smd = {
>         skip <- readSMDHeader(fullname)$NHeaderRecords
>         obj <- read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", fill = TRUE, ...)
>         nspots <- nrow(obj)
>     }, bluefuse = {
>     # MODIF #
>         skip <- my.readBlueFuseHeader( fullname )$NHeaderRecords
>     obj <<-  read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", fill = TRUE, ...)
>            nspots <- nrow(obj)
>     },{
>         skip <- grep(protectMetachar(columns$Rf), readLines(fullname,
>             n = 80)) - 1
>         if (length(skip) == 0)
>             stop("Cannot find column heading in image output file")
>         else skip <- skip[1]
>         obj <- read.table(fullname, skip = skip, header = TRUE,
>             sep = sep, quote = quote, as.is = TRUE, check.names = FALSE,
>             comment.char = "", fill = TRUE, ...)
>         nspots <- nrow(obj)
>     })
>     Y <- matrix(0, nspots, nslides)
>     colnames(Y) <- names
>     RG <- list(R = Y, G = Y, Rb = Y, Gb = Y)
>     if (!is.null(wt.fun))
>         RG$weights <- Y
>     RG$targets <- data.frame(FileName = I(files), row.names = names)
> 
>     if (is.null(annotation))
>         annotation <- switch(source,
>             agilent = c("Row", "Col",
>             "Start", "Sequence", "SwissProt", "GenBank", "Primate",
>             "GenPept", "ProbeUID", "ControlType", "ProbeName",
>             "GeneName", "SystematicName", "Description"),
>        
>         genepix = c("Block",
>             "Row", "Column", "ID", "Name"),
>        
>         smd = c("Spot", "Clone ID",
>             "Gene Symbol", "Gene Name", "Cluster ID", "Accession",
>             "Preferred name", "Locuslink ID", "Name", "Sequence Type",
>             "X Grid Coordinate (within sector)", "Y Grid Coordinate 
> (within sector)",
>             "Sector", "Failed", "Plate Number", "Plate Row",
>             "Plate Column", "Clone Source", "Is Verified", "Is 
> Contaminated",
>             "Luid"),
>        
>         smd.old = c("SPOT", "NAME", "Clone ID",
>             "Gene Symbol", "Gene Name", "Cluster ID", "Accession",
>             "Preferred name", "SUID"),
>     # modif #
>         bluefuse = c( "BLOCK", "ID", "SUBGRIDROW", "SUBGRIDCOL", "NAME", 
> "FLAG" ),   
>        
>         NULL)
> 
>     if (!is.null(annotation)) {
>         j <- match(annotation, colnames(obj), 0)
>  
>      objGlob <<- obj
>  
>       # modif #
>       if( source == "bluefuse" )
>           names(obj)[c(3,4,6,7,8,10)] <- c( "Row", "Column", "Block", 
> "Name", "ID", "Flags")
> 
>    
>       if (any(j > 0))
>            RG$genes <- data.frame(obj[, j, drop = FALSE])
>     }
> 
>     if (source == "agilent") {
>         if (!is.null(RG$genes$Row) && !is.null(RG$genes$Col)) {
>             nr <- length(unique(RG$genes$Row))
>             nc <- length(unique(RG$genes$Col))
>             if (nspots == nr * nc)
>                 RG$printer <- list(ngrid.r = 1, ngrid.c = 1,
>                   nspot.r = nr, nspot.c = nc)
>         }
>     }
> 
>    
>    
>     if (!is.null(other.columns)) {
>         other.columns <- as.character(other.columns)
>         j <- match(other.columns, colnames(obj), 0)
>         if (any(j > 0)) {
>             other.columns <- colnames(obj)[j]
>             RG$other <- list()
>             for (j in other.columns) RG$other[[j]] <- Y
>         }
>         else {
>             other.columns <- NULL
>         }
>     }
>     
>     for (i in 1:nslides) {
> 
>           if (i > 1) {
>                 fullname <- slides[i]
>                 if (!is.null(path))
>                     fullname <- file.path(path, fullname)
>                 if (source == "genepix")
>                     skip <- readGPRHeader(fullname)$NHeaderRecords
> 
>             # modif #
>             if ( source == "bluefuse" ){
>                 skip <- my.readBlueFuseHeader(fullname)$NHeaderRecords
>                     obj <- read.table(fullname, skip = skip, header = TRUE,
>                         sep = sep, as.is = TRUE, quote = quote, 
> check.names = FALSE,
>                             comment.char = "", fill = TRUE, nrows = nspots,
>                             ...)
>                 names(obj)[c(3,4,6,7,8,10)] <- c( "Row", "Column", 
> "Block", "Name", "ID", "Flags")
>             }
> 
>             }
> 
>         RG$R[, i] <- obj[, columns$Rf]
>         RG$G[, i] <- obj[, columns$Gf]
>         RG$Rb[, i] <- obj[, columns$Rb]
>         RG$Gb[, i] <- obj[, columns$Gb]
>     
>         if (!is.null(wt.fun))
>                 RG$weights[, i] <- wt.fun(obj)
>        
>         if (!is.null(other.columns))
>               for (j in other.columns){
>                     RG$other[[j]][, i] <- obj[, j]
>                 }
>             if (verbose)
>                 cat(paste("Read", fullname, "\n"))
>     }
>     new("RGList", RG)
> }
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
>



More information about the Bioconductor mailing list