Measuring stress level objectively is crucial for personalized health monitoring. While traditional methods require a clinical setting, wearables provide a valuable alternative. In this paper, we approach stress assessment as a regression task, focusing on stress exposure, and evaluate Functional Data Analysis (FDA) to extract richer information from physiological signals. We apply scalar-on-function regression and functional clustering to WESAD, a public dataset which contains signals from wearables and psychometric questionnaires that we use as a ground truth for stress. We compare the results obtained by applying FDA with those achieved by methods using features extracted from signals rather than the signals themselves. The comparison reveals that FDA excels in capturing signal variations and their association with stress, offering new insights into how this association changes with different stressful activities. While non-functional techniques suffice for some analyses, FDA is key to capture overtime patterns linked to stress levels.

We propose a novel methodology to validate software product line (PL) models by integrating Statistical Model Checking (SMC) with Process Mining (PM). We consider the feature-oriented language QFLan from the PL engineering domain. QFLan allows to model PL equipped with rich cross-tree and quantitative constraints, as well as aspects of dynamic PLs such as the staged configurations. This richness allows us to easily obtain models with infinite state-space, calling for simulation-based analysis techniques, like SMC. For example, we use a running example with infinite state space. SMC is a family of analysis techniques based on the generation of samples of the dynamics of a system. SMC aims at estimating properties of a system like the probability of a given event (e.g., installing a feature), or the expected value of quantities in it (e.g., the average price of products from the studied family). Instead, PM is a family of data-driven techniques that uses logs collected on the execution of an information system to identify and reason about its underlying execution process. This often regards identifying and reasoning about process patterns, bottlenecks, and possibilities for improvement. In this paper, to the best of our knowledge, we propose, for the first time, the application of Process Mining (PM) techniques to the byproducts of Statistical Model Checking (SMC) simulations. This aims to enhance the utility of SMC analyses. Typically, if SMC gives unexpected results, the modeler has to discover whether these come from actual characteristics of the system, or from bugs in the model. This is done in a black-box manner, only based on the obtained numerical values. We improve on this by using PM to get a white-box perspective on the dynamics of the system observed by SMC. Roughly speaking, we feed the samples generated by SMC to PM tools, obtaining a compact graphical representation of the observed dynamics. This mined PM model is then transformed into a mined QFLan model, making it accessible to PL engineers. Using two well-known PL models, we show that our methodology is effective (helps in pinpointing issues in models, and in suggesting fixes), and that it scales to complex models. We also show that it is general, by applying it to the security domain.

Feature Selection and Functional Data Analysis are two dynamic areas of research, with important applications in the analysis of large and complex data sets. Straddling these two areas, we propose a new highly efficient algorithm to perform Group Elastic Net with application to function-on-scalar feature selection, where a functional response is modeled against a very large number of potential scalar predictors. First, we introduce a new algorithm to solve Group Elastic Net in ultra-high dimensional settings, which exploits the sparsity structure of the Augmented Lagrangian to greatly reduce the computational burden. Next, taking advantage of the properties of Functional Principal Components, we extend our algorithm to the function-on-scalar regression framework. We use simulations to demonstrate the CPU time gains afforded by our approach compared to its best existing competitors, and present an application to data from a Genome-Wide Association Study on childhood obesity.

The H-score (or Mean Squared Residue score) underlies Cheng and Church’s (2000) biclustering algorithm, one of the best-known and most widely employed algorithms in bioinformatics and computational biology, and many subsequent algorithms (e.g. FLOC, Yang et al., 2005 and CBEB, Huang et al., 2012). Cheng and Church’s algorithm has ∼2600 citations to date, 650 since 2015 and 230 in 2018–2019 alone. It was the first to be applied to gene microarray data, and it is one of the main tools available in biclustering packages (e.g. the ‘biclust’ R library) and in gene expression data analysis packages (e.g. IRIS-EDA, Monier et al., 2019). In addition, it is widely used as a benchmark: almost all published biclustering algorithms include a comparison with it. Squared residue measures such as H-scores have a double role in biclustering methods. On the one hand, they are employed by many algorithms as merit functions to guide the discovery of biclusters (see e.g. the reviews in Madeira and Oliveira, 2004; Pontes et al., 2015). On the other hand, they are used to assess solutions—in particular, H-scores are used to assess the ‘homogeneity’ of the discovered biclusters. Both uses involve the comparisons of biclusters which may have different numbers of rows and columns. Our findings document a bias that can distort biclustering results. We prove, both analytically and by simulation, that the average H-score increases with the number of rows/columns in a bicluster—even in the ‘ideal’ (and simplest) case of a single bicluster generated by a constant model plus a white noise. This biases the H-score, and hence all H-score based algorithms, toward small biclusters. Importantly, our analytical proof provides a straightforward way to correct this bias. An optimal bicluster is a submatrix with the lowest possible H-score H(I,J)=0 (Madeira and Oliveira, 2004). This corresponds to a bicluster defined by the additive model aij=μ+αi+βj⁠, where μ is the mean of the bicluster, and αi and βj are additive adjustments for rows and columns, respectively. As an example, in the case of gene expression data, such a bicluster is a group of genes with expression levels that tend to fluctuate in unison across a group of conditions. In general, an additive error term εij is also present, leading to H(I,J)>0⁠. The algorithm proposed by Cheng and Church (2000) starts from the entire matrix A, iteratively deletes rows and/or columns which contribute to the H-score the most, and stops when the current submatrix has H(I,J)<δ—a given threshold. This identifies a so-called δ-bicluster. To find additional δ-biclusters, the procedure is repeated after replacing the entries of the prior δ-bicluster(s) with random numbers. Following Cheng and Church (2000), many other biclustering algorithms based on the H-score have appeared in the literature. For example, Yang et al. (2005) proposed FLOC (FLexible Overlapped biClustering), a probabilistic algorithm that simultaneously identifies a set of k (possibly overlapping) biclusters with low H-score. The procedure iteratively reduces the H-scores of k randomly initialized biclusters, until the overall biclustering quality stops improving. Angiulli et al. (2008) proposed an algorithm based on a greedy technique combined with a local search strategy to escape poor local minima. This also employs the H-score, together with the row (gene) variance and the size of the bicluster. The Reactive GRASP Biclustering (RGRASP-B, Dharan and Nair, 2009) also uses the H-score to evaluate bicluster quality, and the algorithm in Bryan et al. (2006) uses a modified version of it. H-scores were recently extended to guide the search of triclusters, i.e. of homogeneous subsets within 3D data (Bhar et al., 2012; Dede and Oğul, 2013). Three-dimensional H-scores have in fact been incorporated in several triclustering algorithms. For instance, Gutiérrez-Avilés and Rubio-Escudero (2014) define triclusters as subsets with H-score below a given coherence strength δ, while Liu et al. (2008) discover triclusters by solving a multiobjective optimization problem where H-scores are one of three merit functions employed. The biclustering and triclustering algorithms cited here are only a small selection of those that rely on H-scores to discover bi- or triclusters and to validate and evaluate results (see e.g. Henriques and Madeira, 2019; Madeira and Oliveira, 2004; Pontes et al., 2015, for comprehensive reviews). (a) H-score as function of the number of rows (columns), in a single bicluster of 10 columns (rows) generated by a constant model plus a white noise aij=μ+ϵij⁠. Red dots indicate average H-scores H¯n⁠. (b) Ratio rn,n+1 between average H-scores, with n=2,…,9 according to simulation (black dot) and Theorem 1 (red cross). (Color version of this figure is available at Bioinformatics online.) The infinite product in (4) converges, and in particular r2,2+m→2⁠. As a consequence, the H-score bias is at most 1, and it becomes small for comparisons between biclusters of large size (Table 1). We also observe that the value H¯n only depends on the variance of the error term εij⁠, and not on the values of μ, αi and βj, nor on the distribution of εij (see Table 1 and Supplementary Figures S1 and S2). Average H-scores and ratios for different models (with/without row and column differential terms), errors (Gaussian, Uniform) and bicluster sizes Average H-scores and ratios for different models (with/without row and column differential terms), errors (Gaussian, Uniform) and bicluster sizes Our results show that employing H-scores to compare biclusters with different numbers of rows or columns could lead to biased results. While this bias is small and likely inconsequential for large biclusters, it can be substantial and rather misleading for small biclusters. Suppose one is comparing biclusters with n and n + m rows; Equation (3) suggests that this bias can be straightforwardly corrected normalizing the H-score ratio by the factor ∏i=nn+m−1i2i2−1 (an identical reasoning holds for columns). Notably, this correction should also be employed to adjust the H-scores and the H-score threshold δ when finding δ-biclusters. In particular, in order to adjust the H-scores of multiple biclusters of different dimensions using (3), one has to select the dimension (n, p) of the reference matrix to be used. The absolute value of the adjusted H-scores will change depending on this choice (see Supplement Section 3 and Supplement Fig. S3) and is thus not very informative per se. However, the ordering of biclusters based on adjusted H-scores is independent of the selected reference dimensions and generally different from the ordering based on unadjusted H-scores (see Supplement Section 3 and Supplement Fig. S4). As a consequence, unadjusted H-scores could lead to erroneous evaluations of biclusters. Delta-based algorithms should adjust all H-scores and define the δ-threshold based on the same reference dimensions (n, p). Although reference matrices can in principle extend their dimensions to the entire dataset, we recommend the use of a reference matrix of dimensions (2, 2) since this choice will work for every dataset, regardless of its size. This choice also directly links the δ-threshold to the level of noise allowed in a bicluster, regardless of the specific dataset analyzed. With this choice, Equation (3) normalizes the H-score of a bicluster with 2+m rows and 2+q columns dividing it by the quantity ∏i=22+m−1i2i2−1∏j=22+q−1j2j2−1⁠. Importantly, both the H-score bias and the proposed correction naturally extend to the triclustering framework. Indeed, in the case of triclusters and 3D H-scores, we obtain the same result as in (3) for all three dimensions (n, p, s). Finally, it is worth noticing that neither H-scores nor adjusted H-scores are measures of the statistical significance of bi- or triclusters, which should be assessed with specific methods and tools, such as BSig (Henriques and Madeira, 2018). Considering the seminal role and ubiquitousness of H-scores in biclustering algorithms, and the importance of biclustering algorithms in bioinformatics and computational biology, we believe the H-score bias should be taken into serious consideration. The correction we propose is simple to implement, and could help shape the conclusions and insights provided by a broad range of applications both in biclustering and in triclustering frameworks. This work was partially funded by the Eberly College of Science, the Institute of CyberScience and the Huck Institutes of the Life Sciences (the Pennsylvania State University); National Science Foundation award DMS-1407639; Tobacco Settlement and CURE funds, PA Department of Health (the Department specifically disclaims responsibility for any analyses, interpretations or conclusions). Conflict of Interest: none declared.

SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Summary: With increased generation of high-resolution sequence-based 'Omics' data, detecting statistically significant effects at different genomic locations and scales has become key to addressing several scientific questions. IWTomics is an R/Bioconductor package (integrated in Galaxy) that, exploiting sophisticated Functional Data Analysis techniques (i.e. statistical techniques that deal with the analysis of curves), allows users to pre-process, visualize and test these data at multiple locations and scales. The package provides a friendly, flexible and complete workflow that can be employed in many genomic and epigenomic applications. Availability and implementation: IWTomics is freely available at the Bioconductor website (http://bioconductor.org/packages/IWTomics) and on the main Galaxy instance (https://usegalaxy.org/). Supplementary information: Supplementary data are available at Bioinformatics online.

Endogenous retroviruses (ERVs), the remnants of retroviral infections in the germ line, occupy ~8% and ~10% of the human and mouse genomes, respectively, and affect their structure, evolution, and function. Yet we still have a limited understanding of how the genomic landscape influences integration and fixation of ERVs. Here we conducted a genome-wide study of the most recently active ERVs in the human and mouse genome. We investigated 826 fixed and 1,065 in vitro HERV-Ks in human, and 1,624 fixed and 242 polymorphic ETns, as well as 3,964 fixed and 1,986 polymorphic IAPs, in mouse. We quantitated >40 human and mouse genomic features (e.g., non-B DNA structure, recombination rates, and histone modifications) in ±32 kb of these ERVs' integration sites and in control regions, and analyzed them using Functional Data Analysis (FDA) methodology. In one of the first applications of FDA in genomics, we identified genomic scales and locations at which these features display their influence, and how they work in concert, to provide signals essential for integration and fixation of ERVs. The investigation of ERVs of different evolutionary ages (young in vitro and polymorphic ERVs, older fixed ERVs) allowed us to disentangle integration vs. fixation preferences. As a result of these analyses, we built a comprehensive model explaining the uneven distribution of ERVs along the genome. We found that ERVs integrate in late-replicating AT-rich regions with abundant microsatellites, mirror repeats, and repressive histone marks. Regions favoring fixation are depleted of genes and evolutionarily conserved elements, and have low recombination rates, reflecting the effects of purifying selection and ectopic recombination removing ERVs from the genome. In addition to providing these biological insights, our study demonstrates the power of exploiting multiple scales and localization with FDA. These powerful techniques are expected to be applicable to many other genomic investigations.

BACKGROUND: Infection with feline immunodeficiency virus (FIV) causes an immunosuppressive disease whose consequences are less severe if cats are co-infected with an attenuated FIV strain (PLV). We use virus diversity measurements, which reflect replication ability and the virus response to various conditions, to test whether diversity of virulent FIV in lymphoid tissues is altered in the presence of PLV. Our data consisted of the 3' half of the FIV genome from three tissues of animals infected with FIV alone, or with FIV and PLV, sequenced by 454 technology. RESULTS: Since rare variants dominate virus populations, we had to carefully distinguish sequence variation from errors due to experimental protocols and sequencing. We considered an exponential-normal convolution model used for background correction of microarray data, and modified it to formulate an error correction approach for minor allele frequencies derived from high-throughput sequencing. Similar to accounting for over-dispersion in counts, this accounts for error-inflated variability in frequencies - and quite effectively reproduces empirically observed distributions. After obtaining error-corrected minor allele frequencies, we applied ANalysis Of VAriance (ANOVA) based on a linear mixed model and found that conserved sites and transition frequencies in FIV genes differ among tissues of dual and single infected cats. Furthermore, analysis of minor allele frequencies at individual FIV genome sites revealed 242 sites significantly affected by infection status (dual vs. single) or infection status by tissue interaction. All together, our results demonstrated a decrease in FIV diversity in bone marrow in the presence of PLV. Importantly, these effects were weakened or undetectable when error correction was performed with other approaches (thresholding of minor allele frequencies; probabilistic clustering of reads). We also queried the data for cytidine deaminase activity on the viral genome, which causes an asymmetric increase in G to A substitutions, but found no evidence for this host defense strategy. CONCLUSIONS: Our error correction approach for minor allele frequencies (more sensitive and computationally efficient than other algorithms) and our statistical treatment of variation (ANOVA) were critical for effective use of high-throughput sequencing data in understanding viral diversity. We found that co-infection with PLV shifts FIV diversity from bone marrow to lymph node and spleen.

Insertions and deletions (indels) cause numerous genetic diseases and lead to pronounced evolutionary differences among genomes. The macaque sequences provide an opportunity to gain insights into the mechanisms generating these mutations on a genome-wide scale by establishing the polarity of indels occurring in the human lineage since its divergence from the chimpanzee. Here we apply novel regression techniques and multiscale analyses to demonstrate an extensive regional indel rate variation stemming from local fluctuations in divergence, GC content, male and female recombination rates, proximity to telomeres, and other genomic factors. We find that both replication and, surprisingly, recombination are significantly associated with the occurrence of small indels. Intriguingly, the relative inputs of replication versus recombination differ between insertions and deletions, thus the two types of mutations are likely guided in part by distinct mechanisms. Namely, insertions are more strongly associated with factors linked to recombination, while deletions are mostly associated with replication-related features. Indel as a term misleadingly groups the two types of mutations together by their effect on a sequence alignment. However, here we establish that the correct identification of a small gap as an insertion or a deletion (by use of an outgroup) is crucial to determining its mechanism of origin. In addition to providing novel insights into insertion and deletion mutagenesis, these results will assist in gap penalty modeling and eventually lead to more reliable genomic alignments.