R/swathData_assignToImage.R

In markdunning/beadarray-devel: Quality assessment and low-level analysis for Illumina BeadArray data

assignToImage <- function(txt, sectionName, inputDir = NULL, twocolour = TRUE, locslist, GrnTiffSuffix1 = "-Swath1_Grn.tif", GrnTiffSuffix2 = "-Swath2_Grn.tif", verbose = FALSE) {
    
    ## Function to determine which swath each bead in the .txt came from
    ##
    ## Arguments:  txt - 4 or 7 column matrix generated by reading bead-level text file
    ##             sectionName - string containing name of section to read
    ##             twocolour - is this a two colour array or not
    ##             locslist - list of locs files associated with the array
    ##             GrnTiffSuffix1 - suffix of tif file from green swath 1 
    ##             GrnTiffSuffix2 - suffix of tif file from green swath 2 
    ##             verbose - how much should be reported to the user?
    ##
    ## Output: 5 or 8column matrix identical to input but with an additional  
    ##         column defining which swath the bead came from
    

#     ## First check that the stated number of channels is consistent with the locs list
#     if(twocolour){
#         if(length(locslist) != 4){stop("Number of locs files does not match number of colours")}
#     }
#     if(!twocolour){
#         if(length(locslist) != 2){stop("Number of locs files does not match number of colours")}
#     }

    ## combine the txt and locs files using BeadDataPackR
        
    #locsg <- rbind(locslist$glocs1,locslist$glocs2)
    locsg <- matrix(ncol = 2, nrow = 0)
    for(i in 1:length(locslist$Grn)) 
        locsg <- rbind(locsg, locslist$Grn[[i]])
    if(twocolour){
        #locsr <- rbind(locslist$rlocs1,locslist$rlocs2)
        locsr <- matrix(ncol = 2, nrow = 0)
        for(i in 1:length(locslist$Red)) 
            locsr <- rbind(locsr, locslist$Red[[i]])
    }

    if(verbose) cat("Combining files... ")
    if(twocolour){
        tmp <- BeadDataPackR:::combineFiles(txt, locsGrn = locsg, locsRed = locsr)
    }
    else {
        tmp <- BeadDataPackR:::combineFiles(txt, locsGrn = locsg)
    }
    if(verbose) cat("Done\n")
    #idx <- NULL
    #for(i in 1:length(locslist$Grn)) 
    #    idx <- c(idx, rep(i, nrow(locslist$Grn[[i]])))
    #idx <- idx[order(order(tmp[, "LocsIdx"]))]
    #tmp[, "LocsIdx"] <- idx
    #colnames(tmp)[colnames(tmp) == "LocsIdx"] <- "SwathIdx";

    ## remove non-decoded beads
    nondecoded <- tmp[which(tmp[,"Code"] == 0),]
    tmp <- tmp[-which(tmp[,"Code"] == 0),]

    ## if the result is smaller than the original .txt we must have some duplicates.
    ## find which they are and try using intensities to determine which image they're from
    unassigned <- NULL;
    if( nrow(tmp) < nrow(txt) ) {

        if(verbose) cat("Matching duplicate spots... ")

        if(length(locslist$Grn) == 2) {
            #unassigned <- duplicateSpots_2swaths(inputDir = inputDir, sectionName = sectionName, txt = txt, tmp = tmp, locsLength = nrow(locslist$Grn[[1]]), verbose = verbose)
            message("Assigning duplicated spot not currently supported;\n ", nrow(txt) - nrow(tmp), " bead(s) cannot be assigned to a swath and will be ignored");
        }
        else if(length(locslist$Grn) == 3) {
            message("Assigning duplicated spot not supported with 3 swaths;\n ", nrow(txt) - nrow(tmp), " beads cannot be assigned to a swath and will be ignored");
        }


    }
            
    tmp <- rbind(tmp, nondecoded, unassigned)
    ## reorder by probeID
    tmp <- tmp[order(tmp[,"Code"]),]
    return(tmp)

}


duplicateSpots_2swaths <- function(inputDir, sectionName, txt, tmp, locsLength, verbose) {
    ## if the tiff files exist we can calculate intesities from both images and compare to the txt file
    tiffFile1 <- file.path(inputDir, paste(sectionName, "-Swath1_Grn.tif", sep = ""))
    tiffFile2 <- file.path(inputDir, paste(sectionName, "-Swath2_Grn.tif", sep = ""))
    if( file.exists(tiffFile1) & file.exists(tiffFile2) ) {

        roundedX <- .Call("roundLocsFileValues", tmp[,3], PACKAGE = "BeadDataPackR");
        roundedY <- .Call("roundLocsFileValues", tmp[,4], PACKAGE = "BeadDataPackR");
        ## find the duplicated coordinates and then sort them into pairs
        txtDups <- which(!paste(txt[,3], txt[,4]) %in% paste(roundedX, roundedY));
        txtDups <- txtDups[order(txt[txtDups,3], txt[txtDups,4])]
        ## create a data.frame for these unassigned beads
        unassigned <- cbind(txt[txtDups,,drop = FALSE], vector(length = length(txtDups)))
        i = 1;
        while(i <= length(txtDups)) {
            g1 <- singleBeadIntensity(tiffFile1, txt[txtDups[i],3:4])
            g2 <- singleBeadIntensity(tiffFile2, txt[txtDups[i],3:4])
            ## assign the bead to the image which gives the closest intensity
            ## we set the "LocsIdx" column to be the largest idx for the appropriate swath
            ## remember this is not really the index of the bead!!!!
            unassigned[i,ncol(unassigned)] <- locsLength * which.min( abs( diff( c(txt[txtDups[i],2], txt[txtDups[i],2], g1, g2), lag = 2 ) ) )
            ## if the pair has been decoded, assign it to the other image
            if(i < length(txtDups) & all(txt[txtDups[i+1],3:4] == txt[txtDups[i],3:4])) {
                unassigned[i+1,ncol(unassigned)] <- locsLength * which.min( abs( diff( c(txt[txtDups[i],2], txt[txtDups[i],2], g1, g2), lag = 2 ) ) )
                i = i + 1;
            }
            i = i + 1;    
        }
        if(verbose) message("Done");
    }
    else {
        message("Unable to find TIFF images; ", nrow(txt) - nrow(tmp), " beads cannot be assigned to a swath and will be ignored");
        unassigned <- matrix(nrow = 0, ncol = ncol(txt))
    }
    return(unassigned)
}