DWI QC

Rmarkdown notebook for examining DWI QC.

We’re starting by using the “eddy_outlier_n_stdev_map” from the EDDY documentation: “The numbers denote how many standard deviations off the mean difference between observation and prediction is.” By “the prediction”, I’m assuming they’re talking about the diffusion model. This is similar (in principle to the residuals of the tensor model). In short, people who are “noisy” tend to have more consistently higher SD values across all volumes rather than a few.

In examining these data, I found that collapsing the data into volumes/directions by taking the standard deviation across slices worked well. This essentially asks what is the mean distance of each slice from every other slice? As can be seen below, people with noisier brains are more similar to each other (smaller distance from each other - see distance matrix).

Data are given per person in a matrix where each column represents a slice, and each row represents a scan. Therefore for the POND data, we have a 69 (scans) x 140 (slices). The following code will need to be re-written to be flexibly applied to all datasets (i.e. scan and volume numbers are hard-coded at the moment).

Dependencies:

library(tidyverse)  # data manipulation
library(cluster)    # clustering algorithms
library(factoextra) # clustering algorithms & visualization
library(reshape2)
library(ggplot2)
library(gplots) 
library(readbulk)
library(ggthemr)
library(ggpubr)
library(doBy)
library(Rmisc)
library(dendextend)
library(dendextendRcpp)
library(kableExtra)
ggthemr("chalk")

#function to prepare SD files for analysis
QC_DWI_Prep <- function(x){
  input <- read_bulk(x, skip = 1, header =FALSE, sep =' ')
  input$V141 <- NULL
  data <- melt(input, id.vars = "File")
  #this function (dcast) combines the data by calculating an SD per volume across slices.  
  #Thus, we have a measure of how variable the variability is per volume.  
  b <- dcast(data, variable ~ File, fun.aggregate = sd, na.rm = TRUE)
  b$variable <- NULL
  d <- t(b)
  ds <- scale(d)
  my_list <- list(d, ds)
return(my_list)
}

First read in the data using “readbulk”:

First just look at the distance matrix

#distance matrix visualization
distance <- get_dist(ds)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

Now try Ward’s method for clustering

# Ward Hierarchical Clustering
d <- dist(d, method = "euclidean") # distance matrix
fit <- hclust(d, method="ward") 
groups <- cutree(fit, k=5) # cut tree into 5 clusters
# draw dendogram with red borders around the 5 clusters 
# using piping to get the dend
dend <- fit %>% as.dendrogram

# plot + color the dend's branches before, based on 3 clusters:


par(mar = c(2,2,2,15))

dend %>% color_branches(k=5) %>% plot(horiz=TRUE)

# add horiz rect
#dend %>% rect.dendrogram(k=5,horiz=TRUE)

# add horiz (well, vertical) line:
text(50, 10, "bad \n participants", col = "#C84D6F")
text(50, 27, "good \n participants", col = "#9C4DCC")
text(50, 17, "?")

Now try a K-means solution

#Now to try a K-means solution 

k2 <- kmeans(ds, centers = 2, nstart = 25)
k3 <- kmeans(ds, centers = 3, nstart = 25)
k4 <- kmeans(ds, centers = 4, nstart = 25)
k5 <- kmeans(ds, centers = 5, nstart = 25)
k10 <- kmeans(ds, centers = 10, nstart = 25)

# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = ds) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = ds) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = ds) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = ds) + ggtitle("k = 5")
#p10 <- fviz_cluster(k10, geom = "point",  data = d) + ggtitle("k = 10")

library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow = 2)

fviz_nbclust(ds, kmeans, method = "silhouette")

fviz_nbclust(ds, kmeans, method = "wss")

gap_stat <- clusGap(ds, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)
fviz_gap_stat(gap_stat)

#Fuzzy K-Means clustering (2 groups) This approximates what Hajer is doing I think (trying to make a binary decision, but providing us with probabities).  I'd look at any person who gets a score of less than .9 for either category.
fannyd <- fanny(ds, k=2, metric = "euclidean", memb.exp = 1.2)
round(fannyd$membership, 2)

FALSE                                 [,1] [,2]
FALSE sub-0880043_outlier_n_stdev_map 1.00 0.00
FALSE sub-0880107_outlier_n_stdev_map 0.81 0.19
FALSE sub-0880138_outlier_n_stdev_map 1.00 0.00
FALSE sub-0880397_outlier_n_stdev_map 0.52 0.48
FALSE sub-0880418_outlier_n_stdev_map 1.00 0.00
FALSE sub-0880464_outlier_n_stdev_map 1.00 0.00
FALSE sub-0880473_outlier_n_stdev_map 0.24 0.76
FALSE sub-0880533_outlier_n_stdev_map 0.98 0.02
FALSE sub-0880558_outlier_n_stdev_map 0.99 0.01
FALSE sub-0880601_outlier_n_stdev_map 0.99 0.01
FALSE sub-0880624_outlier_n_stdev_map 1.00 0.00
FALSE sub-0880703_outlier_n_stdev_map 0.36 0.64
FALSE sub-1050015_outlier_n_stdev_map 0.32 0.68
FALSE sub-1050019_outlier_n_stdev_map 0.09 0.91
FALSE sub-1050032_outlier_n_stdev_map 0.26 0.74
FALSE sub-1050054_outlier_n_stdev_map 0.15 0.85
FALSE sub-1050065_outlier_n_stdev_map 1.00 0.00
FALSE sub-1050090_outlier_n_stdev_map 1.00 0.00
FALSE sub-1050105_outlier_n_stdev_map 0.98 0.02
FALSE sub-1050131_outlier_n_stdev_map 0.21 0.79
FALSE sub-1050179_outlier_n_stdev_map 0.04 0.96
FALSE sub-1050235_outlier_n_stdev_map 0.99 0.01
FALSE sub-1050250_outlier_n_stdev_map 0.20 0.80
FALSE sub-1050253_outlier_n_stdev_map 0.03 0.97
FALSE sub-1050349_outlier_n_stdev_map 0.79 0.21
FALSE sub-1050353_outlier_n_stdev_map 0.02 0.98
FALSE sub-1050378_outlier_n_stdev_map 0.02 0.98
FALSE sub-1050402_outlier_n_stdev_map 0.99 0.01
FALSE sub-1050429_outlier_n_stdev_map 0.08 0.92
FALSE sub-1050431_outlier_n_stdev_map 1.00 0.00

#Clustering cluster sets with Jaccard index - how similar are the cluster solutions?
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/clusterIndex.R") 
clVlist <- lapply(2:12, function(x) clara(ds, k=x)$clustering); names(clVlist) <- paste("k", "=", 2:12)
e <- sapply(names(clVlist), function(x) sapply(names(clVlist), function(ds) cindex(clV1=clVlist[[ds]], clV2=clVlist[[x]], method="jaccard")[[3]]))
hv <- hclust(as.dist(1-e))
plot(as.dendrogram(hv), edgePar=list(col=3, lwd=4), horiz=T, main="Similarities of 10 Clara Clustering Results for k: 2-12") 

So, it looks like 2-3 clusters are optimal

Let’s plot the 3-groups from the hclust (Ward’s solution) - we can try other versions later

#should likely re-write this to take the clusters from the other solutions rather than hard-coding it to be subdirectories (use the %in% function from dplyr)

group3 <-"/Users/John/Dropbox/qc_solution/group3"
group2 <- "/Users/John/Dropbox/qc_solution/group2"
group1 <-"/Users/John/Dropbox/qc_solution/group1"

QC_DWI_Plot <- function(x){
input <- read_bulk(x, skip = 1, header =FALSE, sep =' ')
input$V141 <- NULL
rows <- as.data.frame(rep(c(0:68), times = length(unique(input$File))))
colnames(rows) <- "rows"
data <- cbind(input, rows)
group_L <- melt(data, id.vars = c("File","rows"))
group_L$variable <- as.numeric(substring(group_L$variable, 2))
return(group_L)
}


group1_L <- QC_DWI_Plot(group1)
group2_L <- QC_DWI_Plot(group2)
group3_L <- QC_DWI_Plot(group3)
all_L <- QC_DWI_Plot("/Users/John/Dropbox/qc_solution/qc_stdev_files")

group1_plot <- ggplot(group1_L, aes(rows, variable)) +
geom_raster(aes(fill = value))+facet_grid(.~File)+ 
  scale_y_continuous(breaks=seq(0,140,10))+
  labs(x = "volume|directions", y = "slice")
group2_plot <- ggplot(group2_L, aes(rows, variable)) +
geom_raster(aes(fill = value))+facet_grid(.~File)+ 
  scale_y_continuous(breaks=seq(0,140,10))+
  labs(x = "volume|directions", y = "slice")
group3_plot <- ggplot(group3_L, aes(rows, variable)) +
geom_raster(aes(fill = value))+facet_grid(.~File)+ 
  scale_y_continuous(breaks=seq(0,140,10))+
  labs(x = "volume|directions", y = "slice")


ggarrange(group1_plot, group2_plot, group3_plot, nrow=3, ncol =1, labels = c("Group1", "Group2","Group3"), common.legend = TRUE, widths = c(1, 0.64, 0.5))

ggplot() + 
  geom_histogram(data = group1_L, aes(x = value), color = "green", fill = "green") +
  geom_histogram(data = group2_L, aes(x = value), color = "blue", fill = "blue") + 
  geom_histogram(data = group3_L, aes(x = value), color = "red", fill = "red")+
  xlim(-3,3)

Trainin a classifier?

library(caret)
library(MASS)

pond_all_training <- read_csv("/Users/John/Dropbox/qc_solution/pond_all_training.csv")

TestRows     <- subset(pond_all_training, Training =="No")
TrainData    <- subset(pond_all_training, Training =="Yes")
TrainClasses <- as.data.frame(TrainData[,142])
TrainData    <- TrainData[,2:141]
TestData     <- TestRows[,2:141]
TestSubs     <- TestRows[,1]
#TestClasses  <- iris[TestRows,5]

TrainClasses <- factor(TrainClasses$Group, labels = c("group_1", "group_2","group_3"))

nnetFit <- train(x=TrainData, y=TrainClasses,
                 method = "nnet",
                 preProcess = "range", 
                 tuneLength = 2,
                 trace = FALSE,
                 maxit = 100)

predicted <-as.data.frame(predict(nnetFit, TestData))
test_preds <- cbind(TestSubs, predicted)
kable(test_preds)

File	predict(nnetFit, TestData)
sub-0880002	group_1
sub-0880019	group_2
sub-0880020	group_2
sub-0880023	group_1
sub-0880044	group_1
sub-0880046	group_1
sub-0880049	group_1
sub-0880050	group_1
sub-0880061	group_1
sub-0880062	group_1
sub-0880079	group_1
sub-0880081	group_2
sub-0880082	group_1
sub-0880085	group_3
sub-0880096	group_1
sub-0880099	group_1
sub-0880100	group_1
sub-0880102	group_1
sub-0880112	group_1
sub-0880114	group_2
sub-0880136	group_1
sub-0880137	group_1
sub-0880139	group_1
sub-0880140	group_2
sub-0880146	group_3
sub-0880152	group_1
sub-0880155	group_1
sub-0880157	group_1
sub-0880167	group_1
sub-0880168	group_1
sub-0880175	group_1
sub-0880182	group_1
sub-0880195	group_1
sub-0880205	group_2
sub-0880208	group_2
sub-0880215	group_1
sub-0880221	group_1
sub-0880229	group_1
sub-0880247	group_2
sub-0880250	group_1
sub-0880252	group_1
sub-0880253	group_1
sub-0880254	group_1
sub-0880255	group_3
sub-0880258	group_1
sub-0880265	group_1
sub-0880276	group_1
sub-0880278	group_1
sub-0880279	group_1
sub-0880284	group_1
sub-0880291	group_2
sub-0880293	group_1
sub-0880303	group_1
sub-0880304	group_1
sub-0880305	group_1
sub-0880308	group_2
sub-0880323	group_3
sub-0880325	group_1
sub-0880328	group_2
sub-0880338	group_1
sub-0880349	group_1
sub-0880351	group_1
sub-0880353	group_1
sub-0880369	group_1
sub-0880370	group_1
sub-0880372	group_1
sub-0880381	group_1
sub-0880386	group_1
sub-0880396	group_1
sub-0880405	group_1
sub-0880409	group_1
sub-0880415	group_1
sub-0880419	group_1
sub-0880428	group_1
sub-0880432	group_2
sub-0880436	group_1
sub-0880447	group_1
sub-0880454	group_1
sub-0880461	group_1
sub-0880462	group_1
sub-0880470	group_1
sub-0880476	group_1
sub-0880481	group_1
sub-0880494	group_1
sub-0880514	group_2
sub-0880515	group_1
sub-0880536	group_1
sub-0880539	group_1
sub-0880552	group_1
sub-0880555	group_1
sub-0880582	group_1
sub-0880584	group_1
sub-0880590	group_1
sub-0880620	group_1
sub-0880642	group_1
sub-0880644	group_1
sub-0880650	group_1
sub-0880663	group_1
sub-0880664	group_1
sub-0880665	group_1
sub-0880669	group_1
sub-0880686	group_1
sub-0880687	group_3
sub-0880688	group_1
sub-0880689	group_1
sub-0880692	group_1
sub-0880693	group_1
sub-0880696	group_2
sub-0880720	group_1
sub-0880738	group_1
sub-0880742	group_1
sub-0880745	group_1
sub-0880759	group_2
sub-0880763	group_1
sub-0880770	group_3
sub-0880785	group_1
sub-0880794	group_3
sub-0880825	group_1
sub-0880856	group_2
sub-0880861	group_2
sub-0880862	group_3
sub-0880863	group_1
sub-0880864	group_1
sub-0880865	group_1
sub-0880869	group_1
sub-0880870	group_1
sub-0880871	group_1
sub-0880872	group_2
sub-0880873	group_3
sub-0880874	group_1
sub-0880881	group_1
sub-0880882	group_1
sub-0880895	group_2
sub-0880896	group_1
sub-0880898	group_1
sub-0880912	group_1
sub-0880913	group_1
sub-0880915	group_2
sub-0880927	group_1
sub-0880928	group_1
sub-0880929	group_1
sub-0880930	group_1
sub-0880931	group_2
sub-0880935	group_1
sub-0880938	group_1
sub-0880941	group_1
sub-0880960	group_2
sub-0880964	group_1
sub-0885004	group_1
sub-0885005	group_3
sub-1050004	group_1
sub-1050007	group_1
sub-1050026	group_1
sub-1050027	group_1
sub-1050029	group_2
sub-1050038	group_2
sub-1050039	group_1
sub-1050040	group_1
sub-1050051	group_2
sub-1050053	group_1
sub-1050055	group_1
sub-1050056	group_1
sub-1050064	group_1
sub-1050070	group_1
sub-1050077	group_2
sub-1050089	group_1
sub-1050092	group_1
sub-1050093	group_1
sub-1050097	group_1
sub-1050098	group_1
sub-1050100	group_1
sub-1050101	group_1
sub-1050102	group_1
sub-1050108	group_1
sub-1050111	group_1
sub-1050112	group_1
sub-1050125	group_3
sub-1050128	group_1
sub-1050132	group_1
sub-1050134	group_1
sub-1050135	group_1
sub-1050136	group_1
sub-1050137	group_1
sub-1050142	group_2
sub-1050145	group_1
sub-1050150	group_1
sub-1050158	group_1
sub-1050165	group_1
sub-1050167	group_2
sub-1050172	group_1
sub-1050173	group_1
sub-1050174	group_1
sub-1050178	group_1
sub-1050180	group_1
sub-1050182	group_1
sub-1050184	group_3
sub-1050188	group_3
sub-1050189	group_1
sub-1050192	group_1
sub-1050194	group_1
sub-1050195	group_1
sub-1050198	group_1
sub-1050207	group_3
sub-1050211	group_1
sub-1050212	group_3
sub-1050213	group_1
sub-1050215	group_2
sub-1050218	group_1
sub-1050219	group_1
sub-1050221	group_1
sub-1050222	group_1
sub-1050224	group_1
sub-1050225	group_1
sub-1050229	group_1
sub-1050230	group_1
sub-1050232	group_2
sub-1050234	group_1
sub-1050239	group_1
sub-1050241	group_1
sub-1050242	group_1
sub-1050247	group_1
sub-1050248	group_2
sub-1050249	group_1
sub-1050254	group_1
sub-1050257	group_1
sub-1050258	group_1
sub-1050261	group_1
sub-1050265	group_1
sub-1050268	group_1
sub-1050273	group_3
sub-1050277	group_1
sub-1050280	group_1
sub-1050281	group_1
sub-1050284	group_1
sub-1050286	group_1
sub-1050287	group_1
sub-1050288	group_1
sub-1050289	group_1
sub-1050292	group_1
sub-1050296	group_2
sub-1050298	group_2
sub-1050307	group_1
sub-1050308	group_1
sub-1050320	group_1
sub-1050321	group_3
sub-1050322	group_3
sub-1050325	group_2
sub-1050331	group_2
sub-1050334	group_1
sub-1050335	group_1
sub-1050336	group_1
sub-1050337	group_1
sub-1050340	group_1
sub-1050341	group_1
sub-1050342	group_1
sub-1050344	group_2
sub-1050345	group_1
sub-1050363	group_1
sub-1050368	group_1
sub-1050369	group_1
sub-1050371	group_2
sub-1050373	group_1
sub-1050375	group_3
sub-1050383	group_1
sub-1050385	group_1
sub-1050389	group_1
sub-1050393	group_2
sub-1050395	group_1
sub-1050398	group_1
sub-1050399	group_1
sub-1050401	group_1
sub-1050403	group_2
sub-1050404	group_3
sub-1050414	group_1
sub-1050420	group_1
sub-1050421	group_1
sub-1050422	group_1
sub-1050427	group_1
sub-1050428	group_1
sub-1050432	group_1
sub-1050438	group_1
sub-1050440	group_1
sub-1050445	group_3
sub-1050446	group_2
sub-1050447	group_1
sub-1050448	group_3
sub-1050451	group_1
sub-1050452	group_1
sub-1050457	group_1
sub-1050460	group_1
sub-1050461	group_1
sub-1050462	group_1
sub-1050467	group_3
sub-1050468	group_3
sub-1050486	group_1
sub-1050487	group_1
sub-1050488	group_1
sub-1050495	group_1
sub-1050496	group_2
sub-1050497	group_3
sub-1050500	group_1
sub-1050502	group_1
sub-1050505	group_1
sub-1050506	group_1
sub-1050509	group_1
sub-1050510	group_1
sub-1050512	group_3
sub-1050513	group_1
sub-1050516	group_1
sub-1050517	group_1
sub-1050518	group_3
sub-1050523	group_1
sub-1050525	group_3
sub-1050526	group_1
sub-1050528	group_2
sub-1050531	group_1
sub-1050532	group_1
sub-1050533	group_1
sub-1050543	group_1
sub-1050544	group_1
sub-1050545	group_3
sub-1050546	group_1
sub-1050547	group_3
sub-1050548	group_3
sub-1050551	group_1
sub-1050552	group_2
sub-1050556	group_1
sub-1050557	group_1
sub-1050561	group_1
sub-1050563	group_3
sub-1050570	group_2
sub-1050577	group_2
sub-1050579	group_3
sub-1050580	group_3
sub-1050582	group_1
sub-1050586	group_3
sub-1050591	group_3
sub-1050593	group_1
sub-1050596	group_1
sub-1050598	group_1
sub-1050600	group_1
sub-1050604	group_2
sub-1050607	group_2
sub-1050608	group_3
sub-1050613	group_1
sub-1050614	group_2
sub-1050615	group_1
sub-1050619	group_1
sub-1050621	group_1
sub-1050627	group_2
sub-1050628	group_1
sub-1050630	group_1
sub-1050633	group_1
sub-1050636	group_1
sub-1050639	group_3
sub-1050643	group_1
sub-1050644	group_1
sub-1050646	group_1
sub-1050647	group_1
sub-1050648	group_1
sub-1050649	group_1
sub-1050650	group_1
sub-1050655	group_2
sub-1050658	group_1
sub-1050659	group_1
sub-1050660	group_2
sub-1050662	group_1
sub-1050665	group_1
sub-1050666	group_3
sub-1050672	group_3
sub-1050680	group_1
sub-1050681	group_1
sub-1050687	group_1
sub-1050689	group_1
sub-1050698	group_1
sub-1050700	group_3
sub-1050701	group_3
sub-1050704	group_2
sub-1130131	group_3
sub-1130160	group_1
sub-1130202	group_1
sub-1130211	group_1
sub-1130217	group_1
sub-1130227	group_1
sub-1140078	group_3
sub-3940103	group_1
sub-3940114	group_1
sub-3940137	group_1
sub-3940147	group_1
sub-3940184	group_1
sub-3940185	group_1
sub-3940186	group_1
sub-3940189	group_1
sub-3940193	group_1
sub-3940195	group_1
sub-3940196	group_1
sub-3940203	group_1

table(TrainClasses, predict(nnetFit)) 

FALSE             
FALSE TrainClasses group_1 group_2 group_3
FALSE      group_1      14       0       0
FALSE      group_2       0       9       0
FALSE      group_3       0       0       7

#table(TestClasses,  predict(nnetFit,TestData))  

#visualization