### R code from vignette source 'chroGPS.Rnw' ################################################### ### code chunk number 1: import1 ################################################### options(width=70) par(mar=c(2,2,2,2)) library(chroGPS) data(s2) # Loading Dmelanogaster S2 modEncode toy example data(toydists) # Loading precomputed distGPS objects s2 ################################################### ### code chunk number 2: mds1 ################################################### # d <- distGPS(s2, metric='avgdist') d mds1 <- mds(d,k=2,type='isoMDS') mds1 mds1.3d <- mds(d,k=3,type='isoMDS') mds1.3d ################################################### ### code chunk number 3: figmds1 ################################################### cols <- as.character(s2names$Color) plot(mds1,drawlabels=TRUE,point.pch=20,point.cex=8,text.cex=.7, point.col=cols,text.col='black',labels=s2names$Factor,font=2) legend('topleft',legend=sprintf('R2=%.3f / stress=%.3f',getR2(mds1),getStress(mds1)), bty='n',cex=1) #plot(mds1.3d,drawlabels=TRUE,type.3d='s',point.pch=20,point.cex=.1,text.cex=.7, #point.col=cols,text.col='black',labels=s2names$Factor) ################################################### ### code chunk number 4: figmds1 ################################################### cols <- as.character(s2names$Color) plot(mds1,drawlabels=TRUE,point.pch=20,point.cex=8,text.cex=.7, point.col=cols,text.col='black',labels=s2names$Factor,font=2) legend('topleft',legend=sprintf('R2=%.3f / stress=%.3f',getR2(mds1),getStress(mds1)), bty='n',cex=1) #plot(mds1.3d,drawlabels=TRUE,type.3d='s',point.pch=20,point.cex=.1,text.cex=.7, #point.col=cols,text.col='black',labels=s2names$Factor) ################################################### ### code chunk number 5: figprocrustes1 ################################################### data(s2Seq) s2Seq # d2 <- distGPS(c(reduce(s2),reduce(s2Seq)),metric='avgdist') mds2 <- mds(d2,k=2,type='isoMDS') cols <- c(as.character(s2names$Color),as.character(s2SeqNames$Color)) sampleid <- c(as.character(s2names$Factor),as.character(s2SeqNames$Factor)) pchs <- rep(c(20,17),c(length(s2),length(s2Seq))) point.cex <- rep(c(8,5),c(length(s2),length(s2Seq))) par(mar=c(2,2,2,2)) plot(mds2,drawlabels=TRUE,point.pch=pchs,point.cex=point.cex,text.cex=.7, point.col=cols,text.col='black',labels=sampleid,font=2) legend('topleft',legend=sprintf('R2=%.3f / stress=%.3f',getR2(mds2),getStress(mds2)), bty='n',cex=1) legend('topright',legend=c('ChIP-Chip','ChIP-Seq'),pch=c(20,17),pt.cex=c(1.5,1)) ################################################### ### code chunk number 6: figprocrustes1 ################################################### data(s2Seq) s2Seq # d2 <- distGPS(c(reduce(s2),reduce(s2Seq)),metric='avgdist') mds2 <- mds(d2,k=2,type='isoMDS') cols <- c(as.character(s2names$Color),as.character(s2SeqNames$Color)) sampleid <- c(as.character(s2names$Factor),as.character(s2SeqNames$Factor)) pchs <- rep(c(20,17),c(length(s2),length(s2Seq))) point.cex <- rep(c(8,5),c(length(s2),length(s2Seq))) par(mar=c(2,2,2,2)) plot(mds2,drawlabels=TRUE,point.pch=pchs,point.cex=point.cex,text.cex=.7, point.col=cols,text.col='black',labels=sampleid,font=2) legend('topleft',legend=sprintf('R2=%.3f / stress=%.3f',getR2(mds2),getStress(mds2)), bty='n',cex=1) legend('topright',legend=c('ChIP-Chip','ChIP-Seq'),pch=c(20,17),pt.cex=c(1.5,1)) ################################################### ### code chunk number 7: figprocrustes2 ################################################### adjust <- rep(c('chip','seq'),c(length(s2),length(s2Seq))) sampleid <- c(as.character(s2names$Factor),as.character(s2SeqNames$Factor)) mds3 <- procrustesAdj(mds2,d2,adjust=adjust,sampleid=sampleid) par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds3,drawlabels=TRUE,point.pch=pchs,point.cex=point.cex,text.cex=.7, point.col=cols,text.col='black',labels=sampleid,font=2) legend('topleft',legend=sprintf('R2=%.3f / stress=%.3f',getR2(mds3),getStress(mds3)), bty='n',cex=1) legend('topright',legend=c('ChIP-Chip','ChIP-Seq'),pch=c(20,17),pt.cex=c(1.5,1)) ################################################### ### code chunk number 8: figpeakwidth1 ################################################### s2.pAdj <- adjustPeaks(c(reduce(s2),reduce(s2Seq)),adjust=adjust,sampleid=sampleid,logscale=TRUE) # d3 <- distGPS(s2.pAdj,metric='avgdist') mds4 <- mds(d3,k=2,type='isoMDS') par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds4,drawlabels=TRUE,point.pch=pchs,point.cex=point.cex,text.cex=.7, point.col=cols,text.col='black',labels=sampleid,font=2) legend('topleft',legend=sprintf('R2=%.3f / s=%.3f',getR2(mds4),getStress(mds4)), bty='n',cex=1) legend('topright',legend=c('ChIP-Chip','ChIP-Seq'),pch=c(20,17),pt.cex=c(1.5,1)) ################################################### ### code chunk number 9: figprocrustes2 ################################################### adjust <- rep(c('chip','seq'),c(length(s2),length(s2Seq))) sampleid <- c(as.character(s2names$Factor),as.character(s2SeqNames$Factor)) mds3 <- procrustesAdj(mds2,d2,adjust=adjust,sampleid=sampleid) par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds3,drawlabels=TRUE,point.pch=pchs,point.cex=point.cex,text.cex=.7, point.col=cols,text.col='black',labels=sampleid,font=2) legend('topleft',legend=sprintf('R2=%.3f / stress=%.3f',getR2(mds3),getStress(mds3)), bty='n',cex=1) legend('topright',legend=c('ChIP-Chip','ChIP-Seq'),pch=c(20,17),pt.cex=c(1.5,1)) ################################################### ### code chunk number 10: figpeakwidth1 ################################################### s2.pAdj <- adjustPeaks(c(reduce(s2),reduce(s2Seq)),adjust=adjust,sampleid=sampleid,logscale=TRUE) # d3 <- distGPS(s2.pAdj,metric='avgdist') mds4 <- mds(d3,k=2,type='isoMDS') par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds4,drawlabels=TRUE,point.pch=pchs,point.cex=point.cex,text.cex=.7, point.col=cols,text.col='black',labels=sampleid,font=2) legend('topleft',legend=sprintf('R2=%.3f / s=%.3f',getR2(mds4),getStress(mds4)), bty='n',cex=1) legend('topright',legend=c('ChIP-Chip','ChIP-Seq'),pch=c(20,17),pt.cex=c(1.5,1)) ################################################### ### code chunk number 11: import2 ################################################### s2.tab[1:10,1:4] ################################################### ### code chunk number 12: mds6 ################################################### d <- distGPS(s2.tab, metric='tanimoto', uniqueRows=TRUE) d mds1 <- mds(d,k=2,type='isoMDS') mds1 mds2 <- mds(d,k=3,type='isoMDS') mds2 ################################################### ### code chunk number 13: figmds6 ################################################### par(mar=c(2,2,2,2)) plot(mds1,point.cex=1.5,point.col=densCols(getPoints(mds1))) #plot(mds2,point.cex=1.5,type.3d='s',point.col=densCols(getPoints(mds2))) ################################################### ### code chunk number 14: uniqueCount ################################################### dim(s2.tab) dim(uniqueCount(s2.tab)) ################################################### ### code chunk number 15: genomeGPS1 ################################################### system.time(mds3 <- mds(d,k=2,type='isoMDS')) mds3 ################################################### ### code chunk number 16: genomeGPS2 ################################################### system.time(mds3 <- mds(d,type='isoMDS',splitMDS=TRUE,split=.5,overlap=.05,mc.cores=1)) mds3 system.time(mds4 <- mds(d,mds3,type='boostMDS',scale=TRUE)) mds4 ################################################### ### code chunk number 17: getxpr ################################################### summary(s2.wt$epigene) summary(s2.wt$gene) ################################################### ### code chunk number 18: figmds7 ################################################### plot(mds1,point.cex=1.5,scalecol=TRUE,scale=s2.wt$epigene, palette=rev(heat.colors(100))) ################################################### ### code chunk number 19: figmds7 ################################################### plot(mds1,point.cex=1.5,scalecol=TRUE,scale=s2.wt$epigene, palette=rev(heat.colors(100))) ################################################### ### code chunk number 20: figclus00 ################################################### h <- hclust(as.dist(as.matrix(d)),method='average') set.seed(149) # Random seed for the MCMC process within density estimation clus <- clusGPS(d,mds1,h,ngrid=1000,densgrid=FALSE,verbose=TRUE, preMerge=TRUE,k=max(cutree(h,h=0.5)),minpoints=20,mc.cores=1) clus ################################################### ### code chunk number 21: figclus0 ################################################### clus clusNames(clus) tabClusters(clus,125) point.col <- rainbow(length(tabClusters(clus,125))) names(point.col) <- names(tabClusters(clus,125)) point.col par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.col=point.col[as.character(clusterID(clus,125))], point.pch=19) ################################################### ### code chunk number 22: figclus1 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.95, 0.50)) plot(clus,type='contours',k=max(cutree(h,h=0.5)),lwd=5,probContour=p, drawlabels=TRUE,labcex=2,font=2) ################################################### ### code chunk number 23: figclus0 ################################################### clus clusNames(clus) tabClusters(clus,125) point.col <- rainbow(length(tabClusters(clus,125))) names(point.col) <- names(tabClusters(clus,125)) point.col par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.col=point.col[as.character(clusterID(clus,125))], point.pch=19) ################################################### ### code chunk number 24: figclus1 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.95, 0.50)) plot(clus,type='contours',k=max(cutree(h,h=0.5)),lwd=5,probContour=p, drawlabels=TRUE,labcex=2,font=2) ################################################### ### code chunk number 25: figclus2 ################################################### plot(clus,type='stats',k=max(cutree(h,h=0.5)),ylim=c(0,1),col=point.col,cex=2,pch=19, lwd=2,ylab='CCR',xlab='Cluster ID',cut=0.75,cut.lty=3,axes=FALSE) axis(1,at=1:length(tabClusters(clus,125)),labels=names(tabClusters(clus,125))); axis(2) box() ################################################### ### code chunk number 26: figclus2 ################################################### plot(clus,type='stats',k=max(cutree(h,h=0.5)),ylim=c(0,1),col=point.col,cex=2,pch=19, lwd=2,ylab='CCR',xlab='Cluster ID',cut=0.75,cut.lty=3,axes=FALSE) axis(1,at=1:length(tabClusters(clus,125)),labels=names(tabClusters(clus,125))); axis(2) box() ################################################### ### code chunk number 27: figloc1 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.5,0.95)) plot(clus,type='contours',k=max(cutree(h,h=0.5)),lwd=5,probContour=p, drawlabels=TRUE,labcex=2,font=2) fgenes <- uniqueCount(s2.tab)[,'HP1a_wa184.S2']==1 set.seed(149) c1 <- contour2dDP(getPoints(mds1)[fgenes,],ngrid=1000,contour.type='none') for (p in seq(0.1,0.9,0.1)) plotContour(c1,probContour=p,col='black') legend('topleft',lwd=1,lty=1,col='black',legend='HP1a contours (10 to 90 percent)',bty='n') ################################################### ### code chunk number 28: figloc1 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.5,0.95)) plot(clus,type='contours',k=max(cutree(h,h=0.5)),lwd=5,probContour=p, drawlabels=TRUE,labcex=2,font=2) fgenes <- uniqueCount(s2.tab)[,'HP1a_wa184.S2']==1 set.seed(149) c1 <- contour2dDP(getPoints(mds1)[fgenes,],ngrid=1000,contour.type='none') for (p in seq(0.1,0.9,0.1)) plotContour(c1,probContour=p,col='black') legend('topleft',lwd=1,lty=1,col='black',legend='HP1a contours (10 to 90 percent)',bty='n') ################################################### ### code chunk number 29: figloc2 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.5,0.95)) plot(clus,type='contours',k=max(cutree(h,h=0.5)),lwd=5,probContour=p, drawlabels=TRUE,labcex=2,font=2) set.seed(149) # Random seed for random gene sampling geneset <- sample(rownames(s2.tab),10,rep=FALSE) mds2 <- geneSetGPS(s2.tab,mds1,geneset,uniqueCount=TRUE) points(getPoints(mds2),col='black',cex=5,lwd=4,pch=20) points(getPoints(mds2),col='white',cex=4,lwd=4,pch=20) text(getPoints(mds2)[,1],getPoints(mds2)[,2],1:nrow(getPoints(mds2)),cex=1.5) legend('bottomright',col='black',legend=paste(1:nrow(getPoints(mds2)), geneset,sep=': '),cex=1,bty='n') ################################################### ### code chunk number 30: figloc2 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.5,0.95)) plot(clus,type='contours',k=max(cutree(h,h=0.5)),lwd=5,probContour=p, drawlabels=TRUE,labcex=2,font=2) set.seed(149) # Random seed for random gene sampling geneset <- sample(rownames(s2.tab),10,rep=FALSE) mds2 <- geneSetGPS(s2.tab,mds1,geneset,uniqueCount=TRUE) points(getPoints(mds2),col='black',cex=5,lwd=4,pch=20) points(getPoints(mds2),col='white',cex=4,lwd=4,pch=20) text(getPoints(mds2)[,1],getPoints(mds2)[,2],1:nrow(getPoints(mds2)),cex=1.5) legend('bottomright',col='black',legend=paste(1:nrow(getPoints(mds2)), geneset,sep=': '),cex=1,bty='n') ################################################### ### code chunk number 31: figmerge0 ################################################### set.seed(149) # Random seed for MCMC within the density estimate process clus2 <- clusGPS(d,mds1,h,ngrid=1000,densgrid=FALSE,verbose=TRUE, preMerge=TRUE,k=max(cutree(h,h=0.2)),minpoints=20,mc.cores=1) ################################################### ### code chunk number 32: figmerge1 ################################################### par(mar=c(2,2,2,2)) clus3 <- mergeClusters(clus2,brake=0,mc.cores=1) clus3 tabClusters(clus3,330) ################################################### ### code chunk number 33: figmerge1 ################################################### par(mar=c(2,2,2,2)) clus3 <- mergeClusters(clus2,brake=0,mc.cores=1) clus3 tabClusters(clus3,330) ################################################### ### code chunk number 34: figmerge21 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.95, 0.50)) plot(clus2,type='contours',k=max(cutree(h,h=0.2)), lwd=5,probContour=p,drawlabels=TRUE,labcex=2,font=2) ################################################### ### code chunk number 35: figmerge22 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.95, 0.50)) plot(clus3,type='contours',k=max(cutree(h,h=0.2)), lwd=5,probContour=p) ################################################### ### code chunk number 36: figmerge21 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.95, 0.50)) plot(clus2,type='contours',k=max(cutree(h,h=0.2)), lwd=5,probContour=p,drawlabels=TRUE,labcex=2,font=2) ################################################### ### code chunk number 37: figmerge22 ################################################### par(mar=c(0,0,0,0),xaxt='n',yaxt='n') plot(mds1,point.cex=1.5,point.col='grey') for (p in c(0.95, 0.50)) plot(clus3,type='contours',k=max(cutree(h,h=0.2)), lwd=5,probContour=p) ################################################### ### code chunk number 38: figmerge31 ################################################### plot(clus2,type='stats',k=max(cutree(h,h=0.2)),ylim=c(0,1),lwd=2, ylab='CCR',xlab='Cluster ID') ################################################### ### code chunk number 39: figmerge32 ################################################### plot(clus3,type='stats',k=max(cutree(h,h=0.2)),ylim=c(0,1),lwd=2, ylab='CCR',xlab='Cluster ID') ################################################### ### code chunk number 40: figmerge31 ################################################### plot(clus2,type='stats',k=max(cutree(h,h=0.2)),ylim=c(0,1),lwd=2, ylab='CCR',xlab='Cluster ID') ################################################### ### code chunk number 41: figmerge32 ################################################### plot(clus3,type='stats',k=max(cutree(h,h=0.2)),ylim=c(0,1),lwd=2, ylab='CCR',xlab='Cluster ID') ################################################### ### code chunk number 42: figprofile1 ################################################### p1 <- profileClusters(s2.tab, uniqueCount = TRUE, clus=clus3, i=max(cutree(h,h=0.2)), log2 = TRUE, plt = FALSE, minpoints=0) # Requires gplots library library(gplots) heatmap.2(p1[,1:20],trace='none',col=bluered(100),margins=c(10,12),symbreaks=TRUE, Rowv=FALSE,Colv=FALSE,dendrogram='none') ################################################### ### code chunk number 43: figprofile1 ################################################### p1 <- profileClusters(s2.tab, uniqueCount = TRUE, clus=clus3, i=max(cutree(h,h=0.2)), log2 = TRUE, plt = FALSE, minpoints=0) # Requires gplots library library(gplots) heatmap.2(p1[,1:20],trace='none',col=bluered(100),margins=c(10,12),symbreaks=TRUE, Rowv=FALSE,Colv=FALSE,dendrogram='none') ################################################### ### code chunk number 44: cyto1 ################################################### # For instance if mds1 contains a valid chroGPS-factors map. # gps2xgmml(mds1, fname='chroGPS_factors.xgmml', fontSize=4, # col=s2names$Color, cex=8) # And use Cytoscape -> File -> Import -> Network (Multiple File Types) # to load the generated .xgmml file