imageHTS-introduction.R
In imageHTS: Analysis of high-throughput microscopy-based screens

### R code from vignette source 'imageHTS-introduction.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: init
###################################################
display = function(...) {invisible()}


###################################################
### code chunk number 3: library
###################################################
library('imageHTS')


###################################################
### code chunk number 4: parseImageConf
###################################################
localPath = tempdir()
serverURL = system.file('submorph', package='imageHTS')
x = parseImageConf('conf/imageconf.txt', localPath=localPath,
                   serverURL=serverURL)


###################################################
### code chunk number 5: configure
###################################################
x = configure(x, 'conf/description.txt', 'conf/plateconf.txt',
                 'conf/screenlog.txt')
x = annotate(x, 'conf/annotation.txt')


###################################################
### code chunk number 6: getUnames
###################################################
unames = setdiff(getUnames(x), getUnames(x, content='empty'))
unames 


###################################################
### code chunk number 7: getWellFeatures
###################################################
getWellFeatures(x, unames[1:3])


###################################################
### code chunk number 8: segmentWells1
###################################################
uname = getUnames(x, content='rluc')[3]
print(uname)
y =  segmentWells(x, uname=uname,
                     segmentationPar='conf/segmentationpar.txt')
display(y$cal)
hseg = highlightSegmentation(0.6*y$cal, y$nseg, y$cseg, thick=TRUE)
display(hseg)


###################################################
### code chunk number 9: segmentWells1h
###################################################
writeImage(y$cal, 'cal.jpeg', quality=80)
writeImage(hseg, 'hseg.jpeg', quality=80)


###################################################
### code chunk number 10: segmentWells2 (eval = FALSE)
###################################################
## unames = setdiff(getUnames(x), getUnames(x, content='empty'))
## segmentWells(x, unames, 'conf/segmentationpar.txt')


###################################################
### code chunk number 11: fileHTS
###################################################
fileHTS(x, type='source', uname=uname, channel=1)
fileHTS(x, type='seg', uname=uname)
fileHTS(x, type='viewfull', uname=uname)        


###################################################
### code chunk number 12: extractFeatures
###################################################
extractFeatures(x, uname, 'conf/featurepar.txt')


###################################################
### code chunk number 13: readFeatures
###################################################
y = readHTS(x, type='ftrs', uname=uname, format='tab')
dim(y)
y[1:10, 1:7]


###################################################
### code chunk number 14: histfig
###################################################
ctub <- y$c.t.b.mean*y$c.s.area
hist(ctub, 20, xlab='Cell tubulin intensity (a.u.)', main='')
abline(v=1600, col=2)
cellid = which(ctub>1600)
print(cellid)


###################################################
### code chunk number 15: traceCell
###################################################
cal = readHTS(x, type='cal', uname=uname, format='rda')
seg = readHTS(x, type='seg', uname=uname, format='rda')
cseg = rmObjects(seg$cseg, setdiff(1:nrow(y), cellid))
hightub = highlightSegmentation(0.6*cal, cseg=cseg, thick=TRUE)
display(hightub)


###################################################
### code chunk number 16: traceCellh
###################################################
writeImage(hightub, 'hightub.jpeg', quality=80)


###################################################
### code chunk number 17: readLearnTS (eval = FALSE)
###################################################
## set.seed(1)
## readLearnTS(x, 'conf/featurepar.txt', 'conf/trainingset.txt')


###################################################
### code chunk number 18: predictCellLabels
###################################################
predictCellLabels(x, uname)


###################################################
### code chunk number 19: predictCellLabelsOutput
###################################################
cellLabels = capture.output(predictCellLabels(x, uname))
cellLabels = unlist(strsplit(cellLabels,"[= ]")) 
cellNumber = function(x) {
  as.integer(cellLabels[which(cellLabels==x)+1])
}


###################################################
### code chunk number 20: drawCellLabels
###################################################
clab = readHTS(x, type='clabels', uname=uname, format='tab')
labid = split(1:nrow(clab), clab$label)
inter = seg$cseg%in%labid$I
mito = seg$cseg%in%labid$M
debris = seg$cseg%in%labid$D
dc = Image(c(inter+mito, inter, debris+inter), colormode='Color',
           dim=c(dim(seg$cseg)[1:2], 3))
dc = highlightSegmentation(0.5*dc+0.2*drop(cal), cseg=seg$cseg,
                           thick=TRUE)
display(dc)


###################################################
### code chunk number 21: drawCellLabelsh
###################################################
writeImage(dc, 'dc.jpeg', quality=80)


###################################################
### code chunk number 22: summarizeWells (eval = FALSE)
###################################################
## summarizeWells(x, unames, 'conf/featurepar.txt')


###################################################
### code chunk number 23: sum
###################################################
profiles = readHTS(x, type='file', filename='data/profiles.tab',
                   format='tab')
wfcontent = 
  factor(as.character(getWellFeatures(x, unames)$controlStatus))
table(wfcontent)
zwf = split(1:nrow(profiles), wfcontent)
ft = c('n', 'med.c.s.area', 'med.c.t.b.mean', 'M', 'D')
avef = do.call(rbind,
               lapply(zwf, function(z) colMeans(profiles[z, ft])))
print(avef)


###################################################
### code chunk number 24: imageconf
###################################################
f = fileHTS(x, 'file', filename='conf/imageconf.txt')
cat(paste(readLines(f), collapse='\n'), '\n')


###################################################
### code chunk number 25: segmentationpar
###################################################
f = fileHTS(x, 'file', filename='conf/segmentationpar.txt')
cat(paste(readLines(f), collapse='\n'), '\n')


###################################################
### code chunk number 26: featurepar
###################################################
f = fileHTS(x, 'file', filename='conf/featurepar.txt')
cat(paste(readLines(f), collapse='\n'), '\n')


###################################################
### code chunk number 27: script (eval = FALSE)
###################################################
## library('imageHTS')
## localPath = tempdir()
## serverURL = system.file('submorph', package='imageHTS')
## x = parseImageConf('conf/imageconf.txt', localPath=localPath,
##                    serverURL=serverURL)
## x = configure(x, 'conf/description.txt', 'conf/plateconf.txt',
##                  'conf/screenlog.txt')
## x = annotate(x, 'conf/annotation.txt')
## unames = setdiff(getUnames(x), getUnames(x, content='empty'))
## segmentWells(x, unames, 'conf/segmentationpar.txt')
## extractFeatures(x, unames, 'conf/featurepar.txt')
## readLearnTS(x, 'conf/featurepar.txt', 'conf/trainingset.txt')
## predictCellLabels(x, unames)
## summarizeWells(x, unames, 'conf/featurepar.txt')


###################################################
### code chunk number 28: init
###################################################
localPath = file.path(tempdir(), 'kimorph')
serverURL = 'http://www.huber.embl.de/cellmorph/kimorph/'
x = parseImageConf('conf/imageconf.txt', localPath=localPath,
                   serverURL=serverURL)
x = configure(x, 'conf/description.txt', 'conf/plateconf.txt',
                 'conf/screenlog.txt')
x = annotate(x, 'conf/annotation.txt')


###################################################
### code chunk number 29: get
###################################################
us = setdiff(getUnames(x), getUnames(x, content='empty'))
wfcontent = getWellFeatures(x, us)$controlStatus
table(wfcontent)
xd = readHTS(x, 'file', filename='data/profiles.tab', format='tab')
xd = xd[match(us, xd$uname),]


###################################################
### code chunk number 30: boxplots
###################################################
colors = c('#ffffff', NA, '#aaffff', '#ffaaff', '#ff44aa', '#aaaaff',
           '#aaffaa', '#ff7777', '#aaaaaa', '#ffff77')
par(mfrow=c(1,2))
boxplot(xd$med.c.g.ss~wfcontent, las=2, col=colors,
        main='Median cell size (a.u.)')
boxplot(xd$med.c.g.ec~wfcontent, las=2, col=colors,
        main='Median cell eccentricity (a.u.)')


###################################################
### code chunk number 31: lda1
###################################################
library("MASS")
z = wfcontent %in% c('rluc', 'ubc', 'trappc3')
ft = 14:50
ld = lda(xd[z, ft], as.character(wfcontent[z]))
py = predict(ld, xd[, ft])
plot(py$x[,1:2])


###################################################
### code chunk number 32: standout
###################################################
unames = us[which(py$x[,1]>500)]
print(unames)
f = fileHTS(x, type='viewunmonted', spot=3, uname=unames[1])
img1 = readImage(f)[1791:2238,1:448,]
display(img1)
f = fileHTS(x, type='viewunmonted', spot=1, uname=unames[2])
img2 = readImage(f)[1:448,1:448,]
display(img2)


###################################################
### code chunk number 33: standouth
###################################################
writeImage(img1, 'img1.jpeg', quality=80)
writeImage(img2, 'img2.jpeg', quality=80)


###################################################
### code chunk number 34: lda2
###################################################
plot(py$x[,1:2], xlim=c(-35,25), ylim=c(-20,20), cex=0.3)
z = wfcontent!='sample'
points(py$x[z,1:2], col=1, bg=colors[wfcontent[z]], pch=21)
col = rep(1, length(levels(wfcontent)))
col[2] = NA
legend('topleft', legend=levels(wfcontent), col=col,
       pt.bg=colors[1:length(wfcontent)], pch=21, ncol=2, cex=0.8) 


###################################################
### code chunk number 35: popCellPicker (eval = FALSE)
###################################################
## popWebQuery(x)
## uname = getUnames(x, content='trappc3')[1]
## popCellPicker(x, uname)


###################################################
### code chunk number 36: sessionInfo
###################################################
toLatex(sessionInfo())