vignettes/GSE30849/alt_workflow.R

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")
axelmuller/geometric2 documentation built on May 25, 2019, 6:26 p.m.