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