with emerging high - throughput sequencing and phenotyping techniques [ 1 , 2 ] , gene functional module analysis has been extensively applied to explore the system - level functionality of gene groups in microorganisms , plants , animals and humans [ 3 - 5 ] . a gene functional module can be considered as a separated substructure of a biological network , i.e. , a group of genes that are related by one or more types of biological interactions such as gene co - expression / co - regulation , protein - protein interaction and functional association . in principle , a functional module is topologically and functionally separable from other modules , suggesting that genes in the same module often have tighter relations among themselves than with genes of other modules . these relationships could be better revealed when the network topology is visualized . zooming in from a whole network to functional modules allows biological researchers to focus on meso - scale gene groups with specific functions and interactions , thus generating testable hypotheses for unknown - function genes more effectively . an organelle is a specialized subcellular compartment that has a specific set of functions , and is usually enclosed within its own lipid bilayer . while some organelles are part of the endomembrane system and can receive and send materials through membranous vesicle trafficking , others rely on proteins and enzymes that are regulated at transcriptional , post - transcriptional and allosteric mechanisms . organelles need to communicate with each other to collaborate and perform complex functions that they can not do individually , and the complex interactions between organelles are often vital for the survival of organisms . however , the underlining mechanism of how organelles coordinate their functions in a cell is still largely unknown [ 13 , 14 ] . a systematic study of functional modules among multiple functionally - associated organelles in a cell may provide key information to understand how organelles coordinate their functions . the existing computational work to study organelle - organelle or organelle - nucleus interactions largely relies on gene co - expression analysis . for instance , to examine functional association between mitochondrial and nuclear proteins in humans , gene expression correlation tests have been performed to identify genes encoding nuclear proteins ( called nuclear genes in the following text ) whose expressions are significantly correlated with genes encoding organelle proteins ( called organelle genes in the following text ) . in such studies , nuclear genes positively - correlated with mitochondrial genes are often thought to be involved in transcriptional regulation , and negatively - correlated genes are more likely to be involved in translation of the mitochondrial genes . while the current organelle interaction models have revealed significant information , they basically ignore the fact that gene expression correlations may vary vastly among proteins that target to the same organelle . for example , under osmotic stress conditions , ribosome genes in arabidopsis thaliana , a model plant , are strongly co - expressed , but genes that encode endoplasmic reticulum ( er ) proteins are weakly co - expressed ( pearson correlation values being 0.89 and 0.64 respectively , ( see table 1 ) for the overview of the pearson correlation values of all organelles under seven abiotic stresses ) . many organelles perform a variety of biochemical functions , so it is not surprising that genes encoding these proteins display different levels of gene expression correlation coefficients . therefore , it makes more sense biologically to identify functional modules within each organelle and then connect these modules to build a network . functional modules need to be studied in context , as high - level biological functions often need the collaboration of multiple organelles to accomplish . for example , exploring functional modules in sensory organelles has successfully revealed the role of cilia - related trafficking in higher organisms . in this article , we introduce a new computational approach to construct networks of organelle functional modules in silico using transcriptomics data . although by definition the nucleus is also a type of organelle , for simplicity here we only refer to subcellular compartments other than the nucleus as organelles . analysis of arabidopsis gene expression data under seven abiotic stresses demonstrates that constructing a summary map of functional modules has significantly improved the interpretability of high - throughput transcriptomic data . methods to infer functional modules from multiple types of biological datasets have been extensively described [ 17 - 20 ] . in general , there are two categories of approaches to discover functional modules . in the first category , functional modules are detected from well - constructed biological networks , similar to detecting communities in a social network . heuristic methods have been proposed to identify tightly interacted nodes in a network by devising a network scoring function and then finding all high - score subnetworks [ 17 - 20 ] . while the aforementioned algorithms detect non - overlapped modules , they fail to detect highly overlapping functional modules . to this end , a soft clustering approach was presented to identify overlapping functional modules based on the idea of iteratively executing markov clustering . this soft clustering method was shown to outperform many approaches in terms of accuracy in functional module inference . in the second category , functional modules are directly detected from high - throughput biological datasets in a matrix format without using biological networks . for example , functional modules are expected to be identified directly from a matrix of gene expression data instead of a network . it is standard to use the pearson correlation coefficient ( pcc ) as a co - expression measurement , and to apply a threshold on pcc such that genes are connected if they have a significant pairwise expression profile association across various environmental perturbations [ 22 - 24 ] . however , this is a hard thresholding strategy that may cause loss of information and is sensitive to the choice of the hard threshold [ 2 , 23 ] . as an alternative approach , algorithm wgcna was proposed for soft thresholding , which weighs each gene pair and then uses average linkage hierarchical clustering to identify gene groups with highly correlated gene expression profiles across multiple conditions . some of the current functional module identification approaches use well - constructed gene networks as input , but for most organisms these networks are unavailable . some other approaches treat every gene equally and disregard the fact that gene correlation coefficients may vary dramatically in different organelles . therefore , traditional approaches are unfit for organelle functional module discoveries . furthermore , none of these algorithms study the relationships among functional modules , because identifying such type of relationships , especially to separate causative events from their effects , is often difficult [ 25 , 26 ] . in - network these studies include mining frequent subgraphs to extract communication patterns in data centers , inferring friendship network structure by using mobile phone data , and identifying spreading pattern of influenza epidemic . the discovered association patterns are tightly associated with node or edge topological properties , such as betweenness and node degree . applying these algorithms to the network of organelle functional modules may discover in - network associations , which can be explored further for a variety of purposes . in this article , we present a new platform to construct and analyze networks of organelle functional modules . our platform has the following advantages : it can identify functional modules in organelles and nucleus using gene expression data . it does not require a biological network as input , which could be difficult to build due to vastly different gene expression correlation levels in different organelles . our soft thresholding algorithm allows us to build a network of functional modules for nucleus and multiple functionally associated organelles . the connections between functional modules of the nucleus and other organelles may reveal regulatory or signal transduction events . through analysis of experimental data obtained from arabidopsis thaliana , we demonstrate the effectiveness of our method over the hard threshold approaches in interpreting arabidopsis gene expression datasets . we propose a new algorithm to discover organelle - to - organelle and organelle - to - nucleus functional associations . second , we use a self - adaptive thresholding approach to connect each functional module with the rest of the genes , resulting in a mixed network . different from the method proposed by amar and shamir 2014 , which connects two functional modules by checking whether their genes are directly connected in an extra dataset ( often being protein - protein interactions or genetic interactions ) , we examine functional module connectivity via the third party , i.e. , genes not in any functional modules . a clear advantage of our method is that it does not rely on extra datasets . using extra dataset may cause information loss and limit the use of the tool only to well - studied organisms , since pair - wise gene interaction datasets usually have much fewer number of genes than gene expression datasets . in order to identify organelle and nuclear functional modules , we first denote the subcellular localization information to every arabidopsis gene by selecting twelve organelle terms in the cellular component category of gene ontology   ( table s1 ) , and associating a gene with an organelle if the gene is contained in the annotation gene set of the corresponding organelle term ( after gene annotation propagation in go ) . we notice that the vast majority of organelle ( such as mitochondrial and plastid ) proteins are encoded in the nucleus , synthesized by cytosolic ribosomes and subsequently imported into the organelles via active protein transport systems . however , despite the profound differences in the organelle import machineries , a certain number of proteins are imported into multiple organelles . specifically , the proportion of the cross - organelle genes in the arabidopsis genome is constantly low in all the abiotic treatment data ( see details in table s2 ) . therefore , for genes that appear in multiple organelles , we make a copy of them for each corresponding organelle . we identify functional modules in each organelle using algorithm wgcna , where each module is a cluster of densely interconnected genes in the same organelle . wgcna firstly measures the similarities between gene expression profiles across all conditions using pcc and normalize them using jackknifed correlation coefficient . the normalized similarities are then transformed into an adjacency matrix . to avoid unnecessary loss of information , wgcna uses two soft adjacency functions instead of a hard threshold . next , an average linkage hierarchical clustering approach is used to group genes with coherent expression profiles . on the hierarchical clustering tree , a height cutoff is chosen to cut branches off the tree , resulting in gene functional modules , i.e. , sets of highly co - expressed genes . the choice of the height cut - off is guided by the topological overlap matrix plot . finally , a go enrichment test using gotermfinder   is applied to all the functional modules identified by wgcna . the background set of the enrichment test is all the significantly expressed genes in the input gene expression dataset . a functional module is considered to be functionally enriched if there exists at least one go biological process term with its adjusted p - value less than 0.05 . all the go enriched functional modules , we connect them to construct a network of organelle functional modules . a straightforward approach is to use the soft adjacency matrix built by wgcna to connect the functional modules . however , a clear drawback of the soft adjacency matrix is that it is not clear what the direct neighbors of a given functional module are . in particular , when neighbors of a functional module need to be explicitly listed , we have to binarize the values in the adjacency matrix . this is equivalent to the standard approach of hard thresholding the co - expression similarities , which are clearly not applicable to the identification of associations of organelle functional modules , since gene correlation coefficients vary greatly across all organelles ( as shown in table 1 ) . therefore , the wgcna soft adjacency matrix can not be directly used . to solve this problem the mixed network has two types of nodes : functional modules and individual genes that do not belong to any functional modules . the latter type of nodes is important because they may indicate the coordination pattern between functional modules . in a mixed network , specifically , we calculate pcc pair - wisely for all the nodes ( for functional modules , their eigengene expressions are used ) , and then compare this pcc with a soft threshold ( see discussion below ) . if a pcc value is greater than its soft threshold , the corresponding entry in the adjacency matrix madj is 1 , otherwise , it is 0 . to efficiently construct the mixed network , we set madj(i , j ) to 0 if the pcc value of node i and j is smaller than a user - provided lower - bound threshold . the lower - bound threshold is set to avoid high computational cost in network construction . it is based on the assumption that an organelle association network , similar to the other biological networks , is a sparse network . the key challenge here is to define a reasonable soft thresholding criterion to determine whether an edge in madj should be preserved or deleted . it has been found that genes that encode proteins functional in the same biological process or pathway are often well co - expressed [ 34 , 35 ] . for example , in a pathway including a functional module m and a gene g , the gene expression correlations between the genes in m is similar to that between g and at least one gene in m. in the context of a mixed network where pathways are not well defined , we assume the probability that both functional module m and gene g are in the same pathway is anti - correlated with the distance between m and g. we construct the mixed network in two steps . given a functional module m , we determine its core threshold ( ct ) using equation 1 : 
 ct(m)=max(sim(m),sim(m , neighbor(m ) ) ) 
  where sim(m ) is the averaged pcc value in functional module m when the genes in m are sparsely connected , meaning that the number of nodes and the number of edges are comparable ( in our experiment , we set the ratio of the two to be 1:1 ) . function neighbor(m ) returns all the direct neighbor genes of m after applying the user - provided lower - bound threshold . the direct neighbor genes of m are genes which pcc scores with m are greater than the lower - bound threshold . function sim(m , neighbor(m ) ) is the averaged pcc value between the eigengene expression profile of m and its neighbors . sim(m ) , the first part of equation 1 is the internal pcc of a functional module , and the second part ( sim(m , neighbor(m ) ) ) is added to avoid high false positives when the functional module internal pcc value is too small . next , we compute a soft threshold for each gene pair ( or gene - module pair ) using gene expression values and the topological character of the mixed network . mathematically , we define soft threshold st for functional module m and gene g as the sum of two components , shown in equation 2 : 
 st(m , g)=ic(m , g)ct(m)+(1ic(m , g))ht 
  where ht is a user - provided hard threshold and ct(m ) is the core threshold of m ( equation 1 ) . we assume that the chance for both functional module m and gene g to be in the same pathway decreases gradually as the distance between them increases . therefore , the influence coefficient function ic(m , g ) is defined as a transformed sigmoid function such that the influence of functional module m decreases from one to zero with the increasing distance from m. mathematically , ic(m , g ) is defined in equation 3 : 
 ic(m , g)=11+edist(m , g)5ct(m ) 
  where dist(m , g ) is the length of the shortest path between m and g ; and ct(m ) is the core threshold of m ( equation 1 ) . here for example , if gene g is 10 steps away from module m , then ic(m , g ) is close to 0 , i.e. st(m , g )  ht . an illustrative example of soft thresholding is shown in ( fig . 2 ) . the soft threshold st changes gradually as the distance between an organelle functional module and a gene increases . for the functional module represented by the green line , its ct is smaller than ht , so the st gradually increases from ct to ht as the distance from m to g increases . for another functional module ( blue line ) , its ct is greater than ht computationally , we adopt the breadth - first search algorithm   to construct the mixed network g(m ) for each functional module m ( see details in function generatemixednetwork in algorithm 1 ) . in the algorithm , we start with g(m ) that only has one node , i.e. , m , and then iteratively add genes and links to the boundary gene set of g(m ) ( see definition 1 ) . in the kth iteration , a mixed network g(m ) includes functional module m and the k gene . for any gene ( g ) not in g(m ) , we compute the edge correlation coefficient between g and all the genes in the boundary gene set of g(m ) . if at least one edge correlation coefficient value is greater than st(m , g ) , we add g and the corresponding links into the mixed network g(m ) . particularly , for edge e between genes g1 and g2 , which are in the boundary gene set , we add e to g(m ) if its edge correlation coefficient value is greater than st(m , g1 ) or st(m , g2 ) . in our experiment , we stop the iteration when the ic(m , g ) decrease to 0.05 . 
 definition 1 . the boundary genes are the genes whose distance to m is greater than any of its neighbors in g(m ) . 3 ) . given a functional module m , the initial mixed network g(m ) has only one node m and the initial boundary gene set is { m}. given two threshold ct(m ) = 0.954 and ht = 0.980 , we iteratively add nodes and edges to g(m ) . in the first iteration , we check all the nodes that are connectable to the boundary gene set and identify two edges < m , a > and < m , b>. since their weights are both larger than st , nodes a , b and edges < m , a > , < m , b > are added to g(m ) . in the next iteration , the boundary gene set is updated to to { a , b}. we check the edges between the boundary genes , i.e. a and b , and all the nodes that are connectable to the boundary gene set . second , we adopt a voting process to integrate all the single - module mixed networks ( see details in function votemixednetworks in algorithm 2 ) . the soft threshold list of each edge has k elements , which can be divided to two groups depending on whether a soft threshold is greater or smaller than the pcc value of the edge . if the summed difference between the pcc and the soft thresholds in the former group is larger than that in the latter group , we add the edge to the final mixed network . after establishing the final mixed network , we build a network of organelle / nuclear functional modules to explore how nuclear genes regulate the expression of organelle genes and how organelle genes are coordinately expressed . for simplicity , we require that in the network each organelle functional module connects to at most one nuclear functional module , while a nuclear functional module can connect to multiple organelle functional modules . subsequently , a nuclear functional module becomes a hub node surrounded by organelle functional modules . in addition to go based sematic similarity [ 37 - 39 ] between all corresponding gene pairs , we adopt three network topological measurements to measure organelle - nuclear functional module connection , i.e. , maximal flow , minimum edge cut , shortest path length   in the mixed network . for fair comparison , we rescale the values in each measurement to the same range ( between 0 and 100 ) and add up all of the four rescaled values . note that we reverse the shortest path length , because the shorter the path between two functional modules , the more likely they have coordinated functions . if a final score is above 75% percentile of all the values , we connect the corresponding organelle and nuclear functional modules . we next determine whether this method can be applied to transcriptomic data obtained from arabidopsis plants under various stress conditions . the responses of the model plant arabidopsis thaliana to abiotic stresses are accompanied by significant changes in transcriptome composition of genes that encode both nuclear and organellar proteins . in order to obtain a comprehensive view of how stresses may coordinately change the expression of organelle and nuclear genes , we analyzed the effects of salt , osmotic , cold , uv , drought , wounding , and heat stresses using publically available microarray data . we implemented the organelle functional module network construction method using python version 2.7.5 ( http://www.python.org/download/releases/2.7.5 ) and package networkx-1.9.1 ( http://networkx.github.io ) . the arabidopsis abiotic gene expression data set was downloaded from geo website under accession number gse 13584 . the gene ontology data set and the arabidopsis gene annotation data set were downloaded from the gene ontology website dated on july 2013 . in the abiotic gene expression data analysis , three biological replicates for abiotic stress lines under either salt , osmotic , cold , heat , uv , wounding , or drought were compared against three biological replicates of mock control lines at the same time points . statistically significant differences in gene expression between treatment and control plants were detected at a fold change of 2 and fdr of 0.01 using the limma ( bioconductor ) package . we denoted the subcellular localization information in gene ontology to every significantly expressed gene under each stress condition , and discarded the genes that lack location information . we also discarded organelles such as glyoxysome and lysosome , since glyoxysome is a type of peroxisome and has only one significantly expressed gene in our datasets , and plants do not have lysosomes . we identified functional modules in each organelle using algorithm wgcna , where each module is a cluster of densely connected genes in the same organelle . the parameters of wgcna that we used are as follows : the parameters in function adjacency : softpower = 6 , type =  signed hybrid  ; hierarchical clustering function : flashclust , method =  average  ; we used the parameters are pamstage = false , deepsplit = 4 , and minclustersize = 5 . for each functional module the go enrichment test shows that 50.0% of the functional modules identified by wgcna are go enriched ( see summary in table 2 and detailed gene lists in supplementary table s3 ) . the lower bound threshold of pcc and ht were set to be 0.95 and 0.98 respectively for all stress conditions . ( table 3 ) shows the number of nodes and edges of each mixed network , and the number of functional modules . the number of significantly expressed genes is the highest under osmotic and the lowest under drought conditions . therefore , the osmotic mixed network is the largest among all the seven networks , whereas the drought network is the smallest . the network of uv has the highest number of nuclear and organellar functional modules ( in total 105 functional modules ) . among all the mixed networks , network topological analysis shows that it is a scale - free network with an r value of 0.84 ( see degree distribution in fig . 5 ) , and three of them are closely connected to a mitochondrial and a chloroplast functional module . gene ontology analysis reveals that the genes in the closely connected functional modules are enriched in signal transduction and immune response ( adjusted p - value<0.05 ) . 6 ) reveals that at least half of the existing edges will be deleted if a hard threshold approach is applied ( ht=0.98 ) , making all the weakly co - expressed functional modules into orphan nodes . with our soft thresholding approach , even using the same threshold ( ht=0.98 ) , we could adaptively connect the weakly co - expressed organelle functional modules to the rest of the network . surprisingly , unlike the network of drought , the mixed networks of cold , osmotic , uv , salt , wounding and heat conditions each have two or three large and similar - sized subnetworks ( see suppl fig . we generated multiple networks under drought treatment using different lt and ht thresholds , in order to test the sensitivity of our algorithm . first , we fixed the ht threshold as 0.98 and varied lt from 0.8 to 0.95 by 0.05 ( see table s4 ) . second , we fixed the lt threshold as 0.95 and varied ht from 0.96 to 0.99 by 0.01 ( see table s5 ) . third , we fixed lt as 0.8 and varied ht from 0.85 to 0.98 ( see table s6 ) . the results show that when one of the thresholds is selected appropriately , our algorithm is not sensitive to the change of the other one ( see table s4 and 5 ) . however , if both thresholds are improperly set , the network can change a lot ( see table s6 ) . we compared our mixed networks with the networks generated by using a hard threshold . to ensure fair comparison , the hard - threshold networks have the same number of genes as our networks . first , we compared the number of chloroplast genes in the hard - threshold network vs. our network . the result shows that our network can remains more organelle genes than the hard - threshold network under all the abiotic stress conditions ( see table s7 ) . second , we counted the number of genes involved in the go enriched functional modules in our networks and compared it with that in the hard - threshold networks . the result shows that our algorithm can constantly identify more go enriched genes than the hard - threshold networks ( see fig . s8 ) . on the mixed network of each stress condition , we measured the connections between the organelle and the nuclear functional modules pair - wisely using the method described in subsection 3.3 . for each stress condition , we generated a hub - like functional module network , in which the center nodes are nuclear functional modules , the surrounding nodes are organelle functional modules , and the edges represent strong functional associations between them . the network of functional modules may lead to new hypotheses on how organelles are co - regulated or signal transduction events between nucleus and organelles . our data analysis results support some of the previous reports on nuclear regulation in stress conditions . for example , the network of functional modules in drought consists of one nuclear functional module and three organelle functional modules ( chloroplast , mitochondria and golgi ) , forming a star graph ( fig . the centered nuclear functional module contains 11 transcription factors ( see the full gene lists in table 4 ) , among them bzip60 ( at1g42990 ) is a transcription factor activated under conditions that induce unfolded protein response ( upr ) , a signaling pathway that up - regulates the expression of er - chaperones . after its activation and translocation into the nucleus , bzip60 regulates the expression of genes that encode components of the upr [ 46 , 47 ] . our model predicts that some of the genes in the golgi , mitochondrial and chloroplast functional modules may be targets of bzip60 and/or other proteins in the nuclear module . the network of functional modules under osmotic stress consists of three major star graphs that each contains more than ten functional modules ( fig . 8a ) are enriched in primary metabolic process and organic substance biosynthetic process , whereas genes in ( fig . 8b ) are enriched in vesicle - mediated transport , golgi vesicle transport , intracellular transport , cellular catabolic process and vacuolar transport . 8c ) are enriched in fatty acid beta - oxidation , lipid oxidation , lipid modification and fatty acid catabolic process . ath1 encodes a transcription factor in chlorophyll biosynthetic and chlorophyll metabolic processes [ 48 , 49 ] , suggesting that some of the genes in the chloroplast modules may be targets for ath1 . there are only one or two major star graphs under the other stress conditions ( see suppl fig . s7 - 12 ) . under all the stress conditions a nuclear functional module always connects to multiple types of organelles , suggesting that these different organelles may need to coordinate their function in various stress conditions
with the rapid accumulation of gene expression data , gene functional module identification has become a widely used approach in functional analysis . 
 however , tools to identify organelle functional modules and analyze their relationships are still missing . 
 we present a soft thresholding approach to construct networks of functional modules using gene expression datasets , in which nodes are strongly co - expressed genes that encode proteins residing in the same subcellular localization , and links represent strong inter - module connections . 
 our algorithm has three steps . 
 first , we identify functional modules by analyzing gene expression data . 
 next , we use a self - adaptive approach to construct a mixed network of functional modules and genes . finally , we link functional modules that are tightly connected in the mixed network . 
 analysis of experimental data from arabidopsis demonstrates that our approach is effective in improving the interpretability of high - throughput transcriptomic data and inferring function of unknown genes .