BACKGROUND: Protein-DNA binding sites of ChIP-seq experiments are identified where the binding affinity is significant based on a given threshold. The choice of the threshold is a trade-off between conservative region identification and discarding weak, but true binding sites. RESULTS: We rescue weak binding sites using MSPC, which efficiently exploits replicates to lower the threshold required to identify a site while keeping a low false-positive rate, and we compare it to IDR, a widely used post-processing method for identifying highly reproducible peaks across replicates. We observe several master transcription regulators (e.g., SP1 and GATA3) and HDAC2-GATA1 regulatory networks on rescued regions in K562 cell line. CONCLUSIONS: We argue the biological relevance of weak binding sites and the information they add when rescued by MSPC. An implementation of the proposed extended MSPC methodology and the scripts to reproduce the performed analysis are freely available at https://genometric.github.io/MSPC/ ; MSPC is distributed as a command-line application and an R package available from Bioconductor ( https://doi.org/doi:10.18129/B9.bioc.rmspc ).

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: ChIP-seq experiments are widely used to detect and study DNA-protein interactions, such as transcription factor binding and chromatin modifications. However, downstream analysis of ChIP-seq data is currently restricted to the evaluation of signal intensity and the detection of enriched regions (peaks) in the genome. Other features of peak shape are almost always neglected, despite the remarkable differences shown by ChIP-seq for different proteins, as well as by distinct regions in a single experiment. RESULTS: We hypothesize that statistically significant differences in peak shape might have a functional role and a biological meaning. Thus, we design five indices able to summarize peak shapes and we employ multivariate clustering techniques to divide peaks into groups according to both their complexity and the intensity of their coverage function. In addition, our novel analysis pipeline employs a range of statistical and bioinformatics techniques to relate the obtained peak shapes to several independent genomic datasets, including other genome-wide protein-DNA maps and gene expression experiments. To clarify the meaning of peak shape, we apply our methodology to the study of the erythroid transcription factor GATA-1 in K562 cell line and in megakaryocytes. CONCLUSIONS: Our study demonstrates that ChIP-seq profiles include information regarding the binding of other proteins beside the one used for precipitation. In particular, peak shape provides new insights into cooperative transcriptional regulation and is correlated to gene expression.