[untitled]

Script

library(affy)

getwd()
#setwd(dirPath)
setwd("/Users/dennisleenheer/Dropbox/DennisLeenheerWorkflows/Data")

pd <- read.AnnotatedDataFrame("PhenotypeData.tsv", header=TRUE, sep="\t", row.names=1)

rawAffyData = ReadAffy(filenames = rownames(pData(pd)), phenoData = pd, verbose = TRUE)

eset <- expresso(rawAffyData, bgcorrect.method="rma", normalize.method="quantiles", pmcorrect.method="pmonly", summary.method="medianpolish")

write.table(exprs(eset),sep="\t", quote=FALSE, col.names = NA, file="WT_vs_PPARa0_6h_Dennis.txt")

R Server

Script

source("http://bioconductor.org/biocLite.R")
biocLite("mouse430a2.db")
library(mouse430a2.db)

#Read the data and create necessary variables
response <- as.factor(c("WT","WT","WT","WT","KO","KO","KO","KO"))
data <- read.table("WT_vs_PPARa0_6h_Dennis.txt", header=T, sep="\t")
phenodata <- as.data.frame(response, ncol=1) 
colnames(phenodata)[1] <- "response"
rownames(phenodata) <- as.character(colnames(data[,2:9]))
genes <- as.data.frame(data[,1], ncol=1)
colnames(genes)[1] <- "probe"
rownames(genes) <- data[,1]
rownames(data) <- data[,1]
phenodata2 <- new("AnnotatedDataFrame", data=phenodata)
genes2 <-  new("AnnotatedDataFrame", data=genes)
exprsdata <- new("ExpressionSet", phenoData=phenodata2, exprs=data[,2:9], featureData=genes2)
eg <- as.list(mouse430a2ENTREZID)

#summarize data at EG level
probe2entrez <- unlist(eg)
matchprobe <- match(featureNames(exprsdata), names(probe2entrez))
egID <- probe2entrez[matchprobe]
list.of.unique.eg <- unique(egID)
list.of.unique.eg <- list.of.unique.eg[!is.na(list.of.unique.eg)]

#average expression values for probes with same EG annotation (warning: this takes a some time!)
data_eg <- array(NA, dim=c(length(list.of.unique.eg),8))
for (i in 1:length(list.of.unique.eg))
{
 temp <- data[which(list.of.unique.eg[i]==egID),2:9]
 if(is.data.frame(temp)) {data_eg[i,] <- apply(temp,2,mean)} else data_eg[i,] <- temp
}

entrezgenes <- as.data.frame(list.of.unique.eg, ncol=1)
colnames(entrezgenes)[1] <- "eg"
rownames(entrezgenes) <- entrezgenes[,1]

rownames(data_eg) <- entrezgenes[,1]
colnames(data_eg) <- colnames(data[,2:9])

write.table(data_eg, "WT_vs_PPARa0_6h_eg_summarized_Dennis.txt", sep="\t",  quote=FALSE, col.names = NA)

R Server

Script

#Load libraries
source("http://bioconductor.org/biocLite.R")
biocLite("globaltest")
library(globaltest)
biocLite("GO.db")
library(GO.db)

#Read the data and create necessary variables
data <- read.table("WT_vs_PPARa0_6h_eg_summarized_Dennis.txt", header=T, sep="\t")
rownames(data) <- data[,1]
response <- as.factor(c("WT","WT","WT","WT","KO","KO","KO","KO"))
phenodata <- as.data.frame(response, ncol=1) 
colnames(phenodata)[1] <- "response"
rownames(phenodata) <- as.character(colnames(data[,2:9]))
genes <- as.data.frame(data[,1], ncol=1)
colnames(genes)[1] <- "eg"
rownames(genes) <- data[,1]
phenodata2 <- new("AnnotatedDataFrame", data=phenodata)
genes2 <-  new("AnnotatedDataFrame", data=genes)
exprsdata <- new("ExpressionSet", phenoData=phenodata2, exprs=data[,2:9], featureData=genes2)

#Perform the standard global test to see if there is enough difference between the wild type and the mutant data
gt.all <- gt( phenoData(exprsdata)$response, exprsdata)
summary(gt.all)

#Perform the global test for Gene Ontology Biological Processes
gt_BP <-  gtGO(response, exprsdata, annotation = "org.Mm.eg", ontology="BP")

#Sort the results and write to output file
result(sort(gt_BP))[1:10,]
write.table(result(gt_BP), "result_BP_standard_GO_Dennis.txt", sep="\t", quote=FALSE, col.names = NA)

R Server

Script

library(globaltest)
library(GO.db)

#Read the data and create necessary variables
data <- read.table("WT_vs_PPARa0_6h_eg_summarized_Dennis.txt", header=T, sep="\t")
rownames(data) <- data[,1]
response <- as.factor(c("WT","WT","WT","WT","KO","KO","KO","KO"))
phenodata <- as.data.frame(response, ncol=1) 
colnames(phenodata)[1] <- "response"
rownames(phenodata) <- as.character(colnames(data[,2:9]))
genes <- as.data.frame(data[,1], ncol=1)
colnames(genes)[1] <- "eg"
rownames(genes) <- data[,1]
phenodata2 <- new("AnnotatedDataFrame", data=phenodata)
genes2 <-  new("AnnotatedDataFrame", data=genes)
exprsdata <- new("ExpressionSet", phenoData=phenodata2, exprs=data[,2:9], featureData=genes2)

 
#Perform the standard global test to see if there is enough difference between the wild type and the mutant data
gt.all <- gt( phenoData(exprsdata)$response, exprsdata)
summary(gt.all)
#Perform the literature weighted global test for Gene Ontology Biological Processes
p2e <- as.list(featureNames(exprsdata))
names(p2e) <- featureNames(exprsdata)
gt_BP <-  gtConcept(response, exprsdata, probe2entrez = p2e, annotation = "org.Mm.eg", conceptmatrix="Fully Formed Anatomical Structure.txt", concept2name="conceptID2name.txt",entrez2concept="entrezGeneToConceptID.txt")

R Server