## 4 Discussion  
The proposed method has several advantages. First, it is a fully linear method that simply applies PCA or HOSVD to matrices and tensors. The only time-consuming process is SD optimization, which typically ends in a few minutes. Second, the proposed method accepts the bed or the bigWig file formats. In contrast to bam/sam files that retain information about how individual short reads are mapped onto the genome, bed or bigWig files do not keep this information and instead retain which regions of the short reads are mapped. Third, it is robust. Independent of samples, species, and platforms, the proposed method achieved similar performance levels. For example, the H3K27 me3 modification is considered three times in Table 5; two are human ChIP-seq, and one is mouse ChIP-seq. Despite the different species, all three are associated with a similar number of gene symbols, 5208, 4994, and 3545, and they are all associated with the same number of enriched  
Figure 10: Histogram of \(1-P_{i}\) computed by the proposed method using the GSE188173 dataset (H3K27 ac). H3K27 me3 profiles in the \"Epigenomics Roadmap HM ChIP-seq\" category of Enrichr. H3K27 ac modifications were considered three times with two distinct protocols, ChIP-seq and CUT&Tag. A similar number of gene symbols, 11,590, 13,061, and 15,548, and the same number of H3K27 ac-enriched profiles, 24, in the \"Epigenomics Roadmap HM ChIP-seq\" category of Enrichr were associated. Usually, although CUT&Tag requires protocols specific to CUT&Tag [19], the proposed method handled both ChIP-seq and CUT&Tag seamlessly. In addition, as seen in Tables 4 and 5, the proposed method successfully deals with all five \"core histone marks\" (H3K27 ac, H3K4 me3, H3K36 me3, H3K27 me3, and H3K9 me3). Unfortunately, H3K4 me1 cannot be properly dealt with using the proposed method. However, this is not a critical problem, since it can process H3K27 ac modifications, which are strongly associated with H3K4 me1 modifications [37]. Thus, the proposed method is still worth investigating. In the above, \(u_{\ell i}\) does not obey a Gaussian distribution, but instead fits a combination of two Gaussian distributions. To confirm these points, the analysis of the Gaussian mixture to \(u_{1i}\) of the GSE24850 dataset was applied (whose histogram \(1-P_{i}\) is shown in Figure 2). Gaussian mixture analysis was performed using the mclust package [38] in R [30]. The mclust function was applied to \(u_{1i}\) for GSE24850. Then, the Bayesian information criterion of the Gaussian mixture drastically improves when the number of clusters increases from one to two and does not improve for more numbers of clusters. Table 6 shows the confusion matrix between clustering and selection of the proposed method when the number of assumed clusters varies from two to four. It is obvious that they are in agreement only when the number of clusters is assumed to be two, as expected. All regions not selected by the proposed method are in one cluster, and the majority of the selected regions belong to another cluster, whereas this does not occur when the number of assumed clusters is either three or four. Thus, as expected, \(u_{\ell i}\) did not obey a Gaussian distribution, but a mixture of two Gaussian distributions. Since the optimization of SD assumes a single Gaussian distribution, it is not expected to work well, but it does work well empirically (Table 5). The limitations of the proposed methods applied to histone modification must be investigated in the future. The proposed method has some weaknesses. It does not work well if it is applied to non-core histone marks. For example, its performance is reduced if it is applied to H3K9 ac modifications. More detailed conditions in which the proposed method may be successfully applied to histone modification need to be investigated. In addition, we do not know how the proposed method outperforms the previous state-of-the-art methods. It is important to note that the proposed method is an unsupervised method.