dataesr/scientific-paragraphs-categorization huggingface dataset:
Each observation was converted to a U-score  (Guyon et al. , 2015) , which has a standard Gaussian distribution under the null hypothesis that we observe no phenotypic difference in the gene knock-down condition. We only consider the negative phenotype and hence a low P-value is associated to a negative siRNA score corresponding to a decrease of chromatin bound OGG1 upon siRNA transfection. For each gene, we converted either the known RefSeq gene ID or gene symbol to Entrez gene ID using bioDBnet  (Mudunuri et al. , 2009)  (latest release bioDBnet 2. 1; May 6, 2015), which was subsequently converted to STRING ID  (Szklarczyk et al. , 2014) , resulting in 4006 genes. Hence our RNAi screening data set is made up of 4006 genes each with 3 scores/P-values (hereafter OGG1 data). The STRING database  (Szklarczyk et al.
Both 'experimental' and 'database' concern interaction information gathered from other PPI databases which are also likely to be affected by an 'annotation bias' among others  (Gillis et al. , 2014). The 'co-expression' evidence concerns genes found to have been co-expressed in a variety of experiments and across a number of different species  (Stuart et al. , 2003; Von Mering et al. , 2005). Conservation of co-expressed genes over multiple species implies a selective advantage and hence that the genes are functionally related   (Xulvi-Brunet and Li, 2010). This is one source in the STRING database that we expect not to be influenced by an annotation bias  (Gillis et al. , 2014)  and with a low unique contribution, this evidence additionally conforms to that from the other sources. Supplementary Figure  S2  shows the vertex degree and the number of functional/GO annotations from DAVID  (Huang et al.
0054). Recall that for the OGG1 data there are three scores and hence three P-values for each gene. However, there is no natural or intrinsic way to order the siRNAs for each gene so that the data can be considered as trivariate. Correspondent trivariate data are achieved by ordering the three P-values so that the first data dimension is the lowest P-value, the second data dimension is the median P-value and the third data dimension is the highest P-value for each gene. Figure  2  presents the density schematic for the OGG1 data. Here, there are four labels of interest, a given number of negative siRNAs, as well as the dominant 'no negative siRNAs' label. Note that the densities in each dimension intercept at increasingly significant Pvalues. This is necessary as due to the ordering of the P-values, the lowest are skewed towards 0 and require a lower intercept so that the 'no negative siRNAs' label remains the dominant label.
The MRF method finds greater pathway enrichment and allows for more detailed determination of hit genes in a lymphoma study We consider data from a study of diffuse large B-cell lymphoma previously analysed using both the BioNet and Knode methods. For each gene, risk association was obtained from a survival analysis and differential expression was measured between the lymphoma subtypes ABC and GCB  (Beisser et al. , 2010; Dittrich et al. , 2008). Hence there are 2 P-values for each gene corresponding to survival analysis (S) and differential expression (T). A method for combining multiple P-values using order statistics was also proposed with BioNet  (Dittrich et al. , 2008; Beisser et al. , 2010). As can be seen in Supplementary Figure  S4 , the low combined P-values are those that have both low S and T P-values and if one of the P-values is high, then the combined P-value is high.
The PPI network for the lymphoma data was sourced from the Human Protein Reference Database  (Prasad et al. , 2009) , which is no longer maintained. We could have considered obtaining another PPI network from a more contemporary source and also attempted to account for the previously identified 'annotation bias'. However, the aim of this section is not to analyse the lymphoma data for its own sake but rather to compare the MRF output to the previously published output from the BioNet and Knode methods. The PPI network is made up of 2034 vertices and 15 512 unweighted edges (density ¼ 0. 0038). We consider the lymphoma data to show a general advantage of the MRF method in that we both find greater pathway enrichment and are able to more generally determine the hit genes as against the BioNet analysis.
S4 ). There is reasonable overlap between the vertex degree hit list and all the other lists except 'S hits only'. This makes sense for the MRF hit lists as vertex degree is explicitly present the MRF score (Equation (  4 )). This is also not surprising for BioNet since vertices with high degree are likely to be included in the maximum-weight connected subgraph since they allow access to the most number of other vertices through their edges. Also shown in Figure  6  are genes associated with the carcinogenic NFjB pathway  (Hoesel and Schmid, 2013) , originally used to evaluate the BioNet output for the lymphoma data  (Dittrich et al. , 2008). We can see that most of the NFjB genes contained in the BioNet list are also present in the 'S and T hits' list, while there are additional such genes in the 'T hits only' list.
A general advantage of the MRF method that it readily allows for more than a single hit list for multiple labels of interest. The additional information available in this way is not possible to determine after the P-values have been combined as was done for input into the BioNet method. The output from Knode (Supplementary Fig. S5 ) generally conforms to the above discussion. Although multivariate P-values do not appear to be a problem in principle for either the BioNet or Knode methods, it is not clear how multiple labels beyond binary 'hit' and 'non-hit' could also be achieved. The output of the NePhe and NEST methods, as well as simply ordering by combined P-value is also provided in the supplementary MATLAB code.
The four labels are: 'S and T hits' (low P-value in both dimensions), 'S hits only' (low P-value in S dimension), 'T hits only' (low P-value in T dimension) and 'no hits' (no low P-value in either dimension) Fig. 6. Venn diagram of the hit lists for the lymphoma data overlaid with genes associated with the NFjB pathway. There are 46 genes in each hit list (black) along with the genes associated with the NFjB pathway (grey and listed at the sides in blue). The NFjB annotation was obtained from KEGG  (Kanehisa et al. , 2016). Below each method is the P-value for the Fisher exact text for enrichment or depletion of the NFjB terms in the hit list. The Venn diagram is based on the layout from http://bioinformatics. psb. ugent. be/webt ools/Venn/ positives in the BioNet hit lists, the number of true positives found is still generally worse than other methods, which was also the case in the original simulation experiment (Supplementary Fig.
Implementation in MATLAB and computational expense The MATLAB code to reproduce the results presented in the paper is given as Supplementary Material. Our MRF based network scoring method was implemented in MATLAB on a mid-2012 MacBook Pro with a 2. 6 GHz Intel Core i7 processor (quad-core) and 16 GB of RAM. The minimum energy labels were found using the a-expansion algorithm  (Boykov et al. , 2001; Boykov and Kolmogorov, 2004; Kolmogorov and Zabih, 2004). Running over n ¼ 1000 values of b took $1. 3 s for the simulated data (1000 vertices, 2 labels), $20 s for the lymphoma data (2034 vertices, 4 labels) and $70 s for the OGG1 data (4006 vertices, 4 labels). Due to the heterogeneous nature of the genomic data and networks, our MRF-based method is currently not suitable for a graphical user interface.
Conclusion We have proposed a network based method to determine hit genes using an MRF and have shown the effectiveness of the method using a broad range of different data sets. The multiple advantages of the MRF method are that it easily allows for multivariate gene scores in addition to multiple labels of interest beyond binary 'hit' and 'nonhit'. The intercept of the underlying densities can be considered as a parameter with an interpretation and for which we have provided general advice on setting. We have provided the MATLAB code and data to reproduce the results presented in the paper, which is freely available for modification. We have shown that our MRF method gave the best results in a simulation experiment originally used to evaluate Knode and BioNet. We have additionally shown that we were able to find greater pathway enrichment as well as further determine hit genes in a lymphoma study concerning survival analysis and gene expression, whereas the multiple observed P-values for each gene were combined for input into the BioNet method.
For the OGG1 screening data we have shown that the MRF method allows for the analysis of the multivariate scores and is able to find relevant functional/GO annotation. For the study of genomic data with an associated network where the 'guilt by association' assumption is appropriate, we have shown that our proposed MRF based method gives an advantageous way to determine network based hit genes as against previously proposed methods. We also noted a major pitfall using functional/GO enrichment for validation when it is strongly correlated with vertex degree in the PPI network and confounded by an 'annotation bias'. We considered a network of only co-expression PPI evidence to account for this issue and also took this into consideration when evaluating our results. This issue does not appear to have been given enough attention in the literature, especially when there are so many papers utilizing both PPI networks as well as functional/GO enrichment analysis.
This tool is also implemented in the open-access workbench OptFlux  [15]. An alternative approach to the same problem that combines information from flux variability profiling and metabolite centrality is RobOKoD, where knockout, overexpression and dampening strategies are suggested based on the profiling of every reaction in the model and its relation to the target production  [16]. Each of these methods is based on different approaches and rationales, therefore leading to distinct solutions. In fact, one key issue is the multitude of obtained results, known to be dependent on the specific algorithms and the software implementation used. In order to address this problem, a consensus-based approach is here proposed and developed. The method is based on running several optimization procedures for knockout prediction and analyzing a posteriori the consensus solutions obtained.
From the available optimization methodologies, OptKnock was chosen, as it is one of the most used methods, as was the corresponding pessimistic prediction strategy, RobustKnock. In addition, OptGene and RobOKoD were also introduced in this procedure, as they represent very different approaches. Naturally, each method has its advantages and disadvantages: OptKnock and RobustKnock are optimization-based methods, and return an optimal solution to the defined mixed integer linear problem (MILP). Notice that there might be numerous solutions with optimal objective value, out of which only one is returned. In contrast, there are heuristic methods, like the Evolutionary Algorithms and Simulated Annealing implemented in the software tool OptFlux  [15]  and OptGene  [14] , respectively. These methods allow nonlinear/non-convex objective functions, and return a set of candidate mutants.
Implementation The pipeline was implemented in MATLAB, based on the COBRA toolbox  [4, 5]  and with the external use of OptFlux  [15]. The full detailed description and documentation of OptPipe can be found together with the software at https://github. com/AndrasHartmann/OptPipe along with all the source code. The OptPipe software consists of three basic directories: (I. ) common_functions contains the functions that are shared between the methods; (II. ) methods contains the implementation of the main methods of the pipeline based on the COBRA toolbox, which is added to the (III. ) external folder together with other libraries and toolboxes that are being used by the pipeline code, such as RobOKoD  [16] , RobustKnock  [27]  and xlwrite, that is a library for xls reading and writing. In addition, the (IV.
Results and discussion In order to test and validate the proposed method and software, OptPipe was applied to the optimization of naringenin production in C. glutamicum, a bacterium that is routinely used in the industrial production of amino acids  [28]. Naringenin is a phenolic compound belonging to the flavonoids family that is naturally produced by some plants but is not endogenous in the C. glutamicum metabolism. This compound has interesting properties for health applications, such as being an antioxidant and chemoprotective  [29, 30] , and it also serves as a precursor for other flavonoids. Naringenin production from p-coumarate requires the expression of three heterologous genes, which code for 4-coumarate: CoA ligase (4CL), chalcone synthase (CHS) and chalcone isomerase (CHI). 4CL activates p-coumarate to its CoA-thioester p-coumaroyl-CoA.
Taking into account the limiting character of the malonyl-CoA flux in naringenin optimization, it is unsurprising that the knockout of this pathway is expected to increase naringenin production. The remaining 11 double mutants had the same q-values, which in some of the cases is due to their equal characteristics and in other cases might reflect a trade-off between higher growth rates and lower distances between the mutant flux distribution and the wild-type. The results obtained with OptGene and RobustKnock were excluded from this list since they did not comply with the specified growth threshold (data not shown). See detailed results and gene annotation in the Additional files 1 and 2. The proposed reaction deletions were categorized into main metabolic subsystems. The enrichment of the overall knockout dataset was computed to understand which sections of the bacterial metabolism were mostly identified as knockout targets when optimizing naringenin production (Fig.
Funding This work was supported by the  European Union   Framework Program 7 "BacHBerry " (www. bachberry. eu),Project no. FP7-613793  and the  Portuguese Foundation for Science & Technology , through IDMEC, under LAETA, project  UID/EMS/50022/2013. SV acknowledges support by  Program Investigador FCT  ( IF/00653/2012 ) from  FCT , co-funded by the  European Social Fund (ESF)  through the  Operational Program Human Potential (POPH).
Availability and requirements Project name: OptPipe Project home page: https://github. com/AndrasHartmann/ OptPipe Operating system(s): Platform independent Programming language: MATLAB (Tested on 2015a and 2016b) Other requirements: MATLAB parallel toolbox COBRA toolbox from the openCOBRA project (release 2. 0), SBML toolbox, Gurobi optimization software installed for MATLAB OptFlux (recommended) License: GPL v3 Any restrictions to use by non-academics: none Install: See detailed description on the project home page Availability of data and materials Software and data available at https://github. com/AndrasHartmann/OptPipe.
Methods The IFM-DFCI cohort. 383 patient samples from the IFM-DFCI 2009 cohort  32  , included from November 2010 to November 2012 (https://clinicaltrials. gov) were initially sequenced in three batches following the protocol described below. All the samples were collected at diagnosis in France and processed at the University Hospital of Nantes. All included patients were newly diagnosed with symptomatic MM based on International Myeloma Working Group 2003 Diagnostic Criteria  33. All patients signed an informed consent form approved by the Toulouse Ethics Committee. Patients were younger than 65 years old, with a mean age of 56. 6 years (median 58) and were 57% males and 43% females. ISS stages were 31% stage I, 47% stage II, and 22% stage III: with a median beta 2 microglobulin level of 3.
We chose 220 based on empirical values obtained by prior alignment using Star 37 (v2. 4. 0). All other parameters were left to default values. We considered a fusion to be found by both algorithms as long as the same partner genes were involved in the same orientation, even in presence of differences in the exact breakpoint predictions. For instance, our definition considers IGH-MMSET and MMSET-IGH to be two different fusions. In this way we identified an initial set of 2501 common candidate fusions. These were further filtered as described next. The pipeline developed for this analysis is available in Supplementary Data 2. Fusion filtering. We used two criteria: (a) we removed fusions which did not have at least three pairs of reads for which one end covers the junction. We implemented this through a script based on the number of reads reported by both algorithms, namely "Number of spanning mate pairs where one end spans a fusion" for TopHat, and the sum of "multiple_paired_read_count" and "unique_paire-d_read_count", i.
e. , the number of multiple and unique mapped reads mapped to the fusion and that are paired with their mates for MapSplice; (b) we removed fusions for which both partners belong to the same gene family, as we are concerned about the potentially high sequence similarity. We implemented this through a script based on the official HGNC gene-family dataset, downloaded from http://www. genenames. org/cgi-bin/genefamilies/download-all/tsv on 08/10/15. The pipeline developed for this filtering is available in Supplementary Data 2. Sequence homology evaluation. For each candidate fusion (including those filtered out according to the previous criteria), we computed sequence homology scores based on TopHat's prediction of breakpoint localization. To this effect, we retrieved the 20 bp sequences before (3ʹ) and after (5ʹ) the breakpoint of each partner gene, and computed both 3ʹ (respectively 5ʹ) sequences' homology using FASTA  38  with default parameters.
Our final score, given in the last column of Supplementary Data 1, is the minimum of both 3ʹ and 5ʹ E-score as obtained by FASTA. Those E-scores can be interpreted as a P-value for sequences to not share similarities: the smaller the E-score, the higher the homology. The code used for this computation is available in Supplementary Data 2. Geneset enrichment. We performed gene-annotation enrichment analysis using The Gene Set Enrichment Analysis tool from the Broad Institute (v5. 1)  39. All fused genes, unique partner genes from Ig fusions and recurrently fused genes (involved in at least two fusions) were used (list provided as Official Gene Symbol identifier under the Gene Identifiers List in GSEA. Annotation was restricted to Homo tissue compendium (Novartis), and summary results limited to Hallmark gene sets.
Acknowledgements This work was supported by  NIH  grants  PO1-155258  and  P50-100707  (NM, H. A. -L. , KA and G. P. );  Department of Veterans Affairs Merit Review Award 1   I01BX001584-01  (NM) and the  CAPTOR program. The CRCT Team 13 is labeled by  ARC. G. P. was supported by  NC1 P30-006516-51. We thank the  Intergroupe Francophone du Myelome  for providing patient samples and clinical data, and  Dr. Keats J. J. for providing RNA-Seq data from MM cell lines.
Data availability. The authors declare that all data supporting the findings of this study are available within the article and its supplementary information files. Data can be accessed from the dataverse portal with doi:10. 7910/DVN/RFL6TK (https:// dataverse. harvard. edu/dataset. xhtml?persistentId=doi:10. 7910/DVN/RFL6TK). Additional information are available from the corresponding author upon reasonable request.
INTRODUCTION The development of ChIP-seq technology  (1)  has enabled the construction of genome-wide maps of protein-DNA interactions. Such maps provide information about transcriptional regulation at the epigenetic level (histone modifications and histone variants) and at the level of transcription factor activity. Recently, thousands of ChIP-seq datasets have been produced by different consortia including EN-CODE  (2)  and the NIH Roadmap Epigenomics Mapping Consortium  (3). The data produced contain histone modification libraries for both normal and cancer cell karyotypes. In cancer, genetic and epigenetic abnormalities cooperate in the process of regulating activities of oncogenes and onco-suppressors  (4). For example, lower levels of trimethylation of lysine 36 of histone H3 (H3K36me3) and trimethylation of lysine 20 of histone H4 (H4K20me3) in proximity of the gene NSD1, contribute to the development of nervous system tumors  (5).
Moreover, some methods can account for experiments with either biological or technical replicates (PePr  (22) , DiffBind  (13)  and csaw (  23 )), while other methods cannot (ODIN  (20) , ChIPDiff (  16 ) and MACS2). In this study, we introduce HMCan-diff, a method for identifying changes in histone modifications from ChIP-seq cancer data. Our method corrects for copy number aberrations, GC-content bias, sequencing depth, mappability, and noise level, thus accounting for different technical artifacts of ChIP-seq data, and utilizes information from replicates to reduce technical variation effects. We compared HMCan-diff with several recent and most commonly-used methods, namely ChIPDiff  (16) , MAnorm  (19) , MEDIPS  (21) , ODIN  (20) , MACS2 (https: //github. com/taoliu/MACS/tree/master/MACS2), DiffBind  (13) , RSEG  (18)  and csaw  (23).
Construction of normalized density profiles HMCan-diff uses density construction and normalization methods implemented in the HMCan algorithm  (10). The normalization steps include normalization for copy number Figure  1. A workflow illustrating HMCan-diff steps. Initially, HMCan-diff constructs a fragment density profile for each provided ChIP-seq or input dataset. Then, it normalizes density profiles of each replicate in each condition for several types of bias, specifically for copy number variation, library size, GC-content bias and noise level. After that, HMCan-diff conducts additional normalization to eliminate further technical variation between conditions. It initializes HMM parameters based on the data. In particular, HMCan-diff defines the HMM emission probability distribution as the joint empirical distribution of normalized density values.
After raw density construction, HMCan-diff reports density values every N bp, where N is a step fixed by the user. If they like, users can choose the step size as minL/2, where minL is the minimum fragment length. Smaller step sizes will not increase the performance but will PAGE 3 OF 9 Nucleic Acids Research, 2017, Vol. 45, No. 8 e58 provide a better profile resolution for further downstream analysis. In this study we set N to be 50 bp. HMCan-diff also proposes to mask densities from genomic regions known to contain false signals. HMCandiff provides blacklisted regions recommended by the EN-CODE consortium (2) as its default set of regions. Correction for copy number alterations. HMCan-diff utilizes the Control-FREEC algorithm  (25, 26)  to learn the copy number profile from input DNA data (control sample).
HMCan-diff learns parameters of the NB distribution in a similar way to the method proposed in  (27). HMCan-diff sets the first option as its default for identifying enriched density regions. GC-content normalization. It has been reported that GCcontent affects the sequencing depth of Illumina reads  (28). Therefore, HMCan-diff implements a correction of ChIPseq density data for GC-content bias. First, HMCan-diff associates each density value with a GC-content value. For each density value, HMCan-diff calculates the GC-content for a window of length equal to twice the median fragment length centered on that value, and assigns it to this density value. Then, HMCan-diff groups GC-values into different strata (groups). For example, the first group will have GC content values from the interval [0, 0.