Hi everybody,
based on Mark's hints I wrote some code to plot FirmaScores with
GenomeGraphs. Mark's way of extracting the data (d <-log2(extractMatrix
(cs,cells=ind,verbose=verbose)) ) is very elegant and effective. But
somehow I didn't manage to pull the firma scores data off in the same
way.
That's why I wrote (a pretty ugly) routine for this:
getPlotDataGenomeGraphs<- function
(id,idname=NULL,CelSet,annotationfile,groups,data="intensities"){
cdf <- getCdf(CelSet)
u<-which(getUnitNames(cdf)==id)
ugcM <- getUnitGroupCellMap(cdf, units=u, retNames=TRUE)
nProbes <- table(ugcM$group)
m <- match(names(nProbes), annotationfile$probeset_id)
probesetdata<-annotationfile[m,]
x <- readUnits(CelSet,units=u)
ind <- 1 # only read the first unit!
this <- x[[ind]]
thisunitname <- names(x)[1] #gene_id
thisgroupnames <- names(this) #probeset_ids
## build intensitymatrix ##
intensities<-matrix(ncol=nbrOfArrays(CelSet),nrow=sum(probesetdata
$probe_count))
lastrow<-0
for (i in 1:nrow(probesetdata)){
thisgroupname <- thisgroupnames [i]
thisnrow <- nrow(x[[thisunitname]][[thisgroupname]][[data]])
thisnprobe <- probesetdata$probe_count[i]
if (thisnrow!=thisnprobe){ #e.g. when working with firma scores
intensities[seq(lastrow+1,lastrow+thisnprobe),]<-matrix(rep(x
[[thisunitname]][[thisgroupname]][[data]],thisnprobe/
thisnrow),ncol=nbrOfArrays(CelSet),nrow=thisnprobe,byrow = TRUE) #no
check here, if (thisnprobe%%thisnrow==0)!!!
}else{
intensities[seq(lastrow+1,lastrow+thisnprobe),]<-x[[thisunitname]]
[[thisgroupname]][[data]]
}
lastrow <- lastrow+thisnprobe
}
intensities <- cbind(intensities ,apply(intensities[,groups[[1]]],
1,mean),apply(intensities [,groups[[2]]],1,mean)) #add means
intensities <- log2(intensities)
## DONE: build intensitymatrix ##
return(list
(id=id,idname=idname,probesetdata=probesetdata,intensities=intensities))
}
The part ## build intensitymatrix ## should be corrected/improved!!!
### SETUP ###
source(paste(getwd(),"/#Scripts_DA/LowLevelAnalysis_ensb.R",sep=""))
#do low-level analysis here and provide variables celsetN,fs,plmEx
library(GenomeGraphs)
mart = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
CelSet<-celsetN #celsetN <- process(QuantileNormalization(celsetBC,
typesToUpdate="pm"))
CelSetFirma<-fs #fs <- getFirmaScores(firma)
CelSetCes<-getChipEffectSet(plmEx) #plmEx <- ExonRmaPlm(celsetN,
mergeGroups=FALSE)
file <- paste(getwd(),"/annotationData/CSV/HuEx.small.csv",sep="")
#HuEx.small.csv is condensed with awk/gawk
psCsv <- read.table(file,skip=1,sep=" ",comment.char="")
colnames(psCsv) <- c("probeset_id", "start", "stop", "probe_count")
groups <- list(group1=c(1:4),group2=c(5:8))
### GENERATE PLOTDATA ###
id<-"ENSG00000030582"; idname<-"PGRN"
plotdata1<-getPlotDataGenomeGraphs
(id=id,idname=idname,CelSet=CelSet,annotationfile=psCsv,groups=groups)
plotdata2<-getPlotDataGenomeGraphs
(id=id,CelSet=CelSetCes,annotationfile=psCsv,groups=groups,data="theta")
plotdata2$intensities<-plotdata2$intensities[,(9:10)] #remove
everything but the means
plotdata3<-getPlotDataGenomeGraphs
(id=id,CelSet=CelSetFirma,annotationfile=psCsv,groups=groups)
plotdata3$intensities<-plotdata3$intensities[,(9:10)] #remove
everything but the means
plotdata3$intensities<- cbind
(plotdata3$intensities,plotdata3$intensities[,1]-plotdata3$intensities
[,2]) #add the diff of means
### MAYBE PLOTDATA ###
file <- paste(getwd(),"/#Scripts_DA/",id,".plotdata",sep="")
if(TRUE){
save(plotdata,plotdata2,plotdata3,file=file)
}else{
load(file)
}
### WRITE PLOTS TO FILE ###
fileformat<-"pdf"
plotpath <- "/#Scripts_DA/"
plotfilename <- paste("GenomeGraphs_",id,sep="")
if (fileformat=="pdf"){pdf(file = paste(getwd
(),plotpath,plotfilename,".pdf", sep=""), width = 12, height = 12,
pointsize = 10, bg = "transparent", family="Times")
}else{png(file = paste(getwd(),plotpath,plotfilename,".png", sep=""),
width = 1600, height = 1200, units = "px", pointsize = 10, bg =
"transparent", res = NA, restoreConsole = TRUE)}
colgroup1="royalblue3"
colgroup2="springgreen4"
coldiff="red"
colgene="orange"
coltranscript="lightskyblue"
col=c(rep("gray70",8),colgroup1,colgroup2)
lwd=c(rep(1,8),rep(2,2))
title = new("Title", title =paste(plotdata1$id,plotdata1$idname,sep="
- "), dp = DisplayPars(color = "darkslategray"))
exon1 = new("ExonArray", intensity = plotdata1$intensities, probeStart
= plotdata1$probesetdata[,2], probeEnd=plotdata1$probesetdata[,3],
probeId = as.character(plotdata1$probesetdata[,1]), nProbes =
plotdata1$probesetdata[,4], dp = DisplayPars(color = col, lwd=lwd,
mapColor = "gray80", plotMap=FALSE), displayProbesets=FALSE)
exon2 = new("ExonArray", intensity = plotdata2$intensities, probeStart
= plotdata2$probesetdata[,2], probeEnd=plotdata2$probesetdata[,3],
probeId = as.character(plotdata2$probesetdata[,1]), nProbes =
plotdata2$probesetdata[,4], dp = DisplayPars(color = c
(colgroup1,colgroup2), lwd=2, mapColor = "gray80", plotMap=FALSE),
displayProbesets=TRUE)
exon3 = new("ExonArray", intensity = plotdata3$intensities, probeStart
= plotdata3$probesetdata[,2], probeEnd=plotdata3$probesetdata[,3],
probeId = as.character(plotdata3$probesetdata[,1]), nProbes =
plotdata3$probesetdata[,4], dp = DisplayPars(color = c
(colgroup1,colgroup2,coldiff), lwd=c(1,1,2), mapColor = "gray80",
plotMap=TRUE), displayProbesets=FALSE)
gene = new("Gene", id = plotdata1$id , biomart = mart,dp = DisplayPars
(color = colgene))
transcript = new("Transcript", id = plotdata1$id , biomart = mart, dp
= DisplayPars(color = coltranscript,plotId=TRUE))
legend = new("Legend", legend = c("mean(group1)","mean(group2)","mean
(group1)-mean(group2)","gene","transcripts"), dp = DisplayPars(color= c
(colgroup1,colgroup2,coldiff,colgene,coltranscript)))
genomeaxis <- new("GenomeAxis", add53 = TRUE)
gdPlot(list(title,exon1,exon2,exon3, gene, genomeaxis, transcript,
legend), minBase = min(ex...@probestart), maxBase=max(ex...@probeend))
#no y-axis labels here
dev.off()
In the section ### SAVE PLOTDATA ### you could save your generated
plotdata. If you want to change your plot later (color, fonts,
highlight-box), you simply can load your data, which is much faster...
Like Mark said, there are a lot customizations possible, but maybe
this is a start to generate your own plots.
Frank
--~--~---------~--~----~------------~-------~--~----~
When reporting problems on aroma.affymetrix, make sure 1) to run the latest
version of the package, 2) to report the output of sessionInfo() and
traceback(), and 3) to post a complete code example.
You received this message because you are subscribed to the Google Groups
"aroma.affymetrix" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/aroma-affymetrix?hl=en
-~----------~----~----~----~------~----~------~--~---
