Notebooks/Qoala-T_Notebook_subset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"../Figures/KoalaFramework-Logo%20copy%202.jpg\" width=\"250\"/>\n",
    "</div>\n",
    "\n",
    "# Qoala-T Notebook for subset-based predictions, June 2020\n",
    "\n",
    "### For larger datasets, it is recommended to use a rated subset of your own data to train the model. This notebook could be used to run Qoala-T once you have rated a subset.\n",
    "\n",
    "Qoala-T is a supervised learning tool that asseses accuracy of manual quality control of T1 imaging scans and their automated neuroanatomical labeling processed in FreeSurfer. It is particularly intended to use in developmental datasets. \n",
    "\n",
    "More information about Qoala-T can be found on GitHub (https://github.com/Qoala-T/QC) and in the accompanying Open Access paper:   \n",
    "Klapwijk, E.T., van de Kamp, F., Meulen, M., Peters, S., Wierenga, L.M. (2019). Qoala-T: A supervised-learning tool for quality control of FreeSurfer segmented MRI data. _NeuroImage, 189_, 116-129. https://doi.org/10.1016/j.neuroimage.2019.01.014\n",
    "\n",
    "\n",
    "The protocol of our in house developed manual QC procedure can be found here (https://github.com/Qoala-T/QC/blob/master/Qoala-T_Manual.pdf).\n",
    "\n",
    "This notebook was created by Eduard Klapwijk (https://github.com/eduardklap) based on code written by Lara Wierenga (https://github.com/larawierenga) and Olga Veth (https://github.com/OPVeth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using this notebook\n",
    "Once you have processed MRI data using FreeSurfer, you can use this notebook to extract the necessary information needed to perform Qoala-T calculations. Next, within this notebook you can run Qoala-T based on the BrainTime model. \n",
    "\n",
    "To run this notebook you need to use Jupyter Notebooks with the R kernel, for example by installing [Anaconda](https://www.anaconda.com/distribution/) and then follow [these instructions](https://docs.anaconda.com/anaconda/navigator/tutorials/r-lang/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: process images in FreeSurfer v6.0\n",
    "To be able to run the Qoala-T model, T1 MRI images should be processed in FreeSurfer V6.0 (https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall). This should be done before using this notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Install packages\n",
    "Run the next block of code to install and load the packages needed to run Qoala-T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "packages <- c(\"caret\", \"corrplot\", \"gbm\", \"plyr\", \"randomForest\", \"e1071\",\n",
    "              \"pROC\", \"DMwR\",\"dplyr\",\"pbkrtest\",\"car\",\"pbkrtest\",\"doParallel\",\"ROSE\",\"repmis\", \"plotly\", \"r2d3\")\n",
    "if (length(setdiff(packages, rownames(installed.packages()))) > 0) {\n",
    "  install.packages(setdiff(packages, rownames(installed.packages())))  \n",
    "}\n",
    "lapply(packages, library, character.only = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Run stats2table\n",
    "- Written by Olga Veth - Leiden University, MSc Computer Science student\n",
    "- Created on 30-09-2019\n",
    "- Most Recent update: 16-10-19\n",
    "- Version 3.0\n",
    "\n",
    "### Two inputs should be provided by the user of this script:\n",
    "1. Directory containing all the directories of the study participants with FreeSurfer output\n",
    "2. Name of the study and/or dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.\n",
    "datasetDir <- \"/path/to/data\" # Change Directory to your data\n",
    "setwd(datasetDir) \n",
    "# 2.\n",
    "dataset_name <- \"EMPA\" # Provide name of your study and/or dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run stats2table R code\n",
    "After you provided the directory and name of your dataset in the previous step, you can run the next block of code without having to make changes:\n",
    "\n",
    "*Note:* be sure there are only subject directories with FreeSurfer output in the root directory (also, the 'fsaverage' file in the root directory will also result in an error during this step)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "readAseg <- function(){\n",
    "  # The Aseg file of a subject is read in\n",
    "  # Volume_mm3 and StructName are selected\n",
    "  aseg_file <- data.frame(read.table(paste(\"./stats/aseg.stats\", sep=\"\"), row.names=1))[,c(3,4)] \n",
    "  asegTable <- t(data.frame(aseg_file[,1], row.names = aseg_file[,2])) # Aseg file - regular\n",
    "  return (asegTable)\n",
    "}\n",
    "\n",
    "readMetaAseg <- function(){\n",
    "  # The Aseg file of a subject is read in\n",
    "  # Its metadata containing 'lhCortex' etc. and their volume are saved\n",
    "  aseg_meta <- readLines(\"./stats/aseg.stats\", n=35)[14:35]\n",
    "  meta1 <- gsub(\"# \", \"\", aseg_meta)\n",
    "  meta <- t(data.frame(strsplit(meta1, \",\")))[,c(2,4)]\n",
    "  metaTable <- t(data.frame(meta[,2]))\n",
    "  colnames(metaTable) <- meta[,1]\n",
    "  return(metaTable)\n",
    "}\n",
    "\n",
    "editCol <- function(side, string, add){\n",
    "  # Change measure areas from 'Areaname'--> 'lh_Areaname_area'\n",
    "  return(paste(side, \"_\", string, add, sep=\"\"))\n",
    "}\n",
    "\n",
    "readAparc <- function(value){\n",
    "  # Aparc files of lh and rh are read in and the  area and thickness values of both files are retreved as well \n",
    "    # as the metadata measurements of both parts\n",
    "  # Areanames are formatted and eventually the data is saved into a data frame\n",
    "  sides <- c(\"lh\", \"rh\")\n",
    "  ifelse((value == \"area\"), pos <- 1, pos <- 2)\n",
    "  \n",
    "  for (x in 1:length(sides)){\n",
    "    areaThickness <- as.data.frame(read.table(paste(\"./stats/\", sides[x], \".aparc.stats\", sep=\"\"), row.names=1))[, c(2,4)]\n",
    "    rowValues <- rownames(areaThickness)\n",
    "    \n",
    "    meta <- readLines(paste(\"./stats/\", sides[x], \".aparc.stats\", sep=\"\"))[c(20, 21)] \n",
    "    meta1 <- gsub(\"# \", \"\", meta)\n",
    "    meta2 <- t(data.frame(strsplit(meta1, \",\")))[, c(2,4)]\n",
    "    meta3 <- data.frame(meta2[pos,2])\n",
    "    value2 <- gsub(\" \", \"\", meta2[pos,1])\n",
    "    \n",
    "    colnames(meta3) <- paste(sides[x], \"_\", value2, \"_\" , value, sep=\"\")\n",
    "    extra <- t(matrix(areaThickness[,pos]))\n",
    "    colnames(extra) <- paste(sides[x], \"_\", rowValues, \"_\", value, sep=\"\")\n",
    "    ifelse(x==1, aparcTable <- cbind(extra, meta3), aparcTable <- cbind(aparcTable, extra, meta3))\n",
    "  }\n",
    "  return(aparcTable)\n",
    "}\n",
    "\n",
    "readFiles <- function(){\n",
    "  # Aseg and Aparc files are read in and all the data is merged, starting with Aseg data followed with Aparc\n",
    "  asegTable <- readAseg()\n",
    "  metaTable <- readMetaAseg()\n",
    "  \n",
    "  areaAparc <- readAparc(\"area\")\n",
    "  thickAparc <- readAparc(\"thickness\")\n",
    "  \n",
    "  subjectTable <- cbind(asegTable, metaTable, areaAparc, thickAparc) # aparcMeta --> WhiteSurface\n",
    "  subjectTable <- data.frame(subjectTable)\n",
    "  return (subjectTable)\n",
    "}\n",
    "\n",
    "preprocTable <- function(subjectTable){\n",
    "  # Columnames are edited or removed from the table\n",
    "  removeCols <- c(\"*.WM-hypointensities$\",\"*.WM.hypointensities$\", \"*pole*\", \"*bankssts*\", \"VentricleChoroidVol\", \"*CerebralWhiteMatterVol\", \"\\\\bSurfaceHoles\\\\b\",\n",
    "                  \"SegVolFile.mri.aseg.mgz.\", \"*CorticalWhiteMatterVol\")\n",
    "  remove <- grep(paste(removeCols, collapse=\"|\"), colnames(subjectTable))\n",
    "  subjectTable <- subjectTable[, -remove]\n",
    "  \n",
    "  colnames(subjectTable) <- gsub(\"^X\\\\.\", \"\", colnames(subjectTable))\n",
    "  colnames(subjectTable) <- gsub(\"_\\\\.\", \"_\", colnames(subjectTable))\n",
    "  colnames(subjectTable) <- gsub(\"-\", \".\", colnames(subjectTable))\n",
    "  colnames(subjectTable) <- gsub(\" \", \"\", colnames(subjectTable))\n",
    "  \n",
    "  colnames(subjectTable)[which(colnames(subjectTable) == \"eTIV\")] <- \"EstimatedTotalIntraCranialVol\"\n",
    "  colnames(subjectTable)[which(colnames(subjectTable) %in% c(\"rd.Ventricle\", \"th.Ventricle\", \n",
    "                                                             \"5th.Ventricle\"))] <- c(\"X4th.Ventricle\", \"X3rd.Ventricle\", \"X5th.Ventricle\") # change to names\n",
    "  \n",
    "  return(subjectTable)\n",
    "}\n",
    "\n",
    "main <- function(){\n",
    "  # It loops through all subjects sub-directories in the given directory\n",
    "  # With every single subject, data is retrieved and written in a row in the final table. \n",
    "  # The result is saved into a .CSV file\n",
    "  subjects <- c()\n",
    "  first <- T\n",
    "  subjectDirs <- unique(list.dirs('.', recursive=FALSE)) # Get all sample subject\n",
    "  for (x in 1:length(subjectDirs)){\n",
    "    setwd(paste(datasetDir, subjectDirs[x], sep=\"\"))\n",
    "    statsDirs <- list.dirs('.', recursive=FALSE)\n",
    "    if (file.exists(\"./stats/aseg.stats\")){\n",
    "      subjectTable <- readFiles()\n",
    "      subjectTable <- preprocTable(subjectTable)\n",
    "      if (first == T){\n",
    "        stats2Table <- subjectTable\n",
    "        subjects <- c(subjects, substring(subjectDirs[x], 3))\n",
    "        first = F\n",
    "      } \n",
    "      else if (ncol(subjectTable) == ncol(stats2Table)&& (first == F)){\n",
    "        stats2Table <- rbind(stats2Table, subjectTable)\n",
    "        subjects <- c(subjects, substring(subjectDirs[x], 3))\n",
    "      } \n",
    "      \n",
    "    }\n",
    "  }\n",
    "  \n",
    "  stats2Table <- data.frame(stats2Table)\n",
    "  stats2Table <- cbind(subjects, stats2Table) \n",
    "  colnames(stats2Table)[1] <- \"participant_id\"  \n",
    "  setwd(datasetDir)\n",
    "  write.csv(stats2Table, paste(\"FS_Output_\", dataset_name,\".csv\", sep=\"\"), row.names = FALSE)\n",
    "    }\n",
    "main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CSV output\n",
    "Now there should be a csv file called \"FreeSurfer_Ouput_(Dataset_name).csv\" in the datasetDir that was declared by the user in the previous step.  \n",
    "Let's have a look at this table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "read.csv(paste(\"FS_Output_\",dataset_name,\".csv\", sep=\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: add rating to a subset of your data\n",
    "\n",
    "For a subset of your data you should do manual quality control. These ratings should be added to the csv file in a new column.\n",
    "\n",
    "Make sure your data format looks like simulated_data_B_model.Rdata:\n",
    "- First column should contain outcome manual quality control --> \"Rating\".\n",
    "- Subset of data is rated, with two factor levels ('Include' and 'Exclude').\n",
    "- Remaining data has no rating (NA)\n",
    "\n",
    "In the next code block the example data (simulated_data_B_model.Rdata) is shown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_data <- get(load(url(\"https://github.com/Qoala-T/QC/blob/master/ExampleData/simulated_data_B_subset.Rdata?raw=true\")))\n",
    "#show head + tail of example data\n",
    "head(example_data)\n",
    "tail(example_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combine FreeSurfer output and Ratings\n",
    "Once you have rated a substantial subset of your data and you have added a \"Rating\" column to your csv file, save the file as \"FS_Output_Ratings_{dataset_name}.csv\" to use in the final step.  \n",
    "\n",
    "Alternatively, if you have saved the Ratings in a separate csv (e.g., called \"Ratings_subset_{dataset_name}.csv\") use the following code to merge the FS_Output and Ratings for your dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load Output\n",
    "FS_data <- read.csv(paste(\"FS_Output_\", dataset_name, \".csv\", sep=\"\"), header=T)\n",
    "# load ratings\n",
    "Ratings <- read.csv(paste(\"Ratings_subset_\", dataset_name, \".csv\", sep=\"\"), header=T)\n",
    "# merge\n",
    "Output_Ratings <- merge(Ratings,FS_data, by = \"participant_id\")\n",
    "# save\n",
    "write.csv(Output_Ratings, paste(\"FS_Output_Ratings_\", dataset_name,\".csv\", sep=\"\"), row.names = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Run the Qoala-T model on your data\n",
    "\n",
    "- written by Lara Wierenga, PhD at Leiden University\n",
    "- Created on March 2018\n",
    "- Most Recent update: 4-11-19\n",
    "- Version 1.2\n",
    "\n",
    "### Run the model and save csv output file\n",
    "Now that you have your FreeSurfer output in the right format you can run the Qoala-T model based on the BrainTime training set. The results will be saved in a csv file in your dataset directory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset <- read.csv(paste(\"FS_Output_Ratings_\", dataset_name, \".csv\", sep=\"\"), header=T, row.names=1)\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# Next match col.names to Qoala_T_model\n",
    "# -----------------------------------------------------------------\n",
    "githubURL <- \"https://github.com/Qoala-T/QC/blob/master/Qoala_T_model.Rdata?raw=true\"\n",
    "rf.tune <- get(load(url(githubURL)))\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# reorder colnames of dataset to match traningset\n",
    "# -----------------------------------------------------------------\n",
    "dataset_names <- c(\"Rating\",names(rf.tune$trainingData)[-ncol(rf.tune$trainingData)])\n",
    "dataset <- dataset[,dataset_names]\n",
    "dataset <- dataset[complete.cases(dataset[-1]),]\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# reorder colnames of dataset to match traningset\n",
    "# -----------------------------------------------------------------\n",
    "dataset_names <- c(\"Rating\",names(rf.tune$trainingData)[-ncol(rf.tune$trainingData)])\n",
    "dataset <- dataset[,dataset_names]\n",
    "dataset <- dataset[complete.cases(dataset[-1]),]\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# Split into training and testing datasets \n",
    "# -----------------------------------------------------------------\n",
    "# select rated data as training data\n",
    "training = dataset[!is.na(dataset$Rating),]\n",
    "training$Rating = as.factor(training$Rating) \n",
    "\n",
    "# select remaining unrated data as testing data\n",
    "testing = dataset[is.na(dataset$Rating),]\n",
    "testing$Rating = as.factor(testing$Rating)\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# Setting up computational nuances of the train function for internal cross validation \n",
    "# -----------------------------------------------------------------\n",
    "ctrl = trainControl(method = 'repeatedcv',  \n",
    "                    number = 2,    \n",
    "                    repeats = 10,         \n",
    "                    summaryFunction=twoClassSummary, \n",
    "                    classProbs=TRUE,        \n",
    "                    allowParallel=FALSE,    \n",
    "                    sampling=\"rose\") # 'rose' is used to oversample the imbalanced data       \n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# Estimate model \n",
    "# -----------------------------------------------------------------\n",
    "rf.tune = train(y=training$Rating,\n",
    "                  x=subset(training, select=-c(Rating)),\n",
    "                  method = \"rf\",\n",
    "                  metric = \"ROC\",\n",
    "                  trControl = ctrl,\n",
    "                  ntree = 501,\n",
    "                  tuneGrid=expand.grid(mtry = c(8)),\n",
    "                  verbose=FALSE)\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# External cross validation on unrated data (1 repetition)\n",
    "# -----------------------------------------------------------------\n",
    "rf.pred <- predict(rf.tune,subset(testing, select=-c(Rating)))\n",
    "rf.probs <- predict(rf.tune,subset(testing, select=-c(Rating)),type=\"prob\") \n",
    "head(rf.probs)\n",
    "\n",
    "# -----------------------------------------------------------------\n",
    "# Saving output\n",
    "# ----------------------------------------------------------------\n",
    "# create empty data frame\n",
    "Qoala_T_predictions_subset_based <- data.frame(matrix(ncol = 4, nrow = nrow(rf.probs)))                   \n",
    "colnames(Qoala_T_predictions_subset_based) = c('participant_id','Scan_QoalaT', 'Recommendation', 'manual_QC_adviced') \n",
    "\n",
    "# fill data frame\n",
    "Qoala_T_predictions_subset_based$participant_id <- row.names(rf.probs)\n",
    "Qoala_T_predictions_subset_based$Scan_QoalaT <- rf.probs$Include*100 \n",
    "Qoala_T_predictions_subset_based$Recommendation <- rf.pred\n",
    "Qoala_T_predictions_subset_based$manual_QC_adviced <- ifelse(Qoala_T_predictions_subset_based$Scan_QoalaT<60&Qoala_T_predictions_subset_based$Scan_QoalaT>40,\"yes\",\"no\")\n",
    "Qoala_T_predictions_subset_based <- Qoala_T_predictions_subset_based[order(Qoala_T_predictions_subset_based$Scan_QoalaT, Qoala_T_predictions_subset_based$participant_id),]\n",
    "  \n",
    "  \n",
    "csv_Qoala_T_predictions_subset_based = paste('Qoala_T_predictions_subset_based_',dataset_name,'.csv', sep = '')\n",
    "write.csv(Qoala_T_predictions_subset_based, file = csv_Qoala_T_predictions_subset_based, row.names=F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CSV output\n",
    "Now there should be a csv file called \"Qoala_T_predictions_subset_based_{Dataset_name}.csv\" in the datasetDir.\n",
    "Let's have a look at this table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "read.csv(paste('Qoala_T_predictions_subset_based_',dataset_name,'.csv', sep = ''))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "excl_rate <- table(Qoala_T_predictions_subset_based$Recommendation)\n",
    "\n",
    "fill_colour <- rev(c(\"#256da8\",\"#CF4A30\"))\n",
    "font_size <- 12\n",
    "text_col <- \"Black\"\n",
    "\n",
    "p <- ggplot(Qoala_T_predictions_subset_based, aes(name = participant_id, x=Scan_QoalaT,y=1,col=Recommendation)) +  \n",
    "  annotate(\"rect\", xmin=30, xmax=70, ymin=1.12, ymax=.88, alpha=0.2, fill=\"#777777\") +\n",
    "  geom_jitter(alpha=.8,height=.1,size=5) +\n",
    "  ggtitle(paste(\"Qoala-T estimation subset based for \",dataset_name,\"\\nMean Qoala-T Score = \",round(mean(Qoala_T_predictions_subset_based$Scan_QoalaT),1),sep=\"\")) + \n",
    "  annotate(\"text\", x=20, y=1.15, label=paste(\"Excluded = \",as.character(round(excl_rate[1])),\" scans\",sep=\"\")) + \n",
    "  annotate(\"text\", x=80, y=1.15, label=paste(\"Included = \",as.character(round(excl_rate[2])),\" scans\",sep=\"\")) + \n",
    "  scale_colour_manual(values=fill_colour) +\n",
    "  theme_bw() +\n",
    "  theme(panel.grid.major = element_blank(), \n",
    "        panel.grid.minor = element_blank(), \n",
    "        panel.border = element_blank(),\n",
    "        axis.text.x = element_text (size = font_size,color=text_col),\n",
    "        axis.text.y = element_blank(),\n",
    "        axis.title.x = element_text (size = font_size,color=text_col), \n",
    "        axis.title.y = element_blank(), \n",
    "        axis.ticks=element_blank(),\n",
    "        plot.title=element_text (size =16,color=text_col,hjust=.5)\n",
    "  )\n",
    "print(p)\n",
    "\n",
    "filename<- paste(\"Figure_Rating_subset_based_\",dataset_name,\".pdf\",sep=\"\")\n",
    "dev.copy(pdf,filename,width=30/2.54, height=20/2.54)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## That's it! We hope these predictions help in the quality control procedure of your study.  \n",
    "  \n",
    "When using Qoala-T please include the following citation:\n",
    "  \n",
    "Klapwijk, E.T., van de Kamp, F., Meulen, M., Peters, S., Wierenga, L.M. (2019). Qoala-T: A supervised-learning tool for quality control of FreeSurfer segmented MRI data. _NeuroImage, 189_, 116-129. https://doi.org/10.1016/j.neuroimage.2019.01.014"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}