R/fs_stats.R

In VertexWiseR: Simplified Vertex-Wise Analyses of Whole-Brain and Hippocampal Surface

#' @title fs_stats()
#'
#' @description Extracts descriptive statistics, for the whole-brain and subcortical region-of-interests (ROI), within a FreeSurfer subjects directory. It reads them from the aseg.stats file, as generated by the default FreeSurfer preprocessing pipeline.
#'
#' @param sdirpath A string object indicating the path to the 'FreeSurfer' subjects directory. Default is the current working directory ("./").
#' @param sublist A string object indicating the path to the subject list generated by SURFvextract as 'sublist.txt' (optional). This allows users to retrieve stats only from a selected list of subjects. The subject list is a list with 1 subject ID per line.
#' @param ROImeasure A string object indicating what summary measure to extract for the subocrtical ROIs. Choices include: 'NVoxels', 'Volume_mm3', 'StructName', 'normMean', 'normStdDev', 'normMin', 'normMax', and 'normRange'. Default is 'Volume_mm3'.
#'
#' @returns A data.frame object with N columns per aseg.stats measures and N row per subjects.
#' @examples
#' fs_stats(sdirpath="freesurfer_subjdir")
#' @importFrom stringr str_split
#'@importFrom utils read.table
#' @export

fs_stats=function(sdirpath="./", sublist, ROImeasure='Volume_mm3') 
{
  
  
  #check if sdirpath contains aseg.stats files 
  dircount = dir(path=sdirpath, recursive=TRUE, pattern="aseg.stats", include.dirs = TRUE, full.names = TRUE)
  if (length(dircount)==0) { return(message('aseg.stats files could not be found in the set sdirpath')) }
  
  #decide subject_IDs to retrieve if sublist argument given
  if(!missing(sublist)) 
  {
    sublist=read.table(sublist)
    
    #get parent directory of stats/aseg.stats
    subdirs = stringr::str_split(dir(path=sdirpath, recursive=TRUE,
                                     pattern="aseg.stats"), pattern='/stats/') 
    subdirs= unlist(lapply(subdirs, function(x) x[1]))
    
    #select only subject IDs match that directory name and use that
    dircount=dircount[which(subdirs %in% sublist[,1])==TRUE]
    
  }  
  
  
  
  #creates main data.frame
  FS_STAT=data.frame(matrix(ncol=1,nrow=length(dircount)))
  colnames(FS_STAT) <- 'subject_ID'
  
  #for each aseg stats, extracts measure values
  for (i in dircount) 
  {
    #read aseg.stats (read.delim to export the first Measures)
    stattable=read.delim(i, sep="\t")
    
    for (l in 1:nrow(stattable)) #for each line
    {  
      #if subjectname in the line store subject ID
      if (grepl("subjectname", stattable[l,])) 
      {
        subject_ID=stringr::str_split(stattable[l,],
                                      pattern="subjectname ")[[1]][2]
        FS_STAT$subject_ID[which(dircount==i)]=subject_ID
      }
      
      #for each Measure, save name and value
      if (grepl("Measure", stattable[l,])) 
      {
        measure=stringr::str_split(stattable[l,],pattern=",")
        meas_name=stringr::str_split(measure[[1]][1], pattern="Measure ")[[1]][2]
        #if no column for measure name create it
        if (meas_name %in% colnames(FS_STAT)==FALSE)
        {  FS_STAT[[meas_name]] <- NA }
        #append measure value
        meas_val=as.numeric(gsub("\\D", "", measure))
        FS_STAT[[which(dircount==i),meas_name]] <- meas_val
      }
    }
    
    #read aseg.stats again (to easily get ROI mean volume)
    stattable2=read.table(i)
    
    for (l2 in 1:nrow(stattable2))
    {
      #get ROI name
      roi_name=stattable2[l2,5]
      #if no column for measure name create it
      if (roi_name %in% colnames(FS_STAT)==FALSE)
      {  FS_STAT[[roi_name]] <- NA }
      
      #select column based on user-specified ROImeasure:
      if(ROImeasure=='NVoxels'){r=3}
      else if (ROImeasure=='Volume_mm3') {r=4}
      else if (ROImeasure=='normMean') {r=6}
      else if (ROImeasure=='normStdDev') {r=7}
      else if (ROImeasure=='normMin') {r=8}
      else if (ROImeasure=='normMax') {r=9}
      else if (ROImeasure=='normRange') {r=10}
      
      #append measure value
      roi_val=stattable2[l2,r]
      FS_STAT[[which(dircount==i),roi_name]] <- roi_val
    }
    
  }
  
  
  return(FS_STAT)  
}