Last active
November 20, 2018 17:40
-
-
Save grimbough/09d150bff2b07ec032de4d47cf201215 to your computer and use it in GitHub Desktop.
Modified example of code provided in https://support.bioconductor.org/p/115224/ to identify error thrown by EBImage
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
WD <- tempdir() | |
setwd(WD) | |
library(EBImage) | |
library(tidyverse) | |
## download example zip file ad unpack | |
download.file('https://www.dropbox.com/s/40e5tfohgdm9wlo/IF015_plate002_test.zip?dl=1', | |
destfile = 'IF015_plate002.zip', mode = "wb") | |
unzip('IF015_plate002.zip') | |
# Parameters | |
IF = "IF015_noSpeck" | |
# List of files######################################################## | |
AllImages <- as.data.frame(list.files(getwd(), recursive = T, pattern = '.tif$'), col.names = 'FileNames') | |
colnames(AllImages)[1] <- 'FileNames' | |
# Extract Metainformation############################################################################## | |
MetaInformation <- function(Files, Path){ | |
Files$Directory <- Path | |
Files$Plate <- substring(Files$FileNames,12,14) ## Plate | |
Files$Plate <- as.numeric(Files$Plate) | |
Files$Stack <- sub(".*test/", "", Files$FileNames) ## Stack | |
Files$DataStack <- substring(Files$Stack,0,35) ## DataStack | |
Files$ID <- sub("--W.*", "", Files$Stack) ## ID | |
Files$Well <- sub(".*--W", "", Files$Stack) ## W # | |
Files$Well <- substring(Files$Well,0,5) ## W # | |
Files$Position <- sub(".*--P", "", Files$Stack) ## Position | |
Files$Position <- substring(Files$Position,0,5) ## Position | |
Files$Time <- sub(".*--T", "", Files$Stack) ## Time | |
Files$Time <- substring(Files$Time,0,5) ## Time | |
Files$Channel <- sub(".*--", "", Files$Stack) ## Channel | |
Files$Channel <- sub(".tif.*", "", Files$Channel) ## Channel | |
return(as.data.frame(Files)) | |
} | |
# Add metainformation (if necessary, more revelant for data analysis)############ | |
IdToMap <- function(FileList){ | |
rN = c(1:16) | |
cN = c(1:24) | |
cID = c(1:24) | |
rID = c('A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P') | |
Row <- data.frame(rID, rN) | |
ID = data.frame(cID, cN) | |
ID$rID <- NA | |
for(i in rID){ | |
IDtemp <- ID[1:24,] | |
IDtemp$rID <- as.character(i) | |
ID <- rbind(ID, IDtemp) | |
} | |
ID <- merge(ID[25:nrow(ID),], Row, by = "rID", all.x = T) | |
ID$ID <- paste(ID$rID, ID$cID, sep = "") | |
FileList <- merge(x = FileList, y = ID, by = c("ID"), all.x = T) | |
} | |
AllImages <- MetaInformation(Files = AllImages, Path = WD) | |
AllImages <- IdToMap(AllImages) | |
unique(AllImages$Plate) | |
str(AllImages) | |
# Read images and extract features##################### | |
ImageProcessing <- function(Image, Plate, Well, Position){ | |
#---------------Rescale pixel intensity---------- | |
Bottom = 2^15/(2^16-1) | |
Top = (2^15+4095)/(2^16-1) | |
Rescale <- function(x){ | |
(x - Bottom) / (Top - Bottom) | |
} | |
Boolean <- AllImages$Plate == Plate & AllImages$ID == Well & AllImages$Position == Position | |
#---------------Load Images------------------------ | |
Img405 = readImage(paste0(AllImages[AllImages$Channel == "DAPI" & Boolean, "Directory"], | |
"/", | |
AllImages[AllImages$Channel == "DAPI" & Boolean, "FileNames"])) | |
Img405 = Rescale(Img405) | |
NImg405 = EBImage::normalize(Img405, inputRange = c(range(Img405)[1], range(Img405)[2])) | |
Img488 = readImage(paste0(AllImages[AllImages$Channel == "GFP" & Boolean, "Directory"], | |
"/", | |
AllImages[AllImages$Channel == "GFP" & Boolean, "FileNames"])) | |
Img488 = Rescale(Img488) | |
NImg488 = EBImage::normalize(Img488) | |
Img568 = readImage(paste0(AllImages[AllImages$Channel == "mCherry" & Boolean, "Directory"], | |
"/", | |
AllImages[AllImages$Channel == "mCherry" & Boolean, "FileNames"])) | |
Img568 = Rescale(Img568) | |
NImg568 = EBImage::normalize(Img568, inputRange = c(range(Img568)[1], range(Img568)[2])) | |
Img647 = readImage(paste0(AllImages[AllImages$Channel == "647" & Boolean, "Directory"], | |
"/", | |
AllImages[AllImages$Channel == "647" & Boolean, "FileNames"])) | |
Img647 = Rescale(Img647) | |
NImg647 = EBImage::normalize(Img647, inputRange = c(range(Img647)[1], range(Img647)[2])) | |
#---------------smooth and threshold nuleus------------------------------ | |
FilterNuc = makeBrush(size = 51, shape = "gaussian", sigma=2) | |
Img405smooth = filter2(Img405, filter = FilterNuc) | |
nucThrManual = thresh(Img405smooth, w = 100, h = 100, offset = 0.0005) | |
nucOpening = opening(nucThrManual, kern = makeBrush(15, shape="disc")) | |
nucSeed = bwlabel(nucOpening) | |
nucFill = fillHull(thresh(Img405smooth, w = 20, h = 20, offset = 0.01)) | |
nucRegions = propagate(Img405smooth, nucSeed, mask = nucFill) | |
nucMask = watershed(distmap(nucOpening), tolerance = 1, ext = 1) | |
#------------------Generate voronoi------------------------------------- | |
VoR = propagate(nucMask, nucMask, lambda = 100000) | |
#------------------Generate donut and bubble for proxy Cytoplasm ------- | |
Bubble = dilate(nucRegions, kern = makeBrush(15, shape = 'disc')) | |
display(colorLabels(Bubble)) | |
BubbleBound = propagate(Bubble, nucMask, Bubble, lambda = 1000) | |
Donut = BubbleBound - nucMask | |
#---------------Features extraction------------------------------------------------------------------------- | |
F405_Nuc = computeFeatures(nucMask, Img405, xname = "Nuc", refnames = "405", haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F488_Nuc = computeFeatures(nucMask, Img488, xname = "Nuc", refnames = "488", haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F488_Donut = computeFeatures(Donut, Img488, xname = "Donut", refnames = "488", haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F488_Bubble = computeFeatures(BubbleBound, Img488, xname = "Bubble", refnames = "488", haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F568_Nuc = computeFeatures(nucMask, Img568, xname = "Nuc", refnames = "568",haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F568_Donut = computeFeatures(Donut, Img568, xname = "Donut", refnames = "568",haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F568_Bubble = computeFeatures(BubbleBound, Img568, xname = "Donut", refnames = "568",haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F647_Nuc = computeFeatures(nucMask, Img647, xname = "Nuc", refnames = "647",haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F647_Donut = computeFeatures(Donut, Img647, xname = "Donut", refnames = "647",haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
F647_Bubble = computeFeatures(BubbleBound, Img647, xname = "Donut", refnames = "647",haralick.nbins = 32, haralick.scales = c(1,2,4,8)) | |
Fc = cbind(F405_Nuc, | |
F488_Nuc, | |
F488_Donut, | |
F488_Bubble, | |
F568_Nuc, | |
F568_Donut, | |
F568_Bubble, | |
F647_Nuc, | |
F647_Donut, | |
F647_Bubble) | |
Fc = as.data.frame(Fc) %>% rownames_to_column | |
#---------------Remove cells at the edge------------------------------------------------------------------------------- | |
#subset for boundary pixels | |
dims = dim(nucMask) | |
border = c(nucMask[1:dims[1],1], | |
nucMask[1:dims[1],dims[2]], | |
nucMask[1,1:dims[2]], | |
nucMask[dims[1],1:dims[2]] | |
) | |
#extract object identifiers at the boundary | |
ids = unique(border[which(border != 0)]) | |
Fc = filter(Fc, !rowname %in% ids) | |
if (nrow(Fc) == 0) { | |
Fc = data.frame(0, nrow = 0, ncol = 141*12+1) | |
} | |
write.csv(Fc, paste0(IF,"_p--",Plate,"_w--",Well,"_pos--",Position,".csv")) | |
} | |
#----------------Run ImageProcessing() to everyfield--------------------------------------------------------------------- | |
for (Plate in unique(AllImages$Plate)) { | |
for(Well in unique(AllImages$ID)) { | |
for(Position in unique(AllImages$Position)) { | |
ImageProcessing(Image = AllImages, Plate = Plate, Well = Well, Position = Position) | |
print(paste0("Plate: ", Plate, "; Well: ", Well, "; Position: ", Position, "; Time:", Sys.time())) | |
} | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment