# Lab 4 for the University of Tulsa's CS-6643 Bioinformatics Course
# Differential Expression
# Professor: Dr. McKinney, Fall 2022
# Noah L. Schrick - 1492657

## Set Working Directory to file directory - RStudio approach
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#### Part A: Preparing Data
load("sense.filtered.cpm.Rdata")
# load phenotype (mdd/hc) data
subject.attrs <- read.csv("Demographic_symptom.csv", 
                          stringsAsFactors = FALSE)

if (!require("dplyr")) install.packages("dplyr")
library(dplyr)
# grab intersecting X (subject ids) and Diag (Diagnosis) from columns
phenos.df <- subject.attrs %>% 
  filter(X %in% colnames(sense.filtered.cpm)) %>%
  dplyr::select(X, Diag)  
mddPheno <- as.factor(phenos.df$Diag)  

# Normalized and transform
if (!require("preprocessCore")) install.packages("preprocessCore")
library(preprocessCore)
mddExprData_quantile <- normalize.quantiles(sense.filtered.cpm)
mddExprData_quantileLog2 <- log2(mddExprData_quantile)
# attach phenotype names and gene names to data 
colnames(mddExprData_quantileLog2) <- mddPheno  
rownames(mddExprData_quantileLog2) <- rownames(sense.filtered.cpm)

length(rownames(sense.filtered.cpm))


#### Part B: Filter noise genes
# coefficient of variation filter sd(x)/abs(mean(x))
CoV_values <- apply(mddExprData_quantileLog2,1,
                    function(x) {sd(x)/abs(mean(x))})
# smaller threshold, the higher the experimental effect relative to the 
# measurement precision
sum(CoV_values<.045)
# there is one gene that has 0 variation -- remove
sd_values <- apply(mddExprData_quantileLog2,1, function(x) {sd(x)})
rownames(mddExprData_quantileLog2)[sd_values==0]  
# filter the data matrix 
GxS.covfilter <- mddExprData_quantileLog2[CoV_values<.045 & sd_values>0,]
dim(GxS.covfilter)


#### Part C: Differential Expression with t-tests
# convert phenotype
pheno.factor <- as.factor(colnames(GxS.covfilter))
pheno.factor
str(pheno.factor)
levels(pheno.factor)

## Run t-tests
myrow <- 2  # first pick a gene row index to test
mygene<-rownames(GxS.covfilter)[myrow]
mygene

# a. traditional R interface 
mdd <- GxS.covfilter[myrow,pheno.factor=="MDD"]
hc <- GxS.covfilter[myrow,pheno.factor=="HC"]
t.result <- t.test(mdd,hc)
t.result  

# b. formula interface ~ saves a step
t.result <- t.test(GxS.covfilter[myrow,] ~ pheno.factor)
t.result
p <- t.result$p.value
t.result$statistic

## Plot the Data
if (!require("ggplot2")) install.packages("ggplot2")
library(ggplot2)
# create data frame for gene
mygene.data.df <- data.frame(gene=GxS.covfilter[myrow,],phenotype=pheno.factor)
# boxplot
p <- ggplot(mygene.data.df, aes(x=phenotype, y=gene, fill=phenotype)) +
          stat_boxplot(geom ='errorbar') + geom_boxplot()
p <- p + xlab("MDD versus HC") + ylab(mygene)
p


#### Part D: t-test for all filtered genes

# Put it all together into a function to run in loop. 
# First write a function that computes t-test for one gene.
# i is the data row for the gene
ttest_fn <- function(i){
  mygene <- rownames(GxS.covfilter)[i]
  t.result <- t.test(GxS.covfilter[i,] ~ pheno.factor)
  tstat <- t.result$statistic
  pval <- t.result$p.value
  # return vector of three things for each gene    
  c(mygene, tstat, pval) 
} 

# Testing on the second gene
ttest_fn(2)

## Testing on the rest:
# initialize an empty matrix to store the results
ttest_allgene.mat <- matrix(0,nrow=nrow(GxS.covfilter), ncol=3)
# run analysis on all gene rows
for (i in 1:nrow(GxS.covfilter)){
  ttest_allgene.mat[i,] <- ttest_fn(i) 
}
# convert matrix to data frame and colnames
ttest_allgene.df <- data.frame(ttest_allgene.mat)
colnames(ttest_allgene.df) <- c("gene ", "t.stat", "p.val")

# sort based on p-value
ttest_allgene.sorted <- ttest_allgene.df                 %>% 
  mutate_at("p.val", as.character) %>%
  mutate_at("p.val", as.numeric)   %>%
  arrange(p.val) # sort 
ttest_allgene.sorted[1:10,] # look at top 10

## Plot the result of the top gene
# create data frame for gene
myrow <- which(ttest_allgene.df$gene==ttest_allgene.sorted[1,1])
mygene<-rownames(GxS.covfilter)[myrow]
mytopgene.data.df <- data.frame(mygene=GxS.covfilter[myrow,],
                                phenotype=pheno.factor)

# boxplot
p <- ggplot(mytopgene.data.df, aes(x=phenotype, y=mygene, fill=phenotype)) +
  stat_boxplot(geom ='errorbar') + geom_boxplot()
p <- p + xlab("MDD versus HC") + ylab(mygene)
p


#### Part E: Interpretation of top genes
top_cutoff <- 200
top_genes <- as.character(ttest_allgene.sorted[1:top_cutoff,1])
write.table(top_genes, sep="\t", file="", quote=F, row.names=F, col.names=F)


#### Optional 1: Identify Outliers
mddCorr<-cor(GxS.covfilter)  # distance based on correlation
d <- sqrt(1-mddCorr)
dim(d)
rownames(d)
mddTree = hclust(as.dist(d))
mddTree$labels <- phenos.df$X
plot(mddTree)


#### Optional 2: Compare MDS and UMAP clustering
if (!require("umap")) install.packages("umap")
library(umap)
# change umap config parameters
custom.config = umap.defaults
custom.config$random_state = 123
custom.config$n_epochs = 500

SxG.df <- data.frame(t(GxS.covfilter))

obs_mds = cmdscale(d, k=2)
#add colors for MDD/HC
colors = rep("black",nrow(SxG.df))
colors[startsWith(rownames(SxG.df),"MDD")] <- "red"
plot(obs_mds, col=colors, 
     main="mds of observations", xlab="mds dim1", ylab="mds dim2")

obs_umap = umap(SxG.df, config=custom.config)
#add colors for MDD/HC
colors = rep("black",nrow(SxG.df))
colors[startsWith(rownames(SxG.df),"MDD")] <- "red"
plot(obs_umap$layout, col=colors, 
     main="umap of observations", xlab="umap dim1", ylab="umap dim2")

#### Optional 3: WGCNA and UMAP of Genes
if (!require("BiocManager")) install.packages("BiocManager")
library(BiocManager)
if (!require("WGCNA")) BiocManager::install("WGCNA")
library(WGCNA)

# change umap config parameters
custom.config = umap.defaults
custom.config$random_state = 123
custom.config$n_epochs = 50
custom.config$n_neighbors=30
custom.config$metric = "pearson"
custom.config$input = "dist"
custom.config$a = 9.5
custom.config$b = 0.5

GxS.df <- data.frame(GxS.covfilter)

obs_umap = umap(GxS.df, config=custom.config)
#add colors for MDD/HC
colors = rep("black",nrow(GxS.df))
colors[startsWith(colnames(GxS.df),"MDD")] <- "red"
plot(obs_umap$layout, col=colors, 
     main="umap of observations", xlab="umap dim1", ylab="umap dim2")

# Plot the dendrogram and colors underneath
num.clust <- 2
mddCuts <- cutree(mddTree,k=num.clust)

sizeGrWindow(8,6)
dynamicMods = cutreeDynamic(dendro = mddTree, distM = d,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = 2, method = "hybrid")
mddColors = labels2colors(dynamicMods)
table(mddColors)
mddColorstable <- table(mddColors,names(mddCuts))
prop.table(mddColorstable, margin = 1)
plotDendroAndColors(mddTree, mddColors, "Dynamic Clusters",
                    dendroLabels = NULL, # hang = -1,
                    addGuide = TRUE, #guideHang = 0.05,
                    main = "Clustering with WGCNA")