## ----knitr, echo=FALSE, results="hide"-----------------------------------
library("knitr")
opts_chunk$set(tidy=FALSE,dev="png",fig.show="as.is",
               fig.width=10,fig.height=6,
               message=FALSE,eval=T,warning=FALSE)

## ----style, eval=TRUE, echo=F, results="asis"----------------------------
BiocStyle::latex()

## ----include=FALSE-------------------------------------------------------
library(knitr)
opts_chunk$set(
concordance=TRUE
)

## ----installExampleData,eval=FALSE---------------------------------------
## source("http://www.bioconductor.org/biocLite.R")
## biocLite(c("beadarrayExampleData", "illuminaHumanv3.db"))
## 

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

require(beadarrayExampleData)

data(exampleSummaryData)

exampleSummaryData



## ----objectPeek----------------------------------------------------------

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



## ----annotations---------------------------------------------------------

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

pData(exampleSummaryData)


## ----subsetting1---------------------------------------------------------

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]


## ----subsetting2---------------------------------------------------------


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


## ----subsetting4---------------------------------------------------------

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

exampleSummaryData[randIDs,]


## ----boxplot1------------------------------------------------------------

boxplot(exampleSummaryData.log2[randIDs,])


## ----boxplot2------------------------------------------------------------

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


## ----boxplot4------------------------------------------------------------

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


## ----boxplot5------------------------------------------------------------

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


## ----addFdata------------------------------------------------------------

annotation(exampleSummaryData)

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

head(fData(exampleSummaryData.log2))

illuminaHumanv3()


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

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

## ----ggplot-layout,fig.height=8,fig.width=8------------------------------
require("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-----------------------------------------------------------------

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]]

## ----normalise1----------------------------------------------------------

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


## ----normalise2----------------------------------------------------------

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


## ----filter--------------------------------------------------------------

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)


## ----deanalysis,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

## ----toGRangeslgr--------------------------------------------------------
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]


## ------------------------------------------------------------------------
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)

## ------------------------------------------------------------------------
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")
## 

## ----eval=FALSE----------------------------------------------------------
## library(GEOquery)
## url <- "ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE33126/"
## 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))
## 

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

## ----eval=FALSE----------------------------------------------------------
## fData(summaryData)$Status <-
##   ifelse(fData(summaryData)$PROBEQUALITY=="No match","negative","regular" )
## 
## Detection(summaryData) <- calculateDetection(summaryData,
##                               status=fData(summaryData)$Status)
## 

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

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

## ----readBeadSummary, 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"))
## 

## ----readIDAT, 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)

## ----options, echo=FALSE, eval=TRUE--------------------------------------
options(width = 80)

## ----sessionInfo---------------------------------------------------------
sessionInfo()