# 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