inst/doc/beadsummary.R

In beadarray: Quality assessment and low-level analysis for Illumina BeadArray data

## ----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-----------------------------------------
#  install.packages("BiocManager")
#  BiocManager::install(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--------------------------------
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-------------------------------------------------------------------

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]


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

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