dataesr/scientific-paragraphs-categorization huggingface dataset:
acnes phylotype I (C. acnes subsp. acnes) rather than II (C. acnes subsp. defendens). Whilst our Marat strain does not cluster with phylotype Ia, the type more commonly associated with skin surface associated acne vulgaris  (Lomholt and Kilian, 2010) , its position, basal to Type Ib strains cannot exclude its involvement in soft or deep tissue infections  (Nazipi et al. , 2017). Delineating contaminants and commensals from plausible pathogens remains challenging from this type of data source, in particular due to the absence of a suitable control. To alleviate this issue, we conducted full taxonomic assignments of two ancient metagenomes generated from historical parchment samples dating to the 17th and 18th centuries (PA1 and PA2 respectively)  (Teasdale et al. , 2015a). Although these samples were obtained from livestock (ruminant) skins whereas we are working with cellulose paper, we anticipate that they may have been used and handled in a comparable way to the newspaper Marat was annotating.
Acknowledgments and funding This work was supported by  Obra Social "La Caixa" and Secretaria d'Universitats i Recerca  ( GRC2017-SGR880 ) (T. M. -B. and C. L. -F. ),  BFU2017-86471-P  and  PGC2018-101927-B-I00  ( MINECO/FEDER, UE)  (T. M. -B) and  PGC2018-095931-B-100  ( MCIU/AEI/FEDER, UE)  (C. L. -F. ). T. M-B. is also supported by a  U01 MH106874  grant and  Howard Hughes International Early Career  and  CERCA Programme del Departament d'Economia i Coneixement de la Generalitat de Catalunya. S. M is funded by a  Welcome Trust post-doctoral fellowship  ( 206478/Z/ 17/Z ). L. v. D and F. B. acknowledge financial support from the  Newton Fund UK-China NSFC initiative  (grant  MR/P007597/1 ) and the  BBSRC  (equipment grant  BB/R01356X/1 ). Full metagenomic reads are available at NCBI under BioProject  ID PRJEB32319.
Background: Integrating genome-wide gene expression patient profiles with regulatory knowledge is a challenging task because of the inherent heterogeneity, noise and incompleteness of biological data. From the computational side, several solvers for logic programs are able to perform extremely well in decision problems for combinatorial search domains. The challenge then is how to process the biological knowledge in order to feed these solvers to gain insights in a biological study. It requires formalizing the biological knowledge to give a precise interpretation of this information; currently, very few pathway databases offer this possibility.
Results: The presented work proposes an automatic pipeline to extract automatically regulatory knowledge from pathway databases and generate novel computational predictions related to the state of expression or activity of biological molecules. We applied it in the context of hepatocellular carcinoma (HCC) progression, and evaluate the precision and the stability of these computational predictions. Our working base is a graph of 3383 nodes and 13,771 edges extracted from the KEGG database, in which we integrate 209 differentially expressed genes between low and high aggressive HCC across 294 patients. Our computational model predicts the shifts of expression of 146 initially non-observed biological components. Our predictions were validated at 88% using a larger experimental dataset and cross-validation techniques.
In particular, we focus on the protein complexes predictions and show for the first time that NFKB1/BCL-3 complexes are activated in aggressive HCC. In spite of the large dimension of the reconstructed models, our analyses over the computational predictions discover a well constrained region where KEGG regulatory knowledge constrains gene expression of several biomolecules. These regions can offer interesting windows to perturb experimentally such complex systems. Conclusion: This new pipeline allows biologists to develop their own predictive models based on a list of genes. It facilitates the identification of new regulatory biomolecules using knowledge graphs and predictive computational methods. Our workflow is implemented in an automatic python pipeline which is publicly available at https://github.
Upregulation of EMT markers has been associated with tumor aggressiveness and bad prognosis  [7, 8]  and associated with inflammatory microenvironment  [9]. However, in vivo monitoring of EMT processes remains difficult, due to the spatio-temporal dynamics of these molecular events and the snap-shot nature of biopsies sampling. Understanding EMT to identify new therapeutic targets require integrative and modeling approaches. To build computational models and integrate experimental data on molecular events, pathway databases can be used. However, despite the fact that numerous publicly available pathway databases currently exist, compiling hundreds of signaling pathways for various biomolecules, very few formal representations linked with automatic inference processes have been proposed so far  [10].
Knowing that signaling cascades are represented by multiple complexes assemblies, this misinterpretation impacts importantly the construction of a mechanistic model when using pathway databases. On the other hand, such tools are very useful to compute topological scores, perform statistical analyses, and to integrate gene expression measurements using enrichment analyses  [17]. They remain, however, limited to extract logical consequences of the representation of the biological mechanisms. The sign-consistency framework proposes a way to automatically confront the logic of large-scale interaction networks and genome-wide experimental measurements, provided that a signed oriented network is given and that the experimental measurements are discretized in 3 expression levels (up-regulated, down-regulated and nochange).
Our case study is composed of two input data which were publicly available. First, gene expression data from patients with HCC was extracted from the International Cancer Genome Consortium (ICGC) database  [21]. Based on the EMT signature from MSigDB  [22] , HCC samples were clustered into either agressive HCCs (high EMT gene expression) or non-agressive HCCs (low EMT gene expression). Second, the up-stream events of the regulatory events of these genes were obtained by querying automatically KEGG to build a causal model from this database. We used Iggy to study what are the regulatory events that explain the differential expression between low and high aggressiveness from the KEGG interaction knowledge (network of 3383 nodes and 13,771 edges). We discovered that 146 nodes were predicted, of them 33 refer to gene expression, 110 were protein activities, and 3 were protein complexes activities.
88% of the predictions were in agreement with the ICGC gene expression measurements. Importantly, we predicted the activation of NFKB1/BCL3 and NFKB2/RELB complexes, two critical regulators of NFKB signalling pathway implicated in tumorigenesis. Finally, we proposed a method to discover sensitive network regions that explains HCC progression. This means network components which were highly constrained by multiple experimental data points that could   2 ) Extract the upstream/downstream signaling pathways for the set of genes from the signed interaction graph using Pathrider, a tool developed in our team to this purpose. Given a list of excluded genes (such as invariant genes), Pathrider filters these genes to reduce the graph size; (3) Check the sign consistency of our datasets to produce signed predictions for unmeasured biomolecules using iggy tool; (4) Validate the predictions made by iggy by computing sub-predictions (prediction 1, 2.
Overview of the pipeline We introduce key-pipeline: a Python package implementing the workflow for identifying key protein complexes associated to tumor progression. The general pipeline implemented is depicted in Fig. 1. It receives as input data: a list of differentially expressed genes, a graph describing signed and directed signaling interactions, and a set of excluded genes to be filtered out from the graph. Our software allows researchers to: 1) construct the set of observations, and the set of gene names from a file of differentially expressed genes (in CSV format with 3 columns: genes names, log 2 (fold-change) and adjusted pvalue), see "Identification of gene differentially expressed between low and high aggressive HCC" section; 2) extract a specific regulatory and signaling network associated with the input genes list from a signed interaction graph (based on KEGG regulatory knowledge, see "Building a signed interaction graph from the KEGG Pathway database" section); 3) apply the Iggy tool to compute the predictions based on the sign-consistency modeling (see "Sign consistency -Iggy" section); 4) perform robustness and stability analyses (see "Computational Validation of the results" section); and finally 5) generate plots of these analyses.
Integration of Gene Expression in Signaling and Regulatory Network A first interaction graph was built from the KEGG Pathway database (see "Building a signed interaction graph from the KEGG Pathway database" section). This graph was composed of 41,546 interactions (gene transcriptions, protein signaling, protein formation and complex formations) and 8861 components (genes, proteins and protein complexes). It is available as input-data with the pipeline 2. Using our pipeline (Fig. 1 , step 2: Pathrider), the 1913 differentially expressed genes between low and high aggressive tumors were used as input to extract a subgraph from the KEGG pathways graph. Only 209 genes from the 1913 were identified and used to extract upstream regulatory events (Additional file 1: Table  S1 ). In this step the biomolecules associated to the 4220 genes whose expression is undetectable were filtered out.
Among the 146 predictions, 143 have a name that matches with a gene name identified in HCC from the ICGC database, but that were not selected in the 1913 genes differentially expressed (< -0. 5 for down-regulated and > 2 for up-regulated in aggressive HCC). If we remove all thresholds and thus consider any positive 3 Laptop computer containing an Intel Core i7-5600U CPU with 4 threads of 2. 60GHz and running Fedora 27 64 bits. fold-change as an over-expression, and any negative foldchange as an under-expression, then 82 components predicted + are coherent with the ICGC data and 8 are not; 44 components predictedare coherent with the ICGC data and 9 are not. This ratio of 88% of matching predictions speaks in favor of our choice of applying Iggy to this specific biological system, with respect to the currently available data in KEGG and ICGC databases.
Activation of NFκB signaling in aggressive HCC Based on the regulatory model (Additional file 1: Figure  S1 ) and differential expression of mRNA between low and high aggressive HCC (see "Identification of gene differentially expressed between low and high aggressive HCC" section), the algorithm Iggy predicts the activation of complexes NFKB1::BCL3 and NFKB2::RELB and the deactivation of complex JUND::NACA. By activation of the complexes we mean that in order to explain the highaggressive versus low-aggressive tumor gene expression datasets, these complexes have to increase their activity. For example, if NFKB1::BCL3 is activated, then we deduce that its effect on gene IL10 (regulated positively by this complex, Additional file 1: Figures S6) is positive, meaning that the level of gene IL10 may increase if it was only regulated by NFKB1::BCL3.
However, HCC results from complex interactions between the tumor cells and the microenvironment involving stromal cells and extracellular matrix. Molecular biological data from tumor tissues recapitulate all this information and we need to build an unique large-scale model without a priori to take into account such complexity. For that purpose, we propose here an original approach aiming at integrating experimental data on a regulatory graph extracted from the KEGG database to predict new markers and regulators of HCC progression. Based on EMT gene expression signature from MSigDB  [22]  we first separated low from high aggressive HCC samples stored in the ICGC database  [21]  and next we sought to predict the regulatory pathways implicated in this transition. For that purpose, we built a model by querying the KEGG database using the KEGG API to extract an initial network.
Identification of gene differentially expressed between low and high aggressive HCC Normalized HTseq counts and clinical data were retrieved from the LIHC-US project  5  (NCI, TCGA-LIHC). These files were downloaded on 2016-07-19, corresponding to release 21. At this date, LIHC-US dataset comprised 294 donors and 345 samples; among them, we selected samples corresponding to solid primary tumors, based on clinical data, by selecting entries containing the expression "Primary tumour -solid tissue" in the specimen table (7 th field). This allowed selecting one sample for each of the 294 donors. Data retrieval and filtering workflow is detailed in Additional file 3: dataset_filtering. sh. From this filtered dataset, we extracted a two-dimensional table of expression values (converted in log 2 ) for 20,502 genes in 294 LIHC samples.
Based on the bimodal distribution of these expression values, we discarded genes whose expression is undetectable (4220 genes), keeping 16,282 genes. Expression values were normalized by the median value in each sample. Based on the established link between epithelial-mesenchymal transition (EMT) and tumor aggressiveness  [41] , we used the MSigDB  [42]  set of 200 genes termed HALLMARK_EPITHELIAL_ MESENCHYMAL_TRANSITION 6 from the Broad Institute as a molecular signature of aggressiveness. From the LIHC dataset, we extracted a table of expression values for 195 entries of this EMT signature for each of the samples (5 genes were undetectable). Based on the expression values of the EMT signature, LIHC samples were classified (hierarchical clustering of euclidean distances) into three groups termed low_EMT (70 samples), medium_EMT (154 samples), and high_EMT (70 samples).
Together, these 1913 differentially expressed genes, listed in Additional file 2: diffexp_filtered. csv were subsequently used to extract a regulatory network, as explained in "Building a signed interaction graph from the KEGG Pathway database" section, and then used as observations for the coloring propagation process, as detailed in 6 Id: M5930, available at http://software. broadinstitute. org/gsea/msigdb/ geneset_page. jsp?geneSetName=HALLMARK_EPITHELIAL_ MESENCHYMAL_TRANSITION 7 Files GSEA_EMThigh_vs_EMTlow_diffexp. csv (differential expression) and LIHC_primary_weakly_expressed_genes. txt (blacklist of weakly expressed genes) at https://github. com/ LokmaneChebouba/key-pipe/tree/master/example. "Sign consistency -Iggy" section. The full list of differentially expressed genes is available in Additional file 2: diffexp_full.
csv. The workflow of data clustering and differential analysis is available in Additional file 3: diffexp_and_clustering. R. To further validate the clinical relevance of the groups of HCC samples identified by the hierarchical clustering method, we compared this classification obtained with the EMT signature to a classification obtained with markers established as predictive of recurrence, which is a major clinical outcome of tumor aggressiveness. For that purpose, we used the Seoul National University recurrence (SNUR) signature  [43]  that previously permitted to classify HCC samples from TCGA database  [44]  and we compared clusters identified by hierarchical clustering method with SNUR groups. Note that LIHC primary tumors correspond to 294 samples but only 183 are annotated with SNUR score.
Building a signed interaction graph from the KEGG Pathway database For this case study, we used a human signaling network derived from the KEGG Pathway database  [13]. 154 human signaling pathways were fetched using the KEGG API and converted to SIF (Simple Interaction Format) in order to provide KEGG's knowledge as a graph representation. This section summarizes how this network was built. A more in-depth description is available in Additional file 1 as a Supplementary Material & Methods. This step of automatic reconstruction of a causal graph from KEGG constitutes one of the novel contributions of the methodological results of this paper.
Signed interaction graph built from KEGG's regulatory knowledge To model the KEGG regulatory knowledge we imposed a distinction between nodes representing genes and nodes representing proteins. In the KEGG Pathway database, this distinction is implicitly embedded in the relation types, particularly PPrel edges (protein-protein relations) and GErel edges (gene expression relations). PPrel edges indicate that both source and target nodes are proteins. GErel edges indicate that source nodes are transcription factors and that target nodes are genes. Therefore, to explicitly differentiate genes and proteins, the source nodes of GErel edges were suffixed with _prot and the target nodes were suffixed with _gen. Concerning PPrel edges, both the source and target nodes were suffixed with _prot. Differently to what proposed by KEGG maps, we modeled protein complexes explicitly by imposing two relations: A_prot PPrel A_prot::B_prot B_prot PPrel A_prot::B_prot where A_prot::B_prot refers to the protein complex formed by proteins A and B.
Altogether, the human signaling network extracted from KEGG was represented as a signed interaction graph composed of protein signaling interactions, complex formations, gene expression regulations and gene-protein relations; accordingly, the nodes of this graph refer to genes, proteins and protein complexes. This decomposed representation allowed us to map on this network data which corresponded only to gene expressions, without assuming that gene expression correlates to protein activity. This KEGG signed and directed graph is available as input data of the pipeline implemented in this work  8.
Modeling our case-study using sign-consistency: inputs and outputs of Iggy Inputs The signed interaction graph G(V , E, σ ) was obtained from the KEGG Pathways database as explained in "Building a signed interaction graph from the KEGG Pathway database" section. The gene expression profile (see "Identification of gene differentially expressed between low and high aggressive HCC" section) is composed of experimental knowledge of 821 over-(sign '+') and 1092 under-expressed (sign '-') genes in in high aggressive tumors compared to low aggressive tumors. From these differentially expressed genes, only 209 were found in the KEGG graph matching nodes with a suffix '_gene'. Thus, the set of observations (S, μ) used for the sign-consistency analysis was composed of 209 elements. Outputs Following the sign consistency analysis, Iggy proposed predictions under minimal MCoS repairs.
Computational Validation of the results Recall that to create the over-and under-expressed genes between low and high aggressive tumors (see "Identification of gene differentially expressed between low and high aggressive HCC" section) we used thresholds of +2 and -0. 5 on the value of log 2 (fold-change). In this section, we aim at checking if these thresholds are justified. To do this, we computed "sub-predictions", that is, predictions on the same extracted graph of "Building a signed interaction graph from the KEGG Pathway database" section but with subsets of observations. To generate these subsets of observations, we considered a range of samplings, from 10% to 95% of the complete observation set, with a step of 5%. For each sampling of x%, 100 experiments were conducted, where an experiment consisted in randomly picking x% of the over-expressed observations (+) and x% of the under-expressed observations (-), and computing the predictions on this subset of observations.
Behavioral analyses All behavioral analyses were conducted using R (http://www. rproject. org/) and relevant packages  (R Core Team, 2014). All the reported P values are two-tailed and P < 0. 05 was considered to be statistically significant. Data visualization was performed via 'ggplot2' package  (Wickham, 2016). Regarding the choice data, we performed a repeated mixedeffect logistic regression on the decision of choosing the 'accept' option by the glmer function in 'lme4' package  (Bates et al. , 2013) , with dilemma (dummy variable; reference level: self-bad cause dilemma; same below) and payoffs for both parties involved in each dilemma (i. e. the monetary gain and the moral cost; meancentered continuous variable; same below) as the fixed-effect predictors. In addition, we included intercepts varying across participants as the random effects.
For the statistical inference on each predictor, we performed the Type II Wald chi-square test on the model fits by using the Anova function in 'car' package  (Fox et al. , 2016) , and reported the odds ratio as relevant effect size. For decision time (DT), we first did a log-transformation due to its non-normal distribution (Anderson-Darling normality test: A = 91. 90, P < 0. 001) and then performed a mixed-effect linear regression on the log-transformed DT by the lmer function in 'lme4' package, with decision (dummy variable; reference level: accept), dilemma, decision × dilemma, as well as payoffs for both beneficiaries as the fixed-effect predictors. Random-effect factors were specified in the same way as above. Similar analyses were also performed on the post-scanning rating except that dilemma was added as the only fixed-effect predictor.
We followed the procedure recommended by  Luke (2017)  to obtain the statistics for each predictor by applying the Satterthwaite approximations on the restricted maximum likelihood model fit via the 'lmerTest' package  (Luke, 2017). In addition, we computed the Cohen's d of each predictor via the 'EMAtools' package  (Kleiman, 2017) , which provided the effect size measure especially for the mixed-effect regressions. For likeness ratings of the selected associations, we compared whether the ratings significantly differed from 0 in each type of selected associations (i. e. charity or morally bad causes), respectively, by the one-sample T-test, and computed the Cohen's d as effect size.
We used the 'hBayesDM' package  (Ahn et al. , 2017)  to fit all aforementioned candidate models using the hierarchical Bayesian analysis (HBA) approach  (Gelman et al. , 2014). The 'hBayesDM' package is developed based on the Stan language (Stan Development Team, 2016), which utilizes a Markov Chain Monte Carlo (MCMC) sampling scheme to perform full Bayesian inference and obtain the actual posterior distribution. We adopted HBA rather than maximum likelihood estimation (MLE) because HBA provides more stable and accurate estimates than MLE  (Ahn et al. , 2011). Following the approach in 'hBayesDM' package, we assumed the individual-level parameters φ were drawn from a group-level normal distribution: φ ∼ Normal (μ φ , σ φ ), where μ φ and σ φ refer to the group-level mean and standard deviation, respectively.
g. AIC). By convention, the lower WAIC score indicates better outof-sample prediction accuracy of the candidate model. A difference score of 10 on the information criterion scale is considered decisive  (Burnham and Anderson, 2004). We selected the model with the lowest WAIC as the winning model for subsequent analysis. In addition, we also implemented a posterior predictive check to further examine the absolute performance of the winning model, i. e. whether the prediction of the winning model could characterize the features of real choices. In specific, we employed each individual's joint posterior MCMC samples (i. e. 4000 times) to generate new choice datasets correspondingly (i. e. 4000 choices per trial per participant), given the actual trial-wise stimuli sequences presented to each participant.