11institutetext: School of Computing Science, Newcastle University

# The Semantic Web takes Wing: Programming Ontologies with Tawny-OWL

Phillip Lord

###### Abstract

The Tawny-OWL library provides a fully-programmatic environment for ontology
building; it enables the use of a rich set of tools for ontology development,
by recasting development as a form of programming. It is built in Clojure – a
modern Lisp dialect, and is backed by the OWL API. Used simply, it has a
similar syntax to OWL Manchester syntax, but it provides arbitrary
extensibility and abstraction. It builds on existing facilities for Clojure,
which provides a rich and modern programming tool chain, for versioning,
distributed development, build, testing and continuous integration. In this
paper, we describe the library, this environment and the its potential
implications for the ontology development process.

##  1 Introduction

Ontology building remains a difficult and demanding task. Partly this is
intrinsic, but partly stems from the tooling. For example, while ontology
editors like Protégé [1] do allow manual ontology development, they are not
ideal for automation or template-driven development; for these reasons
languages such as OPPL[2] have been developed; these allow a slightly higher-
level of abstraction over the base OWL axiomatisation. However, they involve a
move away from OWL syntax, which in turn requires integration into which ever
environment the developers are using. There has also been significant interest
in collaborative development of ontologies, either using collaborative
development tools such as Web-Protege[3], or through copy-modify-merge
versioning[4].

In this work, we111Plurals are used throughout, and do not indicate multiple
authorship. take an alternative approach. Instead of developing tools for
ontology development, many of which are similar or follow on from software
development tools, we attempt to recast ontology development as a software
engineering problem, and then just use the standard tools that exist for
software engineering. We have achieved this through development of a library,
named Tawny OWL, that at its simplest operates as a domain specific language
for OWL, while still retaining the full capabilities of a modern programming
language with all this entails. We demonstrate the application of this library
to a standard exemplar - namely the Pizza Ontology[5], as well as several
other scenarios. Finally, we consider the implications of this approach for
enabling collaborative and more agile forms of ontology development.

##  2 Requirements

Interaction between OWL and a programming API is not a new idea. For example,
OWL2Perl[6] allows generation of Perl classes from an OWL Ontology, while the
OWL API allows OWL ontology development in Java[7]. The OWL API, however, is
rather unwieldy for direct ontology development; for example, it has a complex
type hierarchy, indirect instantiation of objects through factories, and a set
of change objects following a command design pattern; while these support one
of its original intended use case – building a GUI – they would make direct
ontology development cumbersome. One response to this is Brain[8, 9], which is
a much lighter weight facade over the OWL API also implemented in Java. Brain
is, effectively, type-less as expressions are generated using Strings; the API
distinguishes between OWL class creation (addClass) and retrieval (getClass),
throwing exceptions to indicate an illegal state. While Brain is useful, it is
not clear how an ontology should be structured in Java’s object paradigm, and
it suffers the major drawback of Java – an elongated compile-test-debug cycle,
something likely to be problematic for interactive development as the ontology
increases in size.

For programmatic ontology development, we wanted a much more interactive and
dynamic environment; something equivalent to the R environment for statistics,
where the ontology could be explored, extended and reworked on-the-fly,
without restarting. For this reason we choose to build in Clojure; this
language is a modern Lisp derivative with many attractive features: persistent
data structures; specialised mechanisms for storing state. It suffers somewhat
from being built on the Java Virtual Machine (JVM) –in particular this gives
it a rather slow start-up time – however, in this case, it was a key reason
for its use. Interoperability with the JVM is integrated deeply into Clojure
which makes building on top of the OWL API both possible and convenient. Like
all lisps, Clojure has three other advantages: first, it is untyped which, in
common with Brain, in this context, we consider to be an advantage222We do not
argue that type systems are bad; just that the are less appropriate in this
environment; second, it is highly dynamic – almost any aspect of the language
can be redefined at any time – and it has a full featured read-eval-print-loop
(REPL); finally, it has very little syntax, so libraries can manipulate the
look of the language very easily. Consider, for example, a simple class
definition as shown in Listing1, taken from a pizza ontology available at
https://github.com/phillord/tawny-pizza. The syntax has been designed after
Manchester syntax[10].

⬇

(defclass Pizza

:label "Pizza"

:comment

"An␣over-baked␣flat␣bread␣with␣toppings,␣originating␣from␣Italy."

)

Listing 1: A basic class definition

A more complex definition shows the generation of restrictions and anonymous
classes.

⬇

(defclass CheesyPizza

:equivalent

(owland Pizza

(owlsome hasTopping CheeseTopping)))

Listing 2: A Cheesy Pizza

These definitions bind a new symbol (Pizza and CheesyPizza) to a OWLAPI Java
object. These symbols resolve as a normal Var does in Clojure. Strictly, this
binding is not necessary (and can be avoided if the user wishes), however this
provides the same semantics as Brain’s addClass and getClass – classes,
properties, etc must be created before use; this is a valuable feature
protecting against typing errors[11].

###  2.1 Lisp Terminology

Here we give a brief introduction to Clojure and its terminology. Like all
lisps, it has a regular syntax consisting of parenthesis delimited (lists),
defining an expression. The first element is usually a function, giving lisps
a prefix notation. Elements can be literals, such as strings e.g. "Pizza",
symbols e.g. defclass or keywords e.g. :equivalent. Symbols resolve to their
values, keywords resolve to themselves, and literals are, well, literal.
Unlike many languages, these constructs are directly manipulable in the
language itself which combined with macros enable extension of the language.

##  3 A Rich Development Environment

There are a dizzying array of ontology development tools available[12].
Probably the most popular is Protégé; while it provides a very rich
environment for viewing and interacting with an ontology, it lacks many things
that are present in most IDEs. For instance, it lacks support for version
control or adding to ChangeLogs; it is not possible to edit documentation
along side the ontology; nor to edit code in other languages when, for
instance, driving a build process, or using an ontology in an application.

We have previously attempted to work around this problem by providing support
for Manchester syntax – OMN – within Emacs through omn-mode[13]; while this
provides a richer general-purpose environment, the ontology environment is
comparatively poor. In particular, only syntactic completion is available,
there is no support for documentation look-up, beyond file navigation.
Finally, we used Protégé (and the OWL API) to check syntax, which required a
complete re-parse of the file, and with relatively poor feedback from Protégé
when errors occurred333This is not a criticism of the Protégé interface; it
was not designed to operate on hand-edited files.

With tawny, using a general purpose programming language, a richer development
environment comes nearly for free. In this paper, we describe the use within
Emacs; however, support for Clojure is also available within Eclipse,
IntelliJ, Netbeans and other environments[14]. Compared with direct editing of
OMN files, this provides immediate advantages. The use of paren delimiters
makes indentation straight-forward, well-defined, and well-supported; advanced
tools like paredit ensures that expressions are always balanced. Clojure
provides a REPL, and interaction within this allows more semantic completion
of symbols even when they are not syntactically present in the buffer444We
follow Emacs terminology here – a buffer is a file being edited, which is
common when using levels of abstraction (Section 4) or external OWL files
(Section 8). Syntax checking is easy, and can be performed on buffer, marked
region or specific expression. New entities can be added or removed from the
ontology on-the-fly without reloading the entire ontology, enabling
progressive development. We have also provided support for documentation look-
up of OWL entities; this is hooked into Clojure’s native documentation
facility, so should function within all development environments. We do not
currently provide a rich environment for browsing ontologies, except at the
code level; however, Protégé works well here, reloading OWL files when they
are changed underneath it. Similarly, omn-mode can be used to view individual
generated OMN files.

##  4 Supporting Higher Levels of Abstraction

Most ontologies include a certain amount of “boilerplate” code, where many
classes follow a similar pattern. Tools such as OPPL were built specifically
to address this issue; with tawny, the use of a full programming language,
makes the use of levels of abstraction above that in OWL straight-forward. We
have used this in many areas of Tawny; at its simplest, by providing
convenience macros. For example, it is common-place to define many subclasses
for a single superclass; using OMN each subclass must describe its superclass.
Within tawny, a dynamically-scoped block can be used as shown in Listing 3. As
shown here, disjoint axioms can also be added[15]; and, not used here,
covering axioms[16]. The equivalent OMN generated by these expressions is also
shown in Listing 4.

⬇

(as-disjoint-subclasses

PizzaBase

(defclass ThinAndCrispyBase

:annotation (label "BaseFinaEQuebradica" "pt"))

(defclass DeepPanBase

:annotation (label "BaseEspessa" "pt")))

Listing 3: Subclass Specification

⬇

Class: piz:ThinAndCrispyBase

Annotations:

rdfs:label "BaseFinaEQuebradica"@pt

SubClassOf:

piz:PizzaBase

DisjointWith:

piz:DeepPanBase

Class: piz:DeepPanBase

Annotations:

rdfs:label "BaseEspessa"@pt,

SubClassOf:

piz:PizzaBase

DisjointWith:

piz:ThinAndCrispyBase

Listing 4: Subclasses in OMN

It is also possible to add suffixes or prefixes to all classes created within
a lexical scope. For example, we can create classes ending in Topping as shown
in Listing 5. While similar functionality could be provided with a GUI, this
has the significant advantage that the developers intent remains present in
the source; so subsequent addition of new toppings are more likely to preserve
the naming scheme.

⬇

(with-suffix Topping

(defclass GoatsCheese)

(defclass Gorgonzola)

(defclass Mozzarella)

(defclass Parmesan))

Listing 5: Adding Suffixes

Tawny also includes initial support for ontology design patterns; in
particular, we have added explicit support for the value partition[17]. This
generates classes, disjoints and properties necessary to fulfil a pattern, but
is represented in Tawny succinctly (Listing 6)

⬇

(p/value-partition

Spiciness

[Mild

Medium

Hot])

Listing 6: A Value Partition

While some abstractions are generally useful, an important advantage of a
full-programmatic language for OWL is that abstractions can be added to any
ontology including those which likely to be useful only within a single
ontology. These can defined as functions or macros in the same file as their
use. For example, within the pizza ontology, Listing 7 generates two pizzas –
in each case the pizza class comes first, followed by constituent parts; a
closure axiom is added to each pizza. As well, as being somewhat more concise
than the equivalent OMN, this approach also has the significant advantage that
it is possible to change the axiomatisation for all the named pizzas by
altering a single function; this is likely to increase the consistency and
maintainability of ontologies.

⬇

(generate-named-pizza

[MargheritaPizza MozzarellaTopping TomatoTopping]

[CajunPizza MozzarellaTopping OnionTopping PeperonataTopping

PrawnsTopping TobascoPepperSauce TomatoTopping]

Listing 7: Generating Named Pizzas

##  5 Separating Concerns for Different Developer Groups

One common requirement in ontology development is a separation of concerns;
different contributors to the ontology may need different editing
environments, as for instance with RightField or Populous[18]. Tawny enables
this approach also; here, we describe how this enables internationalisation.
Originally, the pizza ontology had identifiers in English and Portuguese but,
ironically, not Italian. While it would be possible to have a translator
operate directly on a tawny source file, this is not ideal as they would need
to need to embed their translations within OWL entity definitions as shown in
Listing 3; this is likely to be particularly troublesome if machine assisted
translation is required due to the non-standard format. We have, therefore
added support with the polyglot library. Labels are stored in a Java
properties file (Listing 8) and are loaded using a single Lisp form (Listing
9). Tawny will generate a skeleton resources file, with no translations, on
demand, and reports missing labels to the REPL on loading.

⬇

AnchoviesTopping=Acciughe Ingredienti

ArtichokeTopping=Carciofi Ingredienti

AsparagusTopping=Asparagi Ingredienti

Listing 8: Italian Resources

⬇

(tawny.polyglot/polyglot-load-label

"pizza/pizzalabel_it.properties" "it")

Listing 9: Loading Multi-Lingual Labels

Currently, only loading labels is supported in this way, but extending this to
comments or other forms of annotation is possible. While, in this case, we are
loading extra-logical aspects of the ontology from file, it would also be
possible to load logical axioms; for instance, named pizzas (Section 4) could
be loaded from text file, spreadsheet or database.

##  6 Collaborative and Distributed Development

Collaborative development is not a new problem; many software engineering
projects involve many developers, geographically separated, in different time
zones, with teams changing over time. Tools for enabling this form of
collaboration are well developed and well supported. Some of these tools are
also available for ontology development; for instance, Web-Protégé enables
online collaborative editing. However, use of this tool requires installation
of a bespoke Tomcat based server, nor does it yet support offline, concurrent
modification[3].

Alternatively, the ContentCVS system does support offline concurrent
modification. It uses the notion of structural equivalence for comparison and
resolution of conflicts[4]; the authors argue that an ontology is a set of
axioms. However, as the named suggests, their versioning system mirrors the
capabilities of CVS – a client-server based system, which is now considered
archaic.

For tawny, the notion of structural equivalence is not appropriate;
critically, it assumes that an ontology is a _set_ of axioms. This is not true
with tawny, for two reasons: first, tawny requires definition of classes
before use, so source code cannot be arbitrarily re-ordered; secondly, even
where this is not the case, only the ontology axioms are a set. Programmer
intent is often represented through non-axiomised sections of the code –
whitespace, indentation and even comments which may drive a “literate”
development approach. A definition of a difference based purely on
axiomatisation cannot account for these differences; the use of a line-
oriented syntactic diff will.

We argue here that by provision of an attractive and well-supported syntax, we
do not need to provide specific collaborative tooling. Tawny itself has been
built using distributed versioning systems (first mercurial and then git).
These are already advanced systems supporting multiple workflows including
tiered development with authorisation, branching, cherry-picking and so on.
While ontology-specific tooling have some advantages, it is unlikely to
replicate the functionality offered by these systems, aside from issues of
developer familiarity.

Later, we also describe support for testing, which can also ease the
difficulty of collaborative working (Section 9).

##  7 Enabling Modularity

Tawny provides explicit support for name spacing and does this by building on
Clojure’s namespace support. It is possible to build a set of ontologies
spread across a number of different files. Normally, each file contains a
single namespace; tawny mirrors this, with each namespace containing a single
ontology, with a defined IRI.

OWL itself does not provide a distribution mechanism for ontologies; the IRI
of an ontology does not need to resolve. In practice, this is often a
distribution mechanism; by default Protégé will check for resolution if other
mechanisms fail; OBO ontologies, for example, are all delivered from their
IRI.

In contrast, Tawny builds on the Clojure environment; most projects are built
using the Leiningen tool which, in turn, uses the repository and dependency
management from Maven. When building the Pizza ontology in Tawny, the build
tool will fetch Tawny itself, the OWL API and HermiT, and their dependencies.
Ontological dependencies can be fetched likewise. Maven builds come with a
defined semantics for versioning, including release and snapshot
differentiation. A key advantage of this system is that multiple versions of a
single ontology can be supported, with different dependency graphs.

##  8 Coping With Semantics Free Identifiers

Importing one ontology from another is straight-forward in tawny. However, not
all ontologies are developed using tawny; we need to be able interact with
external ontologies only accessible through an OWL file. Tawny provides
facilities for this use-case: the library reads the OWL file, creates symbols
for all entities555It is possible to choose a subset, then associates the
relevant Java object with this symbol. This approach is reasonably scalable;
tawny can import the Gene Ontology within a minute on a desktop machine.
Clojure is a highly-dynamic language and allows this form of programmatic
creation of variables as a first-class part of the language; so an ontology
read in this way functions in every sense like a tawny native ontology.
Ontology classes can be queried for their documentation, auto-completion works
and so forth.

However, there is a significant problem with this import mechanism. Tawny must
create a symbol for each OWL entity in the source ontology. By default, tawny
uses the IRI fragment for this purpose; while Clojure symbols have a
restricted character set which is not the same as that of the IRI fragment, in
practice this works well. However, this is unusable for ontologies built
according to the OBO ontology standard, which uses semantics-free, numeric
identifiers such as OBI_0000107. While this is a valid Clojure symbol, it is
unreadable for a developer. This issue also causes significant difficulties
for ontology development in any syntax; OMN is relatively human-readable but
ceases to be so when all identifiers become numeric. We have previously
suggested a number of solutions to this problem either through the use of
comments or specialised denormalisations[19], or through the addition of an
Alias directive providing a mapping between numeric and readable
identifier[20]. However, all of these require changes to the specification and
tooling updates, potentially in several syntaxes.

For tawny, we have worked around this problem by enabling an arbitrary mapping
between the OWL entity and symbol name [21]. For OBO ontologies, a syntactic
transformation of the rdfs:label works well. Thus, OBI_0000107 can be referred
to as provides_service_consumer_with, while BFO_0000051 becomes the rather
more prosaic has_part.

While this solves the usability problem, it causes another issue for ontology
evolution; the label is open to change, independently of any changes in
semantics; unfortunately, any dependent ontology built with tawny will break,
as the relevant symbol will no longer exist. This problem does not exist for
GUI editors such as Protégé because, ironically, they are not WYSIWYG – the
ontology stores an IRI, while the user sees the label; changes to labels
percolate when reloading the dependent ontology. Tawny provides a solution to
this; it is possible to memorise mappings between symbols and IRIs at one
point in time[22]. If the dependency changes its label, while keeping the same
IRI, Tawny will recognise this situation, and generate a deprecated symbol;
dependent ontologies will still work, but will signal warnings stating that a
label has changed and suggesting appropriate updates. Currently these must be
performed manually, although this could be automated.

##  9 Enabling Unit Testing and Continuous Integration

Unit testing is a key additions to the software development process which has
enabled more agile development. Adapting this process for ontology development
has previously been suggested[23], and implemented as a plugin to Protégé
[24]. To add this capability to tawny, we have integrated reasoning; at the
time of writing, only ELK[25] is available as a maven resource in the Maven
Central repository, therefore we have developed a secondary maven build for
HermiT which allows use of this reasoner also[26]666Available at
http://homepages.cs.ncl.ac.uk/phillip.lord/maven/, or on Github, so both these
reasoners are available for use; others can be added trivially as they are
mavenised. A number of test frameworks exist in Clojure; here we use
clojure.test. As shown in Listing 10, we check that various inferences have
occurred (or not as appropriate), using the isuperclass? predicate. We have
also added support for “probe” classes. In our second test, we use the with-
probe-entities macro; this adds a subclass of VegetarianPizza and CajunPizza –
as the latter contains meat, this should result in an incoherent ontology if
both definitions are correct; probe entities are automatically removed by the
macro, returning the ontology to its previous state, ensuring independence of
tests.

⬇

(deftest CheesyShort

(is (r/isuperclass? p/FourCheesePizza p/CheesyPizza))

(is (r/isuperclass? p/MargheritaPizza p/CheesyPizza))

(is

(not (r/isuperclass? p/MargheritaPizza p/FourCheesePizza))))

(deftest VegetarianPizza

(is

(r/isuperclass? p/MargheritaPizza p/VegetarianPizza))

(is

(not

(o/with-probe-entities

[c (o/owlclass "probe"

:subclass p/VegetarianPizza p/CajunPizza)]

(r/coherent?)))))

Listing 10: Unit Testing a Pizza Ontology

The use of Unit testing in this way has implications beyond simple ontology
development; it also allows a richer form of continuous integration where
dependent ontologies can be developed by independent developers, but
continuously checked for breaking changes. The tawny pizza ontology, for
example, is currently being tested using Travis777http://travis-ci.org.
Unlike, other ontology CI systems[27], this requires no installation,
integrates directly with the DVCS in use. It is also useful for integration
with software that operates on the ontology; for example, both our version of
Hermit, the OWL API and tawny-owl are built and tested using this tool.

##  10 Discussion

In this paper, we have described Tawny, a library which enables the user to
develop ontologies, using the tooling and environments that have long been
available to programmers. Although they both involve producing artifacts with
strong computational properties ontology development and software engineering
have long been disjoint. This has significant negative impacts; there are far
more programmers than knowledge engineers, and as a result the tooling that
they use is far better developed. Tawny seeks to address this, not by
providing richer tools for ontology development, but by recasting ontology
development as a form of programming.

By allowing knowledge engineers to use any level of abstraction that they
choose, tawny can also improve current knowledge engineering process
significantly. It can help to remove duplication, for example, in class names.
It can clearly delineate disjoint classes protecting against future additions;
this helps to address a common ontological error[28]. It is also possible to
model directly using common ontology design patterns generating many axioms in
a succinct syntax. Bespoke templates can be built for a specific ontology;
this mirrors functionality of tools like OPPL[2], but uses the power of a full
programming language and environment. Trivially, for example, tawny can log
its activity and comes with debugger support.

Of course, direct use of a programmatic library like tawny is not suitable for
all users; however, even for these users a library like tawny could be useful.
It is straight-forward to integrate ontologies developed directly with tawny
as a DSL with knowledge stored in other formalisms or locations. In this
paper, we described loading multi-lingual labels from properties files,
isolating the translator from the ontology, and interacting with OWL files
generated by another tool. It would also be possible to load axioms from a
database or spreadsheet, using existing JVM libraries.

While with tawny, we have provided a programmatic alternative to many
facilities that exist in other tools, we also seek to provide tooling for a
more agile and reactive form of ontology development. Current waterfall
methodologies, exemplified by BFO-style realism lack agility, failing to meet
the requirement for regular releases to address short-comings, as has been
seen with both BFO 1.1[29] and BFO 2.0[30]. Likewise, the OBO foundry places
great emphasis on a review process which is, unfortunately, backlogged[31] –
in short, as with waterfall software methodologies, the centralised aspects of
this development model appear to scale poorly.

Tawny uses many ready-made and well tested software engineering facilities:
amenability to modern DVCS, a versioning and release semantics, a test
framework and continuous integration. The provision of a test environment is
particularly important; while ontology developers may benefit from testing
their own ontologies, the ability to contribute tests to their ontological
dependencies is even more valuable. They can provide upstream developers
precise and executable descriptions of the facilities which they depend on;
this gives upstream developers more confidence that their changes will not
have unexpected consequences. When this does happen, downstream developers can
track against older versions of their dependencies, obviating the need for co-
ordination of updates; when they do decide to update, the re-factoring
necessary to cope with external changes will be supported by their own test
sets; finally, continuous integration will provide early warning if their own
changes impact others. In short, tawny provides the tools for a more pragmatic
and agile form of ontology development which is more suited to fulfilling the
changing and varied requirements found in the real world[32].

## References

  * [1] Various:  The protégé; ontology editor and knowledge acquisition system.  http://protege.stanford.edu/ [Online. last-accessed: 2012-10-07 18:08:04]
  * [2] Egana Aranguren, M., Stevens, R., Antezana, E.:  Transforming the axiomisation of ontologies: The ontology pre-processor language.  Nature Precedings (Dec 2009)
  * [3] Tudorache, T., Nyulas, C., Noy, N.F., Musen, M.A.:  WebProtégé: a collaborative ontology editor and knowledge acquisition tool for the web.  Semantic Web
  * [4] Jiminez Ruiz, E., Grau, B.C., Horrocks, I., Berlanga, R.:  Supporting concurrent ontology development: Framework, algorithms and tool.  Data & Knowledge Engineering 70(1) (Jan 2011) 146–164
  * [5] Stevens, R.:  Why the pizza ontology tutorial?  http://robertdavidstevens.wordpress.com/2010/01/22/why-the-pizza-ontolo%gy-tutorial/ [Online. last-accessed: 2012-11-09 22:37:14] (2010)
  * [6] Kawas, E., Wilkinson, M.D.:  Owl2perl: creating perl modules from owl class definitions.  Bioinformatics 26(18) (Sep 2010) 2357–2358
  * [7] Bechhofer, S., Volz, R., Lord, P.:  Cooking the semantic web with the OWL API.  In: Internaional Semantic Web Conference. (2003) 659 – 675
  * [8] Croset, S., Overington, J., Rebholz-Schuhman, D.:  Brain: Biomedical knowledge manipulation.  Bioinformatics (2013) Submitted.
  * [9] loopasam:  Brain.  https://github.com/loopasam/Brain
  * [10] Horridge, M., Patel-Schneider, P.:  Owl 2 web ontology language manchester syntax.  http://www.w3.org/TR/owl2-manchester-syntax/ (2012)
  * [11] Lord, P.:  Owl concepts as lisp atoms.  http://www.russet.org.uk/blog/2254 [Online. last-accessed: 2012-10-25 01:36:03] (2012)
  * [12] Bergman, M.:  The sweet compendium of ontology building tools.  http://www.mkbergman.com/862/the-sweet-compendium-of-ontology-building-%tools/ (2010)
  * [13] Lord, P.:  Ontology building with emacs.  http://www.russet.org.uk/blog/2161 [Online. last-accessed: 2012-07-26 09:28:46] (2012)
  * [14] Various:  Getting started - clojure documentation - clojure development.  http://dev.clojure.org/display/doc/Getting+Started [Online. last-accessed: 2013-01-29 08:36:13]
  * [15] Lord, P.:  Disjoints in clojure-owl.  http://www.russet.org.uk/blog/2275 [Online. last-accessed: 2013-02-11 09:34:50] (2012)
  * [16] Stevens, R.:  Closing down the open world: Covering axioms and closure axioms.  http://ontogenesis.knowledgeblog.org/1001 [Online. last-accessed: 2012-06-19 16:13:39] (2011)
  * [17] Rector, A.:  Representing specified values in owl: “value partitions” and “value sets”.  W3C Working Group Note (2005)
  * [18] Jupp, S., Horridge, M., Iannone, L., Klein, J., Owen, S., Schanstra, J., Wolstencroft, K., Stevens, R.:  Populous: a tool for building owl ontologies from templates.  BMC Bioinformatics 13(Suppl 1) (2011) S5
  * [19] Lord, P.:  Obo format and manchester syntax.  http://www.russet.org.uk/blog/1470 [Online. last-accessed: 2012-06-19 16:32:49] (2009)
  * [20] Lord, P.:  Semantics-free ontologies.  http://www.russet.org.uk/blog/2040 [Online. last-accessed: 2012-06-19 16:32:22] (2012)
  * [21] Lord, P.:  Clojure owl 0.2.  http://www.russet.org.uk/blog/2303 [Online. last-accessed: 2012-12-03 16:28:51] (2012)
  * [22] Lord, P.:  Remembering the world as it used to be.  http://www.russet.org.uk/blog/2316 [Online. last-accessed: 2013-01-11 23:00:11] (2013)
  * [23] Vrandec̆ić, D., Gangemi, A.:  Unit tests for ontologies.  In: In OTM Workshops (2
  * [24] Drummond, N.:  Co-ode & downloads & the owl unit test framework.  http://www.co-ode.org/downloads/owlunittest/ [Online. last-accessed: 2013-01-28 15:22:03]
  * [25] Kazakov, Y., Krötzsch, M., Simancik, F.:  Elk reasoner: Architecture and evaluation.  In: Proceedings of the 1st International Workshop on OWL Reasoner Evaluation (ORE-2012). (2012)
  * [26] Various:  Hermit reasoner: Home.  http://hermit-reasoner.com/
  * [27] Mungall, C., Dietze, H., Carbon, S., Ireland, A., Bauer, S., Lewis, S.:  Continuous integration of open biological ontology libraries.  http://bio-ontologies.knowledgeblog.org/405 (2012)
  * [28] Rector, A., Drummond, N., Horridge, M., Rogers, J., Knublauch, H., Stevens, R., Wang, H., Wroe, C.:  OWL pizzas: Practical experience of teaching OWL-DL: common errors & common patterns.  Engineering Knowledge in the Age of the Semantic Web (2004) 63–81
  * [29] Various:  New version of bfo 1.1 available.  https://groups.google.com/d/topic/bfo-discuss/HQSnudUUM4E/discussion
  * [30] Various:  Proposal for an official bfo 1.2 release.  https://groups.google.com/d/topic/bfo-discuss/iKBlfDPv5GM/discussion
  * [31] OBO Foundry Outreach Working Group:  New obo foundry tracker for feedback, requests, and other issues.  http://sourceforge.net/mailarchive/message.php?msg_id=30391720
  * [32] Lord, P., Stevens, R.:  Adding a little reality to building ontologies for biology.  PLoS One (2010)


# Fast computation of approximate weak common intervals in multiple
indeterminate strings

Daniel Doerr  Institute for Medical Biometry and Bioinformatics, Heinrich
Heine University, Düsseldorf, Germany  Bernard M.E. Moret  School of Computer
and Communication Sciences, EPFL, CH-1015, Lausanne, Switzerland

###### Abstract

In ongoing work to define a principled method for syntenic block discovery and
structuring, work based on homology-derived constraints and a generalization
of common intervals, we faced a fundamental computational problem: how to
determine quickly, among a set of indeterminate strings (strings whose
elements consist of subsets of characters), contiguous intervals that would
share a vast majority of their elements, but allow for sharing subsets of
characters subsumed by others, and also for certain elements to be missing
from certain genomes. An algorithm for this problem in the special case of
determinate strings (where each element is a single character of the alphabet,
i.e., “normal” strings) was described by Doerr et al., but its running time
would explode if generalized to indeterminate strings. In this paper, we
describe an algorithm for computing these special common intervals in time
close to that of the simpler algorithm of Doerr et al. and show that can
compute these intervals in just a couple of hours for large collections (tens
to hundreds) of bacterial genomes.

##  1 Introduction

The rapidly increasing number of whole-genome sequences in public repositories
is creating a demand for an effective means of comparing multiple whole
genomes. Such comparisons not only enrich our knowledge about present-day
organisms, they also open a window into their evolutionary past. Many genomes
of interest (such as vertebrate genomes) are much too large for base-pair by
base-pair comparisons and such comparisons are made all the harder when
multiple genomes are being considered. Thus large-scale whole-genome
comparative studies are limited by two factors: (i) the size of the genomes
(typically in the billions of base pairs for vertebrates) and (ii) the number
of genomes one wishes to compare. The standard approach to the first problem
over the last 15 years has been to decompose the genomes into _syntenic
blocks_. In their most common form, syntenic blocks are contiguous blocks
ranging from $10^{4}$ to $10^{7}$ base pairs that form clearly conserved units
across the genomes—as attested by the conservation of closely related sets of
_genomic markers_. Rather than comparing genomes base-pair by base-pair, or
even marker by marker, one compares them syntenic block by syntenic block. In
other words, genomes with billions of base pairs are represented by several
thousands of syntenic blocks. As most of the work on such decompositions was
carried out using genes as genomic markers, one can view the ordering of
syntenic blocks along each genome as a generalization of the ordering of genes
along each genome: just as genes in each genome are homologous to genes in
other genomes and form gene families, so do syntenic blocks. (The other
problem, of having to compare multiple genomes, is perhaps best carried out in
a phylogenetic context, where direct comparisons, the number of which grows
quadratically with the number of genomes, are eschewed in favor of comparisons
between descendants and common ancestors [7], but it is not our focus here.)

In ongoing work to define a principled method for syntenic block discovery and
structuring, work based on homology-derived constraints [4] and a
generalization of _common intervals_ [5], we faced a fundamental computational
problem: how to determine quickly, among a set of indeterminate strings
(strings whose elements consist of subsets of characters), contiguous
intervals that would share a vast majority of their elements, but allow for
sharing subsets of characters subsumed by others, and also for certain
elements to be missing for certain genomes. An algorithm for this problem in
the special case of determinate strings (where each element is a single
character of the alphabet, i.e., “normal” strings) was described by Doerr _et
al._ [3], but its running time would explode if generalized in the “obvious”
way to indeterminate strings. In this paper, we describe an algorithm for
computing these special common intervals in time close to that of the simpler
algorithm of Doerr _et al._ [3] that makes it possible to compute these
intervals at the scale of complete vertebrate genomes.

##  2 Background

We begin with a quick sketch of our approach to the discovery of syntenic
blocks in multiple genomes, to put in proper perspective the computational
problem that is the focus of this paper.

The identification of large conserved regions across genomes, i.e., syntenic
blocks, is in itself useful, as it allows us to circumscribe regions of
interest and focus our studies on these. However, syntenic blocks can be also
be used to describe whole genomes at a higher level of abstraction, enabling
further analyses to be conducted at this new level. Indeed, the number of
levels of decomposition (or abstraction) need not be limited to two or three:
at one extreme, each position in a genomic sequence can be viewed as a
syntenic block while, at the other extreme, each whole genome can be viewed as
a single syntenic block. Neither extreme is of much interest and the study of
short subsequences is well advanced, so our design starts at the level of
_genomic markers_ , that is, short sequences that are both short enough to be
found in (exactly or much) the same form in other genomes and long enough that
their conservation is statistically significant. Such markers can be obtained
through genome segmentation [6] or, more commonly, using already available
annotations of the genomic sequences that mark conserved regions such as
noncoding and protein-coding genes.

Our method relies on predefined _homologies_ between these genomic markers.
(Recall that homology is an equivalence relation that denotes shared
ancestry.) Based on homologies between genomic markers, we derive homologies
between syntenic blocks using the formal criteria first enunciated by
Ghiurcuta and Moret [4]. Informally, according to their criterion, two
syntenic blocks are homologous if each of their constituent genomic markers is
homologous to a genomic marker in the opposite block.

From a conceptual point of view, shared origin is a common feature of genomic
markers; in fact, a perfect oracle for homology would not provide very useful
information—it would lack differentiability and thus be ineffective for
synteny analysis. It is thus fortunate, in a way, that the divergence of
sequences through mutational changes renders inference of homology impossible
after a certain amount of time. Unfortunately, this same divergence renders
the process of homology inference error-prone and unlikely to result in an
equivalence relation. Most current algorithms for homology inference rely on
some type of clustering, often based almost entirely on sequence similarity,
but homology should also take into account relative placements within each
genome. In other words, homology assignments between markers are required for
synteny analysis, yet their embedding in syntenic blocks provides a powerful
indicator of true homology. We break this circularity by not requiring that
the inferred homologies describe an equivalence relation: we drop the
requirement of transitivity. This relaxation leads to a decisive advantage in
practice: genomic markers no longer need to be clustered into families, so
that the synteny analysis can be directly performed on pairwise sequence
similarities.

We address the problem of identifying syntenic blocks under weak homology by
identifying sets of _approximate weak common intervals_ (AWCIs) in
indeterminate strings—and this is the computational problem we address in this
paper. Informally, indeterminate strings are sequences in which each position
is associated with a (nonempty) subset of alphabet characters. Two intervals
in two indeterminate strings are _weak common intervals_ if each position of
both intervals contains at least one character that is also present in the
other interval. If a character encoding for two indeterminate strings
represents (non-transitive) homologies between markers of two genome
sequences, weak common intervals are equivalent to homologous syntenic blocks.
Using AWCIs, we identify syntenic blocks that mildly violate the synteny
criterion by allowing for a limited number of missing homology statements. We
do this to account for (i) false negative errors in the inferred homology
assignment and (ii) insertions or deletions of one or a few genomic markers
that are a result of genome evolution.

##  3 A First Algorithm

###  3.1 Basic definitions

We recall some basic definitions of common intervals in indeterminate strings
from Doerr  _et al._ [3]. An indeterminate string $S$ with $n$ index positions
is a string over the power set $\mathcal{P}(\Sigma)\setminus\emptyset$ of an
alphabet $\Sigma$. That is, for each $i$, $1\leq i\leq n$, we have
$S[i]\subseteq\Sigma$ and $S[i]\neq\emptyset$, where $S[i]$ denotes the
character set associated with the $i$-th position in $S$. We denote the
_length_ of an indeterminate string $S$ with $n$ index positions by $|S|\equiv
n$ and its _cardinality_ , i.e., the number of _all_ elements in $S$, by
$\|S\|\equiv\sum_{i=1}^{n}|S[i]|$. Two positions $i$ and $j$, $1\leq i\leq
j\leq|S|$, induce the (indeterminate) _substring_ $S[i,j]\equiv
S[i]\,S[i+1]\,\dots\,S[j]$. To distinguish intervals in different
indeterminate strings, we indicate the affiliation of an interval $[i,j]$ to
indeterminate string $S$ by the subscript notation $[i,j]_{S}$.

Common intervals were first defined on permutations [5] and subsequently
extended to ordinary, then indeterminate, strings [2, 3]. The idea behind
common intervals is to compare substrings based on their character sets. The
_character set_ of an indeterminate string $S$ is defined as
$\mathcal{C}(S)\equiv\bigcup_{i=1}^{n}S[i]$. Given two indeterminate strings
$S$ and $T$, two intervals, $[i,j]$ in $S$ and $[k,l]$ in $T$, are _weak
common intervals_ with _common character set_
$C=\mathcal{C}(S[i,j])\cap\mathcal{C}(T[k,l])$ if for each $i^{\prime}$,
$i\leq i^{\prime}\leq j$, we have $C\cap S[i^{\prime}]\neq\emptyset$ and for
each $k^{\prime}$, $k\leq k^{\prime}\leq l$, we have $C\cap
T[k^{\prime}]\neq\emptyset$. We study a variant of weak common intervals that
tolerates a limited number of insertions and deletions. Given some threshold
$\delta\geq 0$, two intervals, $[i,j]$ in indeterminate string $S$ and $[k,l]$
in indeterminate string $T$, are ($\delta$-)_approximate weak common
intervals_ with _common character set_
$C=\mathcal{C}(S[i,j])\cap\mathcal{C}(T[k,l])$ if the number of _positions_
with no intersection with $C$ is limited by $\delta$, i.e., if we can write

$|\\{x\mid i\leq x\leq j\colon S[x]\cap C=\emptyset\\}|+|\\{y~{}|~{}k\leq
y\leq l\colon T[y]\cap C=\emptyset\\}|\leq\delta$

We shall call these positions _indels_. A set of intervals is a set of
_approximate weak common intervals_ if every pair of members satisfies the
property.

Finally, we extend the property of pairs of _mutually closed_ intervals of [3]
to sets of two or more intervals:

######  Definition 1 (closed set of intervals)

A set $\mathcal{I}$ of intervals is _closed_ if each interval
$[i,j]_{S}\in\mathcal{I}$ has neither immediate left nor immediate right
neighboring position $p$ with $S[p]\cap\mathcal{C}(S^{\prime}[k,l])$ for each
$[k,l]_{S^{\prime}}\in\mathcal{S}\setminus\\{[i,j]_{S}\\}$.

Note that closedness is a non-hereditary property of sets of intervals: a
subset of a closed set of intervals need not be closed itself. Figure 1 shows
an example of a closed set of $1$-approximate weak common intervals.

$\mathbf{S_{1}}$gbpxnpdosazewfvlhuzjrk$\mathbf{S_{2}}$ckfnpwbdxclmagrawxpfz$\mathbf{S_{3}}$dgbapsnabfmwewkjuhcrz
Figure 1: Example of three indeterminate strings $S_{1}$, $S_{2}$, $S_{3}$.
The black underlined intervals form a closed set of $1$-approximate weak
common intervals.

##  4 Identifying Maximal Closed Sets of AWCIs

A set is _maximal_ if it cannot be extended by including further valid
members.

######  Problem 1 (maximal closed sets)

Given $m$ indeterminate strings $S_{1},\ldots,S_{m}$ and indel threshold
$\delta\geq 0$, discover all maximal closed sets of $\delta$-approximate weak
common intervals that have members in at least $q$ out of the $m$
indeterminate strings.

Because closedness of sets of intervals does not extend to subsets, each set
of AWCIs that is identified in the course of the enumeration needs to be
tested individually for this property. For the same reason, the closed pairs
of AWCIs returned by the discovery algorithms described in [3] cannot be used
for the construction of larger sets of AWCIs in multiple sequences, because
non-closed pairs are omitted that could contribute to larger closed sets of
AWCIs.

Since there can be exponentially many sets of AWCIs in a given set of
indeterminate strings, we do not enumerate them explicitly, but instead
construct a compact representation from which they can be extracted. Given a
set of indeterminate strings, $\\{S_{1},\ldots,S_{m}\\}$, we construct a graph
$G(S_{1},\ldots,S_{m})$, where each vertex corresponds to an interval in one
of the indeterminate strings and any two vertices are connected by an
undirected edge if and only if their corresponding intervals form a pair of
AWCIs. Identifying maximal closed sets of AWCIs then reduces to finding within
the graph maximal cliques of vertices whose corresponding AWCIs are closed. In
constructing the graph, we already omit intervals that cannot participate in
sets of AWCIs of minimum size $q$. We test each vertex for possible
participation in a closed set of AWCIs. That is, if the neighbors of a vertex
$v$ are all connected to a vertex $u$ such that the interval $I_{v}$
associated to vertex $v$ is a proper subinterval of the interval $I_{u}$
corresponding to vertex $u$, and if interval $I_{u}$ extends to at most one
position left and/or right of $I_{v}$ that shares characters with neighbors of
$I_{v}$, then vertex $v$ can be discarded. The resulting graph is enriched
with maximal closed cliques of AWCIs.

###  4.1 MACSI – An extension of ACSI

We now describe MACSI, a simple approach for enumerating AWCIs in
indeterminate strings—see Algorithm 1.

Algorithm 1 MACSI enumerates all paris of approximate weak common intervals
that are shared by at least $q$ out of $m$ indeterminate strings.

0: Indeterminate strings $S_{1},\ldots,S_{m}$, indel threshold $\delta$,
quorum $q$.

0: All approximate weak common intervals of $S_{1},\ldots,S_{m}$ with at most
$\delta$ indels that are conserved in at least $q$ strings

1: Construct character set intersection tables $\textsc{Pos}_{xy}$ for all
sequence pairs $\\{x,y\\}\subseteq[1,m]$

2: for $x\leftarrow 1~{}\textbf{to}{}~{}m-1$ do

3: for $i\leftarrow 1~{}\textbf{to}{}~{}|S_{x}|$ do

4: Determine set $J$ of potential right bounds of AWCI pairs with left bound
$i$ in $S_{x}$

5: for each  $j$ in $J$ do

6: Initialize empty list Ints

7: for $y\leftarrow x+1~{}\textbf{to}{}~{}m$ do

8: construct sorted set $P_{y}$ of positions
$\textsc{Pos}_{xy}[i]\cup\cdots\cup\textsc{Pos}_{xy}[i+\delta]$

9: $p_{\mathit{prev}}\leftarrow-1$

10: for each  $p$ in $P_{y}$ do

11: Append to list Ints all intervals $[k,l]$ in $S_{y}$,
$p_{\mathit{prev}}<k\leq p\leq l$, that are approximate weak common intervals
with $[i,j]$ in $S_{x}$

12: $p_{\mathit{prev}}\leftarrow p$

13: end for

14: end for

15: if Ints contains intervals from at least $q-y+1$ sequences then

16: report
$\\{([i,j]_{S_{x}},[k,l]_{S_{y}})~{}|~{}[k,l]_{S_{y}}\in\textsc{Ints}\\}$

17: end if

18: end for

19: end for

20: end for

MACSI is a simple extension of Algorithm ACSI described by Doerr  _et al._
[3]; however, unlike ACSI, it enumerates any pairs of AWCIs rather than only
closed ones.

Instead of repeatedly determining intersections of character sets across
positions of distinct indeterminate strings, the positions of intersecting
character sets are identified and stored at the beginning. Given two
indeterminate strings $S_{x}$ and $S_{y}$, each position $j$, $1\leq
j\leq|S_{y}|$ in intersection table $\textsc{Pos}_{xy}$ is associated with a
sorted list of positions in $S_{y}$ that have a nonempty intersection with
character set $S_{x}[j]$.

Given a set of sequences $S_{1},\ldots,S_{m}$, each of its members acts once
as _reference_ (lines 2-19), over which the algorithm iterates, thereby fixing
a left bound $i=1,\ldots,|S_{x}|$ (lines 3-20). The algorithm then determines
candidates for right bounds $J\subseteq\\{i+1,\ldots,|S_{x}|\\}$, where each
candidate interval $[i,j]_{S_{x}}$, $j\in J$, is potentially part of a set of
AWCIs conserved in at least $q$ out of the $m$ sequences. Once the interval
$[i,j]$ in reference string $S_{x}$ is fixed, the algorithm iterates over the
remaining $m-x$ sequences to enumerate all such candidate intervals (lines
7-14).

The construction of set $J$ representing potential right interval bounds in
$S_{x}$ (line 4) is easily done in time linear in $|S_{x}|$. For each fixed
interval $[i,j]$ in $S_{x}$, the enumeration of all intervals in $S_{y}$ that
are AWCIs with $[i,j]_{S_{x}}$ can be achieved in
$\mathcal{O}(|S_{y}|\cdot\|S_{y}\|)$ amortized running time. The algorithm
first collects and sorts all positions $P_{y}\in S_{y}$ around which such
intervals may lie (line 8). For each such position $p\in P_{y}$, all intervals
with left bound larger than the previously processed position are identified
(lines 10-13). Testing for the approximate common interval property of an
interval pair $[i,j]_{S_{x}},[k,l]_{S_{y}}$, $k\leq p\leq l$, can be done in
$\mathcal{O}(1)$ by keeping track of observed indices of $[i,j]_{S_{x}}$,

$C_{kl}=\\{i^{\prime}\in\textsc{Pos}_{yx}[k]\cup\textsc{Pos}_{yx}[l]~{}|~{}i\leq
i^{\prime}\leq j\\}$

and of the number of positions $d_{kl}$ in $[k,l]_{S_{y}}$ whose character
sets do not intersect with $[i,j]_{S_{x}}$. Then $[i,j]_{S_{x}},[k,l]_{S_{y}}$
are AWCIs iff $j-i+1-|C_{kl}|+d_{kl}\leq\delta$ holds and the enumeration of
all intervals that are AWCIs with $[i,j]_{S_{x}}$ can be done in
$\mathcal{O}(|S_{y}|\cdot\|S_{y}\|)$ time.

The overall running time of Algorithm 1 is in $\mathcal{O}(m^{2}\cdot
n^{2}\cdot N^{2})$, with $n=\max_{x=1..m}|S_{x}|$ and
$N=\max_{x=1..m}\|S_{x}\|$. This is close to optimal, since the number of
pairwise AWCIs can be as large as the total number of interval pairs, itself
in $\mathcal{O}(m^{2}\cdot n^{4})$. Hereby we make the reasonable assumption
that the cardinality $\|S\|$ of an inteterminate string $S$ is in the order of
$a\cdot|S|$ with $a\ll\Sigma$.

##  5 Speeding Up the Enumeration of AWCIs for Synteny Analysis

We now study in more detail the construction of set $J$ (line 4), which
contains potential right interval bounds in $S_{x}$, and describe an approach
that significantly improves the running time by taking advantage of two
characteristics of genomes. First, there exist many characters shared only
between a subset of all sequences; and second, each sequence in the dataset is
associated with a single genome representing a concatenation of one or more
chromosomes or contigs, yet AWCIs are prohibited from spanning more than one.
As a result, our sequences exhibit a non-negligible number of indels that can
be trivially identified and used to preempt or terminate unsuccessful searches
of sets of AWCIs.

Recall that we are interested in identifying sets of AWCIs with members in a
quorum of at least $q\leq m$ strings. Sets of AWCIs enforce transitivity, but
the AWCI property itself is not transitive. Therefore, even if we identify a
set of intervals $I_{x+1..m}$ with members in at least $q-1$ strings and which
form AWCI pairs with a given interval $[i,j]_{S_{x}}$, we cannot guarantee
that the intervals $I_{x+1..m}$ form AWCI pairs among themselves. Thus, for
each interval $[i,j]_{S_{x}}$, we want to assess whether it could participate
in at least one set of AWCIs that satisfies the quorum parameter prior to
constructing set $I_{x+1..m}$.

###### Observation 1

Given sequences $S_{1},\ldots,S_{m}$, an indel threshold $\delta$, and a
quorum $q$, an interval $[i,j]$ in $S_{x}$, for $1\leq x\leq m$, is _not_ a
candidate for an AWCI pair if we have
$|\\{|\mathcal{C}(S_{x}[i,j])\setminus\mathcal{C}(S_{y})|\leq\delta~{}|~{}y\in[1,m]\setminus
x\\}|<q-1$.

To exploit this observation, we construct two types of tables,
$\textsc{Ridge}_{xy}^{c}$ and $\textsc{Ridge}_{xy}^{t}$, for
$\\{x,y\\}\subseteq[1,m]$, where superscripts $c$ and $t$ stand for _cis_ and
_trans_ , respectively. The first, $\textsc{Ridge}^{c}_{xy}$, holds
information within intersecton table $\textsc{Pos}_{xy}$ to decide whether
$[i,j]_{S_{x}}$ putatively participates in an AWCI pair, whereas
$\textsc{Ridge}^{t}_{xy}$ holds information from $\textsc{Ridge}_{yx}^{c}$,
allowing to test whether interval $[i,j]_{S_{x}}$ is associated to intervals
in $S_{y}$, with whom it potentially forms AWCI pairs. Each position $0\leq
p<|S_{x}|$ in table $\textsc{Ridge}_{xy}^{c}$ represents a count of the number
of positions $0\leq p^{\prime}\leq p$ for which $S_{x}[p^{\prime}]$ contains
no characters of $\mathcal{C}(S_{y})$. Thus interval $[i,j]$ in string $S_{x}$
cannot participate in any AWCI pair if
$\textsc{Ridge}_{xy}^{c}[j]-\textsc{Ridge}_{xy}^{c}[i-1]>\delta$.

We now construct a similar test for trans pairs of our reference interval
$[i,j]_{S_{x}}$ with the help of tables $\textsc{Ridge}_{xy}^{t}$. We say that
two positions $i,j$ of string $S_{x}$, for $0\leq i<j\leq|S_{x}|$, reside on
the same _ridge_ iff
$\textsc{Ridge}_{xy}^{c}[j]-\textsc{Ridge}_{xy}^{c}[i-1]\leq\delta$ holds. For
a trans string $S_{y}$, $x\neq y$, observe that the characters of
$\mathcal{C}(S_{x}[i,j])$ can occur in many different positions in $S_{y}$ and
thus can be associated with different intervals that form an AWCI pair with
$[i,j]_{S_{x}}$. For each candidate interval $[k,l]$ in $S_{y}$ we must have
$\textsc{Ridge}_{yx}^{c}[k]-\textsc{Ridge}_{yx}^{c}[l-1]\leq\delta$. To test
this property quickly, we construct binary vectors to identify ridges in
$S_{y}$ containing characters of $\mathcal{C}(S_{x}[i,j])$; an example is
shown in Figure 2.

Figure 2: The data structure $\textsc{Ridge}_{xy}^{t}$; lines indicate ridges,
i.e., 1s across binary vectors $\textsc{Ridge}_{xy}^{t}[i^{\prime}]$ for some
$i^{\prime}\in[0,|S_{x}|]$.

Each index in any vector associated with positions $i..j$ in $S_{x}$
corresponds to a distinct ridge in $S_{y}$. In order to test whether any two
positions $i^{\prime}$, $j^{\prime}$, $i\leq i^{\prime}<j^{\prime}\leq j$,
contain characters that are located on common ridges in $S_{y}$, it then
suffices to perform a simple bitwise AND operation between its associated
binary vectors. Moreover, using a bitwise AND of all vectors associated with
positions in $i..j$ in $S_{y}$ we can decide whether there exists at least one
ridge in $S_{y}$ that contains characters from all character sets
$S_{x}[i^{\prime}]$, $i\leq i^{\prime}\leq j$.

Data structure $\textsc{Ridge}_{xy}^{t}$ is of size
$\mathcal{O}((\delta+1)\cdot|S_{x}|\cdot\|S_{y}\|)$ and its use increases the
overall running time of MACSI by a factor of $(\delta+1)$. In practice, its
size is far from its asymptotic bound and MACSI gains a significant speed-up.
To minimize the size of these binary vectors in $\textsc{Ridge}_{xy}^{t}$, we
relax our requirement that each ridge be assigned a unique index across all
binary vectors. We can do this without damaging the consistency of our data
structure, since we do not compare arbitrary binary vectors against each
other, but only those whose positions in $S_{x}$ are at most $\delta$ trivial
indels afar. Therefore we reuse indices of ridges in the binary vectors.
Specifically, when constructing the data structure $\textsc{Ridge}_{xy}^{t}$
from left to right, at a given position $j$ we assign a newly observed ridge
$\textsc{Ridge}_{yx}^{c}[k]$ (for some $k\in\textsc{Pos}_{xy}[j]$) to a
previously used index in the binary vector if the ridge associated to this
index has not been observed in any of the previous positions in $S_{x}$ that
are no more than $\delta$ indels afar.

In Algorithm 2,

Algorithm 2 Testing a position $j$ in string $S_{x}$ for the existence of some
interval $[i^{\prime},j^{\prime}]_{S_{x}}$, $i^{\prime}\leq i\leq j\leq
j^{\prime}$—, i.e., AWCIs with intervals in at least $q-1$ other indeterminate
strings.

0: start $i$ and current position $j$ in string $S_{x}$, tables
$\textsc{Ridge}_{xy}^{c}$ and $\textsc{Ridge}_{xy}^{t}$, temporary data
structures $\textsc{Active}^{t}$ and $\Delta^{t}$, indel treshold $\delta$,
quorum $q$, reference $x$.

0: Returns True if $S_{x}[i,j]$ is a candidate for an AWCI set in at least $q$
strings with less than $\delta$ indels and False otherwise.

1: $\mathit{candidates}\leftarrow 0$

2: for $y\textbf{ in }[1,k]\setminus x$ do

// Determine the location of currently masked (trans-) ridges in
$\Delta^{t}[y]$

3:
$d\leftarrow\max(0,\delta+1+\textsc{Ridge}_{xy}^{c}[i]-\textsc{Ridge}_{xy}^{c}[j])$

// Ridges that have been observed before but are not present at this position

4:
$a_{\mathit{off}}\leftarrow\textsc{Active}^{t}[y]\land\lnot(\textsc{Ridge}_{xy}^{t}[j]\lor\Delta^{t}[y][d])$

// Ridges that have not been observed before but are present in
$\textsc{Ridge}_{xy}^{t}[j]$

5:
$a_{\mathit{new}}\leftarrow\textsc{Ridge}_{xy}^{t}[j]\land\lnot(\textsc{Active}^{t}[y]\lor\Delta^{t}[y][d])$

// Update $\Delta^{t}[y]$ for all ridges in $a_{\mathit{off}}$ and
$a_{\mathit{new}}$

6: for $d^{\prime}\leftarrow 0~{}\textbf{to}{}~{}d$ do

7: $tmp\leftarrow\Delta^{t}[y][d^{\prime}]\land a_{\mathit{off}}$

8: $\Delta^{t}[y][d^{\prime}]\leftarrow\Delta^{t}[y][d^{\prime}]\lor
a_{\mathit{off}}$

// Carry ridges of $a_{\mathit{off}}$ already marked in
$\Delta^{t}[y][d^{\prime}]$ over to $\Delta^{t}[y][d^{\prime}+1]$

9: $a_{\mathit{off}}\leftarrow tmp$

10: if $d^{\prime}<j-i$ then

11: $\Delta^{t}[y][d^{\prime}]\leftarrow\Delta^{t}[y][d^{\prime}]\lor
a_{\mathit{new}}$

12: end if

13: end for

// Add ridges of current position to $\textsc{Active}^{t}[y]$

14:
$\textsc{Active}^{t}[y]\leftarrow\textsc{Active}^{t}[y]\lor\textsc{Ridge}_{xy}^{t}[j]$

// Test if $S_{x}[i,j]$ has at least one ridge $S_{y}$ with less than $\delta$
indels

15: if $\textsc{Active}^{t}[y]\land\lnot\Delta^{t}[y][d]\neq 0$ then

16: $\mathit{candidates}\leftarrow\mathit{candidates}+1$

17: end if

18: end for

19: return $\mathit{candidates}\geq q-1$

we develop this approach to allow for a total of $\delta$ indels in
$S_{x}[i,j]$ or any ridge in $S_{y}$, $[1,k]\setminus x$. In doing so, we make
use of two data structures, $\textsc{Active}^{t}$ and $\Delta^{t}$, which are
updated and queried in each iteration for testing interval $i..j$ in $S_{x}$
with fixed left bound $i$ and increasing right bound $j$.
$\textsc{Active}^{t}[y]$ is a binary vector that captures all ridges in
$S_{y}$ containing characters from $\mathcal{C}(S_{x}[i,j])$. $\Delta^{t}[y]$
is a data structure made of $\delta+1$ vectors that act as counters for indels
in ridges of $S_{y}$—$\Delta^{t}[y][d^{\prime}]$ records all ridges in $S_{y}$
that accumulated $0<d^{\prime}\leq\delta$ indels. Clearly, if
$d^{*}=\textsc{Ridge}_{xy}^{c}[j]-\textsc{Ridge}_{xy}^{c}[i]$ trivial indels
have been observed in $S_{x}[i,j]$ (line 3), then only $\delta-d^{*}$ indels
are allowed to occur in any ridge in $S_{y}$. Algorithm 2 uses the
$(\delta+1-d^{*})^{\textit{th}}$ vector in $\Delta^{t}[y]$ as a _mask_ to
ignore all ridges that exceed $\delta$ indels in subsequent calculations.

At each iteration $j$, ridges of $S_{y}$ can be classified in three
categories: (i) ridges in $\textsc{Ridge}^{t}_{xy}[j]$ previously observed in
any of the $i..j-1$ iterations; (ii) previously observed ridges not in
$\textsc{Ridge}^{t}_{xy}[j]$; (iii) new ridges of
$\textsc{Ridge}^{t}_{xy}[j]$. The last two categories require an update of
$\Delta^{t}$ and hence are assigned to vectors $a_{\mathit{off}}$ and
$a_{\mathit{new}}$, respectively (lines 4, 5). The counter for indels of all
ridges of $a_{\mathit{off}}$ gets increased by one, whereas for
$a_{\mathit{new}}$ the counter gets increased by $\min(j-i,\delta)$; both
updates occur in the for loop on lines 6-13. After updating
$\textsc{Active}^{t}[y]$ vector (line 14), Algorithm 2 can test whether there
exists at least one ridge in $S_{y}$ that satisfies the indel constraint with
$S_{x}[i,j]$ (lines 15-17). Algorithm 2 returns True only if $q-1$ strings
fulfill the indel constraint; otherwise, the current position $j$ in reference
string $S_{x}$ gives rise to set $J=[i,j]$ (Algorithm 1, line 4).

Further improvements in speed are made by individually refining the set of
characters $J$ for each position of set $P_{y}$, which is constructed in line
8 in Algorithm 1. In doing so, we improve on a technique already described in
[3]. Note that $P_{y}$ is not dependent on the candidate set of right bounds
$J$ of intervals in $S_{x}$, but rather on the left candidate-interval bound
$i$ and therefore can be computed prior to iterating through set $J$. For each
position $p\in P_{y}$ of indeterminate string $S_{y}$, we collect indices
$J_{p}\subseteq J$ in all sets
$\textsc{Pos}_{yx}[k^{*}]..\textsc{Pos}_{yx}[l^{*}]$ such that $k^{*}$ and
$l^{*}$, $k^{*}\leq p\leq l^{*}$, are at most $\delta$ trivial indels from $p$
afar, i.e.
$J_{p}=(\textsc{Pos}_{yx}[k^{*}]\cup\cdots\cup\textsc{Pos}_{yx}[l^{*}])\cap
J$. Thus we use Algorithm 2 on each position $k^{\prime}$ of indeterminate
string $S_{y}$ that is visited while walking from $p$ towards $k^{*}$. In
doing so, we test for the possible existence of AWCI sets that have members in
at least $q-2$ strings other than $S_{y}$ and $S_{x}$ and contain some
interval $[k,l]_{S_{y}}$, $k\leq k^{\prime}\leq p\leq l$. We stop progressing
towards $k^{*}$ if Algorithm 2 returns False. This process also applies to
positions $l^{\prime}$, $p\leq l^{\prime}\leq l^{*}$.

Just as described in [3], sets $J_{p}$, $p\in P_{y}$, can then be used to
refine the candidate set $J$ of right interval bounds of $S_{x}$. Unlike in
[3], however, the refinement is performed across multiple sequences. For any
set of positions $P_{y}$, we identify the rightmost interval bound found in
the neighborhood of any $p\in P_{y}$ that satisfies the indel treshold in
$S_{x}$, i.e.,

$j^{\star}_{y}:=\max_{p\in
P_{y}}\\{j_{p}~{}|~{}|J_{p}\cap[i,j_{p}]|\leq\delta\\}$

Then the new candidate set of right interval bounds of the next iteration is

$J=\big{[}i,\max_{\begin{subarray}{c}Z\subseteq[1,m]\setminus x\\\ |Z|\geq
q-1\end{subarray}}\big{(}\min_{z\in Z}j^{\star}_{z}\big{)}\big{]}$

We repeat this process to refine each of the sets $J_{p}$. In practice, fewer
than 3 iterations suffice to converge to stable sets $J$ and $J_{p}$.

##  6 Results and Discussion

We demonstrate the applicability of MACSI in comparing genomic sequences on a
dataset of 93 bacterial genomes that was also used in evaluating the
algorithms for discovering pairs of AWCIs presented in Doerr  _et al._ [3] and
that were originally obtained from Ciccarelli  _et al._ [1]. Doerr  _et al._
transformed the genomic sequences into indeterminate strings, by using
annotions of protein-coding genes and by computing pairwise BLASTP hits
between genes. The resulting indeterminate strings have an average length of
$2726$ elements and each pair of indeterminate strings contains $6499$
positions with intersecting character sets.

MACSI is implemented in Python; our implementation allows for parallel
computation, but for the purpose of assessing the runtime of our method, all
computations were performed in a single thread. In addition to indel parameter
$\delta$ and quorum, our implementation allows to set a minimum size for
reported AWCIs, where size is measured between pairs of AWCIs by the number of
positions that contain shared characters. Throughout our experiments, we set
the minimum size to $10$.

(a) . (b) . (c) .

Figure 3: (a) Runtimes of MACSI and (b) bit vector sizes (per position) in
data structure $\textsc{Ridge}_{xy}^{t}$ on genome sets of different sizes
with indel thresholds $\delta\in\\{0,2\\}$; (c) Runtimes of MACSI with indel
threshold $\delta=0$ on sampled sets of 20 genomes with varying quorum
settings.

Experiments were run on a machine with sixty-four 2.3 GHz cores. To assess the
running time of our implementation we sampled genome sets of different sizes
from the pool of 93 bacterial genomes and then ran MACSI using different
parameter settings. The sampling was repeated 10-fold in all experiments.
Figure 3 shows the results. Since our speed-up relies on the frequency of
indels, the indel threshold $\delta$ has the most dramatic effect on the
running time of the algorithm; in comparison, lowering the quorum has only a
mild effect. Comparing 40 full bacterial genomes took on average just around
3.5 hours of computation.

The number of reported AWCI pairs as well as the data structure
$\textsc{Ridge}_{xy}^{t}$ have the most decisive effect on the space
requirements of our algorithm. Because the former is less prone to algorithmic
improvements, we monitored only the latter in our experiments. Figure 3 (b)
illustrates the sizes of bit vectors in $\textsc{Ridge}_{xy}^{t}$ in our
experiments and shows that fewer than 500 bits were needed in runs with
$\delta=0$ and fewer thatn 1400 bits with $\delta=2$.

##  7 Conclusion

We presented a fast algorithm for the discovery of approximate weak common
interval in multiple indeterminate strings. Our main contributions are data
structures that allow for efficient filtering of intervals that cannot be
members of AWCI sets of a given minimal size. We exploited genomic properties,
in particular the fact that our indeterminate strings will exhibit many
trivial indels, yielding easily identifiable boundaries for AWCIs. How to
speed up the general algorithm in the absence of trivial indels remains an
open question.

## Competing Interests

The authors declare that they have no competing interests.

## Author’s Contributions

BM initiated and directed the research project, DD designed and implemented
the herein presented algorithms and ran the experiments; both authors wrote
the manuscript, and read and approved its final version.

## References

  * [1] F. D. Ciccarelli, T. Doerks, C. von Mering, C. J. Creevey, B. Snel, and P. Bork.  Toward automatic reconstruction of a highly resolved tree of life.  Science, 311(5765):1283–1287, 2006.
  * [2] G. Didier, T. Schmidt, J. Stoye, and D. Tsur.  Character sets of strings.  J. Discr. Alg., 5(2):330–340, 2007.
  * [3] D. Doerr, J. Stoye, S. Böcker, and K. Jahn.  Identifying gene clusters by discovering common intervals in indeterminate strings.  In Proc. 12th RECOMB Satellite Workshop on Comparative Genomics RECOMB-CG’14, volume 15 (Suppl. 6) of BMC Genomics, page S2, 2014.
  * [4] C. Ghiurcuta and B. Moret.  Evaluating synteny for improved comparative studies.  In Proc. 22nd Symp. on Intelligent Systems for Mol. Bio. ISMB’14, volume 30(12) of Bioinformatics, pages i9–18, 2014.
  * [5] T. Uno and M. Yagiura.  Fast algorithms to enumerate all common intervals of two permutations.  Algorithmica, 26(2):290–309, 2000.
  * [6] M. Visnovská, T. Vinař, and B. Brejová.  Dna sequence segmentation based on local similarity.  In Proc. 13th Conf. on Info. Technologies—Applications and Theory ITAT’13, pages 36–43, 2013.
  * [7] X. Zhang and B. Moret.  Refining regulatory networks through phylogenetic transfer of information.  IEEE/ACM Trans. on Computational Biology and Bioinformatics, 9(4):1032–1045, 2012.