library(data.table)
library(limma)
library(GEOquery)
library(stringr)
#library(affycoretools)
#(geometric2)
# load biosamps data from package
#data("biosamps")
#setwd("/home/axel/projects/nursa/submission/")
gse <- "GSE11202"
#dir.create(paste0(gse, "_alt"))
my_path <- paste0("/home/axel/projects/nursa/submission/", gse)
setwd(my_path)
eset<-getGEO(gse, destdir = "GEOtemp")[[1]]
#set up targets
targets<-pData(eset)
colnames(targets)
#
#
#targets<-targets[,c(1,2,10)]
targets
sample_review3 <- fread("sample_review3.txt")
names(sample_review3)
sample_review3$title
# The sample titles are of the form "'D' ODN for 9 hours - repeat 3 - mAdbID:103438"
# remove everything from - repeat onwards
treatment <- str_trim(str_remove(sample_review3$title, "- *repeat.*$"))
treatment
targets$rep <- sample_review3$Group
targets$treat <- gsub(" ",".", treatment)
#targets$treat <- gsub(" ",".",sample_review3$title)
targets
#####################################
# design matrix according to Matt, Ana & userguide
exprs <- exprs(eset)
SlideNumber <- rownames(targets)
Cy3 <- make.names(targets$source_name_ch1)
# the previous command failed for GSE30849 due to mistakes in the GEO entry,
# hence, Cy3 has to be created manually.
Cy3 <- rep(make.names("Universal Mouse Reference RNA"), length(Cy3))
#Cy5 <- make.names(targets$source_name_ch2)
Cy5 <- make.names(targets$rep)
#Cy3 <- gsub(" ",".", targets$source_name_ch1)
#Cy5 <- gsub(" ", ".", targets$source_name_ch2)
targets_edited <- data.frame(SlideNumber, Cy3, Cy5)
targets_matrix <- as.matrix(targets_edited)
design <- modelMatrix(targets_matrix, ref="Universal.Mouse.Reference.RNA")
fit1 <- lmFit(exprs, design)
# in this entry different contiditions of each group are stored in source_name_ch2
#contrast_values <- gsub(" ", ".", unique(sample_review3$source_name_ch2))
# alternatively the group names can be used, that might facilitate integration of contrasts
#colnames(design) <- unique(sample_review3$Group)
#colnames(design) <- unique(sample_review3$Group)
contrasts <- fread("contrasts.txt")
contrasts$Contrast
contrasts$Formula
# the following works try to use the contrast file to construct this matrix:
contrast.matrix <- makeContrasts("ODN1555;ODN1466 vs | 0.5h BALB/c" = "a-A",
"ODN1555;ODN1466 vs | 1h BALB/c" = "a-B",
"ODN1555;ODN1466 vs | 24h BALB/c" = "a-C",
"ODN1555;ODN1466 vs | 3h BALB/c" = "a-D",
"ODN1612;ODN1471 vs | 3h BALB/c" = "a-E",
"ODN1555;ODN1466 vs | 72h BALB/c" = "a-F",
"ODN1555;ODN1466 vs | 9h BALB/c" = "a-G",
"ODN1555;ODN1466 vs | 3h" = "h-I",
"ODN1555;ODN1466 vs | 3h TNF-" = "j-K",
levels=c("a", "A", "B", "C", "D", "E", "F", "G", "h", "I", "j", "K"))
contrast.matrix <- makeContrasts("ODN1555;ODN1466 vs | 0.5h" = "a-A",
"ODN1555;ODN1466 vs | 1h" = "a-B",
"ODN1555;ODN1466 vs | 24h" = "a-C",
"ODN1555;ODN1466 vs | 3h" = "a-D",
"ODN1555;ODN1466 vs | 3h" = "a-E",
"ODN1555;ODN1466 vs | 72h" = "a-F",
"ODN1555;ODN1466 vs | 9h" = "a-G",
"ODN1555;ODN1466 vs | 3h" = "h-I",
"ODN1555;ODN1466 vs | 3h" = "j-K",
levels=c("a", "A", "B", "C", "D", "E", "F", "G", "h", "I", "j", "K"))
*** Step 5
contrasts$Contrast
contrasts$Formula
# contrast.matrix.test <- makeContrasts(contrasts$Contrast[1] = contrasts$Formula[1],
# contrasts$Contrast[2] = contrasts$Formula[2],
# levels=c("a", "A", "B"))
contrast.matrix.test <- makeContrasts(contrasts = contrasts$Formula, levels = fit1$design)
# contrast.matrix <- makeContrasts("ODN1555;ODN1466 vs | 14d" = "Untreated BALB/c WT" - "CpG ODN - 14 days",
# "ODN1555;ODN1466 vs | 24h" = "Untreated BALB/c WT" - "CpG ODN - 24 hrs BALB/c WT",
# "ODN1555;ODN1466 vs | 3h" = "Untreated BALB/c WT" - "CpG ODN - 3 hrs BALB/c WT",
# "ODN1612;ODN1471 vs | 3h" = "Untreated BALB/c WT" - "Control ODN - 3 hrs BALB/c WT",
# "ODN1555;ODN1466 vs | 5d" = "Untreated BALB/c WT" - "CpG ODN - 5 days",
# "ODN1555;ODN1466 vs | 72h" = "Untreated BALB/c WT" - "CpG ODN - 72 hrs BALB/c WT",
# "ODN1555;ODN1466 vs | 7d" = "Untreated BALB/c WT" - "CpG ODN - 7 days",
# "ODN1555;ODN1466 vs | 9d" = "Untreated BALB/c WT" - "CpG ODN - 9 days",
# "ODN1555;ODN1466 vs | 9h" = "Untreated BALB/c WT" - "CpG ODN - 9 hrs BALB/c WT",
# levels = design)
#### this yields an error: The levels must by syntactically valid names in R, see help(make.names).
# # use make.names function for Cy3 and Cy5 (see above)
# contrast.matrix <- makeContrasts("ODN1555;ODN1466 vs | 14d" = "Untreated.BALB.c.WT" - "CpG.ODN...14.days",
# "ODN1555;ODN1466 vs | 24h" = "Untreated.BALB.c.WT" - "CpG.ODN...24.hrs.BALB.c.WT",
# "ODN1555;ODN1466 vs | 3h" = "Untreated.BALB.c.WT" - "CpG.ODN...3.hrs.BALB.c.WT",
# "ODN1612;ODN1471 vs | 3h" = "Untreated.BALB.c.WT" - "Control.ODN...3.hrs.BALB.c.WT",
# "ODN1555;ODN1466 vs | 5d" = "Untreated.BALB.c.WT" - "CpG.ODN...5.days",
# "ODN1555;ODN1466 vs | 72h" = "Untreated.BALB.c.WT" - "CpG.ODN...72.hrs.BALB.c.WT",
# "ODN1555;ODN1466 vs | 7d" = "Untreated.BALB.c.WT" - "CpG.ODN...7.days",
# "ODN1555;ODN1466 vs | 9d" = "Untreated.BALB.c.WT" - "CpG.ODN...9.days",
# "ODN1555;ODN1466 vs | 9h" = "Untreated.BALB.c.WT" - "CpG.ODN...9.hrs.BALB.c.WT",
# levels = design)
# # use make.names function for Cy3 and Cy5 (see above)
# contrast.matrix <- makeContrasts("Untreated.BALB.c.WT" - "CpG.ODN...14.days",
# "Untreated.BALB.c.WT" - "CpG.ODN...24.hrs.BALB.c.WT",
# "ODN1555;ODN1466 vs | 3h" = "Untreated.BALB.c.WT" - "CpG.ODN...3.hrs.BALB.c.WT",
# "ODN1612;ODN1471 vs | 3h" = "Untreated.BALB.c.WT" - "Control.ODN...3.hrs.BALB.c.WT",
# "ODN1555;ODN1466 vs | 5d" = "Untreated.BALB.c.WT" - "CpG.ODN...5.days",
# "ODN1555;ODN1466 vs | 72h" = "Untreated.BALB.c.WT" - "CpG.ODN...72.hrs.BALB.c.WT",
# "ODN1555;ODN1466 vs | 7d" = "Untreated.BALB.c.WT" - "CpG.ODN...7.days",
# "ODN1555;ODN1466 vs | 9d" = "Untreated.BALB.c.WT" - "CpG.ODN...9.days",
# "ODN1555;ODN1466 vs | 9h" = "Untreated.BALB.c.WT" - "CpG.ODN...9.hrs.BALB.c.WT",
# levels = design)
#contrast_assignments <- contrasts$Contrast[1]
# contrast.matrix <- makeContrasts("ODN1555;ODN1466 vs | 14d"=a-A,
# "ODN1555;ODN1466 vs | 24h"=a-B,
# levels=c("a", "A", "B"))
#x <- c("B-A","C-B","C-A", "a-A", "a-B")
#makeContrasts(contrasts=x,levels=c("A","B","C", "a"))
fit2 <- contrasts.fit(fit1, contrast.matrix)
fit3 <- eBayes(fit2)
write.fit(fit3, results=NULL, file="fit3.txt", digits=10, adjust="fdr", method="separate",F.adjust="none", sep="\t")
#####################################
library(geometric2)
#step5 gene annotation
fitfiles <- get_fitfiles("Results/")
## get reannotation package
path2GEOmetadb <- "/home/axel/projects/nursa/GEOmetric/GEOmetadb.sqlite"
DSinfo <- get_DSinfo(gse, path2GEOmetadb)
#source("https://bioconductor.org/biocLite.R")
#biocLite(DSinfo$Reannotation)
#library(illuminaMousev2.db)
#process_fitoutput(gse, fitfile = fitfiles[1], illuminaMousev2.db)
# if process_ fitoutput fails try:
process_fitoutput_with_gpl(gse, fitfile = fitfiles, "GPL4200")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.