text:
NTAP: for NimbleGen tiling array ChIP-chip data analysis

Summary:NTAP is designed to analyze ChIP-chip data generated by the NimbleGen tiling array platform and to accomplish various pattern recognition tasks that are useful especially for epigenetic studies. The modular design of NTAP makes the data processing highly customizable. Users can either use NTAP to perform the full process of NimbleGen tiling array data analysis, or choose post-processing modules in NTAP to analyze pre-processed epigenetic data generated by other platforms. The output of NTAP can be saved in standard GFF format files and visualized in GBrowse. Availability and Implementation:The source code of NTAP is freely available at http://ntap.cbi.pku.edu.cn/. It is implemented in Perl and R and can be used on Linux, Mac and Windows platforms. Contact: ntap@mail.cbi.pku.edu.cn; luojc@pku.edu.cn; hekun78@gmail.com

INTRODUCTION
Genome-level high-density tiling arrays are becoming more accessible for genome-wide profiling studies including transcriptome identification (Bertone et al., 2004), transcription factor binding site identification (Lee et al., 2007), histone modification profiling (Gendrel et al., 2005;Li et al., 2008), DNA methylation profiling (Hayashi et al., 2007) and comparative genome hybridization. Specific analysis methods and tools are required for each type of study because the strategies behind different tiling array applications vary extensively. As a result, several models have been proposed and software tools have been developed for the analysis of different types of tiling array data (Chung et al., 2007;Ji et al., 2008;Li et al., 2005;Wang et al., 2006;Zhang et al., 2007). However, there is still room to improve for data analysis of epigenetic features including histone modifications and DNA methylation. The recognition of distribution patterns of modifications at both the local (gene) and global (chromosome) levels are usually required to infer biologically meaningful conclusions (Hayashi et al., 2007;Li et al., 2008). * To whom correspondence should be addressed.
Here, we present a NimbleGen Tiling array data Analysis Package (NTAP) designed for histone modification profiling analysis (Li et al., 2008) that can also be applied to other ChIP-chip data (Lee et al., 2007). The advantage of our package is its ability to generate reports for various pattern recognition questions instead of focusing only on identifying significantly enriched oligos or genomic regions.

FUNCTIONS AND FEATURES
NTAP was developed using the R statistical language to take advantage of the powerful statistical functions of other open source packages especially those from the Bioconductor project (http://www.bioconductor.org/). It contains five main steps for data analysis: importing, normalization, feature identification, oligos mapping and post-processing for pattern recognition.

Data importing
We implemented an R function similar to the 'read.maimages' function in the limma package (Smyth 2004) to import NimbleGen raw data into limma data object formats for normalization.

Data normalization
Users can apply various microarray normalization methods to the imported datasets through the limma package functions 'normalizeBetweenArrays' and 'normalizeWithinArrays'. Unlike the expression profiling arrays whose log transformed ratio distributions are usually symmetric around zero, the distribution of the ChIP-chip result tends to skew to the ChIP channel. Because only the protein-bound DNA fragments will be pulled down by a specific antibody, more positive log transformed ChIP/Input ratios are expected. Thus, the rank-invariant set scheme (Buck and Lieb, 2004) was incorporated for better data normalization.

Feature identification
Tiling arrays usually contain several oligos per single gene rather than one oligo per gene. For example, the traditional whole-genome array for expression studies in Arabidopsis thaliana usually contains only 23k oligos, while a customized whole-genome tiling array © 2009 The Author(s) This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/ by-nc/2.0/uk/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. tiled at ∼250 bp resolution may contain ∼400k oligos. The much larger number of oligos on a single array makes the traditional methods for feature identification unfeasible. For tiling array data, expressed mRNA or pulled-down DNA fragments can cause the signal of a group of neighboring oligos to increase simultaneously. Therefore, our package implements the non-parametric Wilcoxon rank-sum method to compare the signal differences between the ChIP channel and the reference channel for a group of oligos using sliding windows. Under certain circumstances, however, the density of some tiling arrays may not be high enough to use the Wilcoxon method. In these cases, we utilize simple comparison linear models implemented in limma (Smyth 2004) to identify single oligos whose signal increased significantly in the ChIP channel. Then, we consider a genomic region as 'positive' if the region contains a single oligo that meets stringent user-defined criteria or the region contains a group of neighboring oligos that meet less stringent criteria.

Mapping oligos to gene models
Genome data are usually kept up-to-date by genome sequencing consortia or curation groups, who usually release their data as standard XML format files that can be parsed to easily obtain coordinates of gene models. A Perl module was implemented to retrieve records of the gene model position information on each chromosome and to determine the relative position of a specific oligo to its nearby gene model(s). Signal distribution patterns among different groups of genes can then be determined based on the stored relative position information.

Post-processing functions
The following questions are frequently asked in epigenomics research. What is the modification distribution pattern relative to genes and does it vary between different organs/tissues? Is there an association between specific histone modification levels and gene sizes, or gene expression levels? To answer these questions, we implemented several R functions to align genes, to calculate the average ChIP/Input intensity ratio of the oligos within sliding windows, and to plot the final results for different groups (Fig. 1).

Result visualization
Quality control is a key step to guarantee the validity of the overall data analysis. An R function was implemented to calculate the raw intensities correlation coefficient between any pair of two replicates. MA-plots of array hybridization results are also generated in order to examine the intensity ratio (M) versus averaged intensity (A) to discover possible non-linear biases that require special normalization methods. After raw data processing, all the oligos are mapped back to the most up-to-date chromosomes and the ChIP/Input ratio value of each oligo can then be plotted along the chromosome. These values can be displayed either by a program within NTAP or they can be exported in the GFF format to be displayed in the Generic Genome Browser GBrowse (Stein et al., 2002).

IMPLEMENTATION
Most of the functions are implemented in the R statistical language (http://www.r-project.org/) and Perl. Users can also choose any other software to pre-process their data before using our post-processing modules.


extfieldsofstudy:
Computer Science
Medicine

text:
From disease ontology to disease-ontology lite: statistical methods to adapt a general-purpose ontology for the test of gene-ontology associations

Subjective methods have been reported to adapt a general-purpose ontology for a specific application. For example, Gene Ontology (GO) Slim was created from GO to generate a highly aggregated report of the human-genome annotation. We propose statistical methods to adapt the general purpose, OBO Foundry Disease Ontology (DO) for the identification of gene-disease associations. Thus, we need a simplified definition of disease categories derived from implicated genes. On the basis of the assumption that the DO terms having similar associated genes are closely related, we group the DO terms based on the similarity of gene-to-DO mapping profiles. Two types of binary distance metrics are defined to measure the overall and subset similarity between DO terms. A compactness-scalable fuzzy clustering method is then applied to group similar DO terms. To reduce false clustering, the semantic similarities between DO terms are also used to constrain clustering results. As such, the DO terms are aggregated and the redundant DO terms are largely removed. Using these methods, we constructed a simplified vocabulary list from the DO called Disease Ontology Lite (DOLite). We demonstrated that DOLite results in more interpretable results than DO for gene-disease association tests. The resultant DOLite has been used in the Functional Disease Ontology (FunDO) Web application at http://www.projects.bioinformatics.northwestern.edu/fundo. Contact: s-lin2@northwestern.edu

INTRODUCTION
A general-purpose open biomedical ontology can be thought of as a knowledge capture or knowledge representation device. According to the Open Biomedical Ontologies (OBO) Foundry principles, the ontology should contain well-defined terms with well-defined relationships (e.g. is_a, part_of) between the terms, and represent as much of the current knowledge in a given domain as possible. Disease Ontology (DO) is an OBO Foundry ontology, organized such that the path to the root is always true. DO organizes disease concepts in a directed acyclic graph (DAG) so that traversing away from the root of DO moves towards progressively more granular and specific terms. The full DO graph is very useful for organizing a wide spectrum of data, but is not necessarily optimal for specific applications such as identifying disease-gene relationships. By analogy with other biological ontologies, a trimmed-down version is very useful for building a high-level functional summary from a gene list. For instance, a simplified version of the Gene Ontology (GO) called 'GO Slim' provides a broad, integrative overview of molecular and cellular biology by combining and removing finegrained terms in the GO. GO Slim has proven critical in comparing annotations across genomes (Adams et al., 2000) and interpreting biological functions of a gene list (Shah and Fedoroff, 2004).
We are especially interested in interpreting a list of genes in the context of disease, which is a critical step in translating molecular findings from microarrays, proteomics and other types of high-throughput screening methods into clinical relevance. To achieve this goal, we were part of the collaboration that has developed the general-purpose OBO Foundry DO (http://www.diseaseontology.sourceforge.net). The DO has been successfully used here at Northwestern in the NUgene project to build a detailed disease phenotype from data available in the electronic medical record. In addition, the DO has been used as a controlled vocabulary to annotate the human genome in terms of diseases (Osborne et al., 2009).
However, the DO is very complex: revision 26 of the DO contains 11 961 terms in the form of a directed acyclic graph, of which 4399 terms are internal nodes with up to 16 levels of hierarchical structures. Such a complex structure creates a special challenge for functional enrichment tests.
Functional enrichment tests were originally developed to interpret biological meanings of a gene list from microarray experiments (Dopazo, 2006). Such statistical tests are also called GO analysis, because the GO was originally used to define the functional categories. Briefly, k functional categories (diseases) are defined. For each category, a hypergeometric test is used to test the enrichment of proportions (per cent of genes belong to this category) in the identified gene list versus the genome (Falcon and Gentleman, 2007). Although more sophisticated statistical procedures are in development (Goeman and Mansmann, 2008), canonical tests assume that each category is independent and identically distributed (iid, i.e. not in a hierarchical structure). A simplified approach is to treat each term in the ontology graph as an independent category (ignoring the graph structure), which not only violates the iid assumption but also creates another problem of multiple comparisons (k = 11 961 for DO). Thus, trimming the DO into © 2009 The Author(s) This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/ by-nc/2.0/uk/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. simplified categories will make the statistical test tractable as well as the results simpler for human comprehension.
A simple method of trimming DO terms is to select the DO terms at a certain level in the DAG. This method is very arbitrary and it is hard to determine a fixed level suitable for different branches in the DAG; moreover, many important DO terms might be missed. To deal with the similar problem for GO, Alterovitz et al. (2007) proposed to partition the GO database based on information content (Shannon information) of individual GO terms. But how to choose and interpret the level of information content is still a problem. Here, we prefer grouping the DO terms using semantic distance in the context of gene annotations. The semantic distances in previous studies are either information based (Jiang and Conrath, 1997;Resnik, 1999), graph-structure based (Wang et al., 2007;Wu et al., 2005) or a hybrid of both (Sheehan et al., 2008). One limitation of graph-structure-based semantic distances is that it does not consider genes that are annotated by DO terms. As the genes annotated by DO terms are direct evidence of disease classification at the molecular level, DO terms with the same or similar mapped genes are closely related and can be combined. Moreover, the combination of DO terms based on similar gene mappings will remove the redundancies from the results of an enrichment test. Therefore, we define the distance metrics of DO terms based on gene-to-DO mapping profiles. The mapping profile-based distance metric has been used to measure the similarity between two genes, but, to our knowledge, it has not been used to measure the ontological similarity. As we are interested in both the overall similarity and the subset similarity (one DO term-associated gene list is the subset of another one) between two DO terms, we separately define two types of binary distance metric. A compactness-scalable fuzzy clustering method (Du et al., 2005) is then applied to group similar DO terms based on defined distance metrics. For the clustering using subset similarity distance metric, the clustering results are further constrained by the semantic similarities between DO terms to reduce the false clustering. Following these steps, a domain expert was assigned to curate the computational results and assign new names to these DO term groups. We named these simplified controlled vocabulary lists of diseases as 'Disease-Ontology Lite' (DOLite). The details of each step are described in the methodology section. Figure 1 shows the framework of creating the DOLite database based on the DO database. Next, we describe each step in detail.

Pre-filtering the DO terms
In the DO database, many of the DO terms are abstract concepts created for the purpose of the ontological reasoning; they are not used in scholarly communications and have very few genes directly associated with them. For instance, 'cell proliferation disease' only appears three times in 18 million MEDLINE abstracts, but it is still listed as a disease category in the DO graph ( Fig. 2). A quick filtering based on the number of genes directly mapped to each DO term can easily identify these abstract DO terms.
Some other DO terms, which are not very well studied, also have very few genes associated. Because of their small sizes of gene lists, they are easy to be significant by chance in the enrichment test. To reduce the false positives and increase the computational efficiency, we pre-filter these DO  terms and map their associated genes to their direct parent DO terms before further calculations.

Gene-to-DO mapping profile matrix
Each gene in the human genome was annotated with DO, as we reported before (Osborne et al., 2009). Figure 3 shows the mapping profile matrix M. Each element M lk represents the mapping between gene G l and ontology term O k . The values of the matrix elements are as follows: for M lk = 1, there exists evidence showing that gene G l is involved in the disease defined by ontology term O k and for M lk = 0, there is a lack of evidence. Lack of evidence does not imply that there is no relationship but rather that there is no information available on a relationship. For our application, each column is a gene-to-DO mapping profile of the corresponding DO term; it is also equivalent to a DO-term-associated gene list.

i64
From disease ontology to disease-ontology lite

Distance metrics between gene-to-DO mapping profiles
To measure the similarity between DO terms, we first need to define the distance metrics. We chose to use the gene-to-DO mapping matrix M based on distance because we are interested in identifying associations between genes and DO terms. Kappa statistics were used to measure the similarities between genes based on their mapping profiles (rows of mapping profile matrix M in Fig. 3) (Dennis et al., 2003). As kappa statistics put equal importance on 0s and 1s, this is very problematic when applied to cluster DO terms. As shown in Figure 3, M lk = 0 is a lack of evidence of involvement, not information that rules out an involvement between gene G l and ontology term O k . Therefore, our distance metrics will focus on 1s. Considering the special situation of the structured ontology, we defined two types of binary distance metrics, as shown in Equations (1) and (2). The distance values are between 0 and 1, with values closer to 0 indicating closer similarities. In Equations (1) and (2), '∧' is the binary AND operator and '∨' is the binary OR operator. dist 1 (A,B) measures the overall similarity between two binary vectors A and B. In our case, A and B represent the columns of mapping profile matrix M, and the distance is 1 minus the percentage of overlapping genes among all genes associated with DO terms A and B. dist 2 (A,B) measures the subset similarity between two binary vectors A and B. In our case, it is 1 minus the percentage of one DO term-associated gene list being the subset of another DO term-associated gene list. This is of interest particularly when these two ontology terms have a close semantic relationship; e.g. for parent-child relationships in DO DAG graph.

Clustering the gene-to-DO mapping profiles
With respect to the molecular basis of disease, significant redundancy and overlap also exist in DO. For example, six terms ('carcinoma of the large intestine', 'adenocarcinoma of large intestine', 'adenocarcinoma of the colon', 'colorectal neoplasms', 'malignant tumour of colon' and 'colonic neoplasms') in a small portion of the DO graph (as shown in Fig. 2) basically describe the same group of diseases centred on the adenocarcinoma of colon. Some of these terms may include other rare cancers affecting the colon, and others may also include benign colonic tumours. However, these nuances might not be differentiable by their associated genes. As such, the enrichment test will likely report all these terms redundantly, making it more difficult to interpret. In order to find the DO terms with very similar associated genes, we perform clustering based on gene-to-DO mapping matrix M. Two types of clustering based on two different binary distance metrics, as shown in Equations (1) and (2), are performed. The first clustering is based on the binary distance metric, dist 1 , which measures the overall similarity of two ontology terms (i.e. their associated gene lists). As one gene can involve in multiple diseases, semantically distant DO terms may have significant biological pathway overlap and we do not want to prune semantically distant terms as these may be biologically informative. Unlike the conventional clear-cut clustering algorithms, fuzzy clustering methods allow the overlapping assignment of the DO terms to different clusters with varying degrees of membership. Therefore, fuzzy clustering is ideally suited for our situation. Another requirement for the clustering algorithm is the ability to control of cluster compactness, which means the level of similarity among cluster members can be modified by the researcher. On the basis of these considerations, we use the compactnessscalable fuzzy K-means clustering algorithm described in Du et al. (2005). This algorithm adds a scalable Gaussian window to the fuzzy membership function to control the compactness of cluster. The σ Gaussian window function controls the scale of the clustering algorithm, with smaller scales having more compact clusters. For more details, please refer to Du et al. (2005). Just as K-means clustering algorithms, the clustering result of compactness-scalable fuzzy K-means algorithm depends on the initial choice of the cluster centres. To better estimate the initial cluster centres, we first perform a hierarchical clustering using a complete linkage method. By limiting the height of the resultant hierarchical tree, we use the complete linkage method to seed the initial cluster centres for the compactness-scalable fuzzy K-means clustering algorithm.
The second clustering is based on the distance metric, dist 2 , which measures the subset similarity between DO terms. For the case of one DO term is a subset (or partial subset) of another, we are only interested in the case with the terms having high semantic similarity. Therefore, the clustering result will be constrained by the DO semantic distance as described in the Section 2.6). Because of the constraint on DO semantic similarity, the results of fuzzy and hierarchical clustering algorithms are very similar. Since hierarchical clustering algorithms are much more efficient, the second cluster was built by hierarchical clustering.

Semantic distances between DO terms
The estimation of semantic distance can be categorized as informationbased, graph-based and hybrid methods. We select to use a graphbased Union-Intersection (UI) method, as implemented in the GOstats Bioconductor package (Falcon and Gentleman, 2007), for its simplicity, ease of interpretation and consistency with other distance metric defined in Equations (1) and (2). The semantic distance of UI method is 1 minus the percentage of overlapping nodes between two induced graphs among the total number of nodes in two graphs. The resulting distance is in the range of 0-1 with values close to 0 having better similarities. The graph corresponding to each DO term is composed of multiple paths from the DO term to the root of the DO terms (disease). An implication of this method is that two deep sister nodes (i.e. nodes that are distant from the root) will have higher levels of similarity than two shallow sister nodes. Given the design of DO, this is consistent and intuitive with the semantic similarity, or lack thereof, for shallow nodes. This also means that two neighbouring, deep nodes, one of which has two or more distinct paths to the root will be quite distinct from nodes that only have one path to the top, or have differing paths to the top. This again leverages the biological knowledge embedded in DO.

Combine clustering results
We treat the clustering result based on dist 1 as the major clusters. The clustering result based on dist 2 is sub-clustered based on the semantic distance between DO terms. Hierarchical clustering using the single linkage method was used for sub-clustering. Sub-clusters with single elements are filtered from further processing. Other sub-clusters are then merged with major clusters whenever there are overlapping elements.

Expert curation of DOLite
After the statistical treatments as described above, a board-certified pathologist manually edited the clustered DO terms, resolved the computational artefacts and assigned a class label for each DOLite term. DOLite is a simplified version of the DO for functional enrichment tests of genes ( Table 2). Instead of using medical terminologies for the label and i65 definitions, we have attempted to use common English names as used in Wikipedia for DOLite.

RESULTS
Following the procedures shown in Figure 1, we created DOLite based on the DO. The DO terms were filtered to include only those with at least five direct mappings. A gene-to-DO mapping profile matrix was created. In this mapping matrix, each DO term also include the mappings of its offspring. Clustering was performed based on this mapping profile matrix. Because all the distance metrics we used are in the range of 0-1 and are related with the percentage of overlapping, we used the same threshold value for different steps, which include the cut-off thresholds used in clustering and the maximum allowed semantic distances of DO terms. The scale of fuzzy clustering is selected as the half of this threshold value. In our implementation, we set this distance threshold value as 0.2 (80% of similarity). As a reference to the expert curator, we also provide the clustering results based on looser threshold values 0.3 and 0.5. The expert curator decided the final mapping between DO terms and DOLite terms. Although the general-purpose DO is designed to differentiate the nuances of different disease categories, for the application of functional enrichment test of gene lists, we decide to aggregate the related terms and assign common English names for the diseases (as used in Wikipedia). An example is shown in Table 1. Table 2 shows the overview of final DOLite database and its comparison with DO database.

DOLite annotation of the human genome
With the DO-to-DOLite mapping, we convert the previous annotation of the human genes by DO into DOLite annotations. Compared with DO, the DOLite annotation of the human genome is more compact, as evidenced by the larger number of genes assigned to individual disease category (Fig. 4B), and the smaller number of disease categories assigned to individual gene (Fig. 4A).

Validation of DOLite using a benchmark microarray data set
To validate the utility of DOLite, we used a benchmark microarray data set of pancreatic cancer study (Antonov et al., 2008), which was previously utilized to test GO-based annotations (Falcon and Gentleman, 2007), for the functional enrichment test. Briefly, a list of 125 genes identified in that study was used for functional analysis. Hypergeometric test was applied to each DO (Osborne et al., 2009) and DOLite term in the database. We show side-by-side the top 12 categories (ranked by Pvalues) of the DO and DOLite results in Table 3 and Figure 5. As expected, the conventional DO analysis returned highly redundant terms such as 'cancer', 'carcinoma' and 'adenocarcinoma', likely because a number of these genes are shared in the biological process of cancer development (Fig. 5A). However, no pancreasspecific disease entities were picked up. Instead, it returned some diseases unrelated to pancreatic adenocarcinoma, such as respiratory tract disease, soft tissue neoplasm and skin disease. These results may represent non-specific hits that contribute little to data interpretation. In contrast, the analysis based on DOLite did successfully return a number of entities closely related to  the pathology of pancreas, in addition to generic cancer-related terms. These entities include diabetes mellitus, hypoglycemia and obesity and disorders frequently associated with the malfunction of pancreatic islet cells. Another entity, primary biliary cirrhosis, also causes pancreatic injury. These results support a strong pancreas conjunction of the tested gene set, validating the effectiveness of using DOLite to interpret microarray data.
Results from DOLite analysis also lead to interesting discoveries that deserve further clinical investigation. We observed that breast cancer was ranked very high among the DOLite results. This could i66  be simply because breast cancer is one of the most extensively studied cancers and therefore has more associated genes. However, the predominant type of breast cancer is adenocarcinoma, which can resemble some pancreatic adenocarcinomas on histology. So, could these two cancers share some common genetic abnormalities as well, other than those abnormalities seen in most cancers (e.g. genes involved in cell cycle control)? This may be an interesting question to explore.

DISCUSSION
Exploring relations between genes and diseases at the molecular level could greatly facilitate our understanding of pathogenesis, and eventually lead to better diagnosis and treatment. DO was constructed as a general-purpose ontology to define diseases. It is very useful to associate genes with diseases (e.g. a particular cancer) or disease-related processes (e.g. tumour metastasis). However, due to the complex structure of DO, the results of functional enrichment test are usually hard to interpret. In this study, we proposed statistical methods to slim the general-purpose disease ontology to a simplified controlled vocabulary list called 'Disease-Ontology Lite' (DOLite) specifically for the functional enrichment test of a gene list. DOLite has several advantages over the conventional DO for functional enrichment test. The complex hierarchical structure of DO is unnecessary for the desired enrichment test. DOLite contains betterdefined disease entities that are easier to interpret by researchers and clinicians. The entities are also enriched with associated genes, leading to increased sensitivity of detection. The enrichment test using the benchmark microarray data set of pancreatic adenocarcinoma suggests that DOLite returned diseases specific to the organ where original experiments were conducted on, which directly confirm the validity of the experimental results.
In terms of contributions to biological databases, we created DOLite, the first simplified version of the DO, and utilized DOLite to annotate the human genome. In terms of methodological contribution, we defined statistical methods to simplify a generalpurpose ontology. In contrast, the previous construction of GO Slim from GO was largely a subjective process based on expert opinion. The major contributions in methodology include computing the ontology similarity based on gene-toontology mapping profiles [derived from geneRIF (Harris et al., 2004;Shah and Fedoroff, 2004)]; defining two types of binary distance metrics to separately measure the overall similarities and subset similarities and a compactness-scalable fuzzy clustering method (clustering results were verified with the constraints i67 of semantic distance between DO terms). The methodology can be easily adapted to slim other ontologies, like GO. Currently, DOLite is utilized in the Functional Disease Ontology (http://www.projects.bioinformatics.northwestern.edu/fundo) Web application to interpret clinical relevancies of a gene list.
Previous studies have suggested that GO Slim, which is a simplified version of GO, greatly facilitates genome-wide computational analysis (Nam and Kim, 2008). Similarly, we expect that DOLite and its annotation of the human genome can be used as a foundation for the development of more computational and statistical methods to analyse the disease relevance of a gene list. In the future, we will further refine the algorithms for creating the DOLite database. For example, we can better integrate the semantic distance of DO terms into the clustering process. We also plan to add more sensitive methods to test the enrichment of each DOLite term in a gene list.
Funding: This project was supported in part by Award Number UL1RR025741 from the National Center for Research Resources. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Center for Research Resources or the National Institutes of Health.


extfieldsofstudy:
Mathematics
Computer Science

text:
Finding total unimodularity in optimization problems solved by linear programs

A popular approach in combinatorial optimization is to model problems as integer linear programs. Ideally, the relaxed linear program would have only integer solutions, which happens for instance when the constraint matrix is totally unimodular. Still, sometimes it is possible to build an integer solution with the same cost from the fractional solution. Examples are two scheduling problems and the single disk prefetching/caching problem. We show that problems such as the three previously mentioned can be separated into two subproblems: (1) finding an optimal feasible set of slots, and (2) assigning the jobs or pages to the slots. It is straigthforward to show that the latter can be solved greedily. We are able to solve the former with a totally unimodular linear program, from which we obtain simple combinatorial algorithms with improved worst case running time.

Introduction
In this work, we propose a specific approach to give simpler solutions to several optimization problems. Herein we considered three such optimization problems: the first two are scheduling problems : the Tall-Small Jobs Problem [4] and the Equal Length Jobs Problem [6]. The last one is about Offline prefetching and caching to minimize stall time [3]. In the Tall-Small Jobs Problem, we have m machines, n unit length jobs, some of which need to execute on all the machines at the same time. In the Equal Length Jobs Problem, jobs have a given equal length p ≥ 1 and each job executes on a single machine. In both problems jobs have given release times and deadlines in between which they need to execute. The goal is to find a feasible schedule, and moreover, for the equal length jobs problem, a feasible schedule that minimizes total completion time of the jobs. The third optimization problem, Offline prefetching and caching to minimize stall time belongs to a different field: we are given a sequence of n page requests and a cache of size k. We can evict a page from the cache and fetch a new page to replace it. This operation cannot be done in parallel and costs F time units. When a page request is served it costs 1 time unit, unless the page is not yet in the cache, then a stall time is generated until the corresponding fetch completes. The goal is to decide when to evict and fetch pages so as to minimize the total stall time.
Though quite different, those three problems were solved in a similar manner. Unlike previous works where the authors transform the solution of a relaxed integer linear program into an integer solution, we used a new technique which simplifies the linear programs, and allows us to get directly optimal integer solutions: our approach is based on the observation that only the structure of the solution matters in the objective function, jobs and pages don't appear namely. Therefore, we completely dissociate the resolution process into two phases. First a simplified linear program can be used to find an optimal skeleton for the solution, and it is only later that we need to worry about assigning jobs or pages to this skeleton: for scheduling problems, the skeleton is a sequence of slots, and the assignment maps jobs to slots; for the cache problem, the skeleton is a sequence of intervals and the assignment associates to every interval a page to evict at the beginning and a page to fetch at the end. Our skeletons are such that the assignment phase just comes down to running a greedy algorithm. Our contribution is that this strategy, where you don't compute the assignment in the linear program, leads to linear programs with very simple constraint matrices, which not only are totally unimodular, but are (the transpose of ) directed vertex adjacency matrices.
This allows us to reduce our scheduling problems into a shortest path problem and to reduce the caching problem into a min cost flow problem. Here is roughly our approach (see figure 1). The original linear programs have the following structure. The polytopes described by the linear programs have fractional vertices, but for the required objective all optimal vertices are integral. Now we project the solution space to a lower dimensional space. This is done by partitioning the variables, and replacing every set S of variables the original linear program by a single new variable that represents the sum over S. The nice thing is that the resulting linear program is now totally unimodular, and corresponds to shortest path or a min cost flow problem.
As a result, the tall/small job scheduling problem and the prefetch/caching problem can be solved in worst case time O(n 3 ) improving over respectively O(n 10 ) [4] and O * (n 18 ) [2]. Implementations are available from the authors home-pages.  We will first introduce our method on a basic scheduling problem. We have n jobs, each of the same length p. Every job j ∈ [1, n] comes with an interval [r j , D j ] consisting of a release time and a strict deadline. The goal is to find a schedule on m parallel machines, such that each job is assigned to an execution slot consisting of a particular machine and a time interval [s j , s j + p) ⊆ [r j , D j ]. In addition, all execution slots assigned to a particular machine must be disjoint. One possible application could be frequency allocation. A network operator has a link with m optical fiber strings. Users ask for allocations of a frequency band of fixed size, inside the large frequency band that the particular user devices can handle. The goal is to find an assignment which satisfies all users. In addition we want to find the solution (if it exists) that minimizes the total completion time of the jobs. In the standard Graham notation, this problem is called P |r j ; p j = p; D j | C j .
Simons [10] give a complicated greedy-backtrack algorithm running in time O(n 3 log log n), and later improved to O(mn 2 ) [11]. Recently Brucker and Kravchenko [6] gave another algorithm for it, using a completely different approach. While their algorithm has worse complexity it is interesting because of a generalization which permits to solve an open problem, namely minimizing the weighted total completion time, where jobs are given priority weights.
A generalization of the feasibility problem is to find a maximal set of jobs, which can all be scheduled between their release times and deadlines. This problem is still open. Even the more general problem, when jobs come with a weight, and the goal is to find a maximal weighted feasible job set, is not known to be NP-hard.

Previous work
First we observe that without loss of generality we can restrict ourselves to schedules where each execution slot starts at some release time plus a multiple of p, simply by shifting each slot as much to the beginning as possible. Let T = {r i + (a − 1)p : 1 ≤ i, a ≤ n} be this set of time points. And finally for a fixed schedule, if we number the execution slots from left to right, we can always reassign the j-th slot to the machine (j mod m) + 1. This way we don't need to take care of which machines the slots are assigned to, as long as there are at most m slots starting in every time interval of size p, which ensures that slots don't overlap on a particular machine. The linear program of [6] has a variable x jt for each job j and time t ∈ T , with the meaning that x jt = 1 if job j is executed in the slot [t, t + p). Then the program is to minimize jt (t + p)x jt subject to ∀j ∈ [1, n] : t∈T x jt = 1 : (every job cmpl.) ∀s ∈ T : s≤t<s+p j∈ [1,n] x jt ≤ m (no overlapping) It is quite clear that there is an integer solution to this linear program if and only if there is a feasible schedule. While this linear program is not totally unimodular, the authors of [6] were still able to round the fractional solution into an integer solution of the same cost.

Relaxing the linear program
The linear program above computes not only the time slots of the schedule, but also the assignment of jobs to slots. However once we are given the skeleton of a schedule, meaning a set of time slots, it is always possible to assign the jobs greedily in EDD fashion: assign to every slot the job with smallest deadline among the available jobs. We release the linear program from the job assignment, in order to obtain a simpler linear program which only computes a feasible skeleton.
We proceed in several steps. First we weaken equation (every job completes) into the inequality t∈T x jt ≥ 1. Then combining this new constraint with (allowed interval ) leads to ∀j ∈ [1, n] : Now for every pair s, t ∈ T , s ≤ t we sum (1) over all jobs j that have [r j , D j − p] ⊆ [s, t], upper-bounding the left hand side we obtain ∀s, t ∈ T , s ≤ t : The constraints are clearly necessary, and we will show later they are also sufficient to get the optimal solutions. We reduce the number of variables and group j x jt by setting y t := s≤t j x jt . Now y t represents the total number of slots up to time t. To simplify notations we introduce an additional time point t 0 < min T , and set T ′ = T ∪ {t 0 }. For any time t > t 0 , we define the functions round(t) := max{s ∈ T ′ : s ≤ t} and prec(t) := max{s ∈ T ′ : s < t}.
The two first inequalities force y t at the first and last time steps. In fact they are not necessary, but simplify the proof. The order inequalities ensure that (y t ) is a non decreasing sequence.
The load inequalities verify that there are never more than m slots overlapping, and the inclusion inequalities, are there to ensure that there is a feasible mapping from jobs to slots, as we show next. Except the equality y t0 = 0, the linear program has in every constraint exactly two variables, and with respective coefficients +1 and −1. So the dual of the constraint matrix is the incidence matrix of a directed graph, which means the constraint matrix is totally unimodular. Now this property is preserved when adding a row with a single +1 entry, corresponding to y t0 = 0. Therefore our linear program's constraints are in the from Ay ≤ b with A totally unimodular and b integer. This means that if the linear program has a solution then there is an optimal integer solution. Let (y t ) be an optimal integer solution to this linear program. It indeed defines the skeleton of a solution: at each time t ∈ T there will be y t − y prec(t) slots available for scheduling jobs. The standard greedy job-slot assignment, is defined as scheduling at every slot the job with the smallest deadline among the jobs that are not yet completed and already released.

Lemma 1 The greedy assignment produces a valid schedule.
Proof: We can notice that according to the second condition and the inclusion condition on [t 0 , maxT ], y maxT = n. We define V to be the multiset of time slots, such that slot [t, t + p] is contained y t − y prec(t) times. Therefore |V | = n. As mentioned in the previous section, by the load inequality, the slots can be assigned to machines without overlapping. So, it only remains to show that there exist assignments of jobs to slots, which respect release times and deadlines, and then that the greedy assignment is one of them.
Let U be the set of n jobs, and G(U, V, E) a bipartite graph where E contains all edges between a job j and a slot [ We have to show that this graph has an injection from U to V , and will use Hall's theorem for this, see [8] or [5].
For a set of jobs S, we denote the neighboring slots ∂S, as the set of all slots t such that there is a job j ∈ S with (j, t) ∈ E. We need to show that for every set S, |S| ≤ |∂S|, which by Hall's theorem, characterizes the existence of an injection. Let S be a set of jobs. Suppose S can be partitioned into S 1 ∪ S 2 such that for any jobs i ∈ S 1 and j ∈ S 2 the intervals [r i , D i − p] and [r j , D j − p] are disjoint. Then clearly ∂S is the disjoint union of ∂S 1 and ∂S 2 . Therefore we can without loss of generality assume that j∈S [r j , D j ] is a unique interval [r i , D j ], for i = argmin i∈S r i and j = argmax j∈S D j . Then |S| ≤ c ij . Also the number of slots in the interval [r i , D j ) is exactly y t − y s for s = prec(r i ), t = round(D j − p). From the inclusion inequality we get the required inequality and we conclude that there exist a valid assignment. Now since |V | = |U | = n, the injection is in fact a bijection, and there exists at least one perfect matching from jobs to slots with respect to release times and deadlines.
Proving that you can permute jobs in any of these matching to get the greedy matching is a quite standard in scheduling: let be two jobs i, j with D i < D j , and i is scheduled at some time t, while j is scheduled at some time s with r i ≤ s < t. Then it is possible to exchange the jobs i, j in their execution slots [s, s + p) and [t, t + p). By the use of a potential function, decreasing at each exchange, it is possible to transform our schedule in a so called earliest due date schedule. We conclude that since there exists at least a valid assignment, the greedy assignment is valid as well. 2 This means that an optimal integer solution can be found with a standard linear program solve. But our linear program describes in fact the dual of a minimum cost flow problem, with uncapacitated arcs, and a single supply node, which corresponds to a shortest path problem and can be solved in time O(N M ), where N is the number of variables and M the number of constraints [12, p.558].
Theorem 1 Our algorithm solves P |r j ; p j = p; D j | C j in worst case time O(n 4 ).
Proof: Given the instance m, p, r 1 , . . . , r n , D 1 , . . . , D n , we construct the set T of O(n 2 ) time points. Then we compute for every pair of jobs i, j the number of jobs c ij which need to be scheduled in [r i , D j ]. A naive algorithm does it in time O(n 3 ), which would be enough for us. However it can be solved in time O(n 2 ) using the following recursive formula. We assume jobs are indexed in order of release times. For convenience we set c n+1,j = 0. Then This permits to construct the graph G and find in time O(n 4 ) the optimal solution to the linear program, if there is one. Finally we do an earliest due date assignment of the jobs to the slots defined by the solution to the linear program in time O(n log n) using a priority queue. 2 Note that in this section we don't beat the best known algorithm for P |r j ; p j = p; D j | C j which is O(mn 2 ) [11]. However, it allows us to introduce our technique that will be used later on.
3 Scheduling tall and small jobs  In a parallel machine environment, sometimes maintenance tasks are to be done which involve all machines at the same time. Think of business meetings or inventory. Formally we are given n jobs of unit length p = 1, each job j comes with an integer release time and a deadline interval [r j , D j ] in which it must be scheduled. We distinguish two kinds of jobs. The first n 1 jobs are small jobs, in the sense that they must be scheduled on one of the m parallel machines, it does not matter which one. The n 2 = n − n 1 remaining jobs are tall jobs, in the sense that they must be scheduled on all the m machines at the same time.
A time slot is an interval [t, t + 1) for an integer boundary t. The goal is to find a feasible schedule, where each tall job is assigned to a different time slot, and each small job is assigned to a different (machine, time slot) pair for the remaining time slots. In addition the time slot to which some job j is assigned must be included in [r j , D j ].
This problem has been solved by Baptiste and Schieber [4], with a linear program using O(n 2 ) variables and O(n 2 ) constraints. The linear program is not totally unimodular, however they manage to show that for the particular objective function it always has an integer solution. We provide a linear program using only O(n) variables but still O(n 2 ) constraints, but whose constraint matrix is the incidence matrix of a directed graph, and can be solved in time O(n 3 ) with a shortest path algorithm.
Baptiste and Schieber showed that we can assume that the time interval ranges only from 1 to n, otherwise the problem could easily be divided into two disjoint subproblems.
In a similar way than before, we will denote by x t the total number of time slots assigned to tall jobs in [1, t + 1]. For convenience we set x 0 = 0. The number of small jobs that must be scheduled in [s, t] is k s,t = |{j : j ≤ n 1 , [r j , D j ] ⊆ [s, t]}| and the same for tall jobs is ℓ s,t = |{j : j > n 1 , [r j , D j ] ⊆ [s, t]}|. Consider the following linear program, which does not have an objective value.
Inequalities (3) make sure that (x t ) is a non decreasing sequence, (4) that only one tall job can be scheduled per unit-length interval, (5) that there are enough slots for the tall jobs and (6) that there are enough remaining slots for the small jobs.
Once again, the transpose of the constraint matrix is the adjacency matrix of an oriented graph, and the constant vector b is integer. As previously, it has optimal integer solutions.

Theorem 2 Fix an instance of the tall/small scheduling problem. There is an integer solution to this linear program if and only if there is a feasible schedule.
Proof: It is quite obvious that fixing (x t ) according to any feasible schedule will satisfy the constraints.
For the hard direction, let (x t ) be a solution to the linear program, we know it is integer. Then x t − x t−1 -which can be 0 or 1 -is the number of slots for tall jobs at time t. We will again use Hall's theorem to show that there is a valid assignment of the n 2 tall jobs to these slots. Inequality (5) for [s, t] = [1, n] forces x n ≥ n 2 . Now let be G(U, V, E) the bipartite graph, where U are the n 2 tall jobs, and V the x n slots. There is an edge between job j and time slot [t, t + 1] if it is included in [r j , D j ]. We have to show that for every subset S ⊆ U , the number of neighboring slots in V is at least |S|. Let s be the smallest release time among S and t be the largest deadline among S. Again it is sufficient to show this claim for connected sets S in the sense that ∪ j∈S [r j , D j ] = [s, t]. Now |S| ≤ ℓ s,t ≤ x t−1 − x s−1 , where the last expression is the number of slots in [s, t]. This completes the claim that there is a valid assignment from tall jobs to the slots.
For the small jobs, note that a s,t := (t − s) − (x t−1 − x s−1 ) is the number of remaining slots in [s, t] which are not assigned to tall jobs, and a s,t · m small jobs can fit in that interval. Again inequality (6) implies k s,t ≤ m · a s,t , and Hall's theorem shows that there is a valid assignment of small jobs to the remaining slots. 2 In the original paper [4] the author gave a linear program which is solved in expected time O(n 4 ) and worst case time O(n 10 ). Using the transformation into a shortest path problem allows us to improve this complexity.

Corollary 1 The tall/small scheduling problem can be solved in worst case time O(n 3 ).
Proof: As in the second section, we have a linear program with O(n) variables and O(n 2 ) constraints which can be produced in time O(n 2 ). We just take an arbitrary objective function in which all the variable coefficients are positive, and build the associated graph as in the previous section. Then we compute the all shortest paths from the source x 0 , in time O(n 3 ). If this computation detects a negative cycle, then the problem has no solution. Otherwise, we get the skeleton of a solution to the problem that minimize the total completion time of the tall jobs. Finally if there is a solution, the standard earliest due date assignment, first of tall jobs, then of small ones, produces a valid schedule in time O(n log n). 2 Here again, a direction that we are still exploring is to find another shortest path algorithm inspired from [11], better fitted for these specific graphs, that could improve this complexity.

Prefetching
Caches are used to improve the memory access times. In this context the memory unit is called a page, and is stored on a slow disk. The cache can store up to k pages. Now if a page request arrives, and the page is already in the cache, it can be served immediately, otherwise it must first be fetched from the disk, and that introduces a stall time of F units. In the latter case the new page replaces some other page currently in the cache. The idea of prefetching is to fetch a page even before it is requested, so as to reduce the stall time: During a fetch which evicts some page y replacing it by some page z, other requests can be served for pages currently in the cache and different from y or z. In the single disk model we consider here, only a single fetch can occur at the same time. The goal is, knowing in advance the complete request sequence, to come up with a prefetch schedule, which minimizes total stall time.  While the real life problem is on-line, and has been extensively studied by Cao et al. [7], the offline problem has first been solved in 1998 [3], by the use of a linear program, for which it was shown that it always has an optimal integer solution, while not being totally unimodular. Later in 2000 [2], a polynomial time algorithm was given modeling the problem as a multi commodity flow with some postprocessing. Formally the problem can be defined as follows.
The Offline Prefetching problem The input is a page request sequence x 1 , . . . , x n , an initial cache set C 1 , and a fetch duration F . Let k = |C 1 | be the cache size. A fetch is a tuple (s, y, e, z), where y, z are pages and s, t ∈ [1, n] are time points with s ≤ e ≤ s + F . The meaning is that at time s, the page y leaves the cache and at time e the page z enters the cache. Its cost, the induced stall time, is F − (e − s). The goal is to come up with a fetch sequence minimizing the total stall time, such that two fetches intersect in at most one time point, and such that every request can be served, i.e. ∀t ∈ [1, n] : x t ∈ C t , where C t is the cache at time t obtained from C t−1 by evicting/fetching all the pages that had to be evicted/fetched at time t. To simplify notation we assume that the request sequence contains at least k distinct pages, that C 1 consists of the first k distinct requests, and that at time 1, no page has left/entered the cache yet.
Albers, Garg and Leonardi defined a linear program with a characteristic variable for every fetch interval [s, e], and two additional characteristic variables for every pair (y, [s, e]) indicating whether page y enters (resp. leaves) the cache at the beginning (resp. the end) of the fetch [s, e]. Finally they show that the linear program has always an integer solution for the considered objective function.
As observed in [3] without loss of generality the page to be evicted at time t from the cache C t−1 is the page, whose next request is furthest in the future or which is never requested again. Also without loss of generality the page to be fetched at time t is the page whose next request starting from t is nearest in the future. Therefore all the information about the fetches is in the time intervals, and we will write a linear program which produces only the time intervals in which evictions/fetches occur. The actual pages have to be assigned in a post processing, in greedy manner as just mentioned. Rather than having a single variable for every interval and every page, we only count how many pages entered and how many left the cache in total since the beginning. This leaves us with O(n) instead of O(n 2 F ) variables. We denote by I t (resp. O t ) the total number of pages which entered (resp. left) the cache up to time t included. We get the following linear program.
Inequalities (7) make sure that (O t ) and (I t ) are non decreasing sequences, (8) that the cache cannot overflow, (9) that two fetches don't overlap in time, (10) that a fetch length is at most F and (11) that there are enough fetches to serve all requests. Note that |{x s , x s+1 , . . . , x t }| denotes the number of distinct page requests in [s, t], so (11) is the part of linear program that depends on the actual problem instance.
The optimal solution of this linear program is always integer, since it is totally unimodular (for the same reason as in previous section: its constraint matrix the transposed incidence matrix of a directed graph.) Theorem 3 Let (I t , O t ) be an optimal integer solution to the linear program. Then there is valid fetch sequence of the same cost, which can be built by greedy assignment.
Proof: First we observe that the cost function makes sure that O n = I n , which ensures that all interval are eventually closed. The solution defines m = O n intervals as follows. For every j = 1 . . . m, let s j be the smallest time such that O sj ≥ j and e j the smallest time such that I ej ≥ j. Then by (8) and (10) we have s j ≤ e j ≤ s j + F . Which means that all intervals [s j , e j ] are well defined and of length smaller or equal than F . Now by (9), e j ≤ s j+1 (otherwise, we would have I ej + 1 ≥ O ej > O sj+1 but I ej = j and O sj+1 = j + 1 by definition.), and this for all j < m, so the intervals do not overlap (but the ending point of one might be the starting point of another). Moreover the objective value of (I t , O t ), equals the total stall time of these intervals, for at each time t, the difference O t -I t is 1 if an interval is currently opened and is 0 otherwise. It remains to prove that the greedy assignement of pages to evict/fetch to each interval is such that all requests are served, i.e. it remains to show that the constraints (11) are sufficient. We denote by C s the cache obtained at time s, after all entrances and evictions that occur at time s. We will show that the following invariant holds in a solution of our linear program for every time s ∈ [1, n], ∀t ∈ [s, n] : I t − I s ≥ |{x s , . . . , x t }\C s |.
The invariant implies that if the number of pages requested in [s, t] but not in the cache at time s is a, then at least a pages must enter the cache somewhere in [s + 1, t]. In particular it means for t = s, that the page requested at time s will be the in the cache at that moment. The proof of (12) is by induction on s.
Basis case s = 1 Let t 0 be the greatest request time such that x t0 is not in C 1 . Then by the assumption that initially the cache contains the first k distinct requests, we have that for t < t 0 , {x 1 , . . . , x t } ⊆ C 1 . So the right hand side of (12) is 0 and (12) holds by (7). For t ≥ t 0 , since the intersection of {x 1 , . . . , x t } and C 1 is exactly k, the invariant holds by (11) and O 1 = I 1 = 0.
Induction case Assume the invariant holds for some s. Let's show that it also holds for s + 1. Several things can happen at time s + 1. Pages can leave the cache and pages can enter the cache. We will do these operations step by step, transforming I s into I s+1 and C s into C s+1 , and show that each step preserves the invariant (12). By induction hypothesis x s ∈ C s , so {x s , . . . , x t }\C s = {x s+1 , . . . , x t }\C s . Therefore, if nothing happens and no page enters or leaves the cache, then C s+1 = C s , I s = I s+1 and the invariant is preserved for s + 1. Now we deal with the case when there is some page movement at time s + 1, that is I s+1 > I s or O s+1 > O s or both. We artificially decompose this page movement in as many times as needed, so that at each time there is only one operation happening: a fetch or an eviction. The page movements at those intermediary times are set so as to alternatively evict and enter pages, among the O s+1 − O s pages to evict and the I s+1 − I s pages to enter. Of course if a fetch is pending at time s, that is |C s | = k − 1, then we start with entering a new page and otherwise if the cache is full, i.e. |C s | = k, we start with evicting a page. Since the number of total entrances and evictions up to some time can differ by at most one, it is possible to do so. Therefore, we need to do the induction case only in the case when a page is entering the cache or when one is leaving the cache but not both.
When page is entering the cache, we have I s+1 = I s + 1. Let z be the page entering the cache, and let t 0 ≥ s + 1 be the next request time of z. Then if t < t 0 , by the choice of z, all requests of x s+1 , . . . , x t must be in C s+1 , so the right hand side of (12) at time s + 1 is 0, and the inequality holds by (7). Now if t ≥ t 0 , since z ∈ C s+1 but z ∈ C s , the left hand side of (12) at time s + 1 has decreased by 1 compared to time s, but at the same time I s+1 = I s + 1, so both sides of the invariant decrease by 1 and by induction the inequality is preserved at time s + 1. Now consider the case when a page leaves the cache. Let y be the leaving page. Then I s = I s+1 and O s+1 = O s + 1. Let t 0 be the next request time of y or let t 0 = n + 1 if y is never requested again. Then if t < t 0 , removing y from C s+1 does not change the right hand side of (12) when replacing s by s + 1. The left hand side does not change either since no page enters the cache, and the inequality is preserved. For t ≥ t 0 however by the choice of the evicted page y, we have that C s+1 ⊆ {x s+1 , . . . , x t }. So the left hand side of (12) at time s + 1 is |{x s+1 , . . . , x t }| − (k − 1), and I s+1 = O s+1 − 1 since we have just evicted a page. Therefore, (12) holds by (11).   O2  O3  O4  O5  O6  O7  O8  O9  O10   I1   −F  +1  +1  +1  +1  +1  +1  +1  +1  +1   I10  I9  I8  I7  I6  I5  I4  I3  I2   0   1   It could be solved in time O(n 3 log n) using [9]. To solve it in O(n 3 ), when F = O(n), we first explode the source of supply F −1 into F −1 vertices of supply +1 and do the same with the sink of demand 1−F . The new graph has only sources of supply +1 and sinks of demand −1. Clearly there is a bijection between the min cost flows of the new and the original graph. Moreover a min cost flow matches sources to sinks such that the flow between a matched source/sink pair uses a shortest path (since the arcs have unbounded capacity) and such that the total distances are minimal. To obtain this flow we first compute the distances in the graph between all source/sink pairs, in time O(n 3 ) using Floyd-Marshall's algorithm. Then we construct the bi-partite sources/sinks graph, where every edge is weighted with the source-sink distance in the original graph. Then a minimum weighted perfect matching can be computed in time O(n 3 ) provided F ∈ O(n), using Edmond's algorithm with adapted data-structures. The optimal flow then is obtained by adding a unit flow on the shortest path between source i and sink j, for every edge of the matching corresponding to source i and sink j.
Finally to get an optimal solution for the primal linear program, we use the standard technique of computing a shortest path tree in the residual graph obtained from the flow [1, chapter 9]. 2

Conclusion
Further work would include trying to find other optimization problems where our technique may apply, and maybe generalize from it a general framework. We are also interested in improving the combinatorial algorithms that arise from the graph structures in the scheduling problems: indeed those graphs have, among others, the property that once the vertices drawn as points on a line, the arcs from left to right have positive weights and the ones from right to left negative. One idea for instance is to try and extract from Simons and Warmuth's algorithm a shortest path algorithm suitable for our class of graphs.
We wish to thank Arthur Chargueraud, Philippe Baptiste, Miki Hermann and Leo Liberti for helpful comments.