# R/cytoFunctions.R In prada: Data analysis for cell-based functional assays

#### Documented in gatePointsthreePanelPlot

```gatePoints <- function(obj,totmin=0,totmax=1023,gatecol="red",smooth=FALSE){

##  check arguments:
if (!is.matrix(obj)) stop("First argument is no matrix!\n")
if (ncol(obj)!=2) stop("Data matrix does not have two columns!\n")
if (is.null(colnames(obj))) colnames(obj) <- c("x","y")
options(locatorBell=FALSE)

##  auxiliary functions:
myLocator <- function(mymin=0,mymax=1023){
mypoint <- locator(1)
if (!is.null(mypoint)){
if (mypoint\$x < mymin) mypoint\$x <- mymin
if (mypoint\$x > mymax) mypoint\$x <- mymax
if (mypoint\$y < mymin) mypoint\$y <- mymin
if (mypoint\$y > mymax) mypoint\$y <- mymax
} #if (!is.null(mypoint))
return(mypoint)
} #myLocator

drawLine <- function(point1,point2,gatecol="red"){
x1 <- round(point1\$x)
x2 <- round(point2\$x)
y1 <- round(point1\$y)
y2 <- round(point2\$y)
coord <- list(x=c(x1,x2),y=c(y1,y2))
lines(coord,lty=1,lwd=2,col=gatecol)
} #drawLine

insidePolygon <- function(data, polygon){
# determine if a point 'mypoint' is inside a polygon given by its vertices:
# if number of intersections between a horizontal ray emanating from
# mypoint and edges of the polyong is odd -> mypoint inside polygon
vertices <- matrix(as.numeric(unlist(polygon)), ncol=2, byrow=TRUE)
data <- cbind(as.numeric(data[,1]), as.numeric(data[,2]))
res <- .Call("inPolygon", data, vertices,
return(as.logical(res))
} # insidePolygon

if (is.null(totmin)) totmin <- round(min(obj))
if (is.null(totmax)) totmax <- round(max(obj))
mylimits <- c(totmin,totmax)

# start plot
if (smooth)
smoothScatter(obj, xlab=colnames(obj)[1], ylab=colnames(obj)[2],
xlim=mylimits, ylim=mylimits)
else
plot(obj,pch=20, cex=0.5, xlab=colnames(obj)[1], ylab=colnames(obj)[2],
xlim=mylimits, ylim=mylimits, col=densCols(obj))

polyVertices <- c()

# get input
cat("Draw gate (left click: set point, right click: end drawing).\n")
startpoint <- myLocator(totmin,totmax)
polyVertices <- c(polyVertices,list(c(startpoint\$x,startpoint\$y)))
temppoint <- "bla"
lastpoint <- startpoint
while (!is.null(temppoint)){
temppoint <- myLocator(totmin,totmax)
if (!is.null(temppoint)){
polyVertices <- c(polyVertices,list(c(temppoint\$x,temppoint\$y)))
drawLine(lastpoint,temppoint,gatecol)
lastpoint <- temppoint
}
} # while
drawLine(lastpoint,startpoint,gatecol)
polyVertices <- c(polyVertices,list(c(startpoint\$x,startpoint\$y)))

# process drawn box
cat("Determining events within gate...")
#eturn(polyVertices)
dataInBox <- insidePolygon(obj,polyVertices)
if (smooth)
points(obj[which(dataInBox),,drop=FALSE],pch=20,cex=0.5,col=gatecol)
else
points(obj[which(dataInBox),,drop=FALSE],pch=1,cex=0.7,col=gatecol)
return(dataInBox)
} # gatePoints

gate.matrix <- function(object,gate.colour="red",use.smoothScatter=FALSE,
data.min=0,data.max=1023,max.observations=20000)
{
### 0. check arguments ###
if (!is.matrix(object))
stop("Function 'gate.matrix' can only be applied to matrices!\n")
n.obs <- nrow(object)
if (n.obs > max.observations){
object <- object[1:max.observations,,drop=FALSE]
warning(paste("gate.matrix: matrix contained more than",max.observations,
"observations! Only first", max.observations, "observations",
"were used."))
} #  if (n.obs > max.observations)
if (is.null(colnames(object)))
colnames(object) <- paste("Variable",1:ncol(object),sep="")

### 1. define auxiliary functions ###
getGateVariables <- function(object){
varnames <- colnames(object)
selectedVar1 = selectedVar2 <- NA
if ((.Platform\$OS.type=="windows") && interactive()){
cat("Select first variable for gating...")
while (selectedVar1 %in% c(NA,""))
selectedVar1 <- select.list(varnames,multiple=FALSE)
cat(selectedVar1,"\n")
cat("Select second variable for gating...")
while (selectedVar2 %in% c(NA,""))
selectedVar2 <- select.list(varnames,multiple=FALSE)
cat(selectedVar2,"\n")
selectedVars <- c(selectedVar1,selectedVar2)
} else { # Unix
cat("\nChoose two of the following variables for gating:\n")
print(varnames)
while(is.na(selectedVar1)){
if (is.na(selectedVar1))
} #while
while(is.na(selectedVar2)){
if (is.na(selectedVar2))
} #while
selectedVars <- c(selectedVar1,selectedVar2)
cat("Selected Variables:")
print(selectedVars)
} # else Unix
return(selectedVars)
} #getGateVariables

### 2. initialize result ###
inGate <- rep(FALSE,nrow(object))# initialize
keepRows <- 1:nrow(object)

while (length(grep("^[fF]",userAnswer))==0){ # if not finished
combRows <- numeric()
keepobject <- object[keepRows,,drop=FALSE]
selectedVars <- getGateVariables(object)
gate1 <- gatePoints(keepobject[,selectedVars,drop=FALSE],totmin=data.min,totmax=data.max,gatecol=gate.colour,smooth=use.smoothScatter)
cat("found",sum(gate1),"\n")
userAnswer <- readline("(R)edo this gating, (C)ombine this gating with another one, (P)roceed or (F)inish? ")
}#while
if (sum(gate1)==0) cat("No events in drawn gate!\n")
keepRows <- keepRows[gate1]
cat("Keeping only events within previous gate...\n")
} else if (length(grep("^[cC]",userAnswer))!=0){ # combination of gates?
combRows <-  keepRows[gate1]
selectedVars <- getGateVariables(object)
gate2 <- gatePoints(keepobject[,selectedVars,drop=FALSE],totmin=data.min,totmax=data.max,gatecol=gate.colour,smooth=use.smoothScatter)
cat("found",sum(gate2),"\n")
userAnswer <- readline("(R)edo this gating, (A)dd another gating to combination, (U)se this combination? ")
combRows <-  c(combRows,keepRows[gate2])
keepRows <- unique(combRows)
cat("Total of",length(keepRows),"events within this gate combination.\n")
} # else combination of gates
}#while
inGate[keepRows] <- TRUE
inPercent <- 100 * sum(inGate)/length(inGate)
cat("\nWithin last gate: ",inPercent,"% of all",length(inGate),"events.\n\n")
class(inGate) <- "matrixGate"
return(inGate)
} # gate.matrix

gate.cytoFrame <- function(object,gate.colour="red",use.smoothScatter=FALSE,
data.min=0,data.max=1023,max.observations=20000)
{
### 0. check arguments ###
if (class(object)!="cytoFrame")
stop("Function 'gate.cytoFrame' can only be applied to objects of class 'cytoFrame'!\n")

### 1. call 'gate.matrix'
result <- gate.matrix(exprs(object),gate.colour=gate.colour,
use.smoothScatter=use.smoothScatter,
data.min=data.min,data.max=data.max,
max.observations=max.observations)
class(result) <- "cytoFrameGate"

return(result)
} # gate.cytoFrame

gate.cytoSet <- function(object,gate.colour="red",use.smoothScatter=FALSE,
data.min=0,data.max=1023,max.observations=20000)
{
### 0. check arguments ###
if (class(object)!="cytoSet")
stop("Function 'gate.cytoSet' can only be applied to objects of class 'cytoSet'!\n")

### 1. prepare result ###
n.cytoFrames <- length(object)
result <- as.list(rep("cytoFrame not gated",n.cytoFrames))
names(result) <- phenoData(object)\$name
if (is.null(names(result)))
names(result) <- paste("frame",1:n.cytoFrames,sep="")
class(result) <- "cytoSetGate"
on.exit(return(result))

## CONTINUE HERE! ##
### 2. call 'gate.matrix'
for (i in 1:n.cytoFrames){
cat("Working on cytoFrame '",names(result)[i],"'...\n",sep="")
result[[i]] <- gate.cytoFrame(object[[i]],gate.colour=gate.colour,
use.smoothScatter=use.smoothScatter,
data.min=data.min,data.max=data.max,
max.observations=max.observations)
} # for (i in 1:n.cytoFrames)
#return(result)

} # gate.cytoSet

plot.matrixGate <- function(x,y,x.panels=c(1,4,5),y.panels=c(2,3,6),
use.smoothScatter=TRUE, limits=c(0,1023),
plotTitle="Gated Data",gate.colour="red",...)
{
### assume x to be the gate and y to be the data used for gating
stopifnot(is.matrix(y),class(x)=="matrixGate",nrow(y)==length(x))
threePanelPlot(y,x.panels=x.panels,y.panels=y.panels,
use.smoothScatter=use.smoothScatter,limits=limits,
# since x is a logical vector, we can use it to index the data matrix
invisible(NULL)
} # plot.matrixGate

plot.cytoFrameGate <- function(x,y,x.panels=c(1,4,5),y.panels=c(2,3,6),
use.smoothScatter=TRUE, limits=c(0,1023),
plotTitle="Gated Data",gate.colour="red",
new.device=TRUE,verbose=TRUE,...)
{
### assume x to be the gate and y to be the data used for gating
stopifnot(class(y)=="cytoFrame",class(x)=="cytoFrameGate")
data <- exprs(y)[1:length(x),,drop=FALSE]
# in the likely case, that observations were dropped from the cytoFrame
#   for gating, drop them here as well

threePanelPlot(data,x.panels=x.panels,y.panels=y.panels,
use.smoothScatter=use.smoothScatter,limits=limits,
plotTitle=plotTitle,
new.device=new.device,verbose=verbose,...)
# since x is a logical vector, we can use it to index the data matrix
invisible(NULL)
} # plot.cytoFrameGate

plot.cytoSetGate <- function(x,y,x.panels=c(1,4,5),y.panels=c(2,3,6),
use.smoothScatter=TRUE, limits=c(0,1023),
plotTitle="Gated Data",gate.colour="red",
verbose=TRUE,...)
{
### assume x to be the gate and y to be the data used for gating
stopifnot(class(y)=="cytoSet",class(x)=="cytoSetGate")

n.cytoFrames <- length(y)
frame.names  <- phenoData(y)\$name
if (is.null(frame.names))
frame.names <- paste("frame",1:n.cytoFrames,sep="")

### 1. prepare plot
par(mfrow=c(n.cytoFrames,3))

### 2. call 'plot.cytoFrameGate'
for (i in 1:n.cytoFrames){
cat("Working on cytoFrame '",frame.names[i],"'...\n",sep="")

thisFrame <- y[[i]]
thisFrameGate <- x[[i]]
class(thisFrameGate) <- "cytoFrameGate"
thisPlotTitle <- paste("Gated Frame",abbreviate(frame.names[i],15))

plot.cytoFrameGate(thisFrameGate,thisFrame,
x.panels=x.panels,y.panels=y.panels,
gate.colour=gate.colour, plotTitle=thisPlotTitle,
use.smoothScatter=use.smoothScatter,
limits = limits, new.device=FALSE, verbose=verbose,...)

} # for (i in 1:n.cytoFrames)

invisible(NULL)
} # plot.cytoFrameGate

threePanelPlot <- function(data,x.panels=c(1,4,5),y.panels=c(2,3,6),
tot.width=15,tot.height=5.4,maxcells=20000,
limits=c(0,1023),remove.extremes=TRUE,
plotTitle="Three-Panel Plot",use.smoothScatter=TRUE,
palette=colorRampPalette(brewer.pal(9,"Blues")),
new.device=TRUE,verbose=TRUE,
{
# check arguments
if (verbose) cat("Checking arguments...\n")
stopifnot(all(c(x.panels,y.panels) %in% 1:ncol(data)))

# prepare data:
if (verbose) cat("Preparing data...\n")
if (remove.extremes){
data.rowmins <- apply(data,1,min)
data.rowmaxs <- apply(data,1,max)
data <- data[(data.rowmins>limits[1])&(data.rowmaxs<limits[2]),,drop=FALSE]
} # if (remove.extremes)
ncells <- nrow(data)
sampled.events <- sample(1:ncells,min(ncells,maxcells),replace=FALSE)
data <- data[sampled.events,,drop=FALSE]

# prepare plot
if (verbose) cat("Plotting...\n")
if (new.device){
do.call(getOption("device"),list(width=tot.width,height=tot.height))
par(mfrow=c(1,3),font.lab=2,mar=c(5,4,5,4))
} # if (new.device)
for (i in 1:3){
if (use.smoothScatter)
smoothScatter(data[,c(x.panels[i],y.panels[i],drop=FALSE)],
colramp=palette,xlim=limits,ylim=limits,...)
else
plot(data[,c(x.panels[i],y.panels[i],drop=FALSE)],xlim=limits,
ylim=limits,...)
if (all(i==2,all(!is.na(plotTitle)),!is.null(plotTitle)))
title(main=plotTitle,cex.main=1.5,font.main=2)