R/swathData_genSwaths.R
In markdunning/beadarray: Quality assessment and low-level analysis for Illumina BeadArray data
genSwaths <- function(txt, sectionName, twocolour = TRUE, GrnTiffSuffix2 = "-Swath2_Grn.tif", RedTiffSuffix2 = "-Swath2_Red.tif", section_height = 326, section_width = 397, locslist, verbose = TRUE) {
## Function to identify beads in the overlapping region between swaths
## For these beads intensities are calculated from both images for comparison
##
## Arguments: txt - 5 column matrix generated by assignToImage()
## sectionName - string containing name of section to read
##
## Output: list containing two 4-column matrices in the same format as bead-level txt files
if(twocolour){
if( length(locslist$Red) != length(locslist$Grn) ) {
stop("There should be an equal number of Grn and Red locs files")
}
}
swath_overlap <- (2 * nrow(locslist$Grn[[1]]) / (9 * section_height)) - section_width
glocs1<-locslist$Grn[[1]]
glocs2<-locslist$Grn[[2]]
rlocs2 <- NULL
locsg <- rbind(glocs1, glocs2)
if(twocolour){
rlocs1<-locslist$Red[[1]]
rlocs2<-locslist$Red[[1]]
locsr <- rbind(rlocs1, rlocs2)
}
## for now, assume any bead in the overlap that was decoded is in Swath1
swath1 <- txt[which(txt[ ,"LocsIdx"] <= nrow(glocs1)),]
swath2 <- txt[which(txt[ ,"LocsIdx"] > nrow(glocs1)),]
## only continue if we can find the TIFF files
if(file.exists(paste(sectionName, GrnTiffSuffix2, sep = ""))) {
## find the decoded beads in swath 1 locs file and get their IDs (this is really slow)
if(verbose) cat("Identifying decoded beads... ")
locsInTxt1 <- match(paste(glocs1[,1], glocs1[,2]), paste(swath1[,3], swath1[,4]));
locsIDs1 <- swath1[locsInTxt1,1];
locsIns1 <- swath1[locsInTxt1,2];
locsInTxt2 <- match(paste(glocs2[,1], glocs2[,2]), paste(swath2[,3], swath2[,4]));
locsIDs2 <- swath2[locsInTxt2,1];
locsIns2 <- swath2[locsInTxt2,2];
if(verbose) cat("Done\n");
if(verbose) cat("Creating locs grid... ")
grid <- t(sapply(1:nrow(glocs1), locsIndicesToGrid, nrow = section_height, ncol = (section_width+swath_overlap)/2))
if(verbose) cat("Done\n");
## indices for beads in the overlapping region
## Swath1
s1idx <- which(grid[,1] > ((section_width + swath_overlap)/2 - swath_overlap))
## Swath2
s2idx <- which(grid[,1] < swath_overlap+1)
if(length(s1idx) != length(s2idx)) stop("Something is wrong!!!!");
## create matrix with two sets of coordinates and probeIDs in the middle
overlap <- cbind(locsIDs1[s1idx],locsIns2[s2idx],glocs2[s2idx,],locsIns2[s2idx],rlocs2[s2idx,]);
# colnames(overlap) <- c("swath1.X", "swath1.Y", "ProbeID", "swath2.X", "swath2.Y")
## remove non-decoded beads
#s1idx<-s1idx[-which(is.na(overlap[,1]))]
#s2idx<-s2idx[-which(is.na(overlap[,1]))]
#overlap <- overlap[-which(is.na(overlap[,1])),]
s1idx<-s1idx[-which(overlap[,1] == 0)]
s2idx<-s2idx[-which(overlap[,1] == 0)]
overlap <- overlap[-which(overlap[,1] == 0),]
## now we need to compute intensities for both swaths
## (maybe we should keep the .txt ones too, but we'll need to keep track of them,
## which we don't do currently)
tiff <- readTIFF(paste(sectionName, GrnTiffSuffix2, sep = ""));
Gbg <- illuminaBackground(tiff, overlap[,3:4]);
tiff <- illuminaSharpen(tiff);
Gfg <- illuminaForeground(tiff, overlap[,3:4]);
overlap[,2]<-floor(Gfg-Gbg)
if(any(is.na(overlap[,2]))){
s1idx<-s1idx[-which(is.na(overlap[,2]))]
s2idx<-s2idx[-which(is.na(overlap[,2]))]
overlap <- overlap[-which(is.na(overlap[,2])),]
}
if(twocolour) {
tiff <- readTIFF(paste(sectionName, RedTiffSuffix2, sep = ""));
Rbg <- illuminaBackground(tiff, overlap[,6:7]);
tiff <- illuminaSharpen(tiff);
Rfg <- illuminaForeground(tiff, overlap[,6:7]);
overlap[,5]<-floor(Rfg-Rbg)
if(any(is.na(overlap[,5]))) {
s1idx<-s1idx[-which(is.na(overlap[,5]))]
s2idx<-s2idx[-which(is.na(overlap[,5]))]
overlap <- overlap[-which(is.na(overlap[,5])),]
}
}
else {
overlap <- overlap[,-5];
}
#overlap<-cbind(overlap,rep(2,dim(overlap)[1]))
#colnames(overlap) <- colnames(swath2)
overlap<-cbind(overlap, ((nrow(glocs1)+1):nrow(locsg))[s2idx] )
#overlap<-cbind(overlap, (1:(dim(overlap)[1])))
#swath2<-cbind(swath2, ((nrow(glocs1)+1):nrow(locsg))[-s2idx])
swath2<-rbind(swath2,overlap)
colnames(swath2)[ncol(swath2)]<-"Overlap"
swath1<-cbind(swath1,rep(0,dim(swath1)[1]))
colnames(swath1)[ncol(swath1)]<-"Overlap"
swath1[locsInTxt1[s1idx],dim(swath1)[2]]<-(1:(dim(overlap)[1]))
}
else {
## if we cant find the TIFFs then add a column to the swath data indicating there are no overlapping beads
swath1 <- cbind(swath1, rep(0, nrow(swath1)));
colnames(swath1)[ncol(swath1)] <- "Overlap";
swath2 <- cbind(swath2, rep(0, nrow(swath2)));
colnames(swath2)[ncol(swath2)] <- "Overlap";
}
## return a list containing the two swaths (Removing the "Swath Index" column since that's implied)
return(list("Swath1" = swath1[,-(ncol(swath1)-1)], "Swath2" = swath2[,-(ncol(swath2)-1)]));
}
genTwoSwaths <- function(txt, sectionName, inputDir = ".", locsList, twocolour = TRUE, GrnTiffSuffix2 = "-Swath2_Grn.tif", RedTiffSuffix2 = "-Swath2_Red.tif", segmentHeight = 326, segmentWidth = 397, verbose = TRUE) {
## Function to identify beads in the overlapping region between swaths
## For these beads intensities are calculated from both images for comparison
##
## Arguments: txt - 5 column matrix generated by assignToImage()
## sectionName - string containing name of section to read
##
## Output: list containing two 4-column matrices in the same format as bead-level txt files
if(twocolour){
if( length(locsList$Red) != length(locsList$Grn) ) {
stop("There should be an equal number of Grn and Red locs files")
}
}
txt <- txt[order(txt[,"LocsIdx"]),]
s1Width <- nrow(locsList$Grn[[1]]) / (9 * segmentHeight)
s2Width <- nrow(locsList$Grn[[2]]) / (9 * segmentHeight)
swathOverlap <- (s1Width + s2Width) - segmentWidth;
swath1 <- txt[which(txt[,ncol(txt)] <= (segmentHeight * s1Width * 9)), ]
swath2 <- txt[which(txt[,ncol(txt)] > (segmentHeight * s1Width * 9)), ]
swath1res <- matrix(ncol = ncol(swath1) + 1, nrow = nrow(swath1))
swath2res <- matrix(ncol = ncol(swath2) + 1, nrow = nrow(swath2))
swath1res[,ncol(swath1res)] <- swath2res[,ncol(swath1res)] <- 0
colnames(swath1res) <- c(colnames(swath1), "Overlap")
colnames(swath2res) <- c(colnames(swath1), "Overlap")
grid1 <- t(sapply(1:(segmentHeight * s1Width), locsIndicesToGrid, nrow = segmentHeight, ncol = s1Width))
grid2 <- t(sapply(1:(segmentHeight * s2Width), locsIndicesToGrid, nrow = segmentHeight, ncol = s2Width))
tiffG2 <- readTIFF(file.path(inputDir, paste(sectionName, GrnTiffSuffix2, sep = "")))
if(twocolour) {
tiffR2 <- readTIFF(file.path(inputDir, paste(sectionName, RedTiffSuffix2, sep = "")))
}
for(seg in 0:8) {
if(verbose)
message("Overlapping region: Segment ", seg);
swath1seg <- swath1[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)),]
swath2seg <- swath2[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)),]
idxToUpdate1 <- which(grid1[,1] > s1Width - swathOverlap)
idxToUpdate2 <- which(grid2[,1] <= swathOverlap)
## remove the nondecoded beads from swath1 - most of these lie outside the image
tmpIdx <- which(swath1seg[idxToUpdate1, 2] == 0)
idxToUpdate1 <- idxToUpdate1[-tmpIdx];
idxToUpdate2 <- idxToUpdate2[-tmpIdx];
# swath2seg[idxToUpdate2, 1] <- swath1seg[idxToUpdate1, 1]
## calculate new intensities - we need to account for which resolution images were produced
if(ncol(tiffG2) > 1024) { ## High Res (iScan)
intenG <- round(illuminaForeground_6x6(tiffG2, swath2seg[idxToUpdate2,3:4]) - illuminaBackground(tiffG2, swath2seg[idxToUpdate2,3:4]))
if(twocolour) {
intenR <- round(illuminaForeground_6x6(tiffR2, swath2seg[idxToUpdate2,6:7]) - illuminaBackground(tiffR2, swath2seg[idxToUpdate2,6:7]))
}
}
else { ## low res (BeadScan)
bg <- illuminaBackground(tiffG2, swath2seg[idxToUpdate2,3:4])
fg <- illuminaForeground(illuminaSharpen(tiffG2), swath2seg[idxToUpdate2,3:4])
intenG <- round(fg - bg)
if(twocolour) {
bg <- illuminaBackground(tiffR2, swath2seg[idxToUpdate2,6:7])
fg <- illuminaForeground(illuminaSharpen(tiffR2), swath2seg[idxToUpdate2,6:7])
intenR <- round(fg - bg)
}
}
## we can get beads with NA intensity, so lets remove those
if(twocolour)
naIdx <- which(is.na(intenG) | is.na(intenR))
else
naIdx <- which(is.na(intenG))
if(length(naIdx)) {
idxToUpdate1 <- idxToUpdate1[-naIdx]
idxToUpdate2 <- idxToUpdate2[-naIdx]
intenG <- intenG[-naIdx]
## beads in swath2 that can't get an intensity should be given ID 0
swath2seg[naIdx,1] <- 0
if(twocolour)
intenR <- intenR[-naIdx]
}
swath2seg[idxToUpdate2,2] <- intenG
if(twocolour)
swath2seg[idxToUpdate2,5] <- intenR
swath2seg[idxToUpdate2, 1] <- swath1seg[idxToUpdate1, 1]
swath1res[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)), 1:ncol(swath1seg)] <- swath1seg
swath1res[((seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)))[idxToUpdate1], ncol(swath1res)] <- 1
swath2res[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)), 1:ncol(swath2seg)] <- swath2seg
swath2res[((seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)))[idxToUpdate2], ncol(swath2res)] <- 1
}
## add a column with an index if it's a bead with multiple intensities
## for now with 3 swath data we are ignoring these beads
#swath1 <- cbind(swath1, rep(0, nrow(swath1)));
#colnames(swath1res)[ncol(swath1res)] <- "Overlap";
#swath2 <- cbind(swath2, rep(0, nrow(swath2)));
#colnames(swath2res)[ncol(swath2res)] <- "Overlap";
#swath3 <- cbind(swath3, rep(0, nrow(swath3)));
#colnames(swath3res)[ncol(swath3res)] <- "Overlap";
## return a list containing the two swaths (Removing the "Swath Index" column since that's implied)
return(list("Swath1" = swath1res[,-(ncol(swath1res)-1)], "Swath2" = swath2res[,-(ncol(swath2res)-1)]));
}
genThreeSwaths <- function(txt, sectionName, inputDir = ".", locsList, twocolour = TRUE, GrnTiffSuffix2 = "-Swath2_Grn.tif", GrnTiffSuffix3 = "-Swath3_Grn.tif", RedTiffSuffix2 = "-Swath2_Red.tif", RedTiffSuffix3 = "-Swath3_Red.tif", segmentHeight = 326, verbose = TRUE) {
## Function to identify beads in the overlapping region between swaths
## For these beads intensities are calculated from both images for comparison
##
## Arguments: txt - 5 column matrix generated by assignToImage()
## sectionName - string containing name of section to read
##
## Output: list containing two 4-column matrices in the same format as bead-level txt files
txt <- txt[order(txt[,"LocsIdx"]),]
s1Width <- nrow(locsList$Grn[[1]]) / (9 * segmentHeight)
s2Width <- nrow(locsList$Grn[[2]]) / (9 * segmentHeight)
swath1 <- txt[which(txt[,ncol(txt)] <= (segmentHeight * s1Width * 9)),]
swath2 <- txt[which(txt[,ncol(txt)] > (segmentHeight * s1Width * 9) & txt[,ncol(txt)] <= ((segmentHeight * s1Width * 9) + (segmentHeight * s2Width * 9))),]
swath3 <- txt[which(txt[,ncol(txt)] > ((segmentHeight * s2Width * 9) + (segmentHeight * s1Width * 9))),]
swath1res <- matrix(ncol = ncol(swath1) + 1, nrow = nrow(swath1))
swath2res <- matrix(ncol = ncol(swath2) + 1, nrow = nrow(swath2))
swath3res <- matrix(ncol = ncol(swath3) + 1, nrow = nrow(swath3))
swath1res[,ncol(swath1res)] <- swath2res[,ncol(swath1res)] <- swath3res[,ncol(swath1res)] <- 0
colnames(swath1res) <- c(colnames(swath1), "Overlap")
colnames(swath2res) <- c(colnames(swath1), "Overlap")
colnames(swath3res) <- c(colnames(swath1), "Overlap")
grid1 <- grid3 <- t(sapply(1:(segmentHeight * s1Width), beadarray:::locsIndicesToGrid, nrow = segmentHeight, ncol = s1Width))
grid2 <- t(sapply(1:(segmentHeight * s2Width), beadarray:::locsIndicesToGrid, nrow = segmentHeight, ncol = s2Width))
tiffG2 <- readTIFF(file.path(inputDir, paste(sectionName, GrnTiffSuffix2, sep = "")))
tiffG3 <- readTIFF(file.path(inputDir, paste(sectionName, GrnTiffSuffix3, sep = "")))
if(twocolour) {
tiffR2 <- readTIFF(file.path(inputDir, paste(sectionName, RedTiffSuffix2, sep = "")))
tiffR3 <- readTIFF(file.path(inputDir, paste(sectionName, RedTiffSuffix3, sep = "")))
}
for(seg in 0:8) {
if(verbose)
message("Overlapping region: Segment ", seg);
swath1seg <- swath1[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)),]
swath2seg <- swath2[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)),]
swath3seg <- swath3[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)),]
decodedColIdx2 <- unlist(lapply(split(swath2seg[,1], grid2[,1]), FUN = function(x) {
if(length(unique(x)) == 1)
return(TRUE)
else
return(FALSE)
}))
decodedColIdx3 <- unlist(lapply(split(swath3seg[,1], grid3[,1]), FUN = function(x) {
if(length(unique(x)) == 1)
return(TRUE)
else
return(FALSE)
}))
idxToUpdate1 <- which(grid1[,1] > s1Width - max(which(decodedColIdx2)))
idxToUpdate2a <- which(grid2[,1] <= max(which(decodedColIdx2)))
swath2seg[idxToUpdate2a,1] <- swath1seg[idxToUpdate1,1]
#swath2seg[idxToUpdate2,2] <- apply(swath2seg[idxToUpdate2,3:4], 1, FUN = singleBeadIntensity_6x6, tiffFile = file.path(paste(sectionName, "-Swath2_Grn.tif", sep = "")))
swath2seg[idxToUpdate2a,2] <- round(illuminaForeground_6x6(tiffG2, swath2seg[idxToUpdate2a,3:4]) - illuminaBackground(tiffG2, swath2seg[idxToUpdate2a,3:4]))
if(twocolour) {
swath2seg[idxToUpdate2a,5] <- round(illuminaForeground_6x6(tiffR2, swath2seg[idxToUpdate2a,6:7]) - illuminaBackground(tiffR2, swath2seg[idxToUpdate2a,6:7]))
}
idxToUpdate2b <- which(grid2[,1] > s2Width - max(which(decodedColIdx3)))
idxToUpdate3 <- which(grid3[,1] <= max(which(decodedColIdx3)))
swath3seg[idxToUpdate3,1] <- swath2seg[idxToUpdate2b,1]
#swath3seg[idxToUpdate3,2] <- apply(swath3seg[idxToUpdate3,3:4], 1, FUN = singleBeadIntensity_6x6, tiffFile = file.path(paste(sectionName, "-Swath3_Grn.tif", sep = "")))
swath3seg[idxToUpdate3,2] <- round(illuminaForeground_6x6(tiffG3, swath3seg[idxToUpdate3,3:4]) - illuminaBackground(tiffG3, swath3seg[idxToUpdate3,3:4]))
if(twocolour) {
swath3seg[idxToUpdate3,5] <- round(illuminaForeground_6x6(tiffR3, swath3seg[idxToUpdate3,6:7]) - illuminaBackground(tiffG3, swath3seg[idxToUpdate3,6:7]))
}
swath1res[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)), 1:ncol(swath1seg)] <- swath1seg
swath1res[((seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)))[idxToUpdate1], ncol(swath1res)] <- 1
swath2res[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)), 1:ncol(swath2seg)] <- swath2seg
swath2res[((seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)))[c(idxToUpdate2a, idxToUpdate2b)], ncol(swath2res)] <- 1
swath3res[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)), 1:ncol(swath3seg)] <- swath3seg
swath3res[((seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)))[idxToUpdate3], ncol(swath3res)] <- 1
}
## add a column with an index if it's a bead with multiple intensities
## for now with 3 swath data we are ignoring these beads
#swath1 <- cbind(swath1, rep(0, nrow(swath1)));
#colnames(swath1res)[ncol(swath1res)] <- "Overlap";
#swath2 <- cbind(swath2, rep(0, nrow(swath2)));
#colnames(swath2res)[ncol(swath2res)] <- "Overlap";
#swath3 <- cbind(swath3, rep(0, nrow(swath3)));
#colnames(swath3res)[ncol(swath3res)] <- "Overlap";
## return a list containing the two swaths (Removing the "Swath Index" column since that's implied)
return(list("Swath1" = swath1res[,-(ncol(swath1res)-1)], "Swath2" = swath2res[,-(ncol(swath2res)-1)], "Swath3" = swath3res[,-(ncol(swath3res)-1)]));
}
markdunning/beadarray documentation built on May 9, 2019, 8:35 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
genSwaths <- function(txt, sectionName, twocolour = TRUE, GrnTiffSuffix2 = "-Swath2_Grn.tif", RedTiffSuffix2 = "-Swath2_Red.tif", section_height = 326, section_width = 397, locslist, verbose = TRUE) {
## Function to identify beads in the overlapping region between swaths
## For these beads intensities are calculated from both images for comparison
##
## Arguments: txt - 5 column matrix generated by assignToImage()
## sectionName - string containing name of section to read
##
## Output: list containing two 4-column matrices in the same format as bead-level txt files
if(twocolour){
if( length(locslist$Red) != length(locslist$Grn) ) {
stop("There should be an equal number of Grn and Red locs files")
}
}
swath_overlap <- (2 * nrow(locslist$Grn[[1]]) / (9 * section_height)) - section_width
glocs1<-locslist$Grn[[1]]
glocs2<-locslist$Grn[[2]]
rlocs2 <- NULL
locsg <- rbind(glocs1, glocs2)
if(twocolour){
rlocs1<-locslist$Red[[1]]
rlocs2<-locslist$Red[[1]]
locsr <- rbind(rlocs1, rlocs2)
}
## for now, assume any bead in the overlap that was decoded is in Swath1
swath1 <- txt[which(txt[ ,"LocsIdx"] <= nrow(glocs1)),]
swath2 <- txt[which(txt[ ,"LocsIdx"] > nrow(glocs1)),]
## only continue if we can find the TIFF files
if(file.exists(paste(sectionName, GrnTiffSuffix2, sep = ""))) {
## find the decoded beads in swath 1 locs file and get their IDs (this is really slow)
if(verbose) cat("Identifying decoded beads... ")
locsInTxt1 <- match(paste(glocs1[,1], glocs1[,2]), paste(swath1[,3], swath1[,4]));
locsIDs1 <- swath1[locsInTxt1,1];
locsIns1 <- swath1[locsInTxt1,2];
locsInTxt2 <- match(paste(glocs2[,1], glocs2[,2]), paste(swath2[,3], swath2[,4]));
locsIDs2 <- swath2[locsInTxt2,1];
locsIns2 <- swath2[locsInTxt2,2];
if(verbose) cat("Done\n");
if(verbose) cat("Creating locs grid... ")
grid <- t(sapply(1:nrow(glocs1), locsIndicesToGrid, nrow = section_height, ncol = (section_width+swath_overlap)/2))
if(verbose) cat("Done\n");
## indices for beads in the overlapping region
## Swath1
s1idx <- which(grid[,1] > ((section_width + swath_overlap)/2 - swath_overlap))
## Swath2
s2idx <- which(grid[,1] < swath_overlap+1)
if(length(s1idx) != length(s2idx)) stop("Something is wrong!!!!");
## create matrix with two sets of coordinates and probeIDs in the middle
overlap <- cbind(locsIDs1[s1idx],locsIns2[s2idx],glocs2[s2idx,],locsIns2[s2idx],rlocs2[s2idx,]);
# colnames(overlap) <- c("swath1.X", "swath1.Y", "ProbeID", "swath2.X", "swath2.Y")
## remove non-decoded beads
#s1idx<-s1idx[-which(is.na(overlap[,1]))]
#s2idx<-s2idx[-which(is.na(overlap[,1]))]
#overlap <- overlap[-which(is.na(overlap[,1])),]
s1idx<-s1idx[-which(overlap[,1] == 0)]
s2idx<-s2idx[-which(overlap[,1] == 0)]
overlap <- overlap[-which(overlap[,1] == 0),]
## now we need to compute intensities for both swaths
## (maybe we should keep the .txt ones too, but we'll need to keep track of them,
## which we don't do currently)
tiff <- readTIFF(paste(sectionName, GrnTiffSuffix2, sep = ""));
Gbg <- illuminaBackground(tiff, overlap[,3:4]);
tiff <- illuminaSharpen(tiff);
Gfg <- illuminaForeground(tiff, overlap[,3:4]);
overlap[,2]<-floor(Gfg-Gbg)
if(any(is.na(overlap[,2]))){
s1idx<-s1idx[-which(is.na(overlap[,2]))]
s2idx<-s2idx[-which(is.na(overlap[,2]))]
overlap <- overlap[-which(is.na(overlap[,2])),]
}
if(twocolour) {
tiff <- readTIFF(paste(sectionName, RedTiffSuffix2, sep = ""));
Rbg <- illuminaBackground(tiff, overlap[,6:7]);
tiff <- illuminaSharpen(tiff);
Rfg <- illuminaForeground(tiff, overlap[,6:7]);
overlap[,5]<-floor(Rfg-Rbg)
if(any(is.na(overlap[,5]))) {
s1idx<-s1idx[-which(is.na(overlap[,5]))]
s2idx<-s2idx[-which(is.na(overlap[,5]))]
overlap <- overlap[-which(is.na(overlap[,5])),]
}
}
else {
overlap <- overlap[,-5];
}
#overlap<-cbind(overlap,rep(2,dim(overlap)[1]))
#colnames(overlap) <- colnames(swath2)
overlap<-cbind(overlap, ((nrow(glocs1)+1):nrow(locsg))[s2idx] )
#overlap<-cbind(overlap, (1:(dim(overlap)[1])))
#swath2<-cbind(swath2, ((nrow(glocs1)+1):nrow(locsg))[-s2idx])
swath2<-rbind(swath2,overlap)
colnames(swath2)[ncol(swath2)]<-"Overlap"
swath1<-cbind(swath1,rep(0,dim(swath1)[1]))
colnames(swath1)[ncol(swath1)]<-"Overlap"
swath1[locsInTxt1[s1idx],dim(swath1)[2]]<-(1:(dim(overlap)[1]))
}
else {
## if we cant find the TIFFs then add a column to the swath data indicating there are no overlapping beads
swath1 <- cbind(swath1, rep(0, nrow(swath1)));
colnames(swath1)[ncol(swath1)] <- "Overlap";
swath2 <- cbind(swath2, rep(0, nrow(swath2)));
colnames(swath2)[ncol(swath2)] <- "Overlap";
}
## return a list containing the two swaths (Removing the "Swath Index" column since that's implied)
return(list("Swath1" = swath1[,-(ncol(swath1)-1)], "Swath2" = swath2[,-(ncol(swath2)-1)]));
}
genTwoSwaths <- function(txt, sectionName, inputDir = ".", locsList, twocolour = TRUE, GrnTiffSuffix2 = "-Swath2_Grn.tif", RedTiffSuffix2 = "-Swath2_Red.tif", segmentHeight = 326, segmentWidth = 397, verbose = TRUE) {
## Function to identify beads in the overlapping region between swaths
## For these beads intensities are calculated from both images for comparison
##
## Arguments: txt - 5 column matrix generated by assignToImage()
## sectionName - string containing name of section to read
##
## Output: list containing two 4-column matrices in the same format as bead-level txt files
if(twocolour){
if( length(locsList$Red) != length(locsList$Grn) ) {
stop("There should be an equal number of Grn and Red locs files")
}
}
txt <- txt[order(txt[,"LocsIdx"]),]
s1Width <- nrow(locsList$Grn[[1]]) / (9 * segmentHeight)
s2Width <- nrow(locsList$Grn[[2]]) / (9 * segmentHeight)
swathOverlap <- (s1Width + s2Width) - segmentWidth;
swath1 <- txt[which(txt[,ncol(txt)] <= (segmentHeight * s1Width * 9)), ]
swath2 <- txt[which(txt[,ncol(txt)] > (segmentHeight * s1Width * 9)), ]
swath1res <- matrix(ncol = ncol(swath1) + 1, nrow = nrow(swath1))
swath2res <- matrix(ncol = ncol(swath2) + 1, nrow = nrow(swath2))
swath1res[,ncol(swath1res)] <- swath2res[,ncol(swath1res)] <- 0
colnames(swath1res) <- c(colnames(swath1), "Overlap")
colnames(swath2res) <- c(colnames(swath1), "Overlap")
grid1 <- t(sapply(1:(segmentHeight * s1Width), locsIndicesToGrid, nrow = segmentHeight, ncol = s1Width))
grid2 <- t(sapply(1:(segmentHeight * s2Width), locsIndicesToGrid, nrow = segmentHeight, ncol = s2Width))
tiffG2 <- readTIFF(file.path(inputDir, paste(sectionName, GrnTiffSuffix2, sep = "")))
if(twocolour) {
tiffR2 <- readTIFF(file.path(inputDir, paste(sectionName, RedTiffSuffix2, sep = "")))
}
for(seg in 0:8) {
if(verbose)
message("Overlapping region: Segment ", seg);
swath1seg <- swath1[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)),]
swath2seg <- swath2[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)),]
idxToUpdate1 <- which(grid1[,1] > s1Width - swathOverlap)
idxToUpdate2 <- which(grid2[,1] <= swathOverlap)
## remove the nondecoded beads from swath1 - most of these lie outside the image
tmpIdx <- which(swath1seg[idxToUpdate1, 2] == 0)
idxToUpdate1 <- idxToUpdate1[-tmpIdx];
idxToUpdate2 <- idxToUpdate2[-tmpIdx];
# swath2seg[idxToUpdate2, 1] <- swath1seg[idxToUpdate1, 1]
## calculate new intensities - we need to account for which resolution images were produced
if(ncol(tiffG2) > 1024) { ## High Res (iScan)
intenG <- round(illuminaForeground_6x6(tiffG2, swath2seg[idxToUpdate2,3:4]) - illuminaBackground(tiffG2, swath2seg[idxToUpdate2,3:4]))
if(twocolour) {
intenR <- round(illuminaForeground_6x6(tiffR2, swath2seg[idxToUpdate2,6:7]) - illuminaBackground(tiffR2, swath2seg[idxToUpdate2,6:7]))
}
}
else { ## low res (BeadScan)
bg <- illuminaBackground(tiffG2, swath2seg[idxToUpdate2,3:4])
fg <- illuminaForeground(illuminaSharpen(tiffG2), swath2seg[idxToUpdate2,3:4])
intenG <- round(fg - bg)
if(twocolour) {
bg <- illuminaBackground(tiffR2, swath2seg[idxToUpdate2,6:7])
fg <- illuminaForeground(illuminaSharpen(tiffR2), swath2seg[idxToUpdate2,6:7])
intenR <- round(fg - bg)
}
}
## we can get beads with NA intensity, so lets remove those
if(twocolour)
naIdx <- which(is.na(intenG) | is.na(intenR))
else
naIdx <- which(is.na(intenG))
if(length(naIdx)) {
idxToUpdate1 <- idxToUpdate1[-naIdx]
idxToUpdate2 <- idxToUpdate2[-naIdx]
intenG <- intenG[-naIdx]
## beads in swath2 that can't get an intensity should be given ID 0
swath2seg[naIdx,1] <- 0
if(twocolour)
intenR <- intenR[-naIdx]
}
swath2seg[idxToUpdate2,2] <- intenG
if(twocolour)
swath2seg[idxToUpdate2,5] <- intenR
swath2seg[idxToUpdate2, 1] <- swath1seg[idxToUpdate1, 1]
swath1res[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)), 1:ncol(swath1seg)] <- swath1seg
swath1res[((seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)))[idxToUpdate1], ncol(swath1res)] <- 1
swath2res[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)), 1:ncol(swath2seg)] <- swath2seg
swath2res[((seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)))[idxToUpdate2], ncol(swath2res)] <- 1
}
## add a column with an index if it's a bead with multiple intensities
## for now with 3 swath data we are ignoring these beads
#swath1 <- cbind(swath1, rep(0, nrow(swath1)));
#colnames(swath1res)[ncol(swath1res)] <- "Overlap";
#swath2 <- cbind(swath2, rep(0, nrow(swath2)));
#colnames(swath2res)[ncol(swath2res)] <- "Overlap";
#swath3 <- cbind(swath3, rep(0, nrow(swath3)));
#colnames(swath3res)[ncol(swath3res)] <- "Overlap";
## return a list containing the two swaths (Removing the "Swath Index" column since that's implied)
return(list("Swath1" = swath1res[,-(ncol(swath1res)-1)], "Swath2" = swath2res[,-(ncol(swath2res)-1)]));
}
genThreeSwaths <- function(txt, sectionName, inputDir = ".", locsList, twocolour = TRUE, GrnTiffSuffix2 = "-Swath2_Grn.tif", GrnTiffSuffix3 = "-Swath3_Grn.tif", RedTiffSuffix2 = "-Swath2_Red.tif", RedTiffSuffix3 = "-Swath3_Red.tif", segmentHeight = 326, verbose = TRUE) {
## Function to identify beads in the overlapping region between swaths
## For these beads intensities are calculated from both images for comparison
##
## Arguments: txt - 5 column matrix generated by assignToImage()
## sectionName - string containing name of section to read
##
## Output: list containing two 4-column matrices in the same format as bead-level txt files
txt <- txt[order(txt[,"LocsIdx"]),]
s1Width <- nrow(locsList$Grn[[1]]) / (9 * segmentHeight)
s2Width <- nrow(locsList$Grn[[2]]) / (9 * segmentHeight)
swath1 <- txt[which(txt[,ncol(txt)] <= (segmentHeight * s1Width * 9)),]
swath2 <- txt[which(txt[,ncol(txt)] > (segmentHeight * s1Width * 9) & txt[,ncol(txt)] <= ((segmentHeight * s1Width * 9) + (segmentHeight * s2Width * 9))),]
swath3 <- txt[which(txt[,ncol(txt)] > ((segmentHeight * s2Width * 9) + (segmentHeight * s1Width * 9))),]
swath1res <- matrix(ncol = ncol(swath1) + 1, nrow = nrow(swath1))
swath2res <- matrix(ncol = ncol(swath2) + 1, nrow = nrow(swath2))
swath3res <- matrix(ncol = ncol(swath3) + 1, nrow = nrow(swath3))
swath1res[,ncol(swath1res)] <- swath2res[,ncol(swath1res)] <- swath3res[,ncol(swath1res)] <- 0
colnames(swath1res) <- c(colnames(swath1), "Overlap")
colnames(swath2res) <- c(colnames(swath1), "Overlap")
colnames(swath3res) <- c(colnames(swath1), "Overlap")
grid1 <- grid3 <- t(sapply(1:(segmentHeight * s1Width), beadarray:::locsIndicesToGrid, nrow = segmentHeight, ncol = s1Width))
grid2 <- t(sapply(1:(segmentHeight * s2Width), beadarray:::locsIndicesToGrid, nrow = segmentHeight, ncol = s2Width))
tiffG2 <- readTIFF(file.path(inputDir, paste(sectionName, GrnTiffSuffix2, sep = "")))
tiffG3 <- readTIFF(file.path(inputDir, paste(sectionName, GrnTiffSuffix3, sep = "")))
if(twocolour) {
tiffR2 <- readTIFF(file.path(inputDir, paste(sectionName, RedTiffSuffix2, sep = "")))
tiffR3 <- readTIFF(file.path(inputDir, paste(sectionName, RedTiffSuffix3, sep = "")))
}
for(seg in 0:8) {
if(verbose)
message("Overlapping region: Segment ", seg);
swath1seg <- swath1[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)),]
swath2seg <- swath2[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)),]
swath3seg <- swath3[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)),]
decodedColIdx2 <- unlist(lapply(split(swath2seg[,1], grid2[,1]), FUN = function(x) {
if(length(unique(x)) == 1)
return(TRUE)
else
return(FALSE)
}))
decodedColIdx3 <- unlist(lapply(split(swath3seg[,1], grid3[,1]), FUN = function(x) {
if(length(unique(x)) == 1)
return(TRUE)
else
return(FALSE)
}))
idxToUpdate1 <- which(grid1[,1] > s1Width - max(which(decodedColIdx2)))
idxToUpdate2a <- which(grid2[,1] <= max(which(decodedColIdx2)))
swath2seg[idxToUpdate2a,1] <- swath1seg[idxToUpdate1,1]
#swath2seg[idxToUpdate2,2] <- apply(swath2seg[idxToUpdate2,3:4], 1, FUN = singleBeadIntensity_6x6, tiffFile = file.path(paste(sectionName, "-Swath2_Grn.tif", sep = "")))
swath2seg[idxToUpdate2a,2] <- round(illuminaForeground_6x6(tiffG2, swath2seg[idxToUpdate2a,3:4]) - illuminaBackground(tiffG2, swath2seg[idxToUpdate2a,3:4]))
if(twocolour) {
swath2seg[idxToUpdate2a,5] <- round(illuminaForeground_6x6(tiffR2, swath2seg[idxToUpdate2a,6:7]) - illuminaBackground(tiffR2, swath2seg[idxToUpdate2a,6:7]))
}
idxToUpdate2b <- which(grid2[,1] > s2Width - max(which(decodedColIdx3)))
idxToUpdate3 <- which(grid3[,1] <= max(which(decodedColIdx3)))
swath3seg[idxToUpdate3,1] <- swath2seg[idxToUpdate2b,1]
#swath3seg[idxToUpdate3,2] <- apply(swath3seg[idxToUpdate3,3:4], 1, FUN = singleBeadIntensity_6x6, tiffFile = file.path(paste(sectionName, "-Swath3_Grn.tif", sep = "")))
swath3seg[idxToUpdate3,2] <- round(illuminaForeground_6x6(tiffG3, swath3seg[idxToUpdate3,3:4]) - illuminaBackground(tiffG3, swath3seg[idxToUpdate3,3:4]))
if(twocolour) {
swath3seg[idxToUpdate3,5] <- round(illuminaForeground_6x6(tiffR3, swath3seg[idxToUpdate3,6:7]) - illuminaBackground(tiffG3, swath3seg[idxToUpdate3,6:7]))
}
swath1res[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)), 1:ncol(swath1seg)] <- swath1seg
swath1res[((seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)))[idxToUpdate1], ncol(swath1res)] <- 1
swath2res[(seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)), 1:ncol(swath2seg)] <- swath2seg
swath2res[((seg * (segmentHeight * s2Width) + 1):((seg + 1) * (segmentHeight * s2Width)))[c(idxToUpdate2a, idxToUpdate2b)], ncol(swath2res)] <- 1
swath3res[(seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)), 1:ncol(swath3seg)] <- swath3seg
swath3res[((seg * (segmentHeight * s1Width) + 1):((seg + 1) * (segmentHeight * s1Width)))[idxToUpdate3], ncol(swath3res)] <- 1
}
## add a column with an index if it's a bead with multiple intensities
## for now with 3 swath data we are ignoring these beads
#swath1 <- cbind(swath1, rep(0, nrow(swath1)));
#colnames(swath1res)[ncol(swath1res)] <- "Overlap";
#swath2 <- cbind(swath2, rep(0, nrow(swath2)));
#colnames(swath2res)[ncol(swath2res)] <- "Overlap";
#swath3 <- cbind(swath3, rep(0, nrow(swath3)));
#colnames(swath3res)[ncol(swath3res)] <- "Overlap";
## return a list containing the two swaths (Removing the "Swath Index" column since that's implied)
return(list("Swath1" = swath1res[,-(ncol(swath1res)-1)], "Swath2" = swath2res[,-(ncol(swath2res)-1)], "Swath3" = swath3res[,-(ncol(swath3res)-1)]));
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.