216 lines
6.8 KiB
R
216 lines
6.8 KiB
R
# 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")
|