### R code from vignette source 'DESeq2_parathyroid.Rnw' ################################################### ### code chunk number 1: options ################################################### options(digits=3, width=80) library("cacheSweave") setCacheDir("cachedir") ################################################### ### code chunk number 2: libraries ################################################### library( "DESeq2" ) library( "parathyroidSE" ) ################################################### ### code chunk number 3: loadEcs ################################################### data("parathyroidGenesSE") ################################################### ### code chunk number 4: countTable ################################################### head( assay( parathyroidGenesSE ) ) ################################################### ### code chunk number 5: nrowSE ################################################### nrow(parathyroidGenesSE) ################################################### ### code chunk number 6: colData ################################################### colData( parathyroidGenesSE ) ################################################### ### code chunk number 7: rowData ################################################### rowData( parathyroidGenesSE ) ################################################### ### code chunk number 8: columnData ################################################### as.data.frame( colData(parathyroidGenesSE)[,c("sample","patient","treatment","time")] ) ################################################### ### code chunk number 9: split ################################################### allColSamples <- colData(parathyroidGenesSE)$sample sp <- split( seq(along=allColSamples), colData(parathyroidGenesSE)$sample ) ################################################### ### code chunk number 10: addSamples ################################################### countdata <- sapply(sp, function(columns) rowSums( assay(parathyroidGenesSE)[,columns,drop=FALSE] ) ) head(countdata) ################################################### ### code chunk number 11: poormanssum ################################################### a <- assay(parathyroidGenesSE) countdata2 <- cbind( a[,1:8], a[,9]+a[,10], a[,11], a[,12]+a[,13], a[,14:22], a[,23]+a[,24], a[,25], a[,26]+a[,27] ) all( countdata == countdata2 ) ################################################### ### code chunk number 12: subsetMetaData ################################################### coldata <- colData(parathyroidGenesSE)[sapply(sp, `[`, 1),] rownames(coldata) <- coldata$sample coldata ################################################### ### code chunk number 13: colDataSubsetCols ################################################### coldata <- coldata[ , c( "patient", "treatment", "time" ) ] head( coldata ) ################################################### ### code chunk number 14: rowdata ################################################### rowdata <- rowData(parathyroidGenesSE) rowdata ################################################### ### code chunk number 15: makeddsfull ################################################### ddsFull <- DESeqDataSetFromMatrix( countData = countdata, colData = coldata, design = ~ patient + treatment, rowData = rowdata) ddsFull ################################################### ### code chunk number 16: testCountSum ################################################### stopifnot(sum(assay(parathyroidGenesSE)) == sum(counts(ddsFull))) ################################################### ### code chunk number 17: subset ################################################### dds <- ddsFull[ , colData(ddsFull)$treatment %in% c("Control","DPN") & colData(ddsFull)$time == "48h" ] ################################################### ### code chunk number 18: refactor ################################################### dds$patient <- factor(dds$patient) dds$treatment <- factor(dds$treatment) ################################################### ### code chunk number 19: relevel ################################################### dds$treatment <- relevel( dds$treatment, "Control" ) ################################################### ### code chunk number 20: multifactorColData ################################################### colData(dds) ################################################### ### code chunk number 21: runDESeq ################################################### dds <- DESeq(dds) ################################################### ### code chunk number 22: extractResults ################################################### res <- results(dds) res ################################################### ### code chunk number 23: resCols ################################################### mcols(res) ################################################### ### code chunk number 24: checkBaseMean ################################################### all.equal(res$baseMean, rowMeans(counts(dds))) all.equal(res$baseMean, rowMeans(counts(dds,normalized=TRUE))) ################################################### ### code chunk number 25: rawpvalue ################################################### sum( res$pvalue < 0.05, na.rm=TRUE ) sum( is.na(res$pvalue) ) ################################################### ### code chunk number 26: adjpvalue ################################################### sum( res$padj < 0.1, na.rm=TRUE ) ################################################### ### code chunk number 27: strongGenes ################################################### resSig <- res[ which(res$padj < 0.1 ), ] head( resSig[ order( resSig$log2FoldChange ), ] ) ################################################### ### code chunk number 28: strongGenesUp ################################################### tail( resSig[ order( resSig$log2FoldChange ), ] ) ################################################### ### code chunk number 29: proportionSigGenes ################################################### table(sign(resSig$log2FoldChange)) ################################################### ### code chunk number 30: plotMA1 ################################################### plotMA(dds, ylim = c( -1.5, 1.5 ) ) ################################################### ### code chunk number 31: dispPlot ################################################### plotDispEsts( dds ) ################################################### ### code chunk number 32: plotMApadjchange ################################################### plotMA(dds, pvalCutoff = 0.5, ylim = c( -1.5, 1.5) ) ################################################### ### code chunk number 33: pvalHist ################################################### hist( res$pvalue, breaks=100 ) ################################################### ### code chunk number 34: filter ################################################### filterThreshold <- 2.0 keep <- rowMeans( counts( dds, normalized=TRUE ) ) > filterThreshold table( keep ) ################################################### ### code chunk number 35: adjFilter ################################################### min( res$padj[!keep], na.rm=TRUE ) ################################################### ### code chunk number 36: padjFilterCmp ################################################### table( p.adjust( res$pvalue, method="BH" ) < .1 ) table( p.adjust( res$pvalue[keep], method="BH" ) < .1 ) ################################################### ### code chunk number 37: pvalHistFilt ################################################### hist( res$pvalue[keep], breaks=100 ) ################################################### ### code chunk number 38: loadOrg ################################################### library( "org.Hs.eg.db" ) ################################################### ### code chunk number 39: keyType ################################################### cols(org.Hs.eg.db) ################################################### ### code chunk number 40: convertIDs ################################################### convertIDs <- function( ids, fromKey, toKey, db, ifMultiple=c( "putNA", "useFirst" ) ) { stopifnot( inherits( db, "AnnotationDb" ) ) ifMultiple <- match.arg( ifMultiple ) suppressWarnings( selRes <- AnnotationDbi::select( db, keys=ids, keytype=fromKey, cols=c(fromKey,toKey) ) ) if( ifMultiple == "putNA" ) { duplicatedIds <- selRes[ duplicated( selRes[,1] ), 1 ] selRes <- selRes[ ! selRes[,1] %in% duplicatedIds, ] } return( selRes[ match( ids, selRes[,1] ), 2 ] ) } ################################################### ### code chunk number 41: addSymbols ################################################### res$symbol <- convertIDs( row.names(res), "ENSEMBL", "SYMBOL", org.Hs.eg.db ) res ################################################### ### code chunk number 42: writeCSV ################################################### write.csv( as.data.frame(res), file="results.csv" ) ################################################### ### code chunk number 43: loadReactome ################################################### library( "reactome.db" ) ################################################### ### code chunk number 44: addEntrez ################################################### res$entrez <- convertIDs( row.names(res), "ENSEMBL", "ENTREZID", org.Hs.eg.db ) ################################################### ### code chunk number 45: subsetTores2 ################################################### res2 <- res[ res$entrez %in% keys( reactome.db, "ENTREZID" ) & !is.na( res$pvalue) , ] head(res2) ################################################### ### code chunk number 46: queryReactome ################################################### reactomeTable <- AnnotationDbi::select( reactome.db, keys=res2$entrez, keytype="ENTREZID", cols=c("ENTREZID","REACTOMEID") ) head(reactomeTable) ################################################### ### code chunk number 47: CreateIncidenceMatrix ################################################### incm <- do.call( rbind, with(reactomeTable, tapply( ENTREZID, factor(REACTOMEID), function(x) res2$entrez %in% x ) )) colnames(incm) <- res2$entrez str(incm) ################################################### ### code chunk number 48: PruneIncm ################################################### incm <- incm[ rowSums(incm) >= 5, ] ################################################### ### code chunk number 49: testFun ################################################### testCategory <- function( reactomeID ) { isMember <- incm[ reactomeID, ] data.frame( reactomeID = reactomeID, numGenes = sum( isMember ), avgLFC = mean( res2$log2FoldChange[isMember] ), strength = sum( res2$log2FoldChange[isMember] ) / sqrt(sum(isMember)), pvalue = t.test( res2$log2FoldChange[ isMember ] )$p.value, reactomeName = reactomePATHID2NAME[[reactomeID]] ) } ################################################### ### code chunk number 50: testTestFun ################################################### testCategory("109581") ################################################### ### code chunk number 51: runGSEA ################################################### reactomeResult <- do.call( rbind, lapply( rownames(incm), testCategory ) ) ################################################### ### code chunk number 52: padjGSEA ################################################### reactomeResult$padjust <- p.adjust( reactomeResult$pvalue, "BH" ) ################################################### ### code chunk number 53: GSEAresult ################################################### reactomeResultSignif <- reactomeResult[ reactomeResult$padjust < 0.05, ] reactomeResultSignif[ order(reactomeResultSignif$strength), ] ################################################### ### code chunk number 54: topDEGene ################################################### deGeneID <- "ENSG00000099194" res[deGeneID,] deGene <- range(rowData(dds[deGeneID,])[[1]]) names(deGene) <- deGeneID deGene ################################################### ### code chunk number 55: defineRange ################################################### as.character(seqnames(deGene)) ucscChrom <- paste0("chr",as.character(seqnames(deGene))) ucscRanges <- ranges(flank(deGene,width=10e6,both=TRUE)) subsetRange <- GRanges(ucscChrom, ucscRanges) subsetRange ################################################### ### code chunk number 56: downloadERBS (eval = FALSE) ################################################### ## ## ## ## Please do not run this code if you do not have an internet connection, ## ## alternatively use the local file import in the next code chunk. ## ## ## library( "rtracklayer" ) ## trackName <- "wgEncodeRegTfbsClusteredV2" ## tableName <- "wgEncodeRegTfbsClusteredV2" ## trFactor <- "ERalpha_a" ## mySession <- browserSession() ## ucscTable <- getTable(ucscTableQuery(mySession, track=trackName, ## range=subsetRange, table=tableName, ## name=trFactor)) ################################################### ### code chunk number 57: localERBS ################################################### ucscTableFile <- "localUcscTable.csv" ucscTable <- read.csv(ucscTableFile, stringsAsFactors=FALSE) ################################################### ### code chunk number 58: ERBS2Peaks ################################################### peaks <- with(ucscTable, GRanges(chrom, IRanges(chromStart, chromEnd), score=score)) seqlevels(peaks) <- gsub("chr(.+)","\\1",seqlevels(peaks)) seqlevels(peaks) <- seqlevels(deGene) ################################################### ### code chunk number 59: nearestPeakDo ################################################### # suppress warning which exists for BioC 2.12 of changes to distance() suppressWarnings({d2nearest <- distanceToNearest(deGene, peaks)}) ################################################### ### code chunk number 60: nearestPeakShow (eval = FALSE) ################################################### ## d2nearest <- distanceToNearest(deGene, peaks) ################################################### ### code chunk number 61: distanceShow (eval = FALSE) ################################################### ## distance(deGene, peaks) ################################################### ### code chunk number 62: distanceDo ################################################### suppressWarnings({distance(deGene, peaks)}) ################################################### ### code chunk number 63: distanceToNearest ################################################### d2nearest ################################################### ### code chunk number 64: nearestPeakOut ################################################### deGene peaks[subjectHits(d2nearest)] ################################################### ### code chunk number 65: plotPeaksAndGene ################################################### plotRange <- start(deGene) + 1e6 * c(-1,1) peakNearest <- ( seq_along(peaks) == subjectHits(d2nearest) ) plot(x=start(peaks), y=ifelse(peakNearest,.3,.2), ylim=c(0,1), xlim=plotRange, pch='p', col=ifelse(peakNearest,"red","grey60"), yaxt="n", ylab="", xlab=paste("2 Mb on chromosome",as.character(seqnames(deGene)))) points(x=start(deGene),y=.8,pch='g') ################################################### ### code chunk number 66: calculateSumForCaption ################################################### plotGRange <- GRanges(seqnames(peaks[1]),IRanges(plotRange[1],plotRange[2])) numPeaksInPlotRange <- sum(peaks %over% plotGRange) ################################################### ### code chunk number 67: peakDists ################################################### peakDists <- diff(sort(start(peaks))) summary(peakDists) mean(peakDists) ################################################### ### code chunk number 68: rld ################################################### rld <- rlogTransformation(dds) head( assay(rld) ) ################################################### ### code chunk number 69: rldPlot ################################################### par( mfrow = c( 1, 2 ) ) plot( log2( 1+counts(dds, normalized=TRUE)[, 1:2] ), col="#00000020", pch=20, cex=0.3 ) plot( assay(rld)[, 1:2], col="#00000020", pch=20, cex=0.3 ) ################################################### ### code chunk number 70: euclDist ################################################### sampleDists <- dist( t( assay(rld) ) ) sampleDists ################################################### ### code chunk number 71: sampleDistHeatmap ################################################### sampleDistMatrix <- as.matrix( sampleDists ) rownames(sampleDistMatrix) <- paste( colData(rld)$treatment, colData(rld)$patient, sep="-" ) library( "gplots" ) heatmap.2( sampleDistMatrix, trace="none" ) ################################################### ### code chunk number 72: betterheatmap ################################################### library("RColorBrewer") colours = colorRampPalette( rev(brewer.pal(9, "Blues")) )(255) heatmap.2( sampleDistMatrix, trace="none", col=colours) ################################################### ### code chunk number 73: samplePCA ################################################### print( plotPCA( rld, intgroup = c( "patient", "treatment") ) ) ################################################### ### code chunk number 74: topVarGenes ################################################### library( "genefilter" ) topVarGenes <- head( order( rowVars( assay(rld) ), decreasing=TRUE ), 35 ) ################################################### ### code chunk number 75: geneHeatmap ################################################### heatmap.2( assay(rld)[ topVarGenes, ], scale="row", trace="none", dendrogram="column", col = colorRampPalette( rev(brewer.pal(9, "RdBu")) )(255)) ################################################### ### code chunk number 76: sessInfo ################################################### sessionInfo() ################################################### ### code chunk number 77: relocateAnswers ################################################### relocateAnswers = function(con="DESeq2_parathyroid.tex"){ txt = readLines(con=con) i1 = grep("\\begin{answer}", txt, fixed=TRUE) i2 = grep("\\end{answer}", txt, fixed=TRUE) if(!((length(i1)==length(i2))&&all(i2>i1))) stop("\\begin{answer}/\\end{answer} macros are unbalanced.") i = unlist(mapply(`:`, i1, i2)) res = txt[-i] sol = txt[i] dest = grep("\\section{Solutions}", res, fixed=TRUE) if(length(dest)!=1) stop("\\section{Solutions} not found exactly once.") res = c(res[1:dest], sol, res[(dest+1):length(res)]) writeLines(res, con=con) }