==== Front PLoS Comput Biol PLoS Comput Biol plos ploscomp PLoS Computational Biology 1553-734X 1553-7358 Public Library of Science San Francisco, CA USA 33253153 PCOMPBIOL-D-20-00774 10.1371/journal.pcbi.1008422 Research Article Physical Sciences Mathematics Applied Mathematics Algorithms Research and Analysis Methods Simulation and Modeling Algorithms Medicine and Health Sciences Oncology Cancers and Neoplasms Hematologic Cancers and Related Disorders Leukemia Medicine and Health Sciences Hematology Hematologic Cancers and Related Disorders Leukemia Research and Analysis Methods Mathematical and Statistical Techniques Cluster Analysis Hierarchical Clustering Biology and Life Sciences Cell Biology Chromosome Biology Chromatin Biology and Life Sciences Genetics Epigenetics Chromatin Biology and Life Sciences Genetics Gene Expression Chromatin Medicine and Health Sciences Oncology Cancers and Neoplasms Hematologic Cancers and Related Disorders Leukemia Myeloid Leukemia Acute Myeloid Leukemia Medicine and Health Sciences Hematology Hematologic Cancers and Related Disorders Leukemia Myeloid Leukemia Acute Myeloid Leukemia Biology and Life Sciences Genetics Genomics Medicine and Health Sciences Oncology Cancers and Neoplasms Hematologic Cancers and Related Disorders Leukemia Lymphoblastic Leukemia Chronic Lymphoblastic Leukemia Medicine and Health Sciences Hematology Hematologic Cancers and Related Disorders Leukemia Lymphoblastic Leukemia Chronic Lymphoblastic Leukemia Biology and Life Sciences Cell Biology Cellular Types Animal Cells Immune Cells Antibody-Producing Cells B Cells Biology and Life Sciences Immunology Immune Cells Antibody-Producing Cells B Cells Medicine and Health Sciences Immunology Immune Cells Antibody-Producing Cells B Cells Biology and Life Sciences Cell Biology Cellular Types Animal Cells Blood Cells White Blood Cells B Cells Biology and Life Sciences Cell Biology Cellular Types Animal Cells Immune Cells White Blood Cells B Cells Biology and Life Sciences Immunology Immune Cells White Blood Cells B Cells Medicine and Health Sciences Immunology Immune Cells White Blood Cells B Cells Systematic clustering algorithm for chromatin accessibility data and its application to hematopoietic cells Systematic clustering algorithm for chromatin accessibility datahttps://orcid.org/0000-0002-0503-4176Tanaka Azusa ConceptualizationData curationFormal analysisFunding acquisitionInvestigationMethodologyProject administrationResourcesSoftwareSupervisionValidationVisualizationWriting – original draftWriting – review & editing12* https://orcid.org/0000-0002-4461-6152Ishitsuka Yasuhiro ConceptualizationFormal analysisMethodologyProject administrationSoftwareSupervisionValidationVisualizationWriting – original draftWriting – review & editing34* https://orcid.org/0000-0002-8815-9037Ohta Hiroki ConceptualizationFormal analysisMethodologyProject administrationSoftwareSupervisionValidationVisualizationWriting – original draftWriting – review & editing35* https://orcid.org/0000-0002-0075-0800Fujimoto Akihiro ConceptualizationProject administrationSupervisionWriting – review & editing1 https://orcid.org/0000-0002-7939-2080Yasunaga Jun-ichirou Funding acquisitionResourcesWriting – review & editing26 https://orcid.org/0000-0002-0473-754XMatsuoka Masao Funding acquisitionResourcesWriting – review & editing26 1 Department of Human Genetics, Graduate School of Medicine, The University of Tokyo, Tokyo, Japan 2 Laboratory of Virus Control, Institute for Frontier Life and Medical Sciences, Kyoto University, Kyoto, Japan 3 Center for Science Adventure and Collaborative Research Advancement, Graduate School of Science, Kyoto University, Kyoto, Japan 4 Department of Mathematics, Graduate School of Science, Kyoto University, Kyoto, Japan 5 Department of Physics, Graduate School of Science, Kyoto University, Kyoto, Japan 6 Department of Hematology, Rheumatology and Infectious Disease, Faculty of Life Sciences, Kumamoto University, Kumamoto, Japan Schlessinger Avner Editor Icahn School of Medicine at Mount Sinai, UNITED STATES The authors declare that they have no conflict of interest. * E-mail: a-tanaka@m.u-tokyo.ac.jp (AT); yasu-ishi@math.kyoto-u.ac.jp (YI); ohta.hiroki.6c@kyoto-u.ac.jp (HO) 11 2020 30 11 2020 16 11 e10084227 5 2020 6 10 2020 © 2020 Tanaka et al2020Tanaka et alThis is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.The huge amount of data acquired by high-throughput sequencing requires data reduction for effective analysis. Here we give a clustering algorithm for genome-wide open chromatin data using a new data reduction method. This method regards the genome as a string of 1s and 0s based on a set of peaks and calculates the Hamming distances between the strings. This algorithm with the systematically optimized set of peaks enables us to quantitatively evaluate differences between samples of hematopoietic cells and classify cell types, potentially leading to a better understanding of leukemia pathogenesis. Author summary High-throughput sequencing provides us huge amounts of data about gene regulation. In order to extract useful information from the data, data reduction is needed. Although RNA-seq data analysis has been extensively studied, where the focus is mainly on genetic loci, tools for epigenetic sequencing data, such as ATAC-seq data which represent chromatin accessibility, are comparatively lacking. Since the binding of transcription factors mainly occurs in open chromatin regions, it is presumably important to understand how chromatin accessibility landscape affects cell phenotype. In this context, we developed a systematic algorithm to select a set of peaks representing the open state of chromatin for a given sample of ATAC-seq data. This algorithm quantifies the difference between samples by regarding the genome as a string of 1s and 0s with Hamming distances and then performs hierarchical clustering. This algorithm has less computational cost and gives a reasonable cell type classification compared to a previous method. In this work, as an application of this algorithm, we present a comparative analysis of leukemia samples with healthy hematopoietic cells and provide new insights about the relationship between chromatin structures, cell surface proteins, and symptoms in leukemia. Japan Society for the Promotion of Science (JP)JP19K16740https://orcid.org/0000-0002-0503-4176Tanaka Azusa http://dx.doi.org/10.13039/501100001691Japan Society for the Promotion of ScienceJP18J40119https://orcid.org/0000-0002-0503-4176Tanaka Azusa Japan Society for the Promotion of Science (JP)JP19H03689https://orcid.org/0000-0002-0473-754XMatsuoka Masao Japan Society for the Promotion of Science (JP)JP20H03514https://orcid.org/0000-0002-7939-2080Yasunaga Jun-ichirou Japan Agency for Medical Research and DevelopmentJP20fk018088h0002https://orcid.org/0000-0002-0473-754XMatsuoka Masao http://dx.doi.org/10.13039/100009619Japan Agency for Medical Research and DevelopmentJP17km0405207h0002https://orcid.org/0000-0002-0075-0800Fujimoto Akihiro http://dx.doi.org/10.13039/100009619Japan Agency for Medical Research and DevelopmentJP18km0405207S0103https://orcid.org/0000-0002-0075-0800Fujimoto Akihiro http://dx.doi.org/10.13039/100007428Naito Foundationhttps://orcid.org/0000-0002-0503-4176Tanaka Azusa This research was supported by JSPS KAKENHI Grant Numbers JP19K16740 (AT), JP18J40119 (AT), JP19H03689 (MM), JP20H03514 (JiY), and by Japan Agency for Medical Research and Development (AMED) Grant Numbers JP20fk0108088h0002 (MM), JP17km0405207h0002 (AF), JP18km0405207S0103 (AF), and by a grant from the Naito Foundation (AT). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. PLOS Publication Stagevor-update-to-uncorrected-proofPublication Update2020-12-10Data AvailabilityAll ATAC-seq and RNA-seq data needed to reproduce this study have been deposited at the DNA Data Bank of Japan (DDBJ) under accession number DRA010939. The source code is available from https://github.com/tanakanishi/findclosest.Data Availability All ATAC-seq and RNA-seq data needed to reproduce this study have been deposited at the DNA Data Bank of Japan (DDBJ) under accession number DRA010939. The source code is available from https://github.com/tanakanishi/findclosest. ==== Body Introduction Cellular phenotypes are governed by epigenetic mechanisms. For example, information about how human DNA is packed and chemically modified in the nucleus plays an important role in understanding the differentiation and regulation of cells [1–4]. Methods such as chromatin immunoprecipitation sequencing (ChIP-seq) and assay for transposase accessible chromatin using sequencing (ATAC-seq) have proven useful for understanding the modification and detection of open chromatin on a genome-wide scale [5–9]. Those epigenetic data analysis methods usually start with data enrichment along the whole genome, also known as “peak calling” [10, 11]. Compared to RNA-seq data analysis, whose target regions are mainly in certain loci or genes across samples, the target regions on epigenetic sequencing data are undetermined. To determine the target regions, peak calling with an appropriate tool is often performed for the entire genome of every sample, and the target regions are defined as merged peaks among all samples. Then the total number of reads or fragments present in each region is counted for each sample, leading to a matrix, X = (xi,j), where xi,j represents the number of reads/fragments from sample i in region j. The matrix elements are normalized by quantile normalization to reduce the biases arising from variations in the data size over samples, followed by downstream processing [7–9]. However, this process raises two concerns. First, we do not fully understand the effect of merging all the peaks from different samples. For example, if two peaks from different samples slightly overlap, those two peaks are considered as one peak after the peak merging step. Therefore, the difference of the two peak positions, which may reflect cell identity, may be unintentionally ignored. The second concern is that we have no justification for applying quantile normalization over samples that are phenotypically different [12, 13]. Thus, the aim of the present study is to avoid these concerns by constructing an algorithm that systematically classifies epigenetic data obtained from high-throughput sequencing. In this analysis, toward cell type classification, we provide a systematic algorithm to select a set of peaks used for the downstream analysis, where the difference between samples are quantified by using the Hamming distance from information theory [14]. This algorithm has less computational cost while still producing reasonable classification compared to a previous method [7]. As an application of the developed algorithm, we use it to obtain new insights on samples of leukemia cells from chronic lymphocytic leukemia (CLL), acute myeloid leukemia (AML), and adult T-cell leukemia (ATL) at the chromatin level. In particular, using this algorithm, we infer the phenotype of a given leukemia sample as output by using only ATAC-seq data of that sample as input. Results ATAC-seq samples In this paper, we mainly focused on 77 ATAC-seq datasets from 13 human primary blood cell types [7] as test data. The 13 cell types are comprised of hematopoietic stem cells (HSC), multipotent progenitor cells (MPP), lymphoid-primed multipotent progenitor cells (LMPP), common myeloid progenitor cells (CMP), megakaryocyte-erythroid progenitor cells (MEP), granulocyte-macrophage progenitor cells (GMP), common lymphoid progenitor cells (CLP), natural killer cells (NK), B cells, CD4+T cells (CD4+T), CD8+T cells (CD8+T), monocytes (Mono) and erythroids (Ery). These cell types are experimentally categorized by immunophenotypes described by the combination of cell surface markers shown in Table 1. 10.1371/journal.pcbi.1008422.t001Table 1 Immunophenotypes of samples. Types of hematopoietic cells and their corresponding cell surface markers in [7]. For example, CD34+ and CD38- for cell type ν means that a cell of type ν expresses CD34 but not CD38 at its surface. Cell type (ν) Number of replicates Immunophenotypes HSC 7 Lin-, CD34+, CD38-, CD10-, CD90+ MPP 6 Lin-, CD34+, CD38-, CD10-, CD90- LMPP 3 Lin-, CD34+, CD38-, CD10-, CD45RA+ CMP 8 Lin-, CD34+, CD38+, CD10-, CD45RA-, CD123+ MEP 7 Lin-, CD34+, CD38+, CD10-, CD45RA-, CD123- GMP 7 Lin-, CD34+, CD38+, CD10-, CD45RA+, CD123+ CLP 5 Lin-, CD34+, CD38+, CD10+, CD45RA+ NK 6 CD56+ B 4 CD19+, CD20+ CD4+T 5 CD3+, CD4+ CD8+T 5 CD3+, CD8+ Mono 6 CD14+ Ery 8 CD71+, GPA+, CD45-low For convenience, T denotes a set of the thirteen cell types; T={B,CD4+T,CD8+T,CLP,CMP,Ery,GMP,HSC,LMPP,MEP,Mono,MPP,NK}. For all 77 samples, we assigned ATAC-seq reads to reference genome hg19 (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/), and among them only those which had high mapping quality values (MQ ≥ 30) were used for the peak calling by MACS2 (see S1 Appendix for details of the preprocessing) [15]. The peak calling results consisted of the location with a peak width and the associated p-value. Concretely, the location of the k-th peak is expressed by gk = (γk, αk, βk), where γk is the chromosome number, αk is the start position, and βk is the end position. Note that we used MACS2 to call all ATAC-seq peaks with the following parameters (--nomodel --nolambda --keep-dup all -p pG), where the number of peaks is affected by the peak calling parameter ‘-p pG’. The parameter pG is larger than any p-values of the peak calling results. (See Materials and methods for details of the peak-calling). Note that the peak position depends on parameter pG of the MACS2 algorithm as shown in Fig 1. For example, the start and end positions of a peak could change and one peak could split into two peaks depending on pG. Thus, we need to take into account the dependence of a set of peaks on different values of pG for careful analysis. 10.1371/journal.pcbi.1008422.g001Fig 1 The number of reads vs genomic positions. The plots show representative data of Mono obtained from SRA with accession number SRR2920475. (A) The number of reads Yx at each position x along chr 1 (γ = 1) and the peak region (αk, βk) as determined by the MACS2 algorithm with peak calling parameter pG = 10−2 (pink shaded regions) is shown. The peak region and its associated p-value ((αk, βk), pk) are (1092756, 1094068, 10−20.36428). (B) The obtained peak regions are ((1092817, 1093330), 10−20.36428) and ((1093480, 1094025), 10−8.19447) for pG = 10−4. Parameterized binarization First we ranked the peak results in the order of ascending p-values and then investigated the relationship between the peak width and the corresponding ranking. We found that as the p-value increased, the width of the ATAC-seq peaks became shorter statistically, which suggested the feasibility of robust data reduction against small noise in the data by selecting peaks with smaller p-values (Fig 2). 10.1371/journal.pcbi.1008422.g002Fig 2 The statistics of peak width. Distribution of peak width (βk − αk) and its corresponding ranking k obtained from the peak calling result of CD4+T cells with peak calling parameter pG = 10−2. The bin size is 400 × 400. The color code indicates the number of data in each bin. Thus, we define Mcut as the threshold such that only peaks with rankings not greater than Mcut are used for the analysis hereafter. Then, for a given set of (Mcut, pG), we introduce B = {hγ,x}, where hγ,x = 1 when position x in chromosome γ is inside a peak and 0 otherwise (Fig 3). The process to obtain the binary sequence from the reads data is illustrated in Fig 4. Note that we do not perform any coarse-grained description for the genome position x but keep 1bp resolution. (See Materials and methods for details of the binarization). 10.1371/journal.pcbi.1008422.g003Fig 3 How to calculate Hamming distance. Schema of the Hamming distance calculation from the peak locations with two samples c1,c2∈S. Each locus is converted to 1 or 0 based on the peak overlapping status. 10.1371/journal.pcbi.1008422.g004Fig 4 Binarizing the number of reads. (A) The number of reads Yx at each position x along chr 3 (γ = 3) and the peak region (αk, βk) as determined by the MACS2 algorithm with peak calling parameter pG = 10−2 (pink shaded regions). This figure shows representative data of NK cells obtained from SRA with accession number SRR2920495. The peak regions and the associated p-values ((αk, βk), pk) in the left and right peaks are ((188271079, 188271985), 10−422.5872) and ((188286401, 188287077), 10−329.52139), respectively. Thus, the width of the peaks (βk − αk) in the left- and right-hand sides are 906 and 676, respectively. (B) Binary sequence (hx) as determined by the peak regions seen in (A) when we chose Mcut satisfying pMcut≥10−329.52139. Quantifying differences between two binary sequences by Hamming distance Let us move onto the situation when one considers a set of samples to evaluate the difference between two binary sequences B. Here our strategy is to find the proper distance that can be measured from the normalized ATAC-seq data of two samples. Using that distance, we try to obtain hierarchical clustering of a set of hematopoietic cell samples to quantitatively characterize the relationship among those samples. Let Ns be the number of samples. We then write the set of samples as S:={1,2,…,Ns}, where Ns = 77 in this study. For sample c∈S, we add index c to related objects as a superscript. For example, we write a binary sequence B associated to sample c as Bc:={hγ,xc}. There are many methods to evaluate the difference between a binary sequence Bc from sample c∈S and Bc′ from sample c′∈S. In this paper, we evaluated the difference between two samples (c, c′) by using the Hamming distance H(Bc,Bc′) between two binary sequences, Bc and Bc′. H(Bc,Bc′) is calculated as the sum of the number of pairs with different values at every position x between Bc and Bc′ (Fig 5). We used the distance as an initial condition for the hierarchical clustering and then used Ward’s method to complete the hierarchical clustering [16]. Examples of hierarchical clustering with (Mcut, pG) = (2000, 10−2) and (80000, 10−2) are shown in Fig 6. (See Materials and methods for details of the Hamming distance and hierarchical clustering). 10.1371/journal.pcbi.1008422.g005Fig 5 Matrix of Hamming distances. Matrix of Hamming distances dij between samples i and j. This matrix is used for the downstream analysis. 10.1371/journal.pcbi.1008422.g006Fig 6 Examples of clustering dendrograms. Hierarchical clustering obtained by Ward’s method with parameters (Mcut, pG) = (2000, 10-2) (A) and (80000, 10−2) (B). Optimization of hierarchical clustering toward cell-type classification By using the methods explained above, we can obtain a clustering dendrogram that depends on (Mcut, pG). We then need to systematically determine the best clustering, which is the clustering closest to the “perfectly classified dendrogram” where each set Sν of all samples with type ν∈T coincides with an offspring set. This condition can be restated as an optimization problem by introducing a cost function “penalty” for the performance of clustering as follows. Concretely, to quantitatively evaluate the obtained dendrogram for each combination of (Mcut, pG), we define type penalty λν for a given cell type ν∈T. Type penalty λν corresponds to the number of samples from different cell types in cluster ν formed when all samples of cell type ν meet together from the bottom of the dendrogram (Fig 7). Additionally, we define global penalty λ:=∑ν∈Tλν as the “cost function” of the optimization. Note that λ ≥ 0, and a “perfectly classified dendrogram” gives λ = 0. (See Materials and methods for details of the penalty). 10.1371/journal.pcbi.1008422.g007Fig 7 Schema of penalty score calculation. Note that this dendrogram is constructed by artificial data to explain how to calculate the penalty, though we use the same labels such as HSC1. This dendrogram has six leaves, and three of them are classified to type HSC. To explain details of this dendrogram, we freely use the symbols and definitions in Materials and methods in this caption. We can see that τ(HSC) = 10. The corresponding node is n10 (displayed by the blue dot), and the corresponding cluster C10 is the set {HSC1, HSC2, HSC3, MPP} (surrounded by the blue dashed line). Among the elements of C10, one leaf, MPP, is not in type HSC, but the three others are. Hence, the type penalty of HSC in this figure is computed as λHSC = 4 − 3 = 1. Determination of the best parameters for the optimization As mentioned above, the optimization problem we have to solve is to find (Mcut*,pG*) that minimizes the cost function λ(Mcut, pG). The schematic workflow in our algorithm is shown in Fig 8. 10.1371/journal.pcbi.1008422.g008Fig 8 Schematic workflow of our algorithm. See Materials and methods for details. First we took into account all the peaks by setting Mcut = ∞ and checked how the dendrograms and λ(∞, pG) depended on pG, as shown in Fig 9. Considering the tendency of the parameter searching, we concluded that 1.5≤-log10pG*≤4. 10.1371/journal.pcbi.1008422.g009Fig 9 Global penalty without cutoff of reads. Global penalty λ(Mcut = ∞, pG) obtained by Ward’s method. We then sought the best parameters to optimize the dendrograms and found that (Mcut*,pG*) was close to (64000, 10−2), which gave the smallest penalty λ in our searching resolution, as shown in Figs 10 and 11. Note that 64000 is the midpoint of (60000, 62000, 64000, 66000, 68000) which give the same minimum penalty in our searching resolution. Hereafter, to investigate the property of the best clustering, we set (Mcut*,pG*) as (64000, 10−2). In our searching resolution, the increment in terms of Mcut was 2000 near Mcut = 64000. Note that more-refined resolutions might give better estimates of the optimized value (Mcut*,pG*), but naturally the computational costs get higher. Even then, the following procedures are operationally unchanged. 10.1371/journal.pcbi.1008422.g010Fig 10 Penalty with cutoff of reads. The distribution of global penalty λ (A) and type penalty λν for each cell type ν (B) along with Mcut with parameter pG = 10−2 by Ward’s method. 10.1371/journal.pcbi.1008422.g011Fig 11 Our best clustering dendrogram. Hierarchical clustering obtained by Ward’s method with (Mcut, pG) = (64000, 10−2). The value of the minimum penalty achieved at (Mcut*,pG*) was 18. This minimum was smaller than the penalty value of 27 for the clustering of the data from GSE74912_ATACseq_All_Counts.txt in [7]. The procedure of the latter clustering was as follows. First we performed a quantile normalization of the reads count in the distal elements (> 1000 bp away from a transcription start site (TSS)). Then we calculated the Pearson coefficients over all samples leading to a distance matrix where each entry is 1-(Pearson coefficient). By using Ward’s method, we finally obtained the clustering dendrogram. Note that for this case, Ward’s method gives penalty λ = 27 and UPGMA gives λ = 29. Computational cost of the algorithm As explained above, after obtaining data of the reads positions, we perform the MACS2 algorithm to get peak regions, and then finally we produce a hierarchical clustering. Here we consider the computational cost of our algorithm after acquiring the data of the reads positions and until acquiring a distance matrix to produce the hierarchical clustering. Note that the computational cost of the MACS2 algorithm is not more than O(Ns), where O() is the Landau notation and Ns is the total number of samples. We consider two situations. (i) One is the case where new samples to analyze are given. (ii) The other is the case where one new sample to analyze is added to the already analyzed samples, for which peak regions and the distance matrix are already calculated. For case (ii), we use the symbol Ns to write the total number of already analyzed samples. We claim that the computational cost of our algorithm is significantly lower than that of a previous method using target regions merged over samples [7] for large values of Ns for case (ii) and, in our case with Ns = 77, that the computational cost of our algorithm is practically lower for case (i). Specifically, in case (i) for our algorithm, the corresponding computational cost is K1McutNs2, which comes solely from the calculation of the Hamming distance. In case (ii), the corresponding computational cost is K2 Mcut Ns, which also comes solely from the calculation of the Hamming distance. Note that K1 and K2 are constants that do not depend on Mcut or Ns. In the context of estimating the best optimization parameter Mcut*, by using Mm different values for Mcut, the computational cost becomes K1McutMmNs2 for case (i) and K2 Mcut MmNs for case (ii), where Mm does not depend on Ns or genome size L and can be adjusted according to the searching resolution of the optimization. Note that K1 and K2 do not depend on Mm. In addition, we optimize pG by Mp different values for pG. Since this optimization can be done for any algorithm, we do not take into account this cost for the comparison of different algorithms. Typically, we set (Mm, Mp) ≃ (30, 10) in our optimization corresponding to case (i). Note that in the section of “Application to leukemic cells” discussed later, corresponding to case (ii), we use the optimized parameters (Mcut,pG)=(Mcut*,pG*), leading to (Mm, Mp) = (1, 1). The previous method using targeted regions merged over samples in [7] includes (a) the merging of reads before peak calling and (b) calculating the distance matrix by the Pearson coefficients which automatically depend on Ns. Thus, for a given number Nnew of unanalyzed samples, the computational cost corresponding to the process of (a) and (b) is at least KrNrNnew+KLL1Ns2, where Nr is the minimum reads number over all samples, and L1 is the number of target regions merged over all samples. The first term comes from counting the reads and the second term comes from calculating the distance matrix. Note that Kr is a constant that does not depend on Nr or Nnew, and KL is a constant that does not depend on L1 or Ns. This form of the computational cost KrNrNnew+KLL1Ns2 is the same for case (i) with Nnew = Ns and case (ii) with Nnew = 1, leading to the conclusion that the computational cost of our algorithm is significantly lower than the previous method, especially for case (ii) with sufficiently large Ns. We do not have the exact estimate of the coefficients K1, K2, Kr, KL, but because Nr=3265006≫Mcut* and L1=590650≫Mcut* in our case, then KrNrNnew+KLL1Ns2 could be costly compared to K1McutNs2. In practice, even in case (i) with Ns = 77, we numerically found that the computational cost of our algorithm is lower due to our algorithm not using the process of merging reads unlike [7]. How to relate the best parameters to genomic context In order to understand why ATAC-seq data under the condition of (Mcut, pG) = (64000, 10−2) was well classified, we analyzed the properties of the peaks with higher rankings. The result of the previous section suggested that peaks of {gk}k=1Mcut* with Mcut*=64000 included key regions for characterizing cell types. Therefore, we investigated which functional genomic regions such as promoters, enhancers, etc. are dominantly related to these top 64000 peaks. Functional annotation of peaks depending on rank In order to investigate functional annotations on the genome overlap with ATAC-seq peaks data, we applied the top 80000 peaks in three cell types (HSC, B cells, and Mono) to the 15-state ChromHMMmodel data. One can obtain data of the biological functions on the genome for HSC, B cells, and Mono from an integrative analysis of 111 reference human epigenome datasets, where we used the data of E032 for B cells, E035 for HSC, and E029 for Mono (https://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/coreMarks/jointModel/final/) [17]. ATAC-seq peaks were ranked according to p-values and divided into groups consisting of 1000 peaks. Then we calculated the average ratio and the standard deviation for each of the 15 states over all samples in each cell type. For an explicit description, let us introduce a set of functional annotations, W:={Wy}y=115, where Wy is the set of regions on the genome, each of which corresponds to functional annotation y. We want to know how many peaks, k, of every 1000 peaks belong to each functional annotation y. For this purpose, we define Exy:={x≤k