R/pipeline_Sun.R

In Syksy/curatedTools: Tools for generating open multi-omics datasets with curated clinical data

####
#
# Teemu Daniel Laajala
# Processing data by Sun et al.
# GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE25136
# Available data types; homogeneous dataset of disease recurrence in GEX, with total 79 samples
# Raw tarball download address: https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE25136&format=file
#
####


# Following GEOquery tutorial from https://bioconductor.org/help/course-materials/2011/BioC2011/LabStuff/publicDataTutorial.pdf
# Downloading the supplementary raw files instead of manual download/moving/etc
library(GEOquery)
# Downloading raw data for Sun et al.
supfiles <- getGEOSuppFiles('GSE25136')
# Download size: 269.5 MB
# Open the tarball
utils::untar(tarfile=rownames(supfiles))
supfiles2 <- list.files()
supfiles2 <- supfiles2[grep(".gz", supfiles2)]
invisible(lapply(supfiles2, FUN=gunzip))
# CEL files are now available in the working directory
# using 'affy' package
library("affy")
Sun <- ReadAffy()
colnames(Sun) <- gsub(".CEL", "", colnames(Sun))
# Requires annotations for: installing the source package 'hgu133acdf'
# > boxplot(Sun)
# -> Requires normalization
# Need to be careful not to mask 'rma from 'affy' by the 'rma' from 'oligo'
Sun2 <- affy::rma(Sun)
# > boxplot(Sun2)
# -> Values look much better
# Removing .CEL and packaging names from the GEO-compatible sample names
colnames(Sun2) <- gsub(".gz", "", colnames(Sun2))

GEX_Sun <- Sun2
#> class(GEX_Sun)
#[1] "ExpressionSet"
#attr(,"package")
#[1] "Biobase"

# Fetch the corresponding database for probe annotations from Bioconductor
BiocManager::install("hgu133a.db", version = "3.8")
library(hgu133a.db)
#head(hgu133aALIAS2PROBE)
#head(hgu133aGENENAME)
keys <- mappedkeys(hgu133aGENENAME)

#> head(as.character(hgu133aALIAS2PROBE))
#      RFC40        RFC2     HSP70B'       HSPA6        PAX8    C6orf131 
#  "1053_at"   "1053_at"    "117_at"    "117_at"    "121_at" "1255_g_at" 
#Warning message:
#In .local(x, ...) : returned vector has duplicated names

nam <- names(as.character(hgu133aALIAS2PROBE)[match(rownames(GEX_Sun), as.character(hgu133aALIAS2PROBE))])
nam[is.na(nam)] <- "NA"

rownames(GEX_Sun) <- make.unique(nam)

save(GEX_Sun, file="GEX_Sun.RData")