71 lines
2.5 KiB
R
71 lines
2.5 KiB
R
# Lab 5 for the University of Tulsa's CS-6643 Bioinformatics Course
|
|
# Gene Expression Statistical Learning
|
|
# Professor: Dr. McKinney, Fall 2022
|
|
# Noah L. Schrick - 1492657
|
|
|
|
## Set Working Directory to file directory - RStudio approach
|
|
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
|
|
|
|
#### 0: Process and filter data
|
|
# load gene expression data
|
|
load("sense.filtered.cpm.Rdata") # setwd!
|
|
|
|
# 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)
|
|
|
|
# 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)
|
|
|
|
# convert phenotype to factor
|
|
pheno.factor <- as.factor(colnames(GxS.covfilter))
|
|
pheno.factor
|
|
str(pheno.factor)
|
|
levels(pheno.factor)
|
|
|
|
|
|
#### Part A: Logistic Regression
|
|
# make sure HC is the reference level
|
|
pheno.factor.relevel <- relevel(pheno.factor,"HC")
|
|
levels(pheno.factor.relevel)
|
|
# also rename levels "0"/"1" from 1/2
|
|
levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="MDD"] <- 1
|
|
levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="HC"] <- 0
|
|
|
|
# Fit logistic model of first gene to phenotype data
|
|
gene.row <- 2
|
|
gene.name <-rownames(GxS.covfilter)[gene.row]
|
|
gene.expr <- GxS.covfilter[gene.row,]
|
|
gene.fit <- glm(pheno.factor.relevel~gene.expr, family=binomial)
|
|
summary(gene.fit)
|
|
coeff.mat <- coef(summary(gene.fit))
|
|
b0 <- coeff.mat[1,1]
|
|
b1 <- coeff.mat[2,1]
|
|
b1.pval <- coeff.mat[2,4] |