generatePACMENinput = function(Exprs,Mut,PPI=NULL,filename="PACMENinput.txt"){ # last updated 29 October 2016 # Exprs is a normalised gene expression matrix (can also be protein) ## Rows are genes and columns are samples # Mut is a mutation matrix, either binary (t-test), or non-negative score (linear regression) (can also be CNV) ## Rows are genes and columns are samples (columns must match with Exprs) ## You may want to preprocess Mut by removing genes with too few mutations before running this code ## by default this only considers genes with at least 3 non-zero entries, but you may want to filter this further. # PPI is network info as an adjacency matrix. If NULL then test all connections # filename to save tab-delimited results to, for uploading to PACMEN if (ncol(Exprs)!=ncol(Mut)) stop("Exprs and Mut must have same number of columns") if (is.null(PPI)) print("running with no interactome, testing all combinations...") require(limma) LMpairs = NULL #testableGenes = names(which(rowSums(Mut)>=3)) testableGenes = names(which(apply(Mut,1,function(x) sum(x!=0)>=3))) if (!is.null(PPI)) { testableGenes = testableGenes[testableGenes %in% rownames(PPI)] } print(length(testableGenes)) for (gene in testableGenes) { print(gene) if (!is.null(PPI)) { exprgenes = sort(c(gene,names(which(PPI[,gene]==1)))) exprgenes = exprgenes[exprgenes %in% rownames(Exprs)] } else { exprgenes = rownames(Exprs) } if (length(exprgenes)==0) next ysmall = Exprs[exprgenes,,drop=FALSE] rownames(ysmall) <- exprgenes if (nrow(ysmall)==0) next design = model.matrix(~Mut[gene,]) if (sum(design[,2])==0) next fit = lmFit(ysmall,design) fit = eBayes(fit) LMvalue = sign(fit$t[,2]) LMp = fit$p.value[,2] #LMpairsShort = cbind(tail=gene,head=names(LMvalue),pval=LMp,dir=LMvalue) LMpairsShort = cbind(tail=gene,head=exprgenes,pval=LMp,dir=LMvalue) if (ncol(LMpairsShort)!=4) LMpairsShort=NULL LMpairs = rbind(LMpairs,LMpairsShort) } res = LMpairs rownames(res) <- NULL #if (!is.null(filename)) { print("saving data to file") write.table(res,file=filename,row.names=FALSE,col.names=TRUE,sep="\t",quote=FALSE) #} #return(res) }