The main goal of this work is to propose a pipeline capable to analyze a data collection of temporally quali ed mutation pro les. A front-end and a post-processor are implemented. The front-end basically transforms the input data retrieved from medical databases into a set of facts. The post-processor allows the user to visualize the computation results as graphs. The reasoning core is based on Answer Set Programming that, allowing to deal with NP or 2P properties, is capable of deducing complex info from data. The system is modular: the reasoning tool can be replaced by any logic programming tool (e.g., Prolog, ILP) for di erent kinds of data analysis.
Cancer is a game of mutation accumulation in the genome of an organism, that is the set of its genes. One of the most delicate aspects of this game is the fact that, while data is produced by today's technology at an unprecedented rate it is, however, extremely noisy data. More precisely, even though any mutation might be triggering signi cant biological processes involved in cancer development, only some of the mutations are classi ed as \drivers" by experts in the eld. The remaining|large|majority of mutations must be classi ed as \passengers", as they are probably accumulating only as a consequence of the (devastating) side-e ects of a more basic biological mechanisms already at work.
Tackling a classi cation problem for the collection of detected mutations in a tissue requires a multi-faceted approach. In order to build reliable lters suitable to de ne the blueprint of the disease under study it is essential to integrate di erent kinds of knowledge with brute force analysis.
The main goal of this work is to propose a pipeline capable to analyze a data collection of temporally quali ed mutation pro les and synthesize a dependency graph. The information inferred by cross comparisons performed \locally" on input data, combined with information coming from external sources, allows to deduce the output graph that should represent a set of putative dependencies.
The idea is to prove the potential of the employed programming technologies in mining data that, when analyzed over realistically long sequences of time steps, rapidly generates (computationally) heavy collections of putative evolution paths.
ASP allows to encode compactly and elegantly problems that belongs to NP or even to 2P . Encodings allows great exibility. In particular, the already implemented strategies are parametrically designed and use thresholds xing conditions on the addition and removal of edges. Moreover, as it is always the case when dealing with biological data, many auxiliary problems|mostly related with I/O|must be faced. The data processing pipeline described here is interfaced with specialized languages and le formats for the eld.
The paper is organized as follows: some related work are discussed in Sec. 2, the main ideas of our approach are presented in Sec. 3. A discussion on the data set used an the applications we have developed for extracting/converting data are presented in Sec. 4, and a brief analysis of the results found is shown in Sec. 5. Finally, some conclusions are drawn in Sec. 6. 2
Cancer data processing is, obviously, a very active and lively eld in many scienti c environments. We just report here on some work using techniques stemming from the logic programming community. We assume the reader has basic knowledge of logic and logic programming.
Inductive Logic Programming (ILP) has been succesfully used to analyze big
data and, in particular, to analyze cancer databases. ILP allows to implement
a \structured" form of machine learning, by inferring a rst-order theory that
\explains" a set of input ground facts (de ning extensionally, possibly partially,
some predicates). Let us brie y recall the main ideas behind ILP (see, e.g., [
Assume to have a set of positive observations O+ and a set of negative observations O . We can assume O+ and O as two disjoint sets of ground atoms. Assume also to have a (possibly empty) background theory (a set of clauses) P , the general inductive problem is the one of nding a set of hypotheses H (clauses) such that P ^ H j= O+ and such that for all ` 2 O it holds that P ^ H 6j= `.
Several strategies have been investigated and implemented by ILP systems to identify sets H made by clauses as much general as possible. Generality in this context can be de ned in terms of logical entailment: if clauses r1 and r2 are such that P ^ r1 j= r2, then r1 is more general of r2 w.r.t. P .1 Of course, this implies that if P ^ r2 j= B (for some B O+), then P ^ r1 j= B. Therefore, r1 is preferable to r2 in the set of hypotheses H (unless it introduces errors, namely, there is ` 2 O s.t. P ^ r1 j= ` and P ^ r2 6j= `).
In presence of large amounts of data that might be a ected from errors since the beginning (e.g., medical data), one might accept a set H that might be not complete, namely it \explains" only a subset G (good) of O+: P ^ H j= G, and that might be not correct, namely such that P ^ H j= ` for ` in a subset W (wrong) of O . In this case, of course, the ILP algorithm aims at maximizing the size of G and at minimizing the size of W .
An extension of ILP is PILP (probablistic inductive logic programming) [
Inductive Logic Programming has been used to analyze cancer data. For
instance in [
Progol have been also used in [
They use data from 438 cases. The classi cation is based on the observation of the lymph node metastatic status and tumour di erentiation status. They use a set of available lab data (e.g., CEA, CA19-9, Glucose, Elastase I, Serum Amylase, . . . ). After having identi ed the most promizing features using standard feature selection criteria data is divided in three groups (low, normal, and high), as usually made by domain experts. Each patient record in the database represents an abnormality of a patient, split into positive and negative. Data is translated into a set of facts, and Progol is executed in order to rank the rules by their capability of explaining the database.
In [
Answer Set Programming has been used in many Bioinformatics applications.
For an overview, see [
There are other approaches to biological data analysis based on ASP. For
instance, in [
The general idea of our approach is simple and can be summarized as an attempt to infer a graph whose nodes are mutations (in genes) and whose arcs represent the causality relation between mutations, that is their most plausible (temporal) sequence. Data consists of variant allele frequencies and thresholds are used to accept/discard a temporal dependency.
The initial steps consist in a deterministic analysis followed by a blind search
entirely left to the inferential engine. A pre-processing phase is necessary in
order to correctly use data. A le MAF (Mutation Annotation Format) [
In this section we will see our ASP code and how it works over some example instances. Once ltered, input data is stored by two predicates, de ned extensionally by facts. Each patient can be involved in more than one experiment; the ternary predicate sample relates the patient ID with its experiment IDs and their times (given the experiments ordering), while the ternary predicate val reports each instance of mutated gene in any experiment and its value:
In the rst ve rows we de ne the values of the samples of three patients pat-A, pat-B and pat-C, with the (temporal) order of samples. Subsequently, starting at row six, we indicate mutated genes values. For the sake of simplicity we named values a,b,c . . . .
A number of \domain" predicates are computed from the input data (needed to limit the grounding), as follows: 1 gene(G) :- val(_,G,_). 2 patient(P) :- sample(P,_,_). 3 value(V) :- val(_,_,V). 4 time(0..T) :- sample(_,_,T), not sample(_,_,T+1).
The predicate gene(G) forces value G to appear as second parameter of val and, therefore, will hold for each gene uploaded in our samples. Analogous considerations for patient(P) and value(V). time(T), instead, is instanced with a list of values that goes from 0 to a maximum temporal value contained in the predicates sample(P,E,T).
The following predicates allow to retrieve information from the data. Here we report some examples. Being ASP capable to deal with the whole NP class, and even the 2P class if disjunctive heads are used, one can capture/infer a wide family of properties. We will be back on this issue in Sec. 6. 1 variationConsec(P,G,T2) :2 sample(P,E1,T1),sample(P,E2,T2),T2=T1+1, 3 val(E1,G,V1),val(E2,G,V2),variation(V1,V2). 4 variation(A,B) :5 value(A),value(B),|A-B|>10. 6 precedes(P,G1,G2) :7 firstVariation(P,G1,T1),firstVariation(P,G2,T2),T1<T2. 8 firstVariation(P,G,T) :9 variationConsec(P,G,T), not varPrec(P,G,T). 10 varPrec(P,G,T) :11 variationConsec(P,G,T1), variationConsec(P,G,T),T1<T. 12 precedeKTimes(G1,G2,K) :13 gene(G1),gene(G2),K=#count{patient(P):precedes(P,G1,G2)}. 14 reachable(G1,G2) :15 maybePrecede(G1,G2). 16 reachable(G1,G2) :17 maybePrecede(G1,G3), reachable(G3,G2). 18 maybePrecede(G1,G2) :19 precedeKTimes(G1,G2,K1),precedeKTimes(G2,G1,K2),K1>K2, 20 not reachable(G2,G1). 21 #show maybePrecede/2.
variationConsec(P,G,T2) holds whenever for patient P there are two consecutive samples (the second one with time T2) with the same gene G with two values V1 e V2 that satisfy variation, namely their di erence in absolute value is greater than 10. This value is the variation sensibility, in this simple example we choose a low parameter in order to take always the consecutive variation of the mutated gene. In our example we can infer that: variationConsec(pat-A,a,2). variationConsec(pat-A,b,4). variationConsec(pat-B,d,2). variationConsec(pat-B,a,4). variationConsec(pat-C,b,2). variationConsec(pat-C,z,4).
variationConsec(pat-A,c,3). variationConsec(pat-A,c,4). variationConsec(pat-B,a,3). variationConsec(pat-B,b,4). variationConsec(pat-C,a,3).
variationConsec(pat-C,d,5).
Before de ning the main predicate of this example, namely, the predicate precedes(P,G1,G2), let us see the rule varPrec(P,G,T) at row 10: it holds if a gene G of a certain patient P has some variations in the sample at time T and also in some previous samples|and therefore it is not the rst variation detected. This predicate is used in the body of firstVariation(P,G,T): the predicate inferring the minimum time T for which we have a variation of gene G in P.
At this point we can illustrate precedes(P,G1,G2) in row 6: it is true if, for a given patient, the rst variation of G1 happens before the rst variation of G2. Notice that it is not necessary to force that the two genes are di erent as it is implicit in the fact that the time of their rst variation is di erent. So from the previous example we can infer: precedes(pat-A,a,c). precedes(pat-B,d,a). precedes(pat-C,b,a). precedes(pat-C,b,d).
precedes(pat-A,a,b). precedes(pat-B,d,b). precedes(pat-C,b,z). precedes(pat-C,a,d).
precedes(pat-A,c,b). precedes(pat-B,a,b). precedes(pat-C,a,z). precedes(pat-C,z,d).
Predicate precedeKTimes(G1,G2,K) stores in K the number of patients for which the gene G1 precedes the gene G2, while maybePrecede(G1,G2) is inferred if G1 precedes G2 more frequently than G2 precedes G1. We report below the inferred instances of precedeKTimes in which the guard K is greater than zero (namely there is at least one case in which G1 precedes G2.
precedeKTimes(b,a,1). precedeKTimes(d,a,1). precedeKTimes(c,b,1). precedeKTimes(d,b,1). precedeKTimes(a,d,1). precedeKTimes(b,d,1). precedeKTimes(a,z,1). precedeKTimes(b,z,1). precedeKTimes(a,b,2). precedeKTimes(a,c,1). precedeKTimes(z,d,1).
For a and b, there are 2 patients that exhibit variation of a before b, while just one exhibits the opposite variation. Hence we deduce maybePrecede(a,b), while in case of a and d we don't infer any most likely order.
The nal result in our example is: maybePrecede(a,b). maybePrecede(b,z).
maybePrecede(a,c). maybePrecede(c,b).
maybePrecede(a,z). maybePrecede(z,d).
In order to guarantee that the direct graph resulting from maybePrecede is acyclic, we check that the number of times the variation of the rst gene precedes the second is strictly greater than the opposite order and we also impose the absence of the predicate reachable that states that already exists a path between two genes in the graph inferred.
The visualisation of the graph obtained is shown in Figure 2. The table
containing the data (the variant allele frequency value has been omitted for
simplicity) where the pairs highlighted are the variations revealed. In Figure 3
(right) it is shown the result automatically produced by Graphviz, after the
Clingo output has been converted in the dot edges format Figure 3 (left).
We tested our ASP program over the|real and publicly available|data set
Breast cancer patient xenografts [
The open source platform cBioPortal provides data from many cancer genomics data-sets (by now there are 150 of them) allowing their download, analysis, and visualisation. It was originally developed at Memorial Sloan Kettering Cancer Center (MSK) and now its software is available under an open source license via GitHub at https://github.com/cBioPortal. The portal is designed to bring cancer researchers near to the complexities of human genome data, allowing a rapid, intuitive, and precise visualization of expression pro les as well as clinical attributes from large-scale cancer genomics projects. A brief tutorial of the tools provided by cBioPortal is available on the web site; further documentation is available at https://cbioportal.readthedocs.io/en/latest/README.html. From the website it is possible to query and/or download data from one or more studies satisfying given criteria (for instance, studies related to speci c gene's values) or it is possible to download the entire archive. From the home page it is also possible to query databases speci ed within studies, select the genetic pro le, the single patient or a group of them and then specify collections of genes of interest. Users can submit even more speci c queries by using the Onco Query Language (OQL), for which a brief guide is available on the website. The same information can be obtained through web API and library in R and Matlab. However, the faster way to visualize all genes from a study, consists in directly downloading the whole study database from https://github.com/cBioPortal/datahub/tree/master/public, containing all complete archives available from the portal. Once unzipped the archive, a series of le is at hand, among which data mutations extended written in MAF that will be our input.
elds used by our program, redirecting the interested reader to the o cial page https://wiki.nci.nih.gov/display/TCGA/Mutation+Annotation+Format+(MAF) +Speci cation for a full- edged description of the format.
Mutations are detected by alignment and comparison between DNA sequences
relative to biopsies and the human genome reference (without mutations) available
at NCBI (National Center for Biotechnology Information) [
The format of a MAF le is tab-delimited columns. The rst row contains the version used, the second one the elds headers, the rst 34 of which are mandatory.
Variant Allele Frequency (VAF) is given by reads depth on cancer's sample, referring to an altered allele t_alt_count, with respect to total reads t_depht. We will use precisely this ratio to describe mutation level, leaving the possibility of more elaborated/precise measurements, open in our code. 4.2
Data from the above mentioned study Breast cancer patient xenografts [
The study shows the clonal dynamics of initial engraftment and subsequent serial propagation of primary and metastatic human breast cancers in immunode cient mice. One of the results showed by the study is that similar clonal expansion patterns can emerge in independent grafts of the same starting tumour population, indicating that genomic aberrations can be reproducible determinants of evolutionary trajectories. So we can see that measurement of genomically de ned clonal population dynamics will be highly informative for functional studies using patient-derived breast cancer xenoengraftment.
The study, employing single-cell sequencing, was able to analyze and provide
data on three di erent tissues: subcutaneous, subrenal capsule, and mammary
fat pad. In [
In order to get the complete archives stored in cBioPortal, automatically analyse data, and visualize the graph of variations found, we built a Java GUI called Mutations Evolution Analyzer (Figure 4) that contains a menu with four items: Download, Create, Run, and Help.
The item Download sends a query to the web interface of the Cancer Genomic Data Server (CGDS) and, if available, returns an updated list of the studies, allowing the download over the site https://github.com/cBioPortal/datahub/tree/master/public of one or more of them in a user speci ed directory|see Figure 5.
Once the archive has been unzipped, we nd data mutations extended le. By clicking on Create and then Selected fields the program allows to get, from a MAF le, the values of the selected elds (see Figure 6). To speed up selection, there is a checkbox called Select standard eld, that automatically selects the elds used for our analysis: Hugo Symbol, T umor Sample Barcode, t alt count, and t depth. In particular, for Breast cancer patient xenografts, the identi cation eld is obtained by the union of the patient ID, the type of the sample (tumor or xenograft), the position in temporal order of the sample, the genotype studied, and whether it is relative to the whole genome or if it is targeted. For instance, sample SA429X1A-targeted is relative to patient SA429, it has been obtained over a xenograft (X means xenograft), is the rst sample of this type (1)|as samples have a sequential number that identi es their temporal order|, applies to the genotype A2, for which it has been analyzed a particular subsequence of the genome (targeted). In order to get this single information directly available (instead of searching them every time from the T umor Sample Barcode) it is su cient to select the checkbox named Extract data from identi er that adds those new elds to the selected ones, with headers named patientID, T ype, T ime, Genotype, and W ide.
Subsequently our program creates a new le containing the selected elds, as well as a new one called V alGene, obtained by normalizing the variant allele frequency|namely the ratio between t alt count e t depth|, and then multiplies it by a reasonable value (we choose 107 in order to consider the rst six digits after decimal point) rounding the result. We perform these additional operations to cast variant allele frequencies to integers, as it is easier to work with them in ASP with respect to the numbers in oating point.
Once the le containing the necessary elds has been obtained, by clicking on Create -> ASP file, a graphical interface asking the number of genes that we want analyze is shown. In Figure 7 the reader can see the default value. On click data are ordered, identifying the most frequent genes, taking the rst k elements, with k as speci ed. An Unknown gene is a gene whose sequence is not known, so our program ignores it, as there is not necessarily a correlation between di erent unknown genes.
For the ASP syntax we use lowercase letters for all literal values, and replace
characters \-" and \." by \ ".
2 For our usage it is su cient to know that di erent letters correspond to di erent
types. Refer to [
Then some Java code creates instances and adds them to the main ASP program under construction. Inessential modi cations to adapt the analysis to the data under the study are performed at this stage in order, for example, to discriminate sample from same patient but with di erent genotypes: for each sample we add a fact for the type, wide, and genotype, in particular we consider the last one. Here there are the predicates deduced from sample SA429X1A-targeted:
As for the other rules we get:
2 sample(P,E1,T1),sample(P,E2,T2),T2=T1+1, 3 val(E1,G,V1),val(E2,G,V2),variation(V1,V2), 4 genotype(E1,GT),genotype(E2,GT). 5 precedes(P,GT,G1,G2) :
Fig. 7. Create ! ASP le 6 firstVariation(P,GT,G1,T1),firstVariation(P,GT,G2,T2),T1<T2. 7 firstVariation(P,GT,G,T) :8 variationConsec(P,GT,G,T), not varPrec(P,GT,G,T). 9 varPrec(P,GT,G,T) :10 variationConsec(P,GT,G,T1), variationConsec(P,GT,G,T),T1<T. 11 precedeKTimes(G1,G2,K) :12 gene(G1),gene(G2), 13 K=#count{patient(P),genotype(GT):precedes(P,GT,G1,G2)}.
The arity of those rst four predicates is increased by one in order to discriminate di erent genotype variations. Moreover, we infer variationConsec if genotypes in the same sample turn out to be equal. Predicate precedeKTimes has the same arity but the function count works with di erent pairs of patientgenotype rather than just with patients.
Finally Run -> Clingo and Graphviz allows the user to select the ASP le created and solve it with Clingo (see Fig. 8). The result is then redirected to a temporary text le, from which the row containing inferred predicates maybePrecede is selected. Those are processed into edges of a direct graph in dot format, on which we run Graphviz to display the result|see Fig. 9. This picture shows the variations detected over 100 genes, with the predicate variation that identi es variations with a value above 1000.
Notice that the height of the graph is barely 1. This is because, despite there are on average seven samples for the fteen patients, they are referred to di erent genotypes, so this increases the number of \sample to analyse" but at most (only) three temporal values are available for each of them, signi cantly reducing the number of relations found.
Finally with the last item of the menu, Help, we can open links to the documentation of the resources used: the home page of cBioPortal, the o cial description of MAF, the web site of Potassco and the description of dot format. 6
In this paper we report on a new system capable of dealing with data retrieved
from carcinome analysis. The system comes equipped with a front-end that
preprocesses/ lters standard-format data and transform them into a logic program,
as well as with a post processor that allows to visualize the results in a graph.
The reasoning core is based on ASP solving | the Clingo system is used. The
predicates implemented for the data analysis reported in the paper are just
simple examples and the expressivity of ASP can be exploited by more complex
queries/analyses. In particular, one of our future work will be trying to reconstruct
phylogenetics of cancer (using the ideas presented in [
A further direction of study consists in an integration with ILP, as a means to replace or collaborate with ASP, as far as the reasoning part is concerned. Acknowledgments. We thank Alessandro Dal Palu, Andrea Formisano, and Enrico Pontelli for the several discussions on this material. Agostino Dovier and Alberto Policriti are partially supported by INdAM-GNCS 2015{2017 and by PRID ENCASE projects.