### R code from vignette source 'motifStack.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: DNAseqLogo ################################################### library(motifStack) pcm <- read.table(file.path(find.package("motifStack"), "extdata", "bin_SOLEXA.pcm")) pcm <- pcm[,3:ncol(pcm)] rownames(pcm) <- c("A","C","G","T") motif <- new("pcm", mat=as.matrix(pcm), name="bin_SOLEXA") ##pfm object #motif <- pcm2pfm(pcm) #motif <- new("pfm", mat=motif, name="bin_SOLEXA") opar<-par(mfrow=c(4,1)) plot(motif) #plot the logo with same height plot(motif, ic.scale=FALSE, ylab="probability") #try a different font plot(motif, font="mono,Courier") #try a different font and a different color group motif@color <- colorset(colorScheme='basepairing') plot(motif,font="Times") par(opar) ################################################### ### code chunk number 3: AAseqLogo ################################################### library(motifStack) protein<-read.table(file.path(find.package("motifStack"),"extdata","cap.txt")) protein<-t(protein[,1:20]) motif<-pcm2pfm(protein) motif<-new("pfm", mat=motif, name="CAP", color=colorset(alphabet="AA",colorScheme="chemistry")) plot(motif) ################################################### ### code chunk number 4: logostack ################################################### library(motifStack) #####Input##### pcms<-readPCM(file.path(find.package("motifStack"), "extdata"),"pcm$") motifs<-lapply(pcms,pcm2pfm) ## plot stacks motifStack(motifs, layout="stack", ncex=1.0) ################################################### ### code chunk number 5: treestack ################################################### ## plot stacks with hierarchical tree motifStack(motifs, layout="tree") ################################################### ### code chunk number 6: radialstack ################################################### ## When the number of motifs is too much to be shown in a vertical stack, ## motifStack can draw them in a radial style. ## random sample from MotifDb library("MotifDb") matrix.fly <- query(MotifDb, "Dmelanogaster") motifs2 <- as.list(matrix.fly) ## use data from FlyFactorSurvey motifs2 <- motifs2[grepl("Dmelanogaster\\-FlyFactorSurvey\\-", names(motifs2))] ## format the names names(motifs2) <- gsub("Dmelanogaster_FlyFactorSurvey_", "", gsub("_FBgn\\d+$", "", gsub("[^a-zA-Z0-9]","_", gsub("(_\\d+)+$", "", names(motifs2))))) motifs2 <- motifs2[unique(names(motifs2))] pfms <- sample(motifs2, 50) ## creat a list of object of pfm motifs2 <- lapply(names(pfms), function(.ele, pfms){new("pfm",mat=pfms[[.ele]], name=.ele)} ,pfms) ## trim the motifs motifs2 <- lapply(motifs2, trimMotif, t=0.4) ## setting colors library(RColorBrewer) color <- brewer.pal(12, "Set3") ## plot logo stack with radial style motifStack(motifs2, layout="radialPhylog", circle=0.3, cleaves = 0.2, clabel.leaves = 0.5, col.bg=rep(color, each=5), col.bg.alpha=0.3, col.leaves=rep(color, each=5), col.inner.label.circle=rep(color, each=5), inner.label.circle.width=0.05, col.outer.label.circle=rep(color, each=5), outer.label.circle.width=0.02, circle.motif=1.2, angle=350) ################################################### ### code chunk number 7: motifCloud ################################################### ## assign groups for motifs groups <- rep(paste("group",1:5,sep=""), each=10) names(groups) <- names(pfms) ## assign group colors group.col <- brewer.pal(5, "Set3") names(group.col)<-paste("group",1:5,sep="") ## use MotIV to calculate the distances of motifs jaspar.scores <- MotIV::readDBScores(file.path(find.package("MotIV"), "extdata", "jaspar2010_PCC_SWU.scores")) d <- MotIV::motifDistances(lapply(pfms, pfm2pwm)) hc <- MotIV::motifHclust(d) ## convert the hclust to phylog object phylog <- hclust2phylog(hc) ## reorder the pfms by the order of hclust leaves <- names(phylog$leaves) pfms <- pfms[leaves] ## create a list of pfm objects pfms <- lapply(names(pfms), function(.ele, pfms){ new("pfm",mat=pfms[[.ele]], name=.ele)} ,pfms) ## extract the motif signatures motifSig <- motifSignature(pfms, phylog, groupDistance=0.01, min.freq=1) ## draw the motifs with a tag-cloud style. motifCloud(motifSig, scale=c(6, .5), layout="rectangles", group.col=group.col, groups=groups, draw.legend=TRUE) ################################################### ### code chunk number 8: motifRadialPhylog ################################################### ## get the signatures from object of motifSignature sig <- signatures(motifSig) ## set the inner-circle color for each signature gpCol <- sigColor(motifSig) ## plot the logo stack with radial style. plotMotifStackWithRadialPhylog(phylog=phylog, pfms=sig, circle=0.4, cleaves = 0.3, clabel.leaves = 0.5, col.bg=rep(color, each=5), col.bg.alpha=0.3, col.leaves=rep(rev(color), each=5), col.inner.label.circle=gpCol, inner.label.circle.width=0.03, angle=350, circle.motif=1.2, motifScale="logarithmic") ################################################### ### code chunk number 9: motifCircos ################################################### ## plot the logo stack with radial style. motifCircos(phylog=phylog, pfms=pfms, pfms2=sig, col.tree.bg=rep(color, each=5), col.tree.bg.alpha=0.3, col.leaves=rep(rev(color), each=5), col.inner.label.circle=gpCol, inner.label.circle.width=0.03, col.outer.label.circle=gpCol, outer.label.circle.width=0.03, r.rings=c(0.02, 0.03, 0.04), col.rings=list(sample(colors(), 50), sample(colors(), 50), sample(colors(), 50)), angle=350, motifScale="logarithmic") ################################################### ### code chunk number 10: motifPiles ################################################### ## plot the logo stack with radial style. motifPiles(phylog=phylog, pfms=pfms, pfms2=sig, col.tree=rep(color, each=5), col.leaves=rep(rev(color), each=5), col.pfms2=gpCol, r.anno=c(0.02, 0.03, 0.04), col.anno=list(sample(colors(), 50), sample(colors(), 50), sample(colors(), 50)), motifScale="logarithmic", plotIndex=TRUE, groupDistance=0.01) ################################################### ### code chunk number 11: sessionInfo ################################################### toLatex(sessionInfo())