##' Compare \code{\link{IcaSet}} objects by computing the correlation between either
##' projection values of common features or genes, or contributions of common
##' samples.
##'
##' The user must carefully choose the object on which the correlation will be
##' computed.
##' If \code{level='samples'}, the correlations are based on the mixing
##' matrices of the ICA decompositions (of dimension samples x components). \code{'A'}
##' will be typically chosen when the ICA decompositions were computed on the
##' same dataset, or on datasets that include the same samples.
##' If \code{level='features'} is
##' chosen, the correlation is calculated between the source matrices (of dimension
##' features x components) of the ICA decompositions. \code{'S'} will be typically
##' used when the ICA decompositions share common features (e.g same microarrays).
##' If \code{level='genes'}, the correlations are calculated
##' on the attributes \code{'SByGene'} which store the
##' projections of the annotated features. \code{'SByGene'} will be typically chosen
##' when ICA were computed on datasets from different technologies, for which
##' comparison is possible only after annotation into a common ID, like genes.
##'
##' \code{cutoff_zval} is only used when \code{level} is one of \code{c('genes','features')}, in
##' order to restrict the correlation to the contributing features or genes.
##'
##' When \code{cutoff_zval} is specified, for each pair of components, genes or features that are
##' included in the circle of center 0 and radius \code{cutoff_zval} are excluded from
##' the computation of the correlation.
##'
##' It must be taken into account by the user that if
##' \code{cutoff_zval} is different from \code{NULL} or \code{0}, the computation will be much
##' slowler since each pair of component is treated individually.
##'
##' @title Comparison of IcaSet objects using correlation
##' @param icaSets list of IcaSet objects, e.g results of ICA decompositions
##' obtained on several datasets.
##' @param labAn vector of names for each icaSet, e.g the the names of the
##' datasets on which were calculated the decompositions.
##' @param type.corr Type of correlation to compute, either
##' \code{'pearson'} or \code{'spearman'}.
##' @param cutoff_zval either NULL or 0 (default) if all genes are used to
##' compute the correlation between the components, or a threshold to compute
##' the correlation on the genes that have at least a scaled projection higher
##' than cutoff_zval. Will be used only when correlations are calculated on S or
##' SByGene.
##' @param level Data level of the \code{IcaSet} objects on which is applied the correlation.
##' It must correspond to a feature shared by the IcaSet objects:
##' \code{'samples'} if they were applied to common samples (correlations are computed between matrix \code{A}), \code{'features'} if they were applied to
##' common features (correlations are computed between matrix \code{S}), \code{'genes'} if they share gene IDs after
##' annotation into genes (correlations are computed between matrix \code{SByGene}).
##' @return A list whose length equals the number of pairs of \code{IcaSet} and whose elements
##' are outputs of function \code{\link{cor2An}}.
##' @export
##' @examples
##'
##' dat1 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat1) <- paste("g", 1:1000, sep="")
##' colnames(dat1) <- paste("s", 1:10, sep="")
##' dat2 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat2) <- paste("g", 1:1000, sep="")
##' colnames(dat2) <- paste("s", 1:10, sep="")
##'
##' ## run ICA
##' resJade1 <- runICA(X=dat1, nbComp=3, method = "JADE")
##' resJade2 <- runICA(X=dat2, nbComp=3, method = "JADE")
##'
##' ## build params
##' params <- buildMineICAParams(resPath="toy/")
##'
##' ## build IcaSet object
##' icaSettoy1 <- buildIcaSet(params=params, A=data.frame(resJade1$A), S=data.frame(resJade1$S),
##' dat=dat1, alreadyAnnot=TRUE)$icaSet
##' icaSettoy2 <- buildIcaSet(params=params, A=data.frame(resJade2$A), S=data.frame(resJade2$S),
##' dat=dat2, alreadyAnnot=TRUE)$icaSet
##'
##' listPairCor <- compareAn(icaSets=list(icaSettoy1,icaSettoy2), labAn=c("toy1","toy2"),
##' type.corr="pearson", level="genes", cutoff_zval=0)
##'
##'
##' \dontrun{
##' #### Comparison of 2 ICA decompositions obtained on 2 different gene expression datasets.
##' ## load the two datasets
##' library(breastCancerMAINZ)
##' library(breastCancerVDX)
##' data(mainz)
##' data(vdx)
##'
##' ## Define a function used to build two examples of IcaSet objects
##' treat <- function(es, annot="hgu133a.db") {
##' es <- selectFeatures_IQR(es,10000)
##' exprs(es) <- t(apply(exprs(es),1,scale,scale=FALSE))
##' colnames(exprs(es)) <- sampleNames(es)
##' resJade <- runICA(X=exprs(es), nbComp=10, method = "JADE", maxit=10000)
##' resBuild <- buildIcaSet(params=buildMineICAParams(), A=data.frame(resJade$A), S=data.frame(resJade$S),
##' dat=exprs(es), pData=pData(es), refSamples=character(0),
##' annotation=annot, typeID= typeIDmainz,
##' chipManu = "affymetrix", mart=mart)
##' icaSet <- resBuild$icaSet
##' }
##' ## Build the two IcaSet objects
##' icaSetMainz <- treat(mainz)
##' icaSetVdx <- treat(vdx)
##'
##' ## The pearson correlation is used as a measure of association between the gene projections
##' # on the different components (type.corr="pearson").
##' listPairCor <- compareAn(icaSets=list(icaSetMainz,icaSetVdx),
##' labAn=c("Mainz","Vdx"), type.corr="pearson", level="genes", cutoff_zval=0)
##'
##' ## Same thing but adding a selection of genes on which the correlation between two components is computed:
##' # when considering pairs of components, only projections whose scaled values are not located within
##' # the circle of radius 1 are used to compute the correlation (cutoff_zval=1).
##' listPairCor <- compareAn(icaSets=list(icaSetMainz,icaSetVdx),
##' labAn=c("Mainz","Vdx"), type.corr="pearson", cutoff_zval=1, level="genes")
##' }
##' @export
##' @author Anne Biton
##' @seealso \code{\link{cor2An}}
compareAn <- function (icaSets,
labAn,
type.corr = c("pearson","spearman"),
cutoff_zval=0,
level = c("samples","features","genes")
) {
if (missing(labAn))
if (is.null(names(icaSets)))
labAn <- c(1:nbAn)
else
labAn <- names(icaSets)
nbAn <- length(icaSets)
if (length(labAn) != nbAn)
stop("labAn must have the same length than the list icaSets")
level <- match.arg(tolower(level), choices = c("features","genes","samples"))
switch(level,
features={ object="S"},
genes={object="SByGene"},
samples={object="A"}
)
type.corr <- match.arg(type.corr)
nbCompByAn <- lapply(icaSets, function(x) ncol(A(x)))
matlist <-
lapply(icaSets,
function(icaSet, object) {
data.frame(icaSet[object],check.names = FALSE, stringsAsFactors = FALSE)
}
, object = object
)
allpairs <- combn(x=1:nbAn,m = 2, simplify = FALSE)
pairnames <- unlist(lapply(allpairs, function(x,labAn) paste(labAn[x[1]],labAn[x[2]],sep="-"), labAn = labAn))
listCorPair <-
lapply(allpairs,
function(pa, matlist, labAn, cutoff_zval, type.corr) {
mat1 <- matlist[[pa[1]]]
mat2 <- matlist[[pa[2]]]
resCor <- cor2An(mat1, mat2, lab = c(labAn[pa[1]],labAn[pa[2]]), type.corr = type.corr, cutoff_zval = cutoff_zval)
}
, matlist = matlist
, labAn = labAn
, type.corr = type.corr
, cutoff_zval = if (missing(cutoff_zval)) NA else cutoff_zval
)
names(listCorPair) <- pairnames
return(listCorPair)
}
##' This function measures the correlation between two matrices containing the results of
##' two decompositions.
##'
##' Before computing the correlations, the components are scaled and restricted
##' to common row names.
##'
##' It must be taken into account by the user that if \code{cutoff_zval} is different
##' from NULL or zero, the computation will be slowler since each pair of
##' component is treated individually.
##'
##' When \code{cutoff_zval} is specified, for each
##' pair of components, genes that are included in the circle of center 0 and
##' radius \code{cutoff_zval} are excluded from the computation of the correlation
##' between the gene projection of the two components.
##' @title Correlation between two matrices
##' @param mat1 matrix of dimension features/genes x number of components, e.g
##' the results of an ICA decomposition
##' @param mat2 matrix of dimension features/genes x number of components, e.g
##' the results of an ICA decomposition
##' @param lab The vector of labels for mat1 and mat2, e.g the the names of the
##' two datasets on which were calculated the two decompositions
##' @param type.corr Type of correlation, either \code{'pearson'} or
##' \code{'spearman'}
##' @param cutoff_zval cutoff_zval: 0 (default) if all genes are
##' used to compute the correlation between the components, or a threshold to
##' compute the correlation on the genes that have at least a scaled projection
##' higher than cutoff_zval.
##' @return This function returns a list consisting of:
##' \item{cor}{a matrix of dimensions '(nbcomp1+nbcomp2) x (nbcomp1*nbcomp2)',
##' containing the correlation values between each pair of components,}
##' \item{pval}{ a matrix of dimension '(nbcomp1+nbcomp2) x (nbcomp1*nbcomp2)',
##' containing the p-value of the correlation tests for each
##' pair of components,}
##' \item{inter}{ the intersection between the features/genes of \code{mat1}
##' and \code{mat2},}
##' \item{labAn}{ the labels of the compared matrices.}
##' @author Anne Biton
##' @export
##' @examples
##' cor2An(mat1=matrix(rnorm(10000),nrow=1000,ncol=10), mat2=matrix(rnorm(10000),nrow=1000,ncol=10),
##' lab=c("An1","An2"), type.corr="pearson")
##'
##' @seealso \code{rcorr}, \code{cor.test}, \code{\link{compareAn}}
cor2An <- function (mat1,
mat2,
lab,
type.corr = c("pearson","spearman"),
cutoff_zval = 0
) {
if (is.null(cutoff_zval))
cutoff_zval <- 0
if (cutoff_zval < 0)
stop("'cutoff_zval' must be positive")
comp1 <- comp2 <- NULL
mat1 <- data.frame(scale(mat1), check.names = FALSE, stringsAsFactors = FALSE)
mat2 <- data.frame(scale(mat2), check.names = FALSE, stringsAsFactors = FALSE)
nbcomp1 <- ncol(mat1)
nbcomp2 <- ncol(mat2)
gg <- intersect(rownames(mat1),rownames(mat2))
if (length(gg)<10)
warning(paste("Less than 10 elements are common between the icaSet objects ", lab[1], " and ", lab[2], ", the comparison of these two decompositions should be reconsidered.", sep = ""))
mat1 <- mat1[gg,]
mat2 <- mat2[gg,]
if (missing(lab))
lab <- paste("an",c(1:2),sep = "")
if (cutoff_zval == 0) {
r <- signif(rcorr(as.matrix(mat1), as.matrix(mat2), type = type.corr)$r, digits = 5)
dimnames(r) <- list(c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")),c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")))
rpval <- NULL
rpval <- rcorr(as.matrix(mat1),as.matrix(mat2),type=type.corr)$P
dimnames(rpval) <- list(c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")),c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")))
}
else {
## transform matrix into lists
l1 <- lapply(as.list(mat1),
function(x,n) {
names(x) <- n
return(x)
}
, n = rownames(mat1)
)
l2 <- lapply(as.list(mat2),
function(x,n) {
names(x) <- n
return(x)
}
, n = rownames(mat2)
)
rpval <- matrix(ncol = nbcomp2+nbcomp1 , nrow =nbcomp2+nbcomp1)
r <- matrix(ncol = nbcomp2+nbcomp1 , nrow =nbcomp2+nbcomp1)
for (i in 1:nbcomp1) {
reslist <-
sapply(1:nbcomp2,
function(j, i, l1, l2, cutoff_zval, gg) {
dc <- data.frame(comp1=l1[[i]],comp2=l2[[j]])
rownames(dc) <- gg
g2del <- rownames(subset(dc, comp1^2 + comp2^2<= cutoff_zval ))
inter <- gg[!(gg %in% g2del)]
res <- cor.test(as.matrix(mat1[inter,i]),as.matrix(mat2[inter,j]), method = type.corr)
return(list(cor=res$estimate,pval=res$p.value))
}
, i = i
, l1 = l1
, l2 = l2
, cutoff_zval = cutoff_zval
, gg = gg)
reslist <- as.data.frame(reslist)
r[i,(nbcomp1+1):(nbcomp1+nbcomp2)] <- r[(nbcomp1+1):(nbcomp1+nbcomp2),i] <- unlist(reslist["cor",])
rpval[i,(nbcomp1+1):(nbcomp1+nbcomp2)] <- rpval[(nbcomp1+1):(nbcomp1+nbcomp2),i] <- unlist(reslist["pval",])
for (j in 1:nbcomp1) {
dc <- data.frame(comp1=l1[[i]],comp2=l1[[j]])
rownames(dc) <- gg
g2del <- rownames(subset(dc, comp1^2 + comp2^2<= cutoff_zval ))
inter <- names(l1[[1]])[!(names(l1[[1]]) %in% g2del)]
res <- cor.test(as.matrix(mat1[inter,i]),as.matrix(mat1[inter,j]), method = type.corr)
rpval[i,j] <- rpval[j,i] <- res$p.value
r[i,j] <- r[j,i] <- res$estimate
}
}
for (i in 1:nbcomp2) {
for (j in 1:nbcomp2) {
dc <- data.frame(comp1=l2[[i]],comp2=l2[[j]])
rownames(dc) <- gg
g2del <- rownames(subset(dc, comp1^2 + comp2^2<= cutoff_zval ))
inter <- names(l2[[1]])[!(names(l2[[1]]) %in% g2del)]
res <- cor.test(as.matrix(mat2[inter,i]),as.matrix(mat2[inter,j]), method = type.corr)
rpval[nbcomp1+i,nbcomp1+j] <- rpval[nbcomp1+j,nbcomp1+i] <- res$p.value
r[nbcomp1+i,nbcomp1+j] <- r[nbcomp1+j,nbcomp1+i] <- res$estimate
}
}
dimnames(rpval) <- list(c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")),c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")))
dimnames(r) <- list(c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")),c(paste("comp",c(1:nbcomp1), "_", lab[1], sep = ""), paste("comp",c(1:nbcomp2), "_", lab[2], sep = "")))
}
return(list(cor=r,pval=rpval,inter = gg, nbComp = c(nbcomp1,nbcomp2), labAn = lab))
}
correl2Comp <- function (
### This function computes the correlation between two components.
comp1
### The first component, a vector of projections or contributions indexed by labels
,comp2
### The second component, a vector of projections or contributions indexed by labels
,type.corr = "pearson"
### name of the type of correlation to compute, either 'pearson' or 'spearman'
,plot = FALSE
### if TRUE the plot of comp1 vs comp2 is drawn
,cutoff_zval = 0
### either NULL or 0 (default) if all genes are used to compute the correlation between the components, or a threshold to compute the correlation on the genes that have at least a scaled projection higher than cutoff_zval.
,test = FALSE
### if TRUE the p-value of correlation is returned instead of the correlation value
,alreadyTreat = FALSE
### if TRUE comp1 and comp2 are considered as being already treated (i.e scaled and restricted to common elements)
) {
##details<< Before computing the correlation, the components are scaled and restricted to common labels.
## When \code{cutoff_zval} is different from 0, the elements that are included in the circle of center 0 and radius \code{cutoff_zval} are not taken into account during the computation of the correlation.
if(!alreadyTreat) {
comp1 <- (comp1-mean(comp1))/sd(comp1)
comp2 <- (comp2-mean(comp2))/sd(comp2)
gg <- intersect(names(comp1),names(comp2))
comp1 <- comp1[gg]
comp2 <- comp2[gg]
}
dc <- data.frame(comp1=comp1,comp2=comp2)
rownames(dc) <- gg
g2del <- rownames(subset(dc, comp1^2 + comp2^2<= cutoff_zval ))
genes.intersect <- gg[!(gg %in% g2del)]
if (length(genes.intersect) < 4)
if (test) return(1) else return(0)
comp1.ord = comp1[genes.intersect]
comp2.ord = comp2[genes.intersect]
if (test) {
r = cor.test(comp1.ord,comp2.ord,method = type.corr, alternative = "two.sided")$p.value
}
else {
r = cor(comp1.ord,comp2.ord,method = type.corr)
r = signif(r, digits = 3)
}
if (plot) plot(comp1.ord,comp2.ord,xlab = paste("comp1"),ylab=paste("comp2"),pch = 16,sub = paste("corr_",type.corr," = ",r,sep=""))
return(r)
### This function returns either the correlation value or the p-value of the correlation test.
}
##' This function builds a data.frame describing for each node of the graph
##' its ID and which analysis/data it comes from.
##'
##' The created file is used in Cytoscape.
##' @title Generate node attributes
##' @param nbAn Number of analyses being considered, i.e number of IcaSet
##' objects
##' @param nbComp Number of components by analysis, if of length 1 then
##' it is assumed that each analysis has the same number of components.
##' @param labAn Labels of the analysis, if missing it will be generated
##' as an1, an2, ...
##' @param labComp List containing the component labels indexed by analysis, if missing
##' will be generated as comp1, comp2, ...
##' @param file File where the description of the node attributes will be
##' written
##' @return A data.frame describing each node/component
##' @export
##' @examples
##' ## 4 datasets, 20 components calculated in each dataset, labAn
##' nodeAttrs(nbAn=4, nbComp=20, labAn=c("tutu","titi","toto","tata"))
##'
##' @author Anne Biton
nodeAttrs <- function(nbAn,
nbComp,
labAn,
labComp,
file
) {
nb <- an <- lab <- comp <- NULL
if (missing(labAn))
labAn <- paste(rep("an",nbAn),1:nbAn,sep = "")
else if (length(labAn) != nbAn)
stop("The length of 'labAn' must equal 'nbAn'.")
if (missing(labComp)) {
if (length(nbComp) == 1) {
labid <-
paste(
rep(paste(rep("comp",nbComp),c(1:nbComp),sep = ""),nbAn),
paste(rep("an",nbComp*nbAn),sapply(c(1:nbAn),rep,nbComp),sep=""),
sep = "_"
)
labComp <-rep(paste(rep("comp",nbComp),c(1:nbComp),sep = ""),nbAn)
nbComp <- rep(nbComp,nbAn)
}
else {
if (length(nbComp) != nbAn)
stop("The provided number of components must equal the number of analyses")
labid <- foreach(nb = nbComp, an = labAn) %do% {paste(paste(rep("comp",nb),c(1:nb),sep = ""),rep(an,nb),sep="_")}
labid <- unlist(labid)
labComp <- unlist(foreach(nb = nbComp, an = 1:nbAn) %do% {paste(rep("comp",sum(nb)),c(1:nb),sep = "")})
}
}
else {
if (!is.list(labComp))
stop("The component labels must be provided using a list.")
if (length(labComp) != nbAn)
stop("The length of the list containing the component labels 'labComp' must have the same length than the list containing the analysis labels 'labAn'.")
names(labComp) <- names(labAn)
foreach(lab = labComp, nb = nbComp) %do% {
if (length(lab) != nb)
stop("The length of the list 'labComp' containing component labels does not map to 'nbComp'.")
}
labid <- unlist(foreach(nb = nbComp, an = labAn, comp = labComp) %do% {paste(comp,rep(an,nb),sep="_")})
labComp <- unlist(labComp)
}
compnb <- unlist(lapply(nbComp,function(x) c(1:x)))
names(nbComp) <- labAn
an <- unlist(lapply(names(nbComp), function(x,nb) rep(x,nb[x]), nb = nbComp))
dataAttr <- data.frame(id = labid, labComp = labComp, indComp = compnb, labAn = an, stringsAsFactors = FALSE, check.names = FALSE)
if (!missing(file))
if (!is.null(file))
write.table(dataAttr,file=file,sep = " ", row.names = FALSE, quote = FALSE)
return(dataAttr)
}
##' compareAn2graphfile
##'
##' This function builds a correlation graph from the outputs of function \code{\link{compareAn}}.
##'
##' When correlations are considered (\code{useVal}="cor"), absolute values
##' are used since the components have no direction.
##'
##' If \code{useMax} is \code{TRUE} each component is linked to the most correlated component of
##' each different \code{IcaSet}.
##'
##' If \code{cutoff} is specified, only
##' correlations exceeding this value are taken into account during the graph construction.
##' For example, if \code{cutoff} is 1, only relationships between components
##' that correspond to a correlation value larger than 1 will be included.
##'
##' When \code{useVal="pval"} and \code{useMax=TRUE}, the minimum value is taken
##' instead of the maximum.
##'
##' @param listPairCor The output of the function \code{\link{compareAn}}, containing the
##' correlation between several pairs of objects of class \code{\link{IcaSet}}.
##' @param useMax If TRUE, the graph is restricted to edges that correspond to
##' maximum score, see details
##' @param cutoff Cutoff used to select pairs that will be included in the
##' graph.
##' @param useVal The value on which is based the graph, either \code{"cor"} for
##' correlation or \code{"pval"} for p-values of correlation tests.
##' @param file File name.
##' @return A data.frame with the graph description, has
##' two columns \code{n1} and \code{n2} filled with node IDs, each row denotes that there is an edge from \code{n1} to \code{n2}. Additional columns quantify the strength of association: correlation (\code{cor}), p-value (\code{pval}), (\code{1-abs(cor)}) (\code{distcor}), log10-pvalue (\code{logpval}).
##' @author Anne Biton
##' @seealso \code{\link{compareAn}}, \code{\link{cor2An}}
##' @export
##' @examples
##'
##' dat1 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat1) <- paste("g", 1:1000, sep="")
##' colnames(dat1) <- paste("s", 1:10, sep="")
##' dat2 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat2) <- paste("g", 1:1000, sep="")
##' colnames(dat2) <- paste("s", 1:10, sep="")
##'
##' ## run ICA
##' resJade1 <- runICA(X=dat1, nbComp=3, method = "JADE")
##' resJade2 <- runICA(X=dat2, nbComp=3, method = "JADE")
##'
##' ## build params
##' params <- buildMineICAParams(resPath="toy/")
##'
##' ## build IcaSet object
##' icaSettoy1 <- buildIcaSet(params=params, A=data.frame(resJade1$A), S=data.frame(resJade1$S),
##' dat=dat1, alreadyAnnot=TRUE)$icaSet
##' icaSettoy2 <- buildIcaSet(params=params, A=data.frame(resJade2$A), S=data.frame(resJade2$S),
##' dat=dat2, alreadyAnnot=TRUE)$icaSet
##'
##' resCompareAn <- compareAn(icaSets=list(icaSettoy1,icaSettoy2), labAn=c("toy1","toy2"),
##' type.corr="pearson", level="genes", cutoff_zval=0)
##'
##' ## Build a graph where edges correspond to maximal correlation value (useVal="cor"),
##' compareAn2graphfile(listPairCor=resCompareAn, useMax=TRUE, useVal="cor", file="myGraph.txt")
##'
##'
##' \dontrun{
##' #### Comparison of 2 ICA decompositions obtained on 2 different gene expression datasets.
##' ## load the two datasets
##' library(breastCancerMAINZ)
##' library(breastCancerVDX)
##' data(mainz)
##' data(vdx)
##'
##' ## Define a function used to build two examples of IcaSet objects
##' treat <- function(es, annot="hgu133a.db") {
##' es <- selectFeatures_IQR(es,10000)
##' exprs(es) <- t(apply(exprs(es),1,scale,scale=FALSE))
##' colnames(exprs(es)) <- sampleNames(es)
##' resJade <- runICA(X=exprs(es), nbComp=10, method = "JADE", maxit=10000)
##' resBuild <- buildIcaSet(params=buildMineICAParams(), A=data.frame(resJade$A), S=data.frame(resJade$S),
##' dat=exprs(es), pData=pData(es), refSamples=character(0),
##' annotation=annot, typeID= typeIDmainz,
##' chipManu = "affymetrix", mart=mart)
##' icaSet <- resBuild$icaSet
##' }
##' ## Build the two IcaSet objects
##' icaSetMainz <- treat(mainz)
##' icaSetVdx <- treat(vdx)
##'
##' ## Compute correlation between every pair of IcaSet objects.
##' resCompareAn <- compareAn(icaSets=list(icaSetMainz,icaSetVdx),
##' labAn=c("Mainz","Vdx"), type.corr="pearson", level="genes", cutoff_zval=0)
##'
##' ## Same thing but adding a selection of genes on which the correlation between two components is computed:
##' # when considering pairs of components, only projections whose scaled values are not located within
##' # the circle of radius 1 are used to compute the correlation (cutoff_zval=1).
##' resCompareAn <- compareAn(icaSets=list(icaSetMainz,icaSetVdx),
##' labAn=c("Mainz","Vdx"), type.corr="pearson", cutoff_zval=1, level="genes")
##'
##' ## Build a graph where edges correspond to maximal correlation value (useVal="cor"),
##' ## i.e, component A of analysis i is linked to component B of analysis j,
##' ## only if component B is the most correlated component to A amongst all component of analysis j.
##' compareAn2graphfile(listPairCor=resCompareAn, useMax=TRUE, useVal="cor", file="myGraph.txt")
##'
##' ## Restrict the graph to correlation values exceeding 0.4
##' compareAn2graphfile(listPairCor=resCompareAn, useMax=FALSE, cutoff=0.4, useVal="cor", file="myGraph.txt")
##'
##' }
compareAn2graphfile <- function (listPairCor,
useMax = TRUE,
cutoff = NULL,
useVal = c("cor","pval"),
file = NULL
) {
useVal <- match.arg(useVal)
dataGraphs <-
lapply(listPairCor,
function(res, useMax, cutoff, useVal) {
dd <- abs(res[[useVal]])
nbcomp1 <- res$nbComp[1]
nbcomp2 <- res$nbComp[2]
attrGraph <- c("cor",if ("pval" %in% names(res)) "pval")
if (useMax) {
if (useVal == "pval")
op <- "which.min"
else
op <- "which.max"
## find correlation max both direction an1 <-> an2 since the max is not necessarily reciprocal
subdd <- dd[(nbcomp1+1):(nbcomp1+nbcomp2), 1:nbcomp1]
ind <- apply(subdd,2,op)
coord <- as.list(data.frame(t(cbind(ind,1:nbcomp1))))
dataGraph <- data.frame(n1=rownames(dd)[1:nbcomp1],
n2=rownames(dd)[(nbcomp1+1):(nbcomp1+nbcomp2)][ind])
dataGraph[[useVal]] <- unlist(lapply(coord, function(x,m) m[x[1],x[2]], m = subdd))
addval <- setdiff(attrGraph,useVal)
if (length(addval)==1)
dataGraph[[addval]] <- unlist(lapply(coord, function(x,m) m[x[1],x[2]], m = res[[addval]][(nbcomp1+1):(nbcomp1+nbcomp2), 1:nbcomp1]))
## an2 -> an1
subdd <- dd[1:nbcomp1,(nbcomp1+1):(nbcomp1+nbcomp2)]
ind <- apply(subdd,2,op)
coord <- as.list(data.frame(t(cbind(ind,1:nbcomp2))))
dataGraph2 <- data.frame(n1=rownames(dd)[(nbcomp1+1):(nbcomp1+nbcomp2)],
n2=rownames(dd)[1:nbcomp1][ind])
dataGraph2[[useVal]] <- unlist(lapply(coord, function(x,m) m[x[1],x[2]], m = subdd))
addval <- setdiff(attrGraph,useVal)
if (length(addval)==1)
dataGraph2[[addval]] <- unlist(lapply(coord, function(x,m) m[x[1],x[2]], m = res[[addval]][1:nbcomp1,(nbcomp1+1):(nbcomp1+nbcomp2)]))
dataGraph <- rbind(dataGraph,dataGraph2)
pval <- NULL
if (!is.null(cutoff)) {
if (useVal == "pval")
dataGraph <- subset(dataGraph, pval <= cutoff)
else
dataGraph <- subset(dataGraph, cor >= cutoff)
}
}
else if (!is.null(cutoff)) {
subdd <- dd[1:nbcomp1, (nbcomp1+1):(nbcomp1+nbcomp2)]
dataGraph <-
lapply(as.list(colnames(subdd)),
function(nn, cutoff, subdd, useVal, attrGraph, res,op, subInd) {
x <- subdd[,nn]
ind <- which(eval(parse(text=paste("x",op,"cutoff"))))
if (length(ind) > 0) {
dataGraph <- data.frame(n1 = rep(nn,length(ind)), n2 = rownames(subdd)[ind])
dataGraph[[useVal]] <- subdd[ind,nn]
addval <- setdiff(attrGraph,useVal)
if (length(addval)==1) {
coord <- as.list(data.frame(t(cbind(dataGraph$n1,dataGraph$n2))))
dataGraph[[addval]] <- unlist(lapply(coord, function(x,m) m[x[1],x[2]], m = res[[addval]][subInd[[1]],subInd[[2]]]))
}
return(dataGraph)
}
}
, cutoff = cutoff
, subdd = subdd
, useVal = useVal
, attrGraph = attrGraph
, res = res
, op = if (useVal == "pval") "<=" else ">="
, subInd = list(1:nbcomp1, (nbcomp1+1):(nbcomp1+nbcomp2))
)
dataGraph <- do.call(rbind,dataGraph)
}
return(dataGraph)
}
, useMax = useMax
, cutoff = if (missing(cutoff)) NULL else cutoff
, useVal = useVal
)
dataGraphs <- do.call(rbind,dataGraphs)
dataGraphs$invcor = 1/(10*as.numeric(dataGraphs$cor))
dataGraphs$distcor = 1-abs(as.numeric(dataGraphs$cor))
dataGraphs[,c("cor","pval","invcor","distcor")] <- apply(dataGraphs[,c("cor","pval","invcor","distcor")],2,as.numeric)
dataGraphs <- annotReciprocal(dataGraph = dataGraphs, keepOnlyReciprocal = FALSE)# file
dataGraphs[,c("cor","pval","invcor","distcor")] <- apply(dataGraphs[,c("cor","pval","invcor","distcor")],2,as.numeric)
if (!missing(file))
if (!is.null(file))
write.table(dataGraphs, file = file, row.names = FALSE, quote = FALSE, sep = "\t")
return(dataGraphs)
}
##' annotReciprocal
##'
##' This function notes edges of a graph as reciprocal or not.
##'
##'
##' @param dataGraph data.frame which contains the graph description, must have
##' two columns \code{n1} and \code{n2} filled with node IDs, each row denoting there is an edge from \code{n1} to \code{n2}.
##' @param file file where the graph description is written
##' @param keepOnlyReciprocal if TRUE \code{dataGraph} is restricted to
##' reciprocal edges, else all edges are kept (default).
##' @return This function returns the argument \code{dataGraph} with an
##' additional column named 'reciprocal' which contains TRUE if the edge
##' described by the row is reciprocal, and FALSE if it is not reciprocal.
##' @examples
##' dg <- data.frame(n1=c("A","B","B","C","C","D","E","F"),n2=c("B","G","A","B","D","C","F","E"))
##' annotReciprocal(dataGraph=dg)
##'
##' @export
##' @author Anne Biton
annotReciprocal <- function (dataGraph,
file,
keepOnlyReciprocal = FALSE
) {
names(dataGraph)[1:2] <- c("n1","n2")
dataGraph$keep <- rep(NA, length=nrow(dataGraph))
dataGraph_keep <- apply( dataGraph, MARGIN = 1 ,
function (r, dataGraph) {
dataN2 = subset(dataGraph, dataGraph$n1 == r["n2"])
if (r["n1"] %in% dataN2$n2) r["keep"] = TRUE
else r["keep"] = FALSE
return(r)
}
, dataGraph = dataGraph
)
dataGraph_keep <- as.data.frame(t(dataGraph_keep), check.names = FALSE, stringsAsFactors = FALSE)
names(dataGraph_keep)[ncol(dataGraph_keep)] <- "reciprocal"
if (keepOnlyReciprocal) {
dataGraph_keep<- subset(dataGraph_keep, dataGraph_keep$reciprocal == TRUE )
dataGraph_keep <- dataGraph_keep[,-ncol(dataGraph_keep)]
}
if (!missing(file))
if(!is.null(file))
write.table(dataGraph_keep, file = file, sep = "\t", quote = FALSE, row.names = FALSE)
return(dataGraph_keep)
}
##' This function plots the
##' correlation graph in an interactive device using function \code{tkplot}.
##'
##' You have to slighly move the nodes to see cliques because strongly related nodes are often superimposed.
##' The \code{edgeWeight} column is used to weight the edges within the
##' fruchterman.reingold layout available in the package \code{igraph}.
##'
##' The argument
##' \code{nodeCol} typically denotes the column containing the names of the datasets.
##' Colors are automatically
##' attributed to the nodes using palette Set3 of package \code{RColorBrewer}. The
##' corresponding colors can be directly specified in the 'col' argument. In
##' that case, 'col' must be a vector of colors indexed by the unique elements
##' contained in \code{nodeCol} column (e.g dataset ids).
##'
##' As for colors, one can define
##' the column of \code{nodeAttrs} that is used to define the node shapes. The
##' corresponding shapes can be directly specified in the \code{shape} argument. In
##' that case, \code{shape} must be one of \code{c("circle","square", " vcsquare", "rectangle", "crectangle", "vrectangle")} and must be
##' indexed by the unique elements of \code{nodeShape} column.
##'
##' Unfortunately, shapes
##' can't be taken into account when tkplot is TRUE (interactive plot).
##'
##' If \code{reciproCol} is not missing, it is used to color the edges, either in grey if the
##' edge is not reciprocal or in black if the edge is reciprocal.
##'
##'
##' @title Plots graph using
##' @param dataGraph A data.frame containing the graph description. It must
##' have two columns \code{n1} and \code{n2}, each row denoting that there is an edge from n1
##' to n2. Node labels in columns \code{n1} and \code{n2} of \code{dataGraph} must correspond to
##' node IDs in column \code{id} of \code{nodeAttrs}.
##' @param edgeWeight The column of dataGraph used to weight edges.
##' @param nodeAttrs A data.frame with node description, see function \code{nodeAttrs}.
##' @param nodeShape Denotes the column of \code{nodeAttrs} used to attribute the node shapes.
##' @param nodeCol Denotes the column of \code{nodeAttrs} used to
##' color the nodes in the graph.
##' @param nodeName Denotes the column of \code{nodeAttrs} used
##' as labels for the nodes in the graph.
##' @param col A vector of colors, for the nodes, indexed by the unique elements of \code{nodeCol}
##' column from \code{nodeAttrs}. If missing, colors will be automatically attributed.
##' @param shape A vector of shapes indexed by the unique elements of
##' column \code{nodeShape} from \code{nodeAttrs}. If missing, shapes will be automatically
##' attributed.
##' @param title Title for the plot
##' @param reciproCol Denotes the column of \code{dataGraph} containing \code{TRUE} if the
##' row defines a reciprocal node, else \code{FALSE}. See \code{\link{annotReciprocal}}.
##' @param tkplot If TRUE, performs interactive plot with function \code{tkplot}, else uses \code{plot.igraph}.
##' @param \dots Additional parameters as required by \code{tkplot}.
##' @return A list consisting of \describe{
##' \item{dataGraph}{a data.frame defining the correlation
##' graph}
##' \item{nodeAttrs}{a data.frame describing the node
##' of the graph}
##' \item{graph}{the graph as an object of class \code{igraph}}
##' \item{graphid}{the id of the graph plotted using \code{tkplot}} }
##' @export
##' @author Anne Biton
##' @seealso \code{\link{compareAn}}, \code{\link{nodeAttrs}}, \code{\link{compareAn2graphfile}}, \code{\link{runCompareIcaSets}}
##' @examples
##'
##' dat1 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat1) <- paste("g", 1:1000, sep="")
##' colnames(dat1) <- paste("s", 1:10, sep="")
##' dat2 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat2) <- paste("g", 1:1000, sep="")
##' colnames(dat2) <- paste("s", 1:10, sep="")
##'
##' ## run ICA
##' resJade1 <- runICA(X=dat1, nbComp=3, method = "JADE")
##' resJade2 <- runICA(X=dat2, nbComp=3, method = "JADE")
##'
##' ## build params
##' params <- buildMineICAParams(resPath="toy/")
##'
##' ## build IcaSet object
##' icaSettoy1 <- buildIcaSet(params=params, A=data.frame(resJade1$A), S=data.frame(resJade1$S),
##' dat=dat1, alreadyAnnot=TRUE)$icaSet
##' icaSettoy2 <- buildIcaSet(params=params, A=data.frame(resJade2$A), S=data.frame(resJade2$S),
##' dat=dat2, alreadyAnnot=TRUE)$icaSet
##' icaSets <- list(icaSettoy1, icaSettoy2)
##'
##' resCompareAn <- compareAn(icaSets=list(icaSettoy1,icaSettoy2), labAn=c("toy1","toy2"),
##' type.corr="pearson", level="genes", cutoff_zval=0)
##'
##' ## Build a graph where edges correspond to maximal correlation value (useVal="cor"),
##' dataGraph <- compareAn2graphfile(listPairCor=resCompareAn, useMax=TRUE, useVal="cor", file="myGraph.txt")
##'
##' ## construction of the data.frame with the node description
##' nbComp <- rep(3,2) #each IcaSet contains 3 components
##' nbAn <- 2 # we are comparing 2 IcaSets
##' # labels of components created as comp*i*
##' labComp <- foreach(icaSet=icaSets, nb=nbComp, an=1:nbAn) %do% {
##' paste(rep("comp",sum(nb)),1:nbComp(icaSet),sep = "")}
##'
##' # creation of the data.frame with the node description
##' nodeDescr <- nodeAttrs(nbAn = nbAn, nbComp = nbComp, labComp = labComp,
##' labAn = c("toy1","toy2"), file = "nodeInfo.txt")
##'
##' ## Plot correlation graph, slightly move the attached nodes to make the cliques visible
##' ## use tkplot=TRUE to have an interactive graph
##' res <- plotCorGraph(title = "Compare toy 1 and 2", dataGraph = dataGraph, nodeName = "indComp", tkplot = FALSE,
##' nodeAttrs = nodeDescr, edgeWeight = "cor", nodeShape = "labAn", reciproCol = "reciprocal")
##'
##'
##' \dontrun{
##' ## load two microarray datasets
##' library(breastCancerMAINZ)
##' library(breastCancerVDX)
##' data(mainz)
##' data(vdx)
##'
##' ## Define a function used to build two examples of IcaSet objects
##' treat <- function(es, annot="hgu133a.db") {
##' es <- selectFeatures_IQR(es,10000)
##' exprs(es) <- t(apply(exprs(es),1,scale,scale=FALSE))
##' colnames(exprs(es)) <- sampleNames(es)
##' resJade <- runICA(X=exprs(es), nbComp=10, method = "JADE", maxit=10000)
##' resBuild <- buildIcaSet(params=buildMineICAParams(), A=data.frame(resJade$A), S=data.frame(resJade$S),
##' dat=exprs(es), pData=pData(es), refSamples=character(0),
##' annotation=annot, typeID= typeIDmainz,
##' chipManu = "affymetrix", mart=mart)
##' icaSet <- resBuild$icaSet
##' }
##' ## Build the two IcaSet objects
##' icaSetMainz <- treat(mainz)
##' icaSetVdx <- treat(vdx)
##'
##' icaSets <- list(icaSetMainz, icaSetVdx)
##' labAn <- c("Mainz", "Vdx")
##'
##' ## correlations between gene projections of each pair of IcaSet
##' resCompareAn <- compareAn(icaSets = icaSets, level = "genes", type.corr= "pearson",
##' labAn = labAn, cutoff_zval=0)
##'
##' ## construction of the correlation graph using previous output
##' dataGraph <- compareAn2graphfile(listPairCor=resCompareAn, useMax=TRUE, file="corGraph.txt")
##'
##' ## construction of the data.frame with the node description
##' nbComp <- rep(10,2) #each IcaSet contains 10 components
##' nbAn <- 2 # we are comparing 2 IcaSets
##' # labels of components created as comp*i*
##' labComp <- foreach(icaSet=icaSets, nb=nbComp, an=1:nbAn) %do% {
##' paste(rep("comp",sum(nb)),1:nbComp(icaSet),sep = "")}
##'
##' # creation of the data.frame with the node description
##' nodeDescr <- nodeAttrs(nbAn = nbAn, nbComp = nbComp, labComp = labComp,
##' labAn = labAn, file = "nodeInfo.txt")
##'
##' ## Plot correlation graph, slightly move the attached nodes to make the cliques visible
##' res <- plotCorGraph(title = "Compare two ICA decomsitions obtained on \n two
##' microarray-based data of breast tumors", dataGraph = dataGraph, nodeName = "indComp",
##' nodeAttrs = nodeDescr, edgeWeight = "cor", nodeShape = "labAn", reciproCol = "reciprocal")
##'
##' }
plotCorGraph <- function(
dataGraph,
edgeWeight = "cor",
nodeAttrs,
nodeShape,
nodeCol = "labAn",
nodeName = "indComp",
col,
shape,
title = "",
reciproCol = "reciprocal",
tkplot = FALSE,
...
) {
if (!(edgeWeight %in% colnames(dataGraph))) {
stop (paste(edgeWeight,"is not available within the columns of dataGraph."))
}
nodeName <- match.arg(nodeName, choices = colnames(nodeAttrs))
if (missing(nodeShape))
nodeAttrs$shape <- rep("circle",nrow(nodeAttrs))
else if (!(nodeShape %in% colnames(nodeAttrs))) {
nodeAttrs$shape <- rep("circle",nrow(nodeAttrs))
warning(paste("The column ", nodeShape, " is not available in 'dataGraph'.",sep = ""))
}
else {
potShapes <- c("circle","square", "csquare", "rectangle", "crectangle", "vrectangle")
autoShape <- TRUE
if (!missing(shape)) {
autoShape <- FALSE
if (length(intersect(names(shape),unique(nodeAttrs[[nodeShape]]))) != length(unique(nodeAttrs[[nodeShape]]))) {
warning("'shape' argument is incorrect, some elements of nodeAttrs$'nodeShape' are not indexed. Shapes will be attributed automatically.")
autoShape <- TRUE
}
if (length(setdiff(shape,potShapes))>0) {
warning("'shape' argument is incorrect, the available shapes are: 'circle','square', 'csquare', 'rectangle', 'crectangle', and 'vrectangle'. Shapes will be attributed automatically.")
autoShape <- TRUE
}
}
if (autoShape) {
levs <- unique(nodeAttrs[[nodeShape]])
if (length(levs) > length(potShapes)) {
warning("R graphs can only handle 6 node shapes at the most.")
shape <- rep(potShapes,6)[1:length(levs)]
}
else
shape <- potShapes[1:length(levs)]
names(shape) <- levs
nodeAttrs$shape <- shape[nodeAttrs[[nodeShape]]]
}
else
nodeAttrs$shape <- shape[nodeAttrs[[nodeShape]]]
}
if (!(nodeCol %in% colnames(nodeAttrs))) {
stop(paste("The column ",nodeCol, " is not available in nodeAttrs.",sep=""))
}
else {
nbAn <- length(unique(nodeAttrs[[nodeCol]]))
if (!missing(col) && !is.null(col)) {
if (length(intersect(names(col),unique(nodeAttrs[[nodeCol]]))) != length(unique(nodeAttrs[[nodeCol]]))) {
warning("'col' argument is incorrect, some elements of nodeAttrs$'nodeCol' are not indexed. Colors will be attributed automatically.")
colAn <- brewer.pal(nbAn,"Set3")
names(colAn) <- unique(nodeAttrs[[nodeCol]])
nodeAttrs$col = colAn[nodeAttrs[[nodeCol]]]
}
else
nodeAttrs$col <- col[nodeAttrs[[nodeCol]]]
}
else {
colAn <- brewer.pal(nbAn,"Set3")
names(colAn) <- unique(nodeAttrs[[nodeCol]])
nodeAttrs$col = colAn[nodeAttrs[[nodeCol]]]
}
}
n1 <- NULL
nodes <- as.character(unique(as.character(dataGraph$n1)))
edges <- llply(nodes,
function(n,data,edgeWeight) {
ll <- list(edges = as.character(subset(data, n1 == n)$n2), weights = subset(data, n1 == n)[[edgeWeight]])
},
data = dataGraph, edgeWeight = edgeWeight)
names(edges) <- nodes
g <- new("graphNEL", nodes = nodes, edgeL = edges, edgemode = "directed")
# build igraph object from graphNEL object
ig <- igraph.from.graphNEL(g, name = TRUE, weight = TRUE)
V(ig)$color <- nodeAttrs$col[match(V(ig)$name,nodeAttrs$id)]
V(ig)$shape <- nodeAttrs$shape[match(V(ig)$name,nodeAttrs$id)]
E(ig)$width <- 10^E(ig)$weight
indHighCor <- which(abs(E(ig)$weight)>0.4)
indLowCor <- which(abs(E(ig)$weight)<=0.4)
E(ig)$weight <- 1/abs(log2(abs(E(ig)$weight)))
E(ig)$weight[indHighCor] <- E(ig)$weight[indHighCor]+6
E(ig)$weight[indLowCor] <- E(ig)$weight[indLowCor]-0.25
if (!missing(reciproCol)) {
if (reciproCol %in% colnames(dataGraph)) {
nodes <- as.character(unique(as.character(dataGraph$n1)))
colEdges <- llply(nodes,
function(n,data,reciproCol) {
c('TRUE' = "black", 'FALSE' = "gray50")[as.character(subset(data, n1 == n)[[reciproCol]])]
},
data = dataGraph, reciproCol = reciproCol)
colEdges <- unlist(colEdges)
E(ig)$color <- colEdges
}
else {
warning(paste("No column of 'dataGraph' has the name",reciproCol))
E(ig)$color <- "black"
}
}
else
E(ig)$color <- "black"
lay <- layout.fruchterman.reingold(ig,niter=500,area=vcount(ig)^2.3,repulserad=vcount(ig)^100, weights = E(ig)$weight)
graph <- ig
if (!capabilities()[["X11"]])
tkplot <- FALSE
if (tkplot) {
if(capabilities()[["X11"]] & tolower(Sys.info()["sysname"])!="windows") {
dimScreen <- system("xdpyinfo | grep dimensions", intern = TRUE)
dimScreen <- as.numeric(strsplit(gsub(gsub(gsub(dimScreen, pattern = " ", replacement = ""), pattern = "[A-Za-z0-9]*:", replacement = ""), pattern = "pixels\\([A-Za-z0-9]*\\)", replacement = ""), split = "x")[[1]])
}
else
dimScreen <- c(450,450)
graphid <- tkplot(ig,
layout = lay,
vertex.label = nodeAttrs[[nodeName]][match(V(ig)$name,nodeAttrs$id)],
vertex.shape = nodeAttrs$shape[match(V(ig)$name,nodeAttrs$id)],
canvas.width=dimScreen[1], canvas.height=dimScreen[2],
...)
tkplot.fit.to.screen(graphid)
}
else {
graph <- plot(ig,
layout = lay,
vertex.label = nodeAttrs[[nodeName]][match(V(ig)$name,nodeAttrs$id)],
vertex.shape = nodeAttrs$shape[match(V(ig)$name,nodeAttrs$id)],
main = title,
...
)
graph <- ig
graphid = ""
}
return(list(dataGraph=dataGraph, nodeAttrs = nodeAttrs, graph = graph, graphid=graphid))
}
##' runCompareIcaSets
##'
##' This function encompasses the comparison of several IcaSet objects using correlations
##' and the plot of the corresponding correlation graph.
##' The IcaSet objects are compared by calculating the correlation between either
##' projection values of common features or genes, or contributions of common
##' samples.
##'
##' This function calls four functions: \code{\link{compareAn}} which computes the
##' correlations, \code{\link{compareAn2graphfile}} which builds the graph,
##' \code{\link{nodeAttrs}} which builds the node description data, and \code{\link{plotCorGraph}}
##' which uses tkplot to plot the graph in an interactive device.
##'
##' If the user wants to see the correlation graph in Cytoscape, he must fill
##' the arguments \code{fileDataGraph} and \code{fileNodeDescr}, in order to
##' import the graph and its node descriptions as a .txt file in Cytoscape.
##'
##' When \code{labAn} is missing, each element i of \code{icaSets} is labeled as
##' 'Ani'.
##'
##' The user must carefully choose the data level used in the comparison:
##' If \code{level='samples'}, the correlations are based on the mixing
##' matrices of the ICA decompositions (of dimension samples x components). \code{'A'}
##' will be typically chosen when the ICA decompositions were computed on the
##' same dataset, or on datasets that include the same samples.
##' If \code{level='features'} is
##' chosen, the correlation is calculated between the source matrices (of dimension
##' features x components) of the ICA decompositions. \code{'S'} will be typically
##' used when the ICA decompositions share common features (e.g same microarrays).
##' If \code{level='genes'}, the correlations are calculated
##' on the attributes \code{'SByGene'} which store the
##' projections of the annotated features. \code{'SByGene'} will be typically chosen
##' when ICA were computed on datasets from different technologies, for which
##' comparison is possible only after annotation into a common ID, like genes.
##'
##' \code{cutoff_zval}
##' is only used when \code{level} is one of \code{c('features','genes')}, in
##' order to restrict the correlation to the contributing features or genes.
##'
##' When \code{cutoff_zval} is specified, for each pair of components, genes or features that are
##' included in the circle of center 0 and radius \code{cutoff_zval} are excluded from
##' the computation of the correlation.
##'
##' It must be taken into account by the user that if cutoff_zval
##' is different from NULL or zero, the computation will be much slowler since
##' each pair of component is treated individually.
##'
##' Edges of the graph are built based on the correlation values between the
##' components. Absolute values of correlations are used since
##' components have no direction.
##'
##' If \code{useMax} is \code{TRUE} each component will be linked to only one component of
##' each other IcaSet that corresponds to the most correlated component among
##' all components of the same IcaSet. If \code{cutoff_graph} is specified, only
##' correlations exceeding this value are taken into account to build the graph.
##' For example, if \code{cutoff} is 1, only relationships between components
##' that correspond to a correlation value higher than 1 will be included.
##' Absolute correlation values are used since the components have no direction.
##'
##' The contents of the returned list are \describe{
##' \item{dataGraph:}{\code{dataGraph} data.frame that describes the correlation
##' graph,} \item{nodeAttrs:}{\code{nodeAttrs} data.frame that describes the node
##' of the graph} \item{graph}{\code{graph} the graph as an igraph-object,}
##' \item{graphid:}{\code{graphid} the id of the graph plotted using tkplot.} }
##'
##' @param icaSets List of \code{\link{IcaSet}} objects, e.g results of ICA decompositions
##' obtained on several datasets.
##' @param labAn Vector of names for each icaSet, e.g the the names of the
##' datasets on which were calculated the decompositions.
##' @param type.corr Type of correlation to compute, either
##' \code{'pearson'} or \code{'spearman'}.
##' @param cutoff_zval Either NULL or 0 (default) if all genes are used to
##' compute the correlation between the components, or a threshold to compute
##' the correlation using the genes that have at least a scaled projection higher
##' than cutoff_zval. Will be used only when \code{level} is one of \code{c("features","genes")}.
##' @param level Data level of the \code{IcaSet} objects on which is applied the correlation.
##' It must correspond to a data level shared by the IcaSet objects:
##' \code{'samples'} if they were applied to common samples (correlations are computed between matrix
##' \code{A}), \code{'features'} if they were applied to common features (correlations are computed between matrix \code{S}), \code{'genes'}
##' if they share gene IDs after
##' annotation into genes (correlations are computed between matrix \code{SByGene}).
##' @param fileNodeDescr File where node descriptions are saved (useful when the
##' user wants to visualize the graph using Cytoscape).
##' @param fileDataGraph File where graph description is saved (useful when the
##' user wants to visualize the graph using Cytoscape).
##' @param plot if \code{TRUE} (default) plot the correlation graph
##' @param title title of the graph
##'
##' @param col vector of colors indexed by elements of labAn; if missing, colors
##' will be automatically attributed
##' @param cutoff_graph the cutoff used to select pairs that will be included in
##' the graph
##' @param useMax if \code{TRUE}, the graph is restricted to edges that correspond to
##' maximum correlation between components, see details
##' @param tkplot If TRUE, performs interactive plot with function \code{tkplot}, else uses \code{plot.igraph}
##' @return A list consisting of \describe{
##' \item{dataGraph:}{a data.frame defining the correlation
##' graph}
##' \item{nodeAttrs:}{a data.frame describing the node
##' of the graph,}
##' \item{graph:}{the graph as an object of class \code{igraph},}
##' \item{graphid}{the id of the graph plotted with \code{tkplot}}. }
##' @export
##' @author Anne Biton
##' @seealso \code{\link{compareAn2graphfile}}, \code{\link{compareAn}}, \code{\link{cor2An}}, \code{\link{plotCorGraph}}
##' @examples
##'
##' dat1 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat1) <- paste("g", 1:1000, sep="")
##' colnames(dat1) <- paste("s", 1:10, sep="")
##' dat2 <- data.frame(matrix(rnorm(10000),ncol=10,nrow=1000))
##' rownames(dat2) <- paste("g", 1:1000, sep="")
##' colnames(dat2) <- paste("s", 1:10, sep="")
##'
##' ## run ICA
##' resJade1 <- runICA(X=dat1, nbComp=3, method = "JADE")
##' resJade2 <- runICA(X=dat2, nbComp=3, method = "JADE")
##'
##' ## build params
##' params <- buildMineICAParams(resPath="toy/")
##'
##' ## build IcaSet objects
##' icaSettoy1 <- buildIcaSet(params=params, A=data.frame(resJade1$A), S=data.frame(resJade1$S),
##' dat=dat1, alreadyAnnot=TRUE)$icaSet
##' icaSettoy2 <- buildIcaSet(params=params, A=data.frame(resJade2$A), S=data.frame(resJade2$S),
##' dat=dat2, alreadyAnnot=TRUE)$icaSet
##'
##' ## compare IcaSet objects
##' ## use tkplot=TRUE to get an interactive graph
##' rescomp <- runCompareIcaSets(icaSets=list(icaSettoy1, icaSettoy2), labAn=c("toy1","toy2"),
##' type.corr="pearson", level="genes", tkplot=FALSE)
##'
##'
##' \dontrun{
##' ## load the microarray-based gene expression datasets
##' ## of breast tumors
##' library(breastCancerMAINZ)
##' library(breastCancerVDX)
##' data(mainz)
##' data(vdx)
##'
##' ## Define a function used to build two examples of IcaSet objects
##' ## and annotate the probe sets into gene Symbols
##' treat <- function(es, annot="hgu133a.db") {
##' es <- selectFeatures_IQR(es,10000)
##' exprs(es) <- t(apply(exprs(es),1,scale,scale=FALSE))
##' colnames(exprs(es)) <- sampleNames(es)
##' resJade <- runICA(X=exprs(es), nbComp=10, method = "JADE", maxit=10000)
##' resBuild <- buildIcaSet(params=buildMineICAParams(), A=data.frame(resJade$A), S=data.frame(resJade$S),
##' dat=exprs(es), pData=pData(es), refSamples=character(0),
##' annotation=annot, typeID= typeIDmainz,
##' chipManu = "affymetrix", mart=mart)
##' icaSet <- resBuild$icaSet
##' }
##' ## Build the two IcaSet objects
##' icaSetMainz <- treat(mainz)
##' icaSetVdx <- treat(vdx)
##'
##' ## compare the IcaSets
##' runCompareIcaSets(icaSets=list(icaSetMainz, icaSetVdx), labAn=c("Mainz","Vdx"), type.corr="pearson", level="genes")
##' }
runCompareIcaSets <- function(icaSets,
labAn,
type.corr = c("pearson","spearman"),
cutoff_zval=0,
level = c("genes","features","samples"),
fileNodeDescr = NULL,
fileDataGraph = NULL,
plot = TRUE,
title = "",
col,
cutoff_graph = NULL,
useMax = TRUE,
tkplot = FALSE
) {
nb <- NULL
if (missing(labAn))
labAn <- paste("An",1:length(icaSets))
else
if (length(labAn) != length(icaSets))
stop("Length of 'labAn' is different from length of 'icaSets'.")
if (!missing(col))
if (length(col) != length(icaSets))
stop("Length of 'col' is different from length of 'icaSets'.")
else
if (is.null(names(col)) | length(intersect(names(col), labAn)) != length(labAn))
names(col) <- labAn
type.corr <- type.corr[1]
resCompareAn <- compareAn(icaSets = icaSets, level = level,type.corr= type.corr, labAn = labAn, cutoff_zval=cutoff_zval)
dataGraph <- compareAn2graphfile(listPairCor=resCompareAn, useMax = useMax, file = fileDataGraph, cutoff = cutoff_graph)
nbComp <- unlist(lapply(icaSets,function(x) length(indComp(x))))
nbAn <- length(icaSets)
icaSet <- NULL
labComp <- (foreach(icaSet = icaSets, nb = nbComp, an = 1:nbAn) %do% {paste(rep("comp",sum(nb)),indComp(icaSet),sep = "")})
nodeDescr <- nodeAttrs(nbAn = nbAn, nbComp = nbComp, labComp = labComp, labAn = labAn, file = fileNodeDescr)
if (plot) {
res <- plotCorGraph(title = title, dataGraph = dataGraph, nodeName = "indComp", nodeAttrs = nodeDescr, edgeWeight = "cor", col = if(missing(col)) NULL else col, nodeShape = "labAn", reciproCol = "reciprocal", tkplot = tkplot)
return(list(dataGraph=res$dataGraph,
nodeAttrs=res$nodeAttrs,
graph=res$graph,
graphid=res$graphid
))
}
else
return(list(dataGraph=dataGraph, nodeAttrs = nodeAttrs))
}
