Ellipsoid calculator function

This is an R function, written to calculate 95% density ellipsoids from the Purkinje cell endpoint distribution, and compare them beween different genotypes. The ‘Coordinates.list’ input is a named list of the finalised X and Z coorinates of Purkinje cell endpoints, after the cell planarity alignment.

Ellipsoid.calculator<-function(Coordinates.list,
                               Genotypes=c("WT","HET","HO"),
                               Output.dir="",
                               Output.prefix="",
                               Colours=c(Ellipse="red",Insidedot="grey",Outsidedot="blue")){
  library("car")
  #https://stackoverflow.com/questions/30825387/getting-the-parameters-of-a-data-ellipse-produced-by-the-car-package-in-r
  #https://math.stackexchange.com/questions/76457/check-if-a-point-is-within-an-ellipse
  ell.info <- cov.wt(Coordinates.list[[Genotypes[1]]])
  eigen.info <- eigen(ell.info$cov)
  axis.lengths <- sqrt(eigen.info$values * 2 * qf(.95, 2, length(Coordinates.list[[Genotypes[1]]]$FinalXcoordinates)-1))
  ell.origin<-ell.info$center
  ell.vertex1<-ell.info$center + axis.lengths[1] * eigen.info$vectors[,1]
  ell.vertex2<-ell.info$center - axis.lengths[1] * eigen.info$vectors[,1]
  ell.covertex1<-ell.info$center + axis.lengths[2] * eigen.info$vectors[,2]
  ell.covertex2<-ell.info$center - axis.lengths[2] * eigen.info$vectors[,2]
  ell.rx<-sqrt((ell.origin[1] - ell.vertex1[1])^2 + (ell.origin[2] - ell.vertex1[2])^2)
  ell.ry<-sqrt((ell.origin[1] - ell.covertex1[1])^2 + (ell.origin[2] - ell.covertex1[2])^2)
  
  #https://stats.stackexchange.com/questions/9898/how-to-plot-an-ellipse-from-eigenvalues-and-eigenvectors-in-r
  alpha <- atan(eigen.info$vectors[ , 1][2] / eigen.info$vectors[ , 1][1])
  theta <- seq(0, 2 * pi, length=(1000))
  x0 <- ell.origin[1]
  y0 <- ell.origin[2]
  x <- x0 + sqrt(eigen.info$values * 2 * qf(.95, 2, length(Coordinates.list[[Genotypes[1]]]$FinalXcoordinates)-1))[1] *
    cos(theta) * cos(alpha) - sqrt(eigen.info$values * 2 * qf(.95, 2, length(
      Coordinates.list[[Genotypes[1]]]$FinalXcoordinates)-1))[2] * sin(theta) * sin(alpha)
  y <- y0 + sqrt(eigen.info$values * 2 * qf(.95, 2, length(Coordinates.list[[Genotypes[1]]]$FinalXcoordinates)-1))[1] *
    cos(theta) * sin(alpha) + sqrt(eigen.info$values * 2 * qf(.95, 2, length(
      Coordinates.list[[Genotypes[1]]]$FinalXcoordinates)-1))[2] * sin(theta) * cos(alpha)
  
  #https://stackoverflow.com/questions/7946187/point-and-ellipse-rotated-position-test-algorithm
  #(((cos(alpha)*(x-h)+sin(alpha)*(y-k))^2)/(rx^2))+(((sin(alpha)*(x-h)+cos(alpha)*(y-k))^2)/(ry^2))<=1
 pdf(file = paste0(Output.dir,"/Overall_endpoints_vs_95pc_WT_ellipse_",Output.prefix,".pdf"),width = 9) 
  for(i in 1:length(Genotypes)){
  plot(Coordinates.list[[Genotypes[i]]][apply(X = Coordinates.list[[Genotypes[i]]],
                                              MARGIN = 1,
                                              FUN = function(x){
                                                (((cos(-alpha)*(x[1]-ell.origin[1])+sin(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.rx^2))+(
                                                  ((sin(-alpha)*(x[1]-ell.origin[1])+cos(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.ry^2))<=1
                                              }),],
       pch=20,
       cex=0.5,
       col=rgb(t(col2rgb(Colours["Insidedot"])),alpha = 150,maxColorValue = 255),
       ylim=c(-50,50),
       xlim=c(-125,125),
       xlab="X",
       ylab="Z",
       main=paste0("Sample:",
                   names(Coordinates.list[Genotypes[i]]),"\n",
                   round(x = (sum(!apply(X = Coordinates.list[[Genotypes[i]]],
                                         MARGIN = 1,
                                         FUN = function(x){
                                           (((cos(-alpha)*(x[1]-ell.origin[1])+sin(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.rx^2))+(
                                             ((sin(-alpha)*(x[1]-ell.origin[1])+cos(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.ry^2))<=1
                                         }))*100)/nrow(Coordinates.list[[Genotypes[i]]]),digits = 2),
                   "% of points are outside of WT 95percentile ellipse"))
  points(Coordinates.list[[Genotypes[i]]][!apply(X = Coordinates.list[[Genotypes[i]]],
                                                 MARGIN = 1,
                                                 FUN = function(x){
                                                   (((cos(-alpha)*(x[1]-ell.origin[1])+sin(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.rx^2))+(
                                                     ((sin(-alpha)*(x[1]-ell.origin[1])+cos(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.ry^2))<=1
                                                 }),],
         pch=20,
         cex=0.5,
         col=rgb(t(col2rgb(Colours["Outsidedot"])),alpha = 150,maxColorValue = 255))
  lines(x, y, type = "l",lwd=2,col=Colours["Ellipse"])
}
  dev.off()
  Cell_by_cell.Values<-list()
  for(i in 1:length(Genotypes)){
    Cell_by_cell.Values[[Genotypes[i]]]<-unlist(lapply(X = split(x = Coordinates.list[[Genotypes[i]]],
                                                 f = unlist(lapply(X = strsplit(x = rownames(Coordinates.list[[Genotypes[i]]]),split = "\\."),
                                                                   FUN = function(x){x[1]}))),
                                       FUN = function(q){
                                         round(x = (sum(!apply(X = q,
                                                               MARGIN = 1,
                                                               FUN = function(x){
                                  (((cos(-alpha)*(x[1]-ell.origin[1])+sin(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.rx^2))+(
                                    ((sin(-alpha)*(x[1]-ell.origin[1])+cos(-alpha)*(x[2]-ell.origin[2]))^2)/(ell.ry^2))<=1
                                                                 }))*100)/nrow(q),digits = 2)
                                         }))
  }
  pdf(file = paste0(Output.dir,"/Overall_cellwise_endpoints_vs_95pc_WT_ellipse_",Output.prefix,"_boxplot.pdf"),width = 9) 
  boxplot(Cell_by_cell.Values,main="% of endpoints outside of WT 95percentile ellipse")
  dev.off()
  write.table(x = data.frame(Percent_outside=do.call(what = c,args = Cell_by_cell.Values)),
              file = paste0(Output.dir,"/Overall_cellwise_endpoints_vs_95pc_WT_ellipse_",Output.prefix,"_values.txt"),
              sep = "\t",
              row.names = T,
              col.names = T,
              quote = F)
}

Plots

Example plots using real data