final_project_mamourian.Rmd

---
title: "BMIN503/EPID600 Project: Structural Brain Imaging Biomarker for Traumatic Brain Injury"
author: "Elizabeth Mamourian"
output: 
  html_document:
    toc: true
    toc_float:
      collapsed: false
    keep_md: true
    depth: 3 
    theme: paper 
    highlight: tango
---
```{r set-options, echo=FALSE, cache=FALSE}
options(width = 400)
```  

# Overview
This project investigates the structural brain change resulting from traumatic brain injury (TBI), 
specifically probing for a reliable measure to quantify the severity of TBI, which is agnostic to 
the type and location of the trauma. The aim is to quantify structural brain volumes, to identify 
regions with disease-related atrophy, and to compare the relative distribution of brain volume in 
TBI stubjects against normal controls.
   
 
https://github.com/mamourie/BMIN503_Final_Project


# Introduction 
Each year, 1.5 million Americans sustain a traumatic brain injury, which is the cause of a 
long-term disability in 80,00-90,000 people [1]. Yet, there is no quantifiable measure that 
can differentiate between those people who recover and those who suffer from a long-term 
disability. For that reason, it is more difficult for those people with a disability to qualify 
for the services that they need to deal with their injury. The tau pathology in chronic 
traumatic encephalopathy is well documented post-mortem [2], but there is currently no in vivo 
method for quantifying the severity of a traumatic brain injury or diagnosing a patient with CTE.

This project takes advantage of computer vision to address a clinical problem. Using feature 
extraction from medical image analysis, it is possible to use a data-driven approach to evaluate 
abnormal regions of interest and to calculate a non-specific measure of disease-related brain change.




# Methods
The data used are volumetric measurements of brain regions of interest, based on a multi atlas 
segmentation method [3]. All results are shown based on the absolute volumetric measure
and the volume normalized by intracranial volume. Z-scores for each ROI are calculated using
and age and sex matched normal control group, and applied to the case subjects, to see how
the disease affects cases in comparison to normal controls. 

The data was collected by the UK Biobank. The UK Biobank is a large-scale, prospective, 
population-based, observational cohort study. Voluntary participants were recruited across the 
United Kingdom, between 2006-2010. Extensive baseline data was collected, including questionnaires 
about health and lifestyle, anthropometric measures, biological samples, and samples for genetic 
analysis. Neuroimaging data was collected using a 3.0-Tesla MRI scanner (Siemens Skyra, Siemens 
Healthcare, Erlangen, Germany) with a 32-channel radio frequency receiver head coil. The T1 imaging 
sequence used a 3D magnetization-prepared rapid acquisition with gradient echo (voxel: 1 x 1 x 1 mm; 
matrix: 208 x 256 x 256; inversion time: 880 msec; repetition time: 2000 msec). This research has 
been conducted using the UK Biobank Resource under Application Number 35148.

      

# Data Preparation

![](/Users/elizabethmamourian/Desktop/Data_Prep.png)

```{r data_prep_1, echo = TRUE, cache = TRUE, message=FALSE}
library(dplyr)
library(moments)
library(ggplot2)

home <- "/Users/elizabethmamourian/Box/mamourie/ISTAGING/TBI/"
# dir <- "/Users/elizabethmamourian/Box/mamourie/project_copies/ClinicalDataConsolidation_201911/"

fromBeginning <- FALSE
if (fromBeginning==TRUE) {
df <- read.csv(paste0(home,"df.csv"), stringsAsFactors = FALSE)
icv.orig <- read.csv(paste0(home,"UKB_DLICV_Vol.csv"),stringsAsFactors = FALSE)
u.batch.orig <- read.csv(paste0(dir,"Protocols/UKBiobank_Consolidated/ukbiobank_batch0-10.csv"), stringsAsFactors=FALSE)
u.batch <- u.batch.orig %>%
  mutate(PTID=as.character(eid)) %>%
  mutate(Visit_Code=case_when(VISIT==2 ~ '2_0')) %>%
  select(-c(BMI,DSST,TMT_A,TMT_B,eid))
}
```

### select TBI subjects

```{r data_prep_2, echo = TRUE, cache = TRUE, message=FALSE}
if (fromBeginning==TRUE) {

# select TBI subjects from self-report f20002 -> (tbi.sr) ####
sr <- u.batch %>%
  select(c(names(u.batch)[startsWith(names(u.batch),"noncancer_illness_code_selfreported_f20002")],"PTID")) %>%
  filter(!is.na(noncancer_illness_code_selfreported_f20002_0_0))            # select columns containing self-report data
emptycols <- sapply(sr, function (k) all(is.na(k)))                         # select empty columns
sr <- sr[!emptycols]                                                        # remove empty columns
temp = sapply(sr[1:17],function(x) grepl("1240", x, ignore.case = TRUE))    # select for TBI self-report (1240)
sr$x1240 = rowSums(temp) > 0
tbi.sr <- sr[sr$x1240==TRUE,"PTID"]                                         # select PTID of subjects self-reporting
rm(temp,sr,emptycols)                                                       # remove temporary data

# select TBI subjects from ICD10 -> (tbi.icd) ####
icd <- u.batch %>%
  select(c(names(u.batch)[startsWith(names(u.batch),"diagnoses_icd10")],    # select columns containing ICD codes
         names(u.batch)[startsWith(names(u.batch),"diagnoses_icd9")],"PTID")) %>%
  filter(!is.na(diagnoses_icd10_f41270_0_0)) %>%                            # remove subjects/timwpoints with no reported ICD-10 codes
  filter(PTID %in% df$PTID)                                                 # selectsubjects included in the ROI dataset
emptycols <- sapply(icd, function (k) all(is.na(k)))                        # select empty columns
icd <- icd[!emptycols]                                                      # remove empty columns
temp = sapply(icd[1:length(icd)-1],function(x) grepl("S06", x, ignore.case = TRUE))
icd$S06 = rowSums(temp) > 0
tbi.icd <- icd[icd$S06==TRUE,"PTID"]
rm(temp,emptycols)                                                          # remove temporary data

# combine list of subjects from self-report and ICD10 -> (tbi.sub) -> (df$TBI) ####
tbi.sub <- c(tbi.sr,tbi.icd)                                                # combine TBI subjects from self-report and ICD-10 definitions
tbi.sub <- tbi.sub[!duplicated(tbi.sub)]                                    # remove duplicates
df <- df %>%
  mutate(TBI=ifelse(PTID%in%tbi.sub,1,0))

# add deep-learning based ICV calculation ####
icv <- icv.orig %>%
  mutate(ID= substr(ID, 1, 7)) %>%
  dplyr::rename(BATCH= Batch)
df <- merge(df, icv)
}
```

### remove ICD-10 confounders

```{r data_prep_3, echo = TRUE, cache = TRUE, message=FALSE}
if (fromBeginning==TRUE) {
##### - select subjects with other pathologies from control group (ICD10 codes beginning with F&G) -> (icd.sub) ####
temp = sapply(icd[1:(length(icd)-2)], function(x) grepl("G", x, ignore.case = TRUE))
icd$G = rowSums(temp) > 0                  # sums all ICD-10 codes containing "G" for each subject
icd.G <- icd[icd$G==TRUE,"PTID"]           # N = 1576 subjects have an ICD-10 code containing "G"
temp = sapply(icd[1:(length(icd)-3)],function(x) grepl("F", x, ignore.case = TRUE))
icd$xF = rowSums(temp) > 0                 # sums all ICD-10 codes containing "F" for each subject
icd.F <- icd[icd$xF==TRUE,"PTID"]          # N = 1229 subjects have an ICD-10 code containing "F"
icd.sub <- c(icd.G,icd.F)                  # combines subjects with "G" and "F" ICD-10 diagnoses
icd.sub <- icd.sub[!duplicated(icd.sub)]   # N = 2568 unique subjects have ICD-10 diagnoses containging "G" or "F"
rm(icd.F,icd.G,tbi.sr,tbi.icd,temp)        # removes temporary datas

# prepare datasets, selecting controls in (ctrl) and cases in (case) ####
ctrl <- df %>% 
  filter(TBI!=1) %>%
  filter(!PTID %in% icd.sub)               # removes subjects with other pathologies from control group
case <- df %>% filter(TBI==1) %>%
  mutate(range.low= Age_At_Visit -1) %>%
  mutate(range.high= Age_At_Visit +1)
}
```

### calculate z-score for each ROI

```{r data_prep_4, echo = TRUE, cache = TRUE, message=FALSE}
if (fromBeginning==TRUE) {
# select age/sex matched controls and calculate z-scores -> (df.z) and (df.z_ctrl) ####
for (row in 1:nrow(case)) {
  sub_data <- case[row,]
  sub_data.muse <- sub_data %>%
    select(c(names(df)[startsWith(names(df),"MUSE")]))
  controls <- ctrl %>%
    filter(Age_At_Visit>=sub_data$range.low) %>%
    filter(Age_At_Visit<=sub_data$range.high) %>%
    filter(Sex== sub_data$Sex)
  sub_controls <- controls %>% 
    select(PTID)
  muse_controls <- controls %>%
    select(c(names(df)[startsWith(names(df),"MUSE")]))
  mean <- sapply(muse_controls, function(k) mean(k))
  sd <- apply(muse_controls, 2, sd)
  nm <- rbind(mean,sd)
  nm.d <- as.data.frame(nm)
  z <- rbind(nm.d,sub_data.muse)
  zz <- sapply(names(z), function(x) (z[3,x] -z[1,x])/z[2,x])
  names(zz) <- paste0(names(z),"_z")
  zz <- as.data.frame(t(zz))
  PTID <- sub_data[,"PTID"]
  age <- sub_data %>% select(Age_At_Visit)
  sex <- sub_data %>% select(Sex)
  N.controls <- nrow(sub_controls)
  if (row==1) {
    df.z <- data.frame(PTID, age, sex, N.controls, zz)
  } else {
    df.z <- rbind(df.z,data.frame(PTID, age, sex, N.controls, zz))
  }
  z_ctrl <- controls %>%
    select(c(PTID, names(df)[startsWith(names(df),"MUSE")]))
  z_ctrl[cols <- c(names(df)[startsWith(names(df),"MUSE")])] <- 
    scale(z_ctrl[cols <- c(names(df)[startsWith(names(df),"MUSE")])])
  z_ctrl$sbj.info <- paste0(sub_data$PTID,"_control")
  names(z_ctrl)[startsWith(names(z_ctrl),"MUSE")] <- paste0(names(z_ctrl)[startsWith(names(z_ctrl),"MUSE")],"_z")
  if (row==1) {
    df.z_ctrl <- data.frame(z_ctrl)
  } else {
    df.z_ctrl <- rbind(df.z_ctrl,z_ctrl)
  }
}
 rm(nm,age,muse_controls,nm.d,sub_data,sub_data.muse,z,zz,PTID,sd,mean,N.controls,row,sub_controls,controls,sex, z_ctrl)
}
```

### ICV-corrected z-scores

```{r data_prep_5, echo = TRUE, cache = TRUE, message=FALSE}
if (fromBeginning==TRUE) {
# select age/sex matched controls and calculate icv corrected z-scores -> (df.z.icv) and (df.z_ctrl.icv) ####
for (row in 1:nrow(case)) {
  sub_data <- case[row,]
  sub_data.muse <- sub_data %>%
    select(c(names(df)[startsWith(names(df),"MUSE")], DLICV))
  sub_data.muse.icv <- as.data.frame(t(sapply(names(sub_data.muse), function(x) sub_data.muse[,x]/sub_data.muse[,"DLICV"])))
  controls <- ctrl %>%
    filter(Age_At_Visit>=sub_data$range.low) %>%
    filter(Age_At_Visit<=sub_data$range.high) %>%
    filter(Sex== sub_data$Sex)
  sub_controls <- controls %>% 
    select(PTID)
  muse_controls <- controls %>%
    select(c(names(df)[startsWith(names(df),"MUSE")], DLICV))
  muse_controls.icv <- as.data.frame(sapply(names(muse_controls), function(x) muse_controls[,x]/muse_controls[,"DLICV"]))
  mean <- sapply(muse_controls.icv, function(k) mean(k))
  sd <- sapply(muse_controls.icv, function(m) sd(m))
  nm <- rbind(mean,sd)
  nm.d <- as.data.frame(nm)
  z <- rbind(nm.d,sub_data.muse.icv)
  zz <- sapply(names(z), function(x) (z[3,x] -z[1,x])/z[2,x])
  names(zz) <- paste0(names(z),"_z")
  zz <- as.data.frame(t(zz))
  PTID <- sub_data[,"PTID"]
  age <- sub_data %>% select(Age_At_Visit)
  sex <- sub_data %>% select(Sex)
  N.controls <- nrow(sub_controls)
  if (row==1) {
    df.z.icv <- data.frame(PTID, age, sex, N.controls, zz)
  } else {
    df.z.icv <- rbind(df.z.icv,data.frame(PTID, age, sex, N.controls, zz))
  }
  z_ctrl <- controls %>%
    select(c(PTID, names(df)[startsWith(names(df),"MUSE")], DLICV))
  z_ctrl.icv_tmp <- as.data.frame(sapply(names(df)[startsWith(names(df),"MUSE")], 
                                                      function(x) z_ctrl[,x]/z_ctrl[,"DLICV"]))
  z_ctrl.icv <- cbind(PTID= z_ctrl$PTID, z_ctrl.icv_tmp)
  z_ctrl.icv[cols <- c(names(df)[startsWith(names(df),"MUSE")])] <- 
    scale(z_ctrl.icv[cols <- c(names(df)[startsWith(names(df),"MUSE")])])
  z_ctrl.icv$sbj.info <- paste0(sub_data$PTID,"_control")
  names(z_ctrl.icv)[startsWith(names(z_ctrl.icv),"MUSE")] <- paste0(names(z_ctrl.icv)[startsWith(names(z_ctrl.icv),"MUSE")],"_z")
    if (row==1) {
    df.z_ctrl.icv <- data.frame(z_ctrl.icv)
  } else {
    df.z_ctrl.icv <- rbind(df.z_ctrl.icv,z_ctrl.icv)
  }
}
rm(nm,age,muse_controls,nm.d,sub_data,sub_data.muse,z,zz,PTID,sd,mean,N.controls,row,sub_controls,controls,sex, z_ctrl)
}
```

### de-identify z-scores

```{r data_prep_6, echo = TRUE, cache = TRUE, message=FALSE}
if (fromBeginning==TRUE) {
# remove any subject-specific identifiers
df.z_503 <- df.z[,c(names(df.z)[startsWith(names(df.z),"MUSE")])]
df.z.icv_503 <- df.z.icv[,c(names(df.z.icv)[startsWith(names(df.z.icv),"MUSE")])]
df.z_ctrl_503 <- df.z_ctrl[,c(names(df.z_ctrl)[startsWith(names(df.z_ctrl),"MUSE")])]
df.z_ctrl.icv_503 <- df.z_ctrl.icv[,c(names(df.z_ctrl.icv)[startsWith(names(df.z_ctrl.icv),"MUSE")])]

# write z-score datasets, to be used in Final Project
write.csv(df.z_503, paste0(home,"df_z_503.csv"), row.names=FALSE)
write.csv(df.z.icv_503, paste0(home,"df_z_icv_503.csv"), row.names=FALSE)
write.csv(df.z_ctrl_503, paste0(home,"df_z_ctrl_503.csv"), row.names=FALSE)
write.csv(df.z_ctrl.icv_503, paste0(home,"df_z_ctrl_icv_503.csv"), row.names=FALSE)
 
} else {
# read z-score datasets, to be used in Final Project
df.z <- read.csv(paste0(home,"df_z_503.csv"), stringsAsFactors = FALSE)
df.z.icv <- read.csv(paste0(home,"df_z_icv_503.csv"), stringsAsFactors = FALSE)
df.z_ctrl <- read.csv(paste0(home,"df_z_ctrl_503.csv"), stringsAsFactors = FALSE)
df.z_ctrl.icv <- read.csv(paste0(home,"df_z_ctrl_icv_503.csv"), stringsAsFactors = FALSE)
}
```

# Analysis

## summarize ROI abnormality

```{r analysis_1, echo = TRUE, cache = TRUE, message=FALSE}
##### select most affected absolute ROIs by looking z-scores less than -2, skewness, kurtosis
df.z.o <- df.z
temp <- colSums(sapply(df.z.o, function(x) x<=-2))
skew <- sapply(df.z.o, moments::skewness)
kurt <- sapply(df.z.o, moments::kurtosis)

##### select most affected ICV-corrected ROIs by looking z-scores less than -2, skewness, kurtosis
df.z.icv.o <- df.z.icv
temp.icv <- colSums(sapply(df.z.icv.o, function(x) x<=-2))
skew.icv <- sapply(df.z.icv.o, moments::skewness)
kurt.icv <- sapply(df.z.icv.o, moments::kurtosis)

#ROI summary measures: N.outliers, skewness, kurtosis
roi_summary <- rbind(temp, skew, kurt, temp.icv, skew.icv, kurt.icv)
```

## name ROIs

```{r analysis_2, echo = TRUE, cache = TRUE, message=FALSE}
# read MUSE ROI dictionary (muse.dict.t) ####
muse.dict <- read.csv(paste0(home,"doc_MUSE_ROI_Dictionary.csv"),stringsAsFactors = FALSE)
muse.dict$ROI_NAME_INDEX <- paste0(muse.dict$ROI_INDEX,"-",muse.dict$ROI_NAME)
muse.dict.t <- as.data.frame(t(muse.dict))
names(muse.dict.t) <- paste0(muse.dict.t["ROI_COL",],"_z")
muse.dict.t <- muse.dict.t %>%
  select(-c(names(muse.dict.t)[!names(muse.dict.t)%in%names(df.z)]))

df.z.dict <- rbind(df.z,muse.dict.t)
df.z.icv.dict <- rbind(df.z.icv,muse.dict.t) 

roi.sel <- df.z.dict[df.z.dict$PTID=="ROI_NAME_INDEX",]
rename <- as.data.frame(c(roi.sel[,1:ncol(roi.sel)]))
df_z <- rbind(df.z, rename)
names(df_z) <- df_z[44,]
df_z <- df_z[-44,] %>%
  sapply(as.numeric) 
df_z.icv <- rbind(df.z.icv, rename)
names(df_z.icv) <- df_z.icv[44,]
df_z.icv <- df_z.icv[-44,] %>%
  sapply(as.numeric)
df_z <- as.data.frame(df_z)
df_z.icv <- as.data.frame(df_z.icv)

roi_summary_ <- rbind(roi_summary,ROI= rename[1:ncol(rename)])
roi.nm <- colnames(roi_summary)
roi_summary_a <- rbind(roi_summary, roi.nm)
```

# Results - ROI specific

The first section of results looks at abnormality in the distribution of normative z-scores within each 
ROI in the TBI case group. Abnormality is defined in terms of outliers (with normative z-score less than negative 2), 
and meausures of gaussianity includeing skewness and kurtosis.

## all ROIs - ICV corrected

Summary of normative z-scores in the TBI case group for all ROIs.

```{r results_4, echo = TRUE, cache = TRUE, message=FALSE}
library(summarytools)
library(skimr)

skim(df.z.icv)
```

## outliers

```{r results_1, echo = TRUE, cache = TRUE, message=FALSE}
library(knitr)
roi <- as.data.frame(t(roi_summary_a)) %>%
  arrange(desc(temp.icv)) %>%
  mutate(ROI= roi.nm, N_outliers= as.numeric(temp.icv), Skew= as.numeric(skew.icv), Kurtosis= as.numeric(kurt.icv)) %>% 
  select(ROI, N_outliers, Skew, Kurtosis)

roi[,3] <- round(roi[,3], 3)  # rounds skew to 3 significant figures
roi[,4] <- round(roi[,4], 3)  # rounds kurtosis to 3 significant figures

muse.dict$ROI <- paste0(muse.dict$ROI_COL,"_z")
roi20 <- roi[1:20,]
muse.dict.names <- muse.dict %>% select(c(ROI_NAME, ROI))
roi20.dict <- merge(roi20,muse.dict.names,by="ROI",all.x=TRUE)
kable(roi20.dict, caption= "Top 20 ROIs with Most Outliers among Case Subjects")
```

### Thalamus (L)

```{r, figures-side_60, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_60_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Left Thalamus")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_60_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Left Thalamus")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

### Frontal Lobe

```{r, figures-side_504, fig.show="hold", out.width="50%"}
ggplot(df.z.icv, aes(MUSE_Volume_504_z)) +
                ylim(0, 15) +
                labs(title = "Histograms: z-scored ICV-corrected ROI volumes in TBI case subjects", x= paste0("z-score of ICV-corrected Frontal Lobe")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_504_z)) + 
                ylim(0, 15) +
                labs(title = "Histograms: z-scored absolute ROI volumes in TBI case subjects", x= paste0("z-score of Frontal Lobe")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

### Parietal Lobe (R)

```{r, figures-side_523, fig.show="hold", out.width="50%"}
ggplot(df.z.icv, aes(MUSE_Volume_523_z)) +
                ylim(0, 15) +
                labs(title = "Histograms: z-scored ICV-corrected ROI volumes in TBI case subjects", x= paste0("z-score of ICV-corrected Right Parietal Lobe")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_523_z)) + 
                ylim(0, 15) +
                labs(title = "Histograms: z-scored absolute ROI volumes in TBI case subjects", x= paste0("z-score of Right Parietal Lobe")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)

```

## skew 

```{r results_2, echo = TRUE, cache = TRUE, message=FALSE}
roi.skew <- roi %>% arrange(Skew)
roi20.skew <- roi.skew[1:20,]
roi20.skew.dict <- merge(roi20.skew,muse.dict.names,by="ROI",all.x=TRUE)
kable(roi20.skew.dict, caption= "Top 20 ROIs with Highest Skewness among Case Subjects")
```

### Angular Gyrus (R)

```{r, figures-side_106, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_106_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Angular Gyrus (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_106_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Angular Gyrus (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

### Fusiform Gyrus (R)

```{r, figures-side_122, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_122_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Fusiform Gyrus (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_122_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Fusiform Gyrus (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

### Internal Capsule, Anterior Limb (R)

```{r, figures-side_91, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_91_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Anterior Limb of Internal Capsule (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_91_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Anterior Limb of Internal Capsule (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

### Ventral Diencephalon (R)

```{r, figures-side_61, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_61_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Ventral Diencephalon (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_61_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Ventral Diencephalon (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

## outliers and skew

```{r results_3, echo = TRUE, cache = TRUE, message=FALSE}
roi.mix <- roi %>%
    arrange(roi$N_outliers*roi$Skew)
roi20.mix <- roi.mix[1:20,]
roi20.mix.dict <- merge(roi20.mix,muse.dict.names,by="ROI",all.x=TRUE)
kable(roi20.mix.dict, caption= "Top 20 ROIs with Most Outliers and Most Negative Skewness among Case Subjects")
```

### Fornix (L)

```{r, figures-side_90, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_90_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Fornix (left hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_90_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Fornix (left hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

### Fornix (R)

```{r, figures-side_89, fig.show="hold", out.width="50%"}
library(ggplot2)
ggplot(df.z.icv, aes(MUSE_Volume_89_z)) +
                ylim(0, 15) +
                labs(title = "ICV-corrected volumes", x= paste0("z-score of ICV-corrected Fornix (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "midnightblue", fill = "midnightblue", alpha = 0.7)

ggplot(df.z, aes(MUSE_Volume_89_z)) + 
                ylim(0, 15) +
                labs(title = "absolute volumes", x= paste0("z-score of Fornix (right hemisphere)")) +
                theme(plot.title = element_text(hjust = 0.5)) +
                geom_histogram(breaks = seq(-3.5, 3.5, .5), color = "blue", fill = "blue", alpha = 0.7)
```

# Results - Whole Brain

```{r results_wholeBrain, echo = FALSE, cache = TRUE, message=FALSE}
# keep original datasets
df.z.o <- df.z
df.z.icv.o <- df.z.icv
df.z_ctrl.o <- df.z_ctrl 
df.z_ctrl.icv.o <- df.z_ctrl.icv 

select_outlier <- 2     #select z-score threshold for abnormality

outlier_function <- function(inp) {
  outliers_negative <- inp[sapply(inp, function(x) x<= -select_outlier)]
  N_outliers_neg <- length(outliers_negative)
  outliers_positive <- inp[sapply(inp, function(x) x>= select_outlier)]
  N_outliers_pos <- length(outliers_positive)
  
  o_pos <- round(sum(outliers_positive),3)
  o_pos_sd <- round(sd(outliers_positive),3)
  o_neg <- round(abs(sum(outliers_negative)),3)
  o_neg_sd <- round(sd(outliers_negative),3)
  o <- o_neg + o_pos
  N_sbj <- nrow(inp)
  p_pos <- round(o_pos/N_sbj,3)
  p_neg <- round(o_neg/N_sbj,3)
  p <- round(o/N_sbj,3)

  out <- format(c(N_sbj, paste0(p_pos," (", o_pos_sd,")"), paste0(p_neg," (", o_neg_sd,")"), p), scientific=FALSE, drop0trailing = TRUE)
  return(out)
}

case <- outlier_function(inp= df.z.o)
case.icv <- outlier_function(inp= df.z.icv.o)
ctrl <- outlier_function(inp= df.z_ctrl.o)
ctrl.icv <- outlier_function(inp= df.z_ctrl.icv.o)

names <- c("N", "z-scores greater than 2, subject average (standard deviation)", "z-scores less than 2, subject average (standard deviation)", "total z-scores outside +/- 2, subject average")
outliers <- cbind(names, case, ctrl)
outliers.icv <- cbind(names, case.icv, ctrl.icv)

kable(outliers.icv,caption=  "Average Outliers in TBI cases vs. Controls, ICV-corrected")
kable(outliers,caption=  "Average Outliers in TBI cases vs. Controls")

```

## Shiny visualization

For more ROI-specific information, use this interactive shiny app to visualize histograms of normative
z-score distribution in TBI cases for all ROIs. The purple histogram on top shows the intracranial 
volume (ICV) corrected z-scores. The blue histogram on bottom shows the z-score calculated from
absolute ROI volumes.

- Page 1: top 20 most outliers (z-score < -2)
- Page 2: top 20 most negative skewness
- Page 3: top 20 lowest kurtosis
- Page 4: top 20 lowest score calculated by N outliers * skewness
- Page 5: includes all ROIs
- Page 6: Whole Brain abnormality

https://himeslab.dev/rsconnect/content/15

# Conclusion

In this project, I was able to address my aims: 
- to quantify structural brain volumes in TBI case group,
- to identify regions with disease-related atrophy, and 
- to compare the relative distribution of brain volume in TBI stubjects against normal controls.

I found that there is a significantly different distribution of brain volume in the TBI case group 
in comparison to the control group, with a subject average of 5 more standard deviations below the 
normative ICV-corrected volumetric mean. This measure is informative, as it quantifies the global
abnormality, or whole brain atrophy resulting from TBI. This measure will be especially interesting 
to look at longitudinally, to see the temporal resolution of global atrophy following a traumatic 
brain injury.

One limitation of this study is the lack of temporal information; I did not have any information
about the time that had elapsed since the TBI event. Since the neurodegeneration resulting from
TBI does not happen immediately, temporal information will be critical to a final analysis.

Another limitation of this study is the low number of TBI case subjects. Since the UK Biobank
is a general population study, and TBI subjects were identified primarily by ICD-10 code,
the number of case subjects is very limited. In this project, the large number of control subjects
(between N= 349 and N= 1,402 per case subject) makes up for the low number of cases, but
further validation will need to be done, especially on the ROI specific measures, to ensure
that the difference is truly a disease-related change and not a subject-specific abnormality.

For future work, I would like to extend this study to the TBI clinic at University of 
Pennsylvania. There are a standard MRI protocol for TBI cases at Penn, including a high resolution T1
MPRAGE sequence, which should allow for the extension of this work into a clinical dataset. 


# Acknowledgements

I would like to thank the following individuals for their guidance and support in this project:
  
Christos Davatzikos, PhD, Wallace T. Miller Sr. Professor of Radiology & Electrical and Systems Engineering at
the Perelman School of Medicine         
Haochang Shou, PhD, Assistant Professor of Biostatistics in Biostatistics and Epidemiology at the Perelman School
of Medicine       
Blanca Himes, PhD, Associate Professor of Informatics in Biostatistics and Epidemiology at the Perelman School
of Medicine

And all of my classmates in Data Science for Biomedical Infomatics, Fall 2020, for their
valuable feedback.

# References

[1] https://www.cdc.gov/traumaticbraininjury/pubs/tbi_report_to_congress.html       
[2] McKee, Ann C., et al. “The spectrum of disease in chronic traumatic encephalopathy.” 
Brain 136.1 (2013): 43-64.        
[3] Doshi, Jimit, et al. "MUSE: MUlti-atlas region Segmentation utilizing Ensembles of 
registration algorithms and parameters, and locally optimal atlas selection." Neuroimage 127 
(2016): 186-195.