CNORfuzzy: Addon to CellNOptR: Fuzzy Logic

#
#  This file is part of the CNO software
#
#  Copyright (c) 2011-2012 - EBI - Massachusetts Institute of Technology
#
#  File author(s): CNO developers (cno-dev@ebi.ac.uk)
#
#  Distributed under the GPLv2 License.
#  See accompanying file LICENSE.txt or copy at
#      http://www.gnu.org/licenses/gpl-2.0.html
#
#  CNO website: http://www.ebi.ac.uk/saezrodriguez/software/cno
#
##############################################################################
##the global function is still called writeNetwork and still needs the same arguments
#(modelComprExpanded, optimResT1, optimResT2, modelOriginal, CNOlist), but it is divided into:
#    -a function writeSNetworkW that does the actual writing to files,
#    -a function that gets the info for the sif and dot files
writeFuzzyNetwork<-function(
    postRefThresh,
    allFinalMSEs,
    allRes,
    tag=NULL,verbose=FALSE){

ReferenceRes = allRes[[1]]
CNOlist = ReferenceRes$paramsList$data
allBits = matrix(NA,length(allRes),length(ReferenceRes$bit))
for (eachRes in 1:length(allRes)){
            currIX = max(which(allFinalMSEs[eachRes,]-allFinalMSEs[eachRes,2]<=postRefThresh))
            allBits[eachRes,] = allRes[[eachRes]]$redRef[[currIX]]$bit
            }
meanBit = apply(allBits,2,mean)

    nwInfo<-getNetworkInfoFuzzy(
        meanBit = meanBit,
        modelOriginal=ReferenceRes$paramsList$model,
        modelComprExpanded=ReferenceRes$processedModel,
        optimResT1=ReferenceRes,
        CNOlist=CNOlist,verbose=verbose)

    writeNetworkW(
        dN=nwInfo$dN,
        dM=nwInfo$dM,
        modelOriginal=ReferenceRes$paramsList$model,
        sifFile=nwInfo$sifFile,
        EAweights=nwInfo$EAweights,
        nodesAttr=nwInfo$nodesAttr,
        tag=tag)

    }


#########################################################################
##############these are the functions used above
#########################################################################
#this is the bit that produces the info that is needed above
getNetworkInfoFuzzy<-function(
    meanBit,
    modelOriginal,
    modelComprExpanded,
    optimResT1,
    CNOlist,
    verbose){

    #MKM: this function is basically equivalent to the one C. Terve wrote for the boolean case, but I took out the time information and, instead of just using the results of one Res, it takes the meanBit which is saved so it will affect the penwidth of resulting dot file.

    #!!! meanBit relates to strings in the compressed and expanded model
    reacs<-modelComprExpanded$reacID
    reacsOriginal<-modelOriginal$reacID

    #check if there are ANDs that were split, and if yes, set the score of the original AND
    #to the max score of its resulting branches (i.e.the max of the branches, if one branch was present
    #at t1 the AND gets 1, same if one was present at t2, if one of the branches was present at T1
    #and another one at t2, then the AND gets a score of 3)
    
    if(names(modelComprExpanded$SplitANDs) != "initialReac"){

        for(i in 1:length(modelComprExpanded$SplitANDs)){

            map<-match(modelComprExpanded$SplitANDs[[i]],reacs)
            mapScore<-meanBit[map]
            #this forces the score to be 3 if a branch was present at t1 and another at t2
            if(any(mapScore == 1) && any(mapScore == 3)) mapScore[1]<-3
            mapScore<-max(mapScore)
            SplitAnd<-names(modelComprExpanded$SplitANDs)[i]
            SplitAnd<-which(reacs == SplitAnd)
            meanBit[SplitAnd]<-mapScore

            }

        }

    #check if there were ands that were created, and if yes, assign to each branch of the original OR
    #the max of the score of the and or of that particular branch

    if(length(modelComprExpanded$newANDs) != 0){

        for(i in 1:length(modelComprExpanded$newANDs)){
            branches<-modelComprExpanded$newANDs[i]
            branchesPos<-match(unlist(branches),reacs)
            meanBit[unlist(branches)]<-pmax(
                meanBit[match(unlist(branches),reacs)],
                meanBit[match(names(modelComprExpanded$newANDs)[i],reacs)])
            }

        }

    #The new meanBit vector now contain weights that have been adapted to propagate weights coming
    #from new ANDs or split ANDs
    #Now create an adjacency list for the compressed/expanded model: col1=in node, col2=out node, col3=weight

    findOutput<-function(x){
        sp<-which(x == 1)
        sp<-modelComprExpanded$namesSpecies[sp]
        }

    findInput<-function(x){
        sp<-which(x == -1)
        sp<-modelComprExpanded$namesSpecies[sp]
        }

    n<-1
    adj<-matrix(NA,nrow=(length(reacs)*3),ncol=3)

    for(i in 1:length(reacs)){

        input<-findInput(modelComprExpanded$interMat[,i])
        output<-findOutput(modelComprExpanded$interMat[,i])

        weight<-meanBit[i]

        for(j in 1:length(input)){
            adj[n,1]<-input[j]
            adj[n,2]<-output
            adj[n,3]<-weight
            n<-n+1
            }

        }

    adj<-adj[-which(is.na(adj[,1])),]

    #remove duplicates that arise from new ANDs and split ANDs
    adj<-cbind(
        adj[order(adj[,3],decreasing=TRUE),1],
        adj[order(adj[,3],decreasing=TRUE),2],
        adj[order(adj[,3],decreasing=TRUE),3])
    adj<-unique(adj)

    for(i in 1:dim(unique(adj[,1:2]))[1]){

        findMatch<-apply(adj,1,function(x) all(x[1:2] == adj[i,1:2]))

        if(sum(findMatch) > 1){
            findMatch[which(findMatch == TRUE)[1]]<-FALSE
            adj<-adj[-which(findMatch==TRUE),]
            }

        }

        #Now do the same for the original model: col1=in node, col2=out node, col3=reaction #

        findOutput<-function(x){
            sp<-which(x == 1)
            sp<-modelOriginal$namesSpecies[sp]
            }

        findInput<-function(x){
            sp<-which(x == -1)
            sp<-modelOriginal$namesSpecies[sp]
            }

        n<-1
        adjOrig<-matrix(NA,nrow=(length(reacsOriginal)*3),ncol=3)

        for(i in 1:length(reacsOriginal)){

            input<-findInput(modelOriginal$interMat[,i])
            output<-findOutput(modelOriginal$interMat[,i])

        for(j in 1:length(input)){
            adjOrig[n,1]<-input[j]
            adjOrig[n,2]<-output
            adjOrig[n,3]<-i
            n<-n+1
            }

        }

        adjOrig<-adjOrig[-which(is.na(adjOrig[,1])),]

        #Go through the adjacency list of the compressed model and find which path is the shortest
        #between those species
        #this will hold the max weight of anything that has been mapped to this row of adjOrig so far
        OrigMap<-rep(0,dim(adjOrig)[1])

        for(i in 1:dim(adj)[1]){
            #If the row of the compressed adjacency list matches a row of the original adjacency list,
            #we just copy the weights
            rowAdj<-adj[i,1:2]
            matchR<-apply(adjOrig[,1:2],1,function(x) all(x == rowAdj))

            if(sum(matchR) > 0){

                OrigMap[which(matchR)]<-max(OrigMap[which(matchR)],adj[i,3])

            }else{

                #Otherwise we need to find the shortest path between the nodes in that row of adj, within the rows
                #of adjOrig, and map the weights of the row of adj to each of the row of adjOrig that are involved
                #in that path
                #This only works for paths of length 2 at most: improve this!!!!!
                ins<-which(adjOrig[,1] == adj[i,1])
                outN<-adjOrig[ins,2]
                target<-FALSE
                n<-1

                while(n <= length(outN) && !target){
                    outs<-which(adjOrig[,1] == outN[n])
                    outN2<-adjOrig[outs,2]

                    if(any(outN2 == adj[i,2])){

                        rOut<-which(outN2 == adj[i,2])
                        rOut<-outs[rOut]
                        rIn<-ins[n]
                        n<-n+1
                        target<-TRUE

                        }else{

                            n<-n+1
                            #this sets those variables to NA if the path couldn't be mapped
                            rIn<-NA
                            rOut<-NA

                        }

                    }
                    #if the path could be mapped, we record the mapping
                    if(!is.na(rIn) && !is.na(rOut)){

                        OrigMap[rIn]<-max(OrigMap[rIn],adj[i,3])
                        OrigMap[rOut]<-max(OrigMap[rOut],adj[i,3])
                    #if the path couldn't be mapped (path too long), we print a warning, nad OrigMap
                    #stays 0
                    }else{

                        if(verbose == TRUE){
                            print("Please be aware that when mapping the scaffold network back to the PKN, compressed paths of length > 2 are ignored.")
                        }

                    }

                }
            }
            meanBitOrig<-cbind(adjOrig[,3],OrigMap)

            for(i in 1:dim(meanBitOrig)[1]){
                meanBitOrig[i,2]<-max(meanBitOrig[i,2],meanBitOrig[which(meanBitOrig[,2] == meanBitOrig[i,1]),2])
            }

            meanBitOrig<-unique(meanBitOrig)
            meanBitOrig<-meanBitOrig[,2]

            #These mini functions are used to find the outputs and inputs of a reaction
    
            findOutput<-function(x){
                sp<-which(x == 1)
                sp<-modelOriginal$namesSpecies[sp]
            }

            reacOutput<-apply(modelOriginal$interMat,2,findOutput)

            findInput<-function(x){
                sp<-which(x == -1)
                sp<-modelOriginal$namesSpecies[sp]
            }

            reacInput<-apply(modelOriginal$interMat,2,findInput)

            #This mini function is used to create a reaction label as used in a cystoscape edge attribute file
    
            createReac<-function(x){
                r<-paste(x[1]," (",x[2],") ",x[3],sep="")
                return(r)
            }

            #if the class of reacInput is not a list, then there are no AND gates

            if(class(reacInput) != "list"){

                isNeg<-function(x){
                isNegI<-any(x == 1)
                return(isNegI)
                }

                inpSign<-apply(modelOriginal$notMat,2,isNeg)
                inpSign<-!inpSign
                inpSign[inpSign]<-1
                inpSign[!inpSign]<--1
                sifFile<-cbind(reacInput,inpSign,reacOutput)

                EApresent<-apply(sifFile,1,createReac)
                EAweights<-cbind(EApresent, meanBitOrig)

                }else{

                    #in this case there are AND gates and so we need to create dummy "and#" nodes
                    sifFile<-matrix(0,nrow=200,ncol=4)
                    nR<-1
                    nANDs<-1

                    for(i in 1:length(reacOutput)){

                        if(length(reacInput[[i]]) == 1){

                            sifFile[nR,1]<-reacInput[[i]]
                            sifFile[nR,3]<-reacOutput[i]
                            sifFile[nR,2]<-ifelse(any(modelOriginal$notMat[,i] == 1),-1,1)
                            sifFile[nR,4]<-meanBitOrig[i]
                            nR<-nR+1

                        }else{

                        for(inp in 1:length(reacInput[[i]])){

                            sifFile[nR,1]<-reacInput[[i]][inp]
                            sifFile[nR,3]<-paste("and",nANDs,sep="")
                            sifFile[nR,2]<-ifelse(modelOriginal$notMat[inp,i] == 1,-1,1)
                            sifFile[nR,4]<-meanBitOrig[i]
                            nR<-nR+1

                            }

                        sifFile[nR,1]<-paste("and",nANDs,sep="")
                        sifFile[nR,3]<-reacOutput[i]
                        sifFile[nR,2]<-1
                        sifFile[nR,4]<-meanBitOrig[i]
                        nANDs<-nANDs+1
                        nR<-nR+1

                        }
                }

            sifFile<-sifFile[1:(nR-1),]
            EApresent<-apply(sifFile[,1:3],1,createReac)
            EAweights<-cbind(EApresent,sifFile[,4])

            }
    #this bit creates a matrix of reaction that will be used for the dot file
    meanBitOrigRound = meanBitOrig
    meanBitOrigRound[meanBitOrig>0] = 1
    dM<-cbind(sifFile,meanBitOrigRound,meanBitOrig)

    #this mini function makes edge attributes in the format e1 (sign) e2 = attr

    makeEA<-function(x){
        ea<-paste(x[1],"=",x[2])
        return(ea)
        }

    EAweights<-apply(EAweights,1,makeEA)

    #write the nodes attributes
    nodesCompr<-modelComprExpanded$speciesCompressed
    indexes<-indexFinder(CNOlist,modelOriginal)
    nodesSig<-modelOriginal$namesSpecies[indexes$signals]
    nodesInh<-modelOriginal$namesSpecies[indexes$inhibited]
    nodesStim<-modelOriginal$namesSpecies[indexes$stimulated]
    nodesNCNO<-findNONC(modelOriginal,indexes)
    nodesAttrNames<-nodesSig
    nodesAttr<-rep("signal",length(nodesSig))

    if(length(nodesInh) != 0){
        nodesAttrNames<-c(nodesAttrNames,nodesInh)
        nodesAttr<-c(nodesAttr,rep("inhibited",length(nodesInh)))
        }

    if(length(nodesStim) != 0){
        nodesAttrNames<-c(nodesAttrNames,nodesStim)
        nodesAttr<-c(nodesAttr,rep("stimulated",length(nodesStim)))
        }

    if(length(nodesNCNO) != 0){
        nodesAttrNames<-c(nodesAttrNames,nodesNCNO)
        nodesAttr<-c(nodesAttr,rep("ncno",length(nodesNCNO)))
        }

    if(length(nodesCompr) != 0){
        nodesAttrNames<-c(nodesAttrNames,nodesCompr)
        nodesAttr<-c(nodesAttr,rep("compressed",length(nodesCompr)))
        }

    nodesAttr<-cbind(nodesAttrNames,nodesAttr)

#this bit creates a matrix of nodes attributes that will be used for the dot file
    dN<-nodesAttr

#this is the node attributes in the format as will be used by cytoscape
    nodesAttr<-apply(nodesAttr,1,makeEA)

#nodes that are neither inh/stim/sig nor NONC will not be in thi snode attribute so I'll add them
    other<-setdiff(modelOriginal$namesSpecies,dN[,1])
    other<-paste(other," = NA",sep="")
    nodesAttr<-c(nodesAttr, other)

#this is the edge attribute matrix that contains, for each edge, whether it is
#absent from the model (0), present at t1(1) or present at t2(2)
    return(list(
        dN=dN,
        dM=dM,
        sifFile=sifFile,
        nodesAttr=nodesAttr,
        EAweights=EAweights))
}


writeNetworkW<-function(dN,dM,modelOriginal,sifFile,EAweights,nodesAttr,
    tag=NULL){

    create_filename<-function(x, tag=NULL){
         if (is.null(tag)){
            return(x)
       }
       else{
           return(paste(c(tag, "_", x), collapse=""))
        }
    }


    writeDot(
        dotNodes=dN,
        dotMatrix=dM,
        model=modelOriginal,
        filename=create_filename("PKN.dot", tag=tag))
    write.table(
        sifFile[,1:3],
        file=create_filename("PKN.sif", tag=tag),
        row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t")
    write.table(
        EAweights,
        file=create_filename("TimesPKN.EA", tag=tag),
        row.names=FALSE,col.names="Times",quote=FALSE,sep="\t")

    write.table(
        nodesAttr,
        file=create_filename("nodesPKN.NA", tag=tag),
        row.names=FALSE,col.names="NodesType",quote=FALSE,sep="\t")

    }