#load and install the fabia package source("http://bioconductor.org/biocLite.R") biocLite("fabia") require(fabia) export.fabia.threshL = function(fabiaRes, threshZ = 0.5) { noL <- as.matrix(L(fabiaRes)) nZ <- as.matrix(Z(fabiaRes)) X <- as.matrix(X(fabiaRes)) n <- nrow(noL) p <- ncol(noL) l <- ncol(nZ) mom <- 0 for (i in 1:p) { mom <- mom + sum(noL[, i]^2) * sum(nZ[i, ]^2) } mom <- mom/(as.double(n) * as.double(l) * as.double(p)) threshL <- sqrt(mom)/threshZ print(threshL) threshL } export.fabia = function(fabiaRes, baseName, clusterAssignments = NULL, clusterNames = NULL, threshZ = 0.5, threshL = NULL) { x = X(fabiaRes) #gene x sample l = L(fabiaRes) #gene x bicluster z = t(Z(fabiaRes)) #sample x bicluster n = nrow(l) p = ncol(l) q <- ncol(z) write.table(x,paste0(baseName,'_X.csv'),quote=FALSE,col.names=NA,dec=".",sep="\t") write.table(l,paste0(baseName,'_L.csv'),quote=FALSE,col.names=NA,dec=".",sep="\t") write.table(z,paste0(baseName,'_Z.csv'),quote=FALSE,col.names=NA,dec=".",sep="\t") #create a threshold table in the form #rows: bicluster ids #col 1 named L : L thresholds #col 2 named Z : Z thresholds threshZ = rep(threshZ,length.out=p) if (is.null(threshL)) { threshL = export.fabia.threshL(fabiaRes,threshZ) } threshL = rep(threshL,length.out=p) thresh = cbind(threshZ,threshL) rownames(thresh) = colnames(l) colnames(thresh) = c('Z','L') write.table(thresh,paste0(baseName,'_thresholds.csv'),quote=FALSE,col.names=NA,dec=".",sep="\t") if (!is.null(clusterAssignments)) { if (!is.null(clusterNames)) assignments = as.matrix(sapply(clusterAssignments,function(x) { clusterNames[x] })) else assignments = as.matrix(clusterAssignments) rownames(assignments) = colnames(x) colnames(assignments) = paste("Chemical Clustering",1:ncol(assignments)) write.table(assignments,paste0(baseName,'_chemicalClusters.csv'),quote=FALSE,col.names=NA,dec=".",sep="\t") } } export.fabiaData = function() { require(fabiaData) data(Breast_A) X <- as.matrix(XBreast) resBreast1 <-fabia(X,5,0.1,400) raBreast1 <- extractBic(resBreast1) exportBiCluster(resBreast1,"breast1",CBreast) data(Multi_A) X <- as.matrix(XMulti) resMulti1 <- fabia(X,5,0.06,300,norm=2) raMulti1 <- extractBic(resMulti1) exportBiCluster(resMulti1,"multi1",CMulti) data(DLBCL_B) X <- as.matrix(XDLBCL) resDLBCL1 <- fabia(X,5,0.1,400,norm=2) raDLBCL1 <- extractBic(resDLBCL1) exportBiCluster(resDLBCL1,"DLBCL1",CDLBCL) } #usage: #export.fabia(resFabia,"ws1") #export.fabia(resFabia,"screenOnco",clusterAssignments,clusterNames)