Automated_Rotating_Beam_Analysis/R_scripts/analysis_functions.R

###1. Preparation of Data###

###1.1 Converting .csv to dataframe.

#'This is a modified function of fread from data.table package by adjusting to the format of DeepLabCut files
#'
#'@param [file]: file path of .csv file from the DLC-analyzed video.
#'@return: data frame, has to be assigned to a variable.


read.dlc = function(file){
  Table = fread(file, skip=2, data.table = FALSE)
  col = colnames(fread(file, skip=1, data.table = FALSE))
  col[1] = "seq"
  for(i in 2:length(col)){
    if(i %in% seq(2, length(col), 3)){
      col[i] = paste(col[i], "_x", sep="")
    }
    else if(i %in% seq(3, length(col), 3)){
      col[i] = paste(col[i], "_y", sep="")
    }
    else if(i %in% seq(4, length(col), 3)){
      col[i] = paste(col[i], "_likeli", sep="")
    }
  }
  colnames(Table) = col
  return(Table)
}

###1.2 Inverting y axis###

#'Inverts the Y axis.
#' 
#'@param [Table]: dataframe after application of read.dlc.
#'@return: dataframe, has to be assigned to a variable.

y_invert<-function(Table){
  for(i in seq(3, ncol(Table)-3, 3)){
    for(j in 1:nrow(Table)){
      Table[j,i]=-(Table[j,i]-Table$height[1])
    }
  }
  return(Table)
}

###1.3 Deleting Low Likelihood###

#'Deletes X and Y coordinates under a given threshold. X and Y coordinates of a data point under the desired threshold turn into NAs.
#'
#'@param [Table]: dataframe after application of read.dlc and y_invert.
#'@param [threshold]: value of desired minimum likelihood.
#'@return: dataframe, has to be assigned to a variable.


lowlikeli = function(Table, threshold = .90){
  likelirows = grep("likeli", names(Table))
  for(i in likelirows){
    for(j in 1:nrow(Table)){
      if(Table[j,i] < threshold){
        Table[j,i-1]<-NA
        Table[j,i-2]<-NA
      }
    }
  }
  return(Table)
}

###1.4 Adding Video Resolution###

#'Adds three columns with the Height, the Width and the FPS of the Video for easy access
#'
#'@param [Table]: dataframe after application of read.dlc, y_invert and lowlikeli.
#'@param [Height]: Height of the analyzed video.
#'@param [Width]: Width of the analyzed video.
#'@param [FPS]: Frames per second of the analyzed video.
#'@return: dataframe, has to be assigned to a variable.

addreso = function(Table, Height, Width, FPS){
  Table$height = Height
  Table$width = Width
  Table$fps = FPS
  return(Table)
}

###1.5 Combines all preparation functions###

#'Combines all preparation functions. First reads DLC file, inverts the Y axis, removes data points under a given threshold and adds the columns Height, Width and FPS.
#'
#'@param [Path]: file path of .csv file from the DLC-analyzed video.
#'@param [FPS]: Frames per second of the analyzed video.
#'@param [Height]: Height of the analyzed video.
#'@param [Width]: Width of the analyzed video.
#'@param [Threshold]: value of desired minimum likelihood.
#'@return: dataframe, has to be assigned to a variable.

prep.dlc<-function(Path, FPS, Height, Width, Threshold=.90){
  Table=read.dlc(Path)
  Table=addreso(Table, Height, Width, FPS)
  Table=y_invert(Table)
  Table=lowlikeli(Table, Threshold)
  print("data process succeeded!")
  return(Table)
}

###2. Conditions###

###2.1 Beam Left###

#'Ignores data points on the left side of the left beam landmark.
#'
#'@param [Table]: dataframe after preparation of data. For preparation of data go to point 1.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: 0 or 1. If 0, data point is on the right side of the left beam landmark. If 1, the data point is on the left side of the left beam landmark.

beam_left<-function(Table, x){
  left=0
  if(x<median(Table$'beam_left_top_x', na.rm=T)){
    left=1
  }
  return(left)
}

###2.2 Beam Right###

#'Ignores data points on the right side of the right beam landmark.
#'
#'@param [Table]: dataframe after preparation of data. For preparation of data go to point 1.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: 0 or 1. If 0, data point is on the left side of the right beam landmark. If 1, the data point is on the right side of the right beam landmark.

beam_right<-function(Table, x){
  right=0
  if(x>median(Table$'beam_right_top_x', na.rm=T)){
    right=1
  }
  return(right)
}

###2.3 Beam Top###

#'Gives the corresponding Y coordinate of the upper beam edge to a X coordinate of a data point.
#'
#'@param [Table]: dataframe after preparation of data. For preparation of data go to point 1.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: Y coordinate of the upper beam edge at the given X coordinate. Will be compared to the Y coordinate of the given X coordinate.

beam_top<-function(Table, x){
  m = (median(Table$'beam_left_top_y', na.rm=T) - median(Table$'beam_right_top_y', na.rm=T))/
    (median(Table$'beam_left_top_x', na.rm=T) - median(Table$'beam_right_top_x', na.rm=T))
  b = median(Table$'beam_right_top_y', na.rm=T) - m*median(Table$'beam_right_top_x', na.rm=T)
  beam_y=m*x+b
  return(beam_y)
}

###2.4 Beam Bottom###

#'Gives the corresponding Y coordinate of the lower beam edge to a X coordinate of a data point.
#'
#'@param [Table]: dataframe after preparation of data. For preparation of data go to point 1.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: Y coordinate of the lower beam edge at the given X coordinate. Will be compared to the Y coordinate of the given X coordinate.

beam_bottom<-function(Table, x){
  m = (median(Table$'beam_left_bottom_y', na.rm=T) - median(Table$'beam_right_bottom_y', na.rm=T))/
    (median(Table$'beam_left_bottom_x', na.rm=T) - median(Table$'beam_right_bottom_x', na.rm=T))
  b = median(Table$'beam_right_bottom_y', na.rm=T) - m*median(Table$'beam_right_bottom_x', na.rm=T)
  beam_y=m*x+b
  return(beam_y)
}

###3. Evaluation###

###3.1 Start###

#'Determines the start point of the evaluation. Start point is set to the first time point, when both forepaw and hindpaw are close to the upper beam edge.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@return: frame number of the start of the evaluation.

start<-function(Table){
  for(i in 1:nrow(Table)){
    if(is.na(Table$forepaw_right_x[i])==F && is.na(Table$forepaw_right_y[i])==F && is.na(Table$hindpaw_right_x[i])==F && is.na(Table$hindpaw_right_y[i])==F && Table$forepaw_right_y[i]<(beam_top(Table, Table$forepaw_right_x[i]+(Table$height[1]*0.01))) && Table$hindpaw_right_y[i]<(beam_top(Table, Table$hindpaw_right_x[i]+(Table$height[1]*0.01)))){
      return(Table$seq[i])
      break()
    }
  }
}

###3.2 End###

###3.2.1 Turn###

#'Determines whether the mouse turns around on the beam at a given X coordinate.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: 0 or frame index of X coordinate. If 0 the mouse don't face in the wrong direction at the given X coordinate.

turn<-function(Table, x){
  turn=0
  if(is.na(Table$snout_x[x]) || is.na(Table$tailbase_x[x])){
    return(turn)
  }
  if(Table$snout_x[x]<Table$tailbase_x[x]){
    turn=Table$seq[x]
  }
  return(turn)
}

###3.2.2 Drop###

#'Determines whether the mouse dropped of the beam at a given X coordinate.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: 0 or frame index of X coordinate. If 0, the mouse is still on the beam at the given X coordinate.

drop<-function(Table, x){
  drop=0
  if(is.na(Table$snout_y[x]) || is.na(Table$tailbase_y[x])){
    return(drop)
  }
  if(Table$snout_y[x]<(beam_bottom(Table, Table$snout_x[x])-Table$height[1]*0.01) && Table$tailbase_y[x]<(beam_bottom(Table, Table$tailbase_x[x])-Table$height[1]*0.01)){
    drop=Table$seq[x]
  }
  return(drop)
}

###3.2.3 End of Beam###

#'Determines whether to mouse reached the last 5 centimeter of the beam at a given X coordinate.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [x]: Frame number of a data point which will be checked.
#'@return: 0 or frame index of X coordinate. If 0, the mouse isn't at the end of the beam at the given X coordinate.
#'
end_beam<-function(Table, x){
  end_beam=0
  if(is.na(Table$forepaw_right_x[x])){
    return(end_beam)
  }
  if(distance(Table, Table$forepaw_right_x[x])>115){
    end_beam=Table$seq[x]
  }
  return(end_beam)
}

###3.2.4 Stop###

#'Determines whether to mouse stops while running on the beam at a given X coordinate.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: 0 or frame index of X coordinate. If 0, the mouse isn't stopping at the given X coordinate. 

stop<-function(Table, x){
  stop=0
  if(is.na(Table$forepaw_right_x[x]) || is.na(Table$hindpaw_right_x[x])){
    return(stop)
  }
  if(Table$seq[x]>start(Table) && Table$seq[i]>Table$fps[1] && Table$forepaw_right_x[x]<(Table$forepaw_right_x[x-Table$fps[1]]+Table$width[1]*0.05) && Table$hindpaw_right_x[x]<(Table$hindpaw_right_x[x-Table$fps[1]]+Table$width[1]*0.05)){
    stop=Table$seq[x]
  }
  return(stop)
}

###3.2.5 End###

#'Determines the end point of the evaluation. End point is set when the mouse either turns, drops, stops, reaches the end of the beam or 60 seconds past.
#'
#'@param [Table]: dataframe after preparation of data, the conditions and determination of the start point. For preparation of data go to point 1. For conditions go to point 2. For determination of start point go to point 3.1.
#'@return: frame index of the end of the evaluation.

end<-function(Table){
  for(i in 1:nrow(Table)){
    if(Table$seq[i]<start(Table)){
      next()
    }
    if(turn(Table, i)!=0){
      return(Table$seq[i])
      break()
    }
    if(drop(Table, i)!=0){
      return(Table$seq[i])
      break()
    }
   if(end_beam(Table, i)!=0){
      return(Table$seq[i])
      break()
    }
  if(stop(Table, i)!=0){
     return(Table$seq[i])
     break()
   }
   if(Table$seq[i]>=(Table$fps[1]*60)){
      return(Table$seq[i])
      break()
    }
  }
}

rfend<-function(Table){
  for(i in 1:nrow(Table)){
    if(Table$seq[i]<start(Table)){
      next()
    }
    if(turn(Table, i)!=0){
      return("turn")
      break()
    }
    if(drop(Table, i)!=0){
      return("drop")
      break()
    }
    if(end_beam(Table, i)!=0){
      return("end of beam")
      break()
    }
    if(stop(Table, i)!=0){
      return("stop")
      break()
    }
    if(Table$seq[i]>=(Table$fps[1]*60)){
      return("time over")
      break()
    }
  }
}

###3.3 Speed###

###3.3.1 Distance###

#'Calculates the traveled distance of the mouse at a given X coordinate.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [x]: X coordinate of a data point which will be checked.
#'@return: traveled distance of the mouse at the given X coordinate in centimeter [cm]. Minimum 0 cm, Maximum 120 cm.

distance<-function(Table, x){
  y<-((120/(median(Table$'beam_right_top_x', na.rm=T)-median(Table$'beam_left_top_x', na.rm=T)))*(x-median(Table$'beam_left_top_x', na.rm=T)))
  if(y>115){
    y=120
  }
  return(y)
}

###3.3.2 Time###

#'Determines the time from the start till the end of the evaluation.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [start]: either frame number of the start of the evaluation or function start with @param [Table], see 3.1.
#'@param [end]: either frame number of the end of the evaluation or function end with @param [Table], see 3.2.5.
#'@param [fps]: Frames per second of the analyzed video.
#'@return: time of evaluation in seconds [s].

time<-function(Table, start, end){
  time=(end-start)/Table$fps[1]
  return(time)
}

###3.3.2 Speed###

#'Determines the speed of the mouse during the evaluation.
#'
#'@param [Table]: dataframe after preparation of data, the conditions and determination of the start and the end point of the evaluation. For preparation of data go to point 1. For conditions go to point 2.
#'@param [end]: either frame number of the end of the evaluation or function end with @param [Table], see 3.2.5.
#'@param [time]: either time in seconds or function time with @param [Table], @param [start], @param [end] and @param [fps], see 3.3.2.
#'@return: speed of the mouse during the evaluation in centimeter per second [cm/s].

speed<-function(Table, end, time){
  speed=distance(Table, Table$forepaw_right_x[end+1])/time
  return(speed)
}

###3.4 Hindlimb Drop### 

#'Detects right and left hindlimb drops at a given X coordinate.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data got to point 1. For conditions go to point 2.
#'@param [x]: X coordinate of the data point which will be checked.
#'@return: frame number of hindlimb drop 

hdrop<-function(Table, Vec){
  for(i in 1:nrow(Table)){
    redropcount=0
    if(length(Vec)>0){
      for(j in (Vec[length(Vec)]+1):i){
        if(is.na(Table$hindpaw_right_x[j])==F && is.na(Table$hindpaw_right_y[j])==F && is.na(Table$hindpaw_left_x[j])==F && is.na(Table$hindpaw_left_y[j])==F && Table$hindpaw_right_y[j]>beam_bottom(Table, Table$hindpaw_right_x[j]) && Table$hindpaw_left_y[j]>beam_bottom(Table, Table$hindpaw_left_x[j])){
          redropcount=redropcount+1
        }
      }    
    }
    if(length(Vec)>0 && redropcount==0){
        next()
    }
    if(is.na(Table$hindpaw_right_x[i])==F && is.na(Table$hindpaw_right_y[i])==F && Table$hindpaw_right_y[i]<beam_bottom(Table, Table$hindpaw_right_x[i]) && drop(Table, i)==0){
      Vec=cbind(Vec,Table$seq[i])
    }
    if(is.na(Table$hindpaw_left_x[i])==F && is.na(Table$hindpaw_left_y[i])==F && Table$hindpaw_left_y[i]<beam_bottom(Table, Table$hindpaw_left_x[i]) && drop(Table, i)==0){
      Vec=cbind(Vec,Table$seq[i])
    }
  }
  return(Vec)
}


###4. Further Analysis###

###4.1 Angle###

#'Calculates the angle between three bodyparts for all frames.
#'
#'@param [Table]: dataframe after preparation of data and the conditions. For preparation of data go to point 1. For conditions go to point 2.
#'@param [Ray1, Vertex, Ray2]: columnnames of the desired bodyparts with quotation marks, e.g. angle(table, "snout", "neck", "tailbase").
#'@param [Vertex]: chosen bodypart indicates the center (or vertex) of the calculated angles.
#'@param [Ray1, Ray2]: chosen bodyparts indicate the outside points of the calculated angles.
#'@return: dataframe with all angles between the given bodyparts, has to be assigned to a variable.

angle<-function(Table, Ray1, Vertex, Ray2){
  vec=data.frame()
  grepRay1=grep(Ray1, names(Table))
  grepRay2=grep(Ray2, names(Table))
  grepVertex=grep(Vertex, names(Table))
  for(i in 1:nrow(Table)){
    if(is.na(Table[i, grepRay1[1]]) || is.na(Table[i, grepRay2[1]]) || is.na(Table[i, grepVertex[1]]) || is.na(Table[i, grepRay1[2]]) || is.na(Table[i, grepRay2[2]]) || is.na(Table[i, grepVertex[2]])){
      next()
    }
    else{
      R1<-c(Table[i, grepRay1[1]], Table[i, grepRay1[2]])
      R2<-c(Table[i, grepRay2[1]], Table[i, grepRay2[2]])
      V<-c(Table[i, grepVertex[1]], Table[i, grepVertex[2]])
      VR1=R1-V
      VR2=R2-V
      skalar=sum(VR1*VR2)
      BVR1=VR1[1]^2+VR1[2]^2
      BVR2=VR2[1]^2+VR2[2]^2
      rad<-acos(skalar/(sqrt(BVR1)*sqrt(BVR2)))
      deg<-180*rad/pi
      vec[Table$seq[i],1]=deg
    }
  }
  return(vec)
}

singleangle<-function(Table, j, Ray1, Vertex, Ray2){
  grepRay1=grep(Ray1, names(Table))
  grepRay2=grep(Ray2, names(Table))
  grepVertex=grep(Vertex, names(Table))
  if(is.na(Table[j, grepRay1[1]]) || is.na(Table[j, grepRay2[1]]) || is.na(Table[j, grepVertex[1]]) || is.na(Table[j, grepRay1[2]]) || is.na(Table[j, grepRay2[2]]) || is.na(Table[j, grepVertex[2]])){
    return(NA)
  }
  R1<-c(Table[j, grepRay1[1]], Table[j, grepRay1[2]])
  R2<-c(Table[j, grepRay2[1]], Table[j, grepRay2[2]])
  V<-c(Table[j, grepVertex[1]], Table[j, grepVertex[2]])
  VR1=R1-V
  VR2=R2-V
  skalar=sum(VR1*VR2)
  BVR1=VR1[1]^2+VR1[2]^2
  BVR2=VR2[1]^2+VR2[2]^2
  rad<-acos(skalar/(sqrt(BVR1)*sqrt(BVR2)))
  deg<-180*rad/pi
  return(deg)
}