### R code from vignette source 'staRank.Rnw' ################################################### ### code chunk number 1: staRank.Rnw:50-51 ################################################### set.seed(42) ################################################### ### code chunk number 2: staRank.Rnw:53-65 ################################################### # sample parameters for data p<-20 # genes n<-4 # replicates trueEffects<-rnorm(p,0,5) # gene effect s<-rgamma(p,shape=2.5,rate=0.5) # sample variance # draw n replicate values of a gene simData <- matrix(0,nr=p,nc=n, dimnames=list(paste("Gene",c(1:p)),paste("Replicate",c(1:n)))) for(i in 1:p){ simData[i,]<-replicate(n,rnorm(1,mean=trueEffects[i],sd=s[i])) } ################################################### ### code chunk number 3: staRank.Rnw:72-76 ################################################### # load the stability ranking package library(staRank) # implemented ranking methods method<-c('mean','median','mwtest','ttest','RSA') ################################################### ### code chunk number 4: staRank.Rnw:80-86 ################################################### stabilityList<-list() stabilityList[['mean']]<-stabilityRanking(simData,method='mean') stabilityList[['median']]<-stabilityRanking(simData,method='median') stabilityList[['mwtest']]<-stabilityRanking(simData,method='mwtest') stabilityList[['ttest']]<-stabilityRanking(simData,method='ttest') stabilityList[['RSA']]<-stabilityRanking(simData,method='RSA') ################################################### ### code chunk number 5: staRank.Rnw:95-108 ################################################### # 1. create a ranking extRanking <- sort(apply(simData,1,mean)) # 2. create 100 sample rankings sampleRankings<-matrix(NA,nr=nrow(simData),nc=100) for(i in 1:100){ sampleRankings[,i]<- apply(simData,1,function(x){mean(sample(x,length(x),replace=TRUE))}) } rownames(sampleRankings)<-rownames(simData) sampleRankings<-sampleRankings[names(extRanking),] # 3. run stabilityRanking stabilityList[['externalMean']]<- stabilityRanking(extRanking,sampleRankings,method='externalMean') ################################################### ### code chunk number 6: staRank.Rnw:116-120 ################################################### #results for the mean ranking baseRank(stabilityList$mean) stabRank(stabilityList$mean) avrgRank(stabilityList$mean) ################################################### ### code chunk number 7: staRank.Rnw:127-128 ################################################### rankCor(stabilityList[['mean']]) ################################################### ### code chunk number 8: staRank.Rnw:134-135 ################################################### stableSetSize(stabilityList[['mean']]) ################################################### ### code chunk number 9: toyStabilityPlot ################################################### # plot the stable genes per cutoff k* for each method mCol<-rainbow(5) name_out<-c('Mean','Median','Rank sum','t-test','RSA') plot(1,type='n',xlim=c(0,p),ylim=c(0,p),xlab='k*',ylab='stable genes') abline(0,1,col='lightgray') for(i in 5:1){ lines(stableSetSize(stabilityList[[i]]),lwd=2,col=mCol[i]) } legend('topleft',name_out,col=mCol,lwd=2) ################################################### ### code chunk number 10: staRank.Rnw:155-156 ################################################### # plot the stable genes per cutoff k* for each method mCol<-rainbow(5) name_out<-c('Mean','Median','Rank sum','t-test','RSA') plot(1,type='n',xlim=c(0,p),ylim=c(0,p),xlab='k*',ylab='stable genes') abline(0,1,col='lightgray') for(i in 5:1){ lines(stableSetSize(stabilityList[[i]]),lwd=2,col=mCol[i]) } legend('topleft',name_out,col=mCol,lwd=2) ################################################### ### code chunk number 11: staRank.Rnw:176-181 ################################################### #load the data set data(salmonellaInfection) data<-as.matrix(salmonellaInfection) down<-rownames(data[which(apply(data,1,mean)<=-1),]) data<-data[down,] ################################################### ### code chunk number 12: staRank.Rnw:183-187 (eval = FALSE) ################################################### ## salmonellaStability<-list() ## salmonellaStability[['mean']]<-stabilityRanking(data,method='mean') ## salmonellaStability[['mwtest']]<-stabilityRanking(data,method='mwtest') ## salmonellaStability[['RSA']]<-stabilityRanking(data,method='RSA') ################################################### ### code chunk number 13: staRank.Rnw:188-189 (eval = FALSE) ################################################### ## save(salmonellaStability,file='data/salmonellaStability.RData') ################################################### ### code chunk number 14: staRank.Rnw:191-192 ################################################### data(salmonellaStability) ################################################### ### code chunk number 15: staRank.Rnw:196-197 ################################################### rankCor(salmonellaStability$mwtest) ################################################### ### code chunk number 16: staRank.Rnw:202-205 ################################################### # gather all rankings in a matrix rankings<-lapply(salmonellaStability,getRankmatrix) rankings<-Reduce(cbind,lapply(rankings,function(x){x[down,]})) ################################################### ### code chunk number 17: dataplot ################################################### # plot data according to rankings par(mfrow=c(1,5),mar=c(4,4,2,2)) range<-range(data) for(m in c(2,3,1,4,8)){ plot(1,type='n',xlim=c(0,length(down)), ylim=range,main='', ylab='infection score', xlab=paste(colnames(rankings)[m],'rank')) orderedData<-data[order(rankings[,m]),] for(j in 1:3){ points(orderedData[,j],pch=20,cex=0.5) } points(apply(orderedData,1,mean),col='red',pch=20,cex=0.5) } ################################################### ### code chunk number 18: staRank.Rnw:227-228 ################################################### # plot data according to rankings par(mfrow=c(1,5),mar=c(4,4,2,2)) range<-range(data) for(m in c(2,3,1,4,8)){ plot(1,type='n',xlim=c(0,length(down)), ylim=range,main='', ylab='infection score', xlab=paste(colnames(rankings)[m],'rank')) orderedData<-data[order(rankings[,m]),] for(j in 1:3){ points(orderedData[,j],pch=20,cex=0.5) } points(apply(orderedData,1,mean),col='red',pch=20,cex=0.5) } ################################################### ### code chunk number 19: infectionStabilityPlot ################################################### # plot the stable genes per cutoff k* for each method mCol<-rainbow(3) name_out<-c('Mean','Rank sum','RSA') plot(1,type='n',xlim=c(0,length(down)),ylim=c(0,length(down)), xlab='k*',ylab='stable genes') abline(0,1,col='lightgray') for(i in 1:3){ lines(stableSetSize(salmonellaStability[[i]]),lwd=2,col=mCol[i]) } legend('topleft',name_out,col=mCol,lwd=2) ################################################### ### code chunk number 20: staRank.Rnw:255-256 ################################################### # plot the stable genes per cutoff k* for each method mCol<-rainbow(3) name_out<-c('Mean','Rank sum','RSA') plot(1,type='n',xlim=c(0,length(down)),ylim=c(0,length(down)), xlab='k*',ylab='stable genes') abline(0,1,col='lightgray') for(i in 1:3){ lines(stableSetSize(salmonellaStability[[i]]),lwd=2,col=mCol[i]) } legend('topleft',name_out,col=mCol,lwd=2) ################################################### ### code chunk number 21: staRank.Rnw:265-266 ################################################### summary(salmonellaStability$mean)