#Suplementary material: #A Flexible Two Stage Procedure for Identifying Gene Sets that are Differentially Expressed #Ruth Heller, Elisabetta Manduchi, Gregory Grant, Warren Ewens. #---R code for the two stage procedure 3.1----------# library(multtest) library(MASS) #--------------------THE FUNCTIONS------------------# ################################################################################ # Author of the following 3 functions: Dan Nettleton # # Date: August 7, 2007 # # Downloaded from: http://www.public.iastate.edu/~dnett/useR/GeneCatTesting.txt# ################################################################################ e.com=function (y) {#Nettleton K=ncol(y) scale=rep(0,K) for(i in 1:K){ D=dist(y[,i]) scale[i]=1/sum(D) } y=sweep(y,2,scale,FUN="*") return(y) } get.delta=function (x,D,n,N,w) {#Nettleton inter.sum=tapply(1:N,x,FUN="inter.dist",D) inter.mean=inter.sum*2/(n*(n-1)) delta=sum(w*inter.mean) return(delta) } inter.dist=function (index,D) {#Nettleton return(sum(as.matrix(D)[index,index])/2) } mrpp=function (x,y,nperm) { #Nettleton's p-value for a gene set. D=dist(y) x=as.factor(x) n=table(x) N=sum(n) w=n/N delta.obs=get.delta(x,D,n,N,w) delta=rep(0,nperm) for(i in 1:nperm){ delta[i]=get.delta(sample(x),D,n,N,w) } return(p=mean(delta.obs>=delta)) } ################################################################################ # End of Dan Nettleton's code # ################################################################################ fdr = function(pv,q){ if (sum(is.na(pv))>0 | sum(pv<0,na.rm=T)>0 | sum(pv>1,na.rm=T)>0) warning("Some p-values are not-valid or missing! These p-values are ignored") end spv = sort(pv[!is.na(pv)]) G = length(spv) temp = spv<=seq(1,G,1)*q/G pID = ifelse(length(spv[temp==T])>0,max(spv[temp==T]),NA) pID } #-----------------------THE DATA----------------------# #Input: #set2geneindex: a matrix of size Sxn, where S is the number of gene sets and n the size of the largest gene set. # entry (i,j) is the gene id of the j-th gene in the i-th gene set. #ysub: a matrix of size Gxm, where G is the number of genes that are in at least one gene set and m is the number of replicates. entry (i,j) is the expression level of gene i in replicate j. #mlabels: a vector of size m, with entries 0 or 1 corresponding to replicates from the first or the second group respectively. #THE CODE FOR THE DATA USED IN THE APPLICATION TO A MICROARRAY STUDY IN SECTION 4 IS AT THE END OF THIS SCRIPT. #----------------------THE ANALYSIS-----------------# #----screening stage---# #compute screening p-values pmrpp = rep(NA,S) for (i in 1: S){#for every gene set print(paste("gene set ",i)) pmrpp[i] = mrpp(mylabels,e.com(t(ysub[set2geneindex[i,1:n[i]],])),nperm=100) # pmrpp[i] = mrpp(type.b,e.com(t(y[set2geneindex[i,1:n[i]],])),nperm=100) }#for i #BH on screening p-values q=0.05 pID = fdr(pmrpp,q) index = seq(1,S,1) index = index[pmrpp<=pID] length(index) #-----testing after screening stage---# pWY = matrix(nrow=S,ncol=500) #the adjusted p-values by WY for (k in 1: length(index)){ i = index[k] #the index of the rejected gene set print(paste("selected gene set ",k)) out = mt.minP(ysub[set2geneindex[i,1:n[i]],],mylabels,test = "wilcoxon", B=100000) # out = mt.minP(ysub[set2geneindex[i,1:n[i]],],mylabels,test = "wilcoxon",side="upper", B=100000) pWY[i,out\$index] = out\$adjp } save(mylabels, mysample,pmrpp,index, pWY, file ="out2sided") }#for Iindex #--------output----# load("out2sided") S = dim(pWY)[1] q=0.05 alpha = length(index)/S*q numGenesPerSet = apply(pWY<=alpha, 1,sum, na.rm=T) cbind(n[numGenesPerSet>0],numGenesPerSet[numGenesPerSet>0], 100*(numGenesPerSet[numGenesPerSet>0]/n[numGenesPerSet>0])) table(numGenesPerSet[numGenesPerSet>0]) #can display the number of genes detected per gene set in a table. temp = set2geneindex[pWY<=alpha] temp= temp[!is.na(temp)] length(unique(temp)) #--------compare to BH on all genes----# G = length(unique(set2geneindex[!is.na(set2geneindex)])) I = length(mysample)/2 mylabels = c(rep(0,I), rep(1,I)) pw = rep(NA,G) ysub = y[unique(set2geneindex[!is.na(set2geneindex)]),] ysub = ysub[,mysample] for (j in 1:dim(ysub)[1]){ mydata = ysub[j,] pw[j] = wilcox.test(mydata[mylabels==0], mydata[mylabels==1])\$p.value # pw[j] = wilcox.test(mydata[mylabels==0], mydata[mylabels==1],alternative="less")\$p.value }#for j q=0.05 pID = fdr(pw,q) if (!is.na(pID)){ sum(pw<=pID) } #---------output table of discoveries---# mysetdiscoveries = (names.goBP2hgu[setID]) mygeneswithinsetdiscoveries = matrix(nrow = S,ncol= dim(set2geneindex)[2]) numgenediscoveries = rep(0,S) for (i in 1:length(index)){ pv = pWY[index[i],] pv = pv[!is.na(pv)] mysetIDs = set2geneID[index[i],1:length(pv)] if (sum(pv<=alpha)){ numgenediscoveries[index[i]] = sum(pv<=alpha) mygeneswithinsetdiscoveries[index[i],1:sum(pv<=alpha)] = mysetIDs[pv<=alpha] } } rbind(numgenediscoveries[index],n[index]) temp = numgenediscoveries[index]/n[index] sub = c(index[temp==1],index[temp==0.75],index[round(temp,3)==0.667]) numgenediscoveries[sub] mysetdiscoveries[sub] for (i in 1: length(sub)){ print(mygeneswithinsetdiscoveries[sub[i],1:numgenediscoveries[sub[i]]]) } ##################################################################### #--THE DATA FOR THE APPLICATION TO A MICROARRAY STUDY IN SECTION 4--# ##################################################################### #---the gene sets---# library(GO) packageDescription('GO', field='Version') library(hgu95av2) packageDescription('hgu95av2', field='Version') go<-as.list(GOTERM) goBP<-character() i<-1 for (j in 1:length(go)) { if (Ontology(go[[j]]) == "BP") { goBP[i]<-GOID(go[[j]]) names(goBP)[i]<-Term(go[[j]]) i<-i+1 } } length(goBP) hgu<-as.list(hgu95av2GO2ALLPROBES) goINhgu<-character() for (i in 1:length(hgu)) { goINhgu[i]<-names(hgu[i]) } length(hgu) length(goINhgu) goBP2hgu<-list() names.goBP2hgu<-character() i<-1 for (j in 1:length(goBP)) { v<-which(goINhgu==goBP[j]) if (length(v)==1) { goBP2hgu[i]<-hgu[v[1]] names.goBP2hgu[i]<-paste(goBP[j],names(goBP)[j],sep="\t") i<-i+1 } } length(goBP2hgu) length(names.goBP2hgu) names.goBP2hgu[1] goBP2hgu[[1]] length(goBP2hgu[])#4364 gene sets nS = rep(0,length(goBP2hgu[]))#the gene set sizes for (i in 1:length(goBP2hgu[])){ nS[i] = length(unique(goBP2hgu[[i]])) } S = length(nS[nS>1 & nS<=500])# 3367 gene sets of size 2 to 500 - consider these for analysis/ n = rep(0,S) set2geneID = matrix(nrow = S, ncol=500) setID = rep(0,S) #the index of the gene set in goBP2hgu k=1 for (i in 1:length(goBP2hgu[])){ if (nS[i]>1 & nS[i]<=500){ set2geneID[k,1:nS[i]] = unique(goBP2hgu[[i]]) setID[k] =i n[k] = nS[i] k=k+1 } } rm("nS") #---the expression data---# source("http://bioconductor.org/biocLite.R") biocLite("ALL") library(ALL) data(ALL) show(ALL) print(summary(pData(ALL))) dim(exprs(ALL)) set2geneindex = matrix(nrow = S, ncol=500) #for a gene set the rows in expression data that correspond to it for (i in 1:S){ set2geneindex[i,1:n[i]] = pmatch(set2geneID[i,1:n[i]],featureNames(ALL)) } set2geneindex[i,] set2geneID[i,] #The B type and T type in ALLtype. phenoData(ALL)\$BT ALLtype=as.factor(substr(phenoData(ALL)\$BT,1,1)) #Define a new factor that has B- or T-cell ALL type designation. type.b=(ALLtype=="B"); # Indicator for B cell type y=exprs(ALL) #----------------------THE ANALYSIS-----------------# Bindex = sample(seq(1,95,1),30) Tindex = sample(seq(96,128,1),30) Ivec = c(5,10,15,20,25,30) for (Iindex in 1:length(Ivec)){ I = Ivec[Iindex] mysample = c(Bindex[1:I], Tindex[1:I]) ysub = y[,mysample] dim(ysub) mylabels = c(rep(0,I), rep(1,I))