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

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

#### Part A: Loading Data
## 1: Loading Gene Expression Data
load("sense.filtered.cpm.Rdata")
dim(sense.filtered.cpm) 
colnames(sense.filtered.cpm)

## 2: Demographic Data
# Loading
subject.attrs <- read.csv("Demographic_symptom.csv", stringsAsFactors = FALSE)
dim(subject.attrs)  # 160 subjects x 40 attributes
colnames(subject.attrs)  # interested in X (sample ids) and Diag (diagnosis)
subject.attrs$X
subject.attrs$Diag

# Matching gene expression samples with their diagnosis
if (!require("dplyr")) install.packages("dplyr")
library(dplyr)
# create a phenotype vector
# grab X (subject ids) and Diag (Diagnosis) from subject.attrs that 
#   intersect %in% with the RNA-Seq data
phenos.df <- subject.attrs %>% 
  filter(X %in% colnames(sense.filtered.cpm)) %>%
  dplyr::select(X, Diag)  
colnames(phenos.df) # $Diag is mdd diagnosis
# grab Diag column and convert character to factor
mddPheno <- as.factor(phenos.df$Diag)  # this is our phenotype/class vector 

summary(mddPheno) # MDD -- major depressive disorder, HC -- healthy control

#### Part B: Normalization
## 1: log2 transformation
# raw cpm boxplots and histogram of one sample
boxplot(sense.filtered.cpm,range=0,
        ylab="raw probe intensity", main="Raw", names=mddPheno)

hist(sense.filtered.cpm[,1], freq=F, ylab="density", 
     xlab="raw probe intensity", main="Raw Data Density for Sample 1")

# log2 transformed boxplots and histogram
boxplot(log2(sense.filtered.cpm), range=0,
        ylab="log2 intensity", main="Log2 Transformed", names=mddPheno)

hist(log2(sense.filtered.cpm[,1]), freq=F, ylab="density", 
     xlab="log2 probe intensity", main="log2 Data Density for Sample 1")

getmode <- function(v) {
  uniqv <-unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

data <- data.frame(Mean = c(mean(sense.filtered.cpm[,1]),
                            mean(log2(sense.filtered.cpm[,1]))),
                   Mode = c(getmode(sense.filtered.cpm[,1]),
                            getmode(log2(sense.filtered.cpm[,1]))),
                   Median = c(median(sense.filtered.cpm[,1]))
                   )
rownames(data) = c("Original", "Log2 Transformed")
data

## 2: Quantile Normalization
# install quantile normalize
#install.packages("BiocManager")
if (!require("BiocManager")) install.packages("BiocManager")
library(BiocManager)
if (!require("preprocessCore"))
  BiocManager::install("preprocessCore")
library(preprocessCore)    # replace with custom function?

# apply quantile normalization
mddExprData_quantile <- normalize.quantiles(sense.filtered.cpm)

boxplot(mddExprData_quantile,range=0,
        ylab="raw intensity", main="Quantile Normalized")

head(mddExprData_quantile)

## 3: Log2 on quantile normalized data
mddExprData_quantileLog2 <- log2(mddExprData_quantile)
# add phenotype names to matrix
colnames(mddExprData_quantileLog2) <- mddPheno  

boxplot(mddExprData_quantileLog2,range=0,
        ylab="log2 intensity", main="Quantile Normalized Log2")

hist(log2(mddExprData_quantileLog2[,1]), freq=F, 
     ylab="density", xlab="log2 probe intensity", 
     main="log2 Quantile Normalized for Sample 1")

## 4: Means
mean(mddExprData_quantileLog2[,1])
colMeans(mddExprData_quantileLog2)