## ----eval=FALSE---------------------------------------------------------------
#  install.packages("BiocManager")
#  BiocManager::install(c("beadarrayExampleData", "illuminaHumanv3.db"))

## -----------------------------------------------------------------------------
library("beadarray")

require(beadarrayExampleData)

data(exampleSummaryData)

exampleSummaryData


## -----------------------------------------------------------------------------
exprs(exampleSummaryData)[1:5,1:5]
se.exprs(exampleSummaryData)[1:5,1:5]


## -----------------------------------------------------------------------------
head(fData(exampleSummaryData))
table(fData(exampleSummaryData)[,"Status"])

pData(exampleSummaryData)

## -----------------------------------------------------------------------------
channelNames(exampleSummaryData)

exampleSummaryData.log2 <- channel(exampleSummaryData, "G")
exampleSummaryData.unlogged <- channel(exampleSummaryData, "G.ul")


sampleNames(exampleSummaryData.log2)
sampleNames(exampleSummaryData.unlogged)

exprs(exampleSummaryData.log2)[1:10,1:3]
exprs(exampleSummaryData.unlogged)[1:10,1:3]

## -----------------------------------------------------------------------------
exampleSummaryData.log2[,1:4]
exampleSummaryData.log2[1:10,]


## -----------------------------------------------------------------------------
randIDs <- sample(featureNames(exampleSummaryData), 1000)

exampleSummaryData[randIDs,]

## -----------------------------------------------------------------------------
boxplot(exampleSummaryData.log2[randIDs,])

## -----------------------------------------------------------------------------
boxplot(exampleSummaryData.log2[randIDs,], what="nObservations")

## -----------------------------------------------------------------------------
boxplot(exampleSummaryData.log2[randIDs,], SampleGroup="SampleFac")

## -----------------------------------------------------------------------------
boxplot(exampleSummaryData.log2[randIDs,], probeFactor = "Status")

## -----------------------------------------------------------------------------
annotation(exampleSummaryData)

exampleSummaryData.log2 <- addFeatureData(exampleSummaryData.log2, 
toAdd = c("SYMBOL", "PROBEQUALITY", "CODINGZONE", "PROBESEQUENCE", "GENOMICLOCATION"))

head(fData(exampleSummaryData.log2))

illuminaHumanv3()

## -----------------------------------------------------------------------------
ids <- which(fData(exampleSummaryData.log2)[,"SYMBOL"] == "ALB")

boxplot(exampleSummaryData.log2[ids,], 
  SampleGroup = "SampleFac", probeFactor = "IlluminaID")

## ----fig.width=12-------------------------------------------------------------
library(ggplot2)
library(gridExtra)
bp1 <- boxplot(exampleSummaryData.log2[ids,], 
SampleGroup = "SampleFac", probeFactor = "IlluminaID") 
bp1 <- bp1+ labs(title = "ALB expression level comparison") + xlab("Illumina Probe") + ylab("Log2 Intensity")

bp2 <- boxplot(exampleSummaryData.log2[randIDs,], probeFactor = "Status") 
bp2 <- bp2 + labs(title = "Control Probe Comparison")

grid.arrange(bp1,bp2)

## -----------------------------------------------------------------------------
bp1$data

## -----------------------------------------------------------------------------
mas <- plotMA(exampleSummaryData.log2,do.log=FALSE)

mas

## -----------------------------------------------------------------------------
##Added lines on the y axis
mas + geom_hline(yintercept=c(-1.5,1.5),col="red",lty=2)
##Added a smoothed line to each plot
mas+ geom_smooth(col="red")
##Changing the color scale
mas + scale_fill_gradient2(low="yellow",mid="orange",high="red")

## -----------------------------------------------------------------------------
mas <- plotMA(exampleSummaryData.log2,do.log=FALSE,SampleGroup="SampleFac")
mas[[1]]


## -----------------------------------------------------------------------------
exampleSummaryData.norm <- normaliseIllumina(exampleSummaryData.log2, 
method="quantile", transform="none")

## -----------------------------------------------------------------------------
exampleSummaryData.norm2 <- normaliseIllumina(channel(exampleSummaryData, "G.ul"), 
  method="neqc", transform="none")

## -----------------------------------------------------------------------------
library(illuminaHumanv3.db)

ids <- as.character(featureNames(exampleSummaryData.norm))

qual <- unlist(mget(ids, illuminaHumanv3PROBEQUALITY, ifnotfound=NA))

table(qual)

rem <- qual == "No match" | qual == "Bad" | is.na(qual)

exampleSummaryData.filt <- exampleSummaryData.norm[!rem,]

dim(exampleSummaryData.filt)

## ----eval=FALSE---------------------------------------------------------------
#  rna <- factor(pData(exampleSummaryData)[,"SampleFac"])
#  
#  design <- model.matrix(~0+rna)
#  colnames(design) <- levels(rna)
#  aw <- arrayWeights(exprs(exampleSummaryData.filt), design)
#  aw
#  fit <- lmFit(exprs(exampleSummaryData.filt), design, weights=aw)
#  contrasts <- makeContrasts(UHRR-Brain, levels=design)
#  contr.fit <- eBayes(contrasts.fit(fit, contrasts))
#  topTable(contr.fit, coef=1)

## -----------------------------------------------------------------------------
limmaRes <- limmaDE(exampleSummaryData.filt, SampleGroup="SampleFac")
limmaRes

DesignMatrix(limmaRes)
ContrastMatrix(limmaRes)
ArrayWeights(limmaRes)
plot(limmaRes)

## -----------------------------------------------------------------------------
gr <- as(exampleSummaryData.filt[,1:5], "GRanges")
gr

## -----------------------------------------------------------------------------
lgr <- as(limmaRes, "GRanges")
lgr

## -----------------------------------------------------------------------------
lgr <- lgr[[1]]
lgr[order(lgr$LogOdds,decreasing=T)]

lgr[p.adjust(lgr$PValue)<0.05]


## -----------------------------------------------------------------------------
library(GenomicRanges)
  HBE1 <- GRanges("chr11", IRanges(5289580,5291373),strand="-")

  lgr[lgr %over% HBE1]

## ----eval=FALSE---------------------------------------------------------------
#  library(ggbio)
#  library(TxDb.Hsapiens.UCSC.hg19.knownGene)
#  tx <- TxDb.Hsapiens.UCSC.hg19.knownGene
#  p1 <- autoplot(tx, which=HBE1)
#  p2 <-   autoplot(lgr[lgr %over% HBE1])
#  tracks(p1,p2)
#  id <- plotIdeogram(genome="hg19", subchr="chr11")
#  tracks(id,p1,p2)

## ----eval=FALSE---------------------------------------------------------------
#  plotGrandLinear(lgr, aes(y = LogFC))

## ----eval=FALSE---------------------------------------------------------------
#  rawdata <- channel(exampleSummaryData, "G")
#  normdata <- normaliseIllumina(rawdata)
#  
#  makeGEOSubmissionFiles(normdata,rawdata)
#  

## ----eval=FALSE---------------------------------------------------------------
#  download.file(
#  "ftp://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL6nnn/GPL6947/annot/GPL6947.annot.gz",
#  destfile="GPL6947.annot.gz"
#  )
#  
#  makeGEOSubmissionFiles(normdata,rawdata,softTemplate="GPL6947.annot.gz")
#  

## -----------------------------------------------------------------------------
library(GEOquery)

url <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE33nnn/GSE33126/matrix/"
filenm <- "GSE33126_series_matrix.txt.gz"
if(!file.exists("GSE33126_series_matrix.txt.gz")) download.file(paste(url, filenm, sep=""), destfile=filenm)
gse <- getGEO(filename=filenm)
head(exprs(gse))

## -----------------------------------------------------------------------------
summaryData <- as(gse, "ExpressionSetIllumina")
summaryData
head(fData(summaryData))

## -----------------------------------------------------------------------------
fData(summaryData)$Status <- 
  ifelse(fData(summaryData)$PROBEQUALITY=="No match","negative","regular" )

Detection(summaryData) <- calculateDetection(summaryData, 
                              status=fData(summaryData)$Status)

## -----------------------------------------------------------------------------
summaryData.norm <- normaliseIllumina(summaryData,method="neqc", 
    status=fData(summaryData)$Status)
boxplot(summaryData.norm)

## -----------------------------------------------------------------------------
limmaResults <- limmaDE(summaryData.norm, "source_name_ch1")
limmaResults

## ----eval=FALSE---------------------------------------------------------------
#  library(beadarray)
#  dataFile = "AsuragenMAQC-probe-raw.txt"
#  qcFile = "AsuragenMAQC-controls.txt"
#  BSData = readBeadSummaryData(dataFile = dataFile,
#  qcFile = qcFile, controlID = "ProbeID",
#  skip = 0, qc.skip = 0, qc.columns = list(exprs = "AVG_Signal",
#  Detection = "Detection Pval"))
#  

## ----eval=FALSE---------------------------------------------------------------
#  library(beadarray)
#  library(GEOquery)
#  downloadDir <- tempdir()
#  getGEOSuppFiles("GSE27073", makeDirectory = FALSE, baseDir = downloadDir)
#  idatFiles <- list.files(path = downloadDir, pattern = ".idat.gz", full.names=TRUE)
#  sapply(idatFiles, gunzip)
#  idatFiles <- list.files(path = downloadDir, pattern = ".idat", full.names=TRUE)
#  BSData <- readIdatFiles(idatFiles)