Artificial Intelligence In Medicine
 
Analysis of Differential Intra-chromosomal Community Interactions: a Case Study of
Breast Cancer
--Manuscript Draft--
 
Manuscript Number: AIIM-D-23-02555
Article Type: Original research contributions
Keywords: network clustering,
intra-chromosomal community interaction,
estrogen-treated tamoxifen-sensitive breast cancer cell line,
Hi-C data,
epigenetic data,
integrative data analysis
Manuscript Region of Origin: Europe
Abstract: It is challenging to analyze dynamics of intra-chromosomal interactions when
considering multiple high-dimensional epigenetic datasets. A computational integrative
approach was proposed here to tackle the problem by integrating Hi-C data with other
genomic data. A novel hyperparameter tuning method facilitating network clustering is
developed to identify valid intra-chromosomal communities at different interaction
resolutions. The approach was first trained on Hi-C data and other epigenetic data in
an untreated and one hour estrogen (E2)-treated tamoxifen-sensitive (TS) breast
cancer cell line, MCF7, and revealed two major types of valid intra-chromosomal
communities - active/repressive compartments, respectively, which resembles the
properties of open/closed chromatin domains. It was further tested on MCF7 and
tamoxifen-resistant (TR) breast cancer cell line, MCF7TR, and identified 515
differential interacting and differentially expressed genes (DIEGs) within intra-
chromosomal community interactions.  Further analysis of these DIEGs revealed that
endocrine resistance is among the top biological pathways, suggesting a looping-
mediated mechanism in regulating breast cancer tamoxifen resistance. This novel
computational approach may have a broader application in many other biological
systems for identifying a biological-specific differential community interaction.
Powered by Editorial Manager® and ProduXion Manager® from Aries Systems CorporationHighlights 
  Use network clustering algorithm to fi nd community structure of intra -chromosomal 
interacti ons  
  Dimensional  reducti on  and d iff erenti al  analysis  of  intra -chromosomal  community 
interacti ons  
  Integrati ve analysis of Hi-C data and epigeneti c data in estrogen-treated tamoxifen-
sensiti ve breast cancer cell line s 
  Two major types of intra-chromosomal communiti es (acti ve/repressive ) resemble the 
properti es of (open/closed) chromati n domains  
 
Highlights  1
Analysis of Differential Intra-chromosomal Community Interactions: a 
Case Study of Breast Cancer  
 
Zhihao Yao 1,5 , Kun Fang 2,3,4 , Jenny Sofie Dragland 6 , Magnar Bjørås 5,7 , Victor X. 
Jin 2,3,4* , Junbai Wang 1 * 
 
1.  Department of Clinical Molecular Biology (EpiGen), Akershus University Hospital 
and University of Oslo, Lørenskog, Norway 
2.  Division of Biostatistics, Medical College of Wisconsin, Milwaukee, WI 53226, USA 
3.  MCW Cancer Center, Medical College of Wisconsin, Milwaukee, WI 53226, USA 
4.  Mellowes Center for Genomic Sciences and Precision Medicine, Medical College of 
Wisconsin, Milwaukee, WI 53226, USA 
5.  Department of Microbiology, Oslo University Hospital and University of Oslo, Oslo, 
Norway  
6.  Department of Pathology, Oslo University Hospital - Norwegian Radium Hospital, 
Oslo, Norway 
7.  Department of Clinical and Molecular Medicine, Norwegian University of Science 
and Technology, Trondheim, Norway 
 
*To whom correspondence should be addressed. Email: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS> or 
<PRESIDIO_ANONYMIZED_EMAIL_ADDRESS> 
 
 
 
 
 
Manuscript File Click here to view linked References  2
 
 
Abstract  
It  is  challenging  to  analyze  dynamics  of  intra-chromosomal  interactions  when  considering 
multiple  high-dimensional  epigenetic  datasets.  A  computational  integrative  approach  was 
proposed here to tackle the problem by integrating Hi-C data with other genomic data. A novel 
hyperparameter tuning method facilitating network clustering is developed to identify valid 
intra-chromosomal communities at different interaction resolutions. The approach was first 
trained on Hi-C data and other epigenetic data in an untreated and one hour estrogen (E2)-
treated tamoxifen-sensitive (TS) breast cancer cell line, MCF7, and revealed two major types 
of valid intra-chromosomal communities - active/repressive compartments, respectively, which 
resembles the properties of open/closed chromatin domains. It was further tested on MCF7 and 
tamoxifen-resistant  (TR)  breast  cancer  cell  line,  MCF7TR,  and  identified  515  differential 
interacting and differentially expressed genes (DIEGs) within intra-chromosomal community 
interactions.  Further analysis of these DIEGs revealed that endocrine resistance is among the 
top biological pathways, suggesting a looping-mediated mechanism in regulating breast cancer 
tamoxifen resistance. This novel computational approach may have a broader application in 
many  other  biological  systems  for  identifying  a  biological-specific  differential  community 
interaction.  
 
 
 
 
 
 
   3
Introduction 
Over the past decade, the advancement of both genomic profiling and computational analysis 
has revolutionized the field of three-dimensional (3D) chromatin architecture[1-7]. The human 
genome is now understood to function as a complex hierarchy of folded high-order chromatin 
units,  such  as  chromosomal  territories,  compartments,  topologically  associating  domains 
(TADs), and chromatin loops[5, 8]. Chromatin or DNA looping is recognized as one of the 
basic principles  of  gene regulation[9, 10]. We and others have demonstrated that although 
chromatin  interaction  domains  are  conserved  at  one  Mb  size  among  different  cell  types, 
hundreds of Kb-size subdomains are usually diverse and dynamically changed[11-15]. These 
highly dynamically changed domains may underlie mechanistic and functional consequences, 
and  computationally  identified  looping  genes  can  be  used  as  prognostic  signatures  for 
measuring the outcome of the endocrine treatment[16]. Current computational methods still 
mainly focus on constructing or inferring global topological domains and detecting long-range 
interacting loops[16-21]. However, many challenges, especially in the aspect of computational 
modeling and analysis remain. For example, how to computationally integrate Hi-C data with 
many  other  epigenetic  modification  data  to  classify  cell  type-specific  intra-chromosomal 
community  interactions  from  multiple  cell  types  or  at  different  time  points  upon  an 
environmental or cellular stimulus? What are the resulted biological functions associated with 
differential intra-chromosomal community interactions between cell types?   
 
Despite  many  computational  tools  for  processing  Hi-C  data,  no  computational  tools  are 
currently available for systematically integrating Hi-C data with epigenetic modification data 
from multiple cell types or at different time points upon an environmental or cellular stimulus. 
In this study, we will develop an integrative data analysis approach for investigating intra-
chromosomal  community  interactions  (Figure  1),  which  includes  1)  constructing  genomic   4
feature matrices by mapping multi-omics datasets to adjacency matrices obtained from Hi-C 
data, 2) identifying valid community structures based on both the network clustering and the 
optimized parameter determination methods, 3) evaluating enrichment of genomic features in 
intra-chromosomal  community  interaction  with  a  permutation  test,  and  4)  screening  for 
significant  gene  loci  and  corresponding  sub-network  by  comparing  genomic  and  network 
topological features through differential intra-chromosomal community interaction analysis. 
Finally, we apply the computational approach on an ERα+ breast cancer model system, i.e., 
untreated tamoxifen-sensitive (TS) cell line MCF7, one hour estrogen (E2)-treated MCF7 cell 
line, as well as tamoxifen-resistant (TR) cell line MCF7TR. Our goal is to identify distinct 
types of TS- and TR-intra-chromosomal community interactions and to gain insights into the 
functional  role  of  these  TR-intra-chromosomal  community  interactions  in  contributing  to 
breast cancer tamoxifen resistance. Moreover, we anticipate that applying and validating our 
novel  bioinformatics  tools  in  an  ERα+  breast  cancer  model  system  will  make  it  broader 
application in many other biological systems for identifying a biological-specific differential 
community interaction.  
 
Material and Methods 
Multi-omics Data 
Hi-C  data  in  untreated  and  one  hour  of  E2-treated  MCF7  cells,  and  MCF7TR  cells  were 
obtained from our previous publication[11]. Preprocessing of raw Hi-C data was performed by 
HiCUP[22] (e.g., map reads to HG19 human reference genome et al). Aligned BAM files of 
Hi-C  experiments  were  further  processed  by  either  in-house  python  scripts  (e.g., 
normalization[23] and transformation of Hi-C interaction data to Z-scores[24]) or HOMER[25] 
for selecting significant chromosomal interactions. Then, intra-chromosomal interactions were 
studied in 500 kb, 100kb, and 50kb resolution, respectively. Microarray expression profiles of   5
genes in zero and one hour of E2 treated MCF7 cells, as well as RNA-Seq data in tamoxifen-
resistant  MCF7TR  cells  were  retrieved  from  earlier  paper[26,  27],  which  were  quantile 
normalized[28] before transforming them to Z-scores. ChIP-Seq experiments of markers for 
enhancer  (H3K27ac  and  H3K4me1),  promoter  (H3K4me3),  repressor  (H3K27me3  and 
H3K9me3), nucleosome density (DNase-seq and ATAC-Seq), and insulator (CTCF) in the 
same MCF7 cells were downloaded from published work[11, 27] and ENCODE project[29]. 
Preprocessing of these ChIP-Seq and ATAC-Seq experiments was done similarly as published 
works[11, 24, 30]: for example, raw reads were mapped to HG19 human reference genome by 
bowtie2[31],  narrow  (e.g.,  in  CTCF,  H3K27ac  and  H3K4me3)  and  broad  (e.g.,  DNase, 
H3K27me3, H3K4me1 and H3K9me3) peaks were called by MACS2[32], and problematic 
regions  for  high  throughput  sequencing  in  human  genome  were  removed  according  to 
ENCODE Black List[33]. Then,  BEDTools[34] were used to estimate the read counts in a 
200 bp window size for called peaks which were quantile normalized across all markers[35] 
before log-transforming and converting them to Z-scores. 
 
Building adjacency matrices with predefined resolution for Hi-C interactions 
If a chromosome length is L (e.g., 22 human chromosomes and an X chromosome), then the 
number of window bins  𝑛 𝑏 𝑖𝑛  of this chromosome based on a predefined window bin size or 
interaction resolution (e.g., bin=500kb) is 
𝑛 𝑏 𝑖𝑛 = ⌈ 𝐿
𝑏 𝑖𝑛⌉,  
where  ⌈𝑥⌉  is a ceiling function, returning an integer greater than or equal to  𝑥 . Then, the i-th 
interval (or window bin) of a chromosome is  𝑤 𝑖 : [(𝑖 − 1) ∙ 𝑏 𝑖𝑛, min(𝑖 ∙ 𝑏 𝑖𝑛, 𝐿)), 𝑖 ∈ {1,2, ⋯ , 𝑛 𝑏 𝑖𝑛 }.    6
For  a  set  of  𝐻 (ℎ 1 , ⋯ , ℎ 𝐾 ) 𝑇  with  K  Hi-C  intra-chromosomal  interactions  (e.g.,  significant 
interactions identified by HOMER), each pair of interactions  ℎ 𝑘 (𝑥 𝑘 , 𝑦 𝑘 )  will be mapped to a 
pair of window bins  ℎ 𝑘 (𝑥 𝑘 , 𝑦 𝑘 ) → 𝑤 𝑖𝑗 (𝑤 𝑖 , 𝑤 𝑗 ),  
where  𝑥 𝑘  and  𝑦 𝑘  are the positions of a pair of Hi-C interactions in intervals (or window bins)  𝑤 𝑖  and  𝑤 𝑗 , respectively, and  𝑤 𝑖𝑗  is an interaction between window  𝑤 𝑖  and  𝑤 𝑗 . Thus, some  𝑤 𝑖𝑗  
might contain several pairs of Hi-C interactions while others might have none. Subsequently, 
an adjacency matrix A for intra-chromosomal interactions can be obtained   𝐴 𝑖𝑗 = Ι {𝑤 𝑖𝑗  ∪ 𝑤 𝑗𝑖 ≠∅} , 𝑖, 𝑗 ∈ {1,2, ⋯ , 𝑛 𝑏 𝑖𝑛 },  
where  𝐴 𝑖𝑗 = 1  if the union ( ∪ ) of  𝑤 𝑖𝑗  and  𝑤 𝑗𝑖  is not empty ( ∅ ), otherwise  𝐴 𝑖𝑗 = 0 . In this way, 
intra-chromosomal  Hi-C  interactions  are  transformed  to  a  binary  adjacency  matrix  A  with  𝐴 𝑖𝑗 = 𝐴 𝑗𝑖 . Here, if different window bin sizes (e.g., bin size = 500kb, 100kb, or 50kb) are used 
to  map  Hi-C  interactions,  then  intra-chromosomal  interaction  networks  with  different 
resolutions will be obtained. 
 
Generating genomic feature matrices based on multi-omics data for Hi-C interactions 
First, Z-scores of intra-chromosomal interactions obtained from Hi-C data are mapped to the 
aforementioned  adjacency  matrix  (e.g.,  bin  size=500kb),  as  the  edge  weights  of  intra-
chromosomal interactions (e.g., 22 human chromosomes and an X chromosome) for untreated, 
one hours E2-treated MCF7 cells, and MCF7TR cells, respectively. Then, genomic features 
(e.g., nucleosome density, histone modifications, and gene expressions) are mapped to the same 
adjacency  matrix  for  generating  genomic  feature  matrices,  respectively.  For  example,  𝐻 (𝑔) (ℎ 1(𝑔) , ⋯ , ℎ 𝑀(𝑔) ) 𝑇  is  genomic  feature  g  with  M  data  points,  where  each  ℎ 𝑚(𝑔) (𝑥 𝑚 , 𝑧 𝑚 )  
contains position information  𝑥 𝑚  and Z-score  𝑧 𝑚  of feature g in the loci. For a set of genomic 
features  𝑊 𝑖  in a window bin  𝑤 𝑖    7
𝑊 𝑖 = {𝑧 𝑚 |𝑥 𝑚 ∩ 𝑤 𝑖 ≠ ∅, 𝑚 ∈ {1,2, ⋯ , 𝑀 }},  
if loci  𝑥 𝑚  intersects with bin  𝑤 𝑖  (e.g.,  𝑥 𝑚 ∩ 𝑤 𝑖 ≠ ∅ ), then feature  𝑧 𝑚  will be included in this 
bin. Otherwise, the feature  𝑧 𝑚  will not be considered because loci  𝑥 𝑚  is not in bin  𝑤 𝑖   (e.g.,  𝑥 𝑚 ∩ 𝑤 𝑖 = ∅ ). Since there may be multiple (or none)  𝑧 𝑚  in the same bin  𝑤 𝑖 , a genomic feature 
Z-score for the bin  𝑤 𝑖  is 
𝑍 𝑖 = {𝑊 𝑖 ̅̅̅, 𝑊𝑖 ≠ ∅ 0, 𝑊𝑖 = ∅ .  
where  𝑊 𝑖 ̅̅̅  is the mean of these  𝑧 𝑚 . Subsequently, a feature matrix  𝐹 (𝑔)  can be built for each 
genomic feature g based on the adjacency matrix A of Hi-C intra-chromosomal interactions. 
For instance, if there is an interaction between the window bins  𝑤 𝑖  and  𝑤 𝑗  (e.g.,  𝐴 𝑖𝑗 = 1 ), then 
an edge weight  𝐹 𝑖𝑗(𝑔)  of the interaction  𝑤 𝑖𝑗  at feature matrix F is assigned by equation 
𝐹 𝑖𝑗(𝑔) = max(0, 𝑍 𝑖 + 𝑍 𝑗 )  
It indicates the edge weight of feature matrix  𝐹 (𝑔)  is determined by the Z-scores of a pair of 
interacting  bins  (e.g.,  𝑤 𝑖  and  𝑤 𝑗 ).  Such  edge  weight  assignment  (e.g.,  a  weighted  edge 
approach) for a network graph is adopted from a previous work [36], which is used here to 
build genomic feature matrices  𝐹 (𝑔)  for adjacency matrix A.  In order words, if there is an 
interaction in A, then the available genomic features will be considered correspondingly for 
building  genomic  features  matrices  𝐹 (𝑔) .  After  converting  multi-omics  datasets  𝐻 (𝑔)  to 
genomic feature matrices  𝐹 (𝑔) , based on the intra-chromosomal interactions A derived from 
Hi-C experiments (e.g., in untreated, one hours E2-treated MCF7 cells, and MCF7TR cells, 
respectively), information from diverse sources is comparable and integrable through the intra-
chromosomal interactions.  
 
Identifying significant intra-chromosomal interactions from Hi-C interactions   8
For  removing  unreliable  or  biased  intra-chromosomal  interactions  from  Hi-C  data,  two 
approaches are used to screen significant chromosomal interactions. One is the filtering of the 
weakest 0, 5, 10, 20, or 60 percentages of intra-chromosomal interactions from observed Hi-C 
interactions. And the other is applying HOMER[25] directly on the aligned BAM files of Hi-
C data, to obtain significant interacting fragments[37] such as using “makeTagDirectory” and 
“analyzeHiC”  with  significant  P-value  <  0.1.  The  screened  significant  intra-chromosomal 
interactions are converted to adjacency matrices A (or intra-chromosomal interaction matrices) 
with  a  predefined  window  bin  size.  Both  filtering  approaches  were  applied  for  screening 
observed intra-chromosomal interactions (e.g., window bin size=500kb) in untreated and one 
hours E2-treated MCF7 cells, respectively. The aim is to evaluate the robustness of predicted 
intra-chromosomal community interactions (e.g., the network modularity scores, the number 
of network edges, and the number of valid communities or network clusters) based on different 
filtering criteria.  
 
Since  identified  significant  chromosomal  interactions  (e.g.,  adjacency  matrices  of  intra-
chromosomal interactions) at different interaction resolutions (e.g., window bin size = 500kb, 
100kb,  or  50kb)  are  frequently  affected  by  the  program  parameters  (e.g.,  HOMER),  it  is 
essential  to  find  a  solution  for  selecting  optimal  parameters  accordingly.  Though  a  higher 
interaction  resolution  (or  a  smaller  window  bin  size)  may  yield  more  accurate  network 
structure, it cannot be increased infinitely. For instance, the significant interactions identified 
by HOMER are dependent on two major parameters: window resolution – bin size, and its 
super  resolution  –  range  of  signal  averaging.  To  find  optimal  parameters  for  identifying 
significant intra-chromosomal interactions by HOMER, a new approach is designed to evaluate 
the associations between the number of identified interactions and the distance of interactions, 
when different window bin sizes are used in analysis. This new statistical method can minimize   9
the  impact  of  window  bin  sizes  (e.g.,  bin=500kb,  100b,  or  50kb)  on  building  adjacency 
matrices,  when  applying  HOMER  on  Hi-C  interactions.  In  the  end,  the  significant  intra-
chromosomal  interactions  identified  by  HOMER  with  optimized  parameters  will  be  used 
further in differential network analysis. 
 
Predicting intra-chromosomal community interactions 
Applying modularity function to find communities in intra-chromosomal interaction network 
An  intra-chromosomal  interaction  matrix  from  Hi-C  experiments  is  represented  by  an 
adjacency  matrix  in  graph  theory,  where  intra-chromosomal  interactions  are  treated  as  an 
undirected network. It is interesting to see whether there are intra-chromosomal community 
interactions in the network or not. For example, a smart local moving (SLM) algorithm [38] 
can be applied on the Hi-C intra-chromosomal interactions (e.g., an adjacency matrix), to find 
community  structure  of  intra-chromosomal  interactions.  The  SLM  method  is  based  on  a 
modularity  function  [39]  for  detecting  communities  in  a  large  network.  The  modularity 
function was proposed by Newman and Girvan[39] for evaluating the whether the  division of 
a network into communities is good or bad. Let  𝐴 𝑖𝑗  be a node in an adjacency matrix A of the 
network, equivalent to a window bin in a Hi-C intra-chromosomal interaction matrix, if  𝐴 𝑖𝑗 =
1  then two nodes (or window bins  𝑤 𝑖  and  𝑤 𝑗 ) are connected, otherwise  𝐴 𝑖𝑗 = 0 . Currently, 
only unweighted networks are considered in building the network communities. The fraction 
of edges that assigned to communities (e.g., a node i belong to community  𝑐 𝑖 ) is 
𝐹𝐸 = ∑ 𝐴𝑖𝑗
(𝑐𝑖, 𝑐𝑗) 𝑖𝑗 ∑ 𝐴𝑖𝑗 𝑖𝑗 = 1
2𝑚 ∑ 𝐴 𝑖𝑗
(𝑐 𝑖 , 𝑐 𝑗 ) 𝑖𝑗  
where 

(𝑐 𝑖 , 𝑐 𝑗 ) = 1  if  𝑐 𝑖 = 𝑐 𝑗  and 0 otherwise, and m is the total number of edges in the graph. 
The expected fraction of edges assigned to communities in a randomized network is 
𝐸  = 1
2𝑚 ∑ 𝑘 𝑖 𝑘 𝑗 2𝑚

(𝑐 𝑖 , 𝑐 𝑗 ) 𝑖𝑗    10
where  𝑘 𝑖 = ∑ 𝐴 𝑖𝑗 𝑗  is the degree of node i in a network. Subsequently, a modularity function  𝑄 = 𝐹𝐸 − 𝐸  [38] to evaluate the division of the network can be obtained. If the FE is close to 
E, except in the case of randomized network, then the modularity score Q will be close to zero. 
In  order  words,  a  higher  value  from  modularity  function  Q  indicates  a  better  division  of 
community structure in a network. Usually, a modularity score Q above 0.3 may suggest of 
significant  community  structure  in  a  network[39].  To  search  for  optimized  community 
structure of a big network based on the modularity function, SLM algorithm[38] (a procedure 
of the local moving heuristic) is used: to repeatedly move one node from one community to 
another that increasing the modularity score Q. After many iterations, the final assignment of 
nodes to communities will be reached if no nodes can be moved further to increase the  Q. 
Please refer to previous publications [38-40] for more detailed descriptions of SLM algorithm 
in large network structure study.  
 
Determining the minimum size of valid communities 
Here, a new method is developed to estimate the minimum size of valid communities following 
a specific window bin size or interaction resolution. For example, a minimum size of valid 
communities  (or  network  clusters)  is  determined  by  the  predicted  communities  of  intra-
chromosomal interactions at a predefined resolution (e.g., window bin size=500kb). First, the 
communities with only one node are removed, and the number of identified communities in a 
chromosome r is defined as  𝐶 𝑟 . Then, the number of edges in a community (or network cluster) 
j for chromosome r is defined as  𝑒 𝑟𝑗 , where  𝑟 ∈ {1,2, ⋯ ,23} , and  𝑗 ∈ {1,2, ⋯ , 𝐶 𝑟 } . Since both 
the number of edges and the size of communities are dependent on the nodes number of a 
network, which are further determined by the resolution of network (e.g., the choice of window 
bin size), the number of edges  𝑒 𝑟𝑗  may not reflect the real-world scale of a network  when 
different network resolution is used (e.g., window bin size = 500kb, 100kb, or 50kb). Thus,   11
𝑒 𝑟𝑗   needs  to  be  normalized  against  a  specific  window  bin  size  (or  network  resolution)  by 
equation 
𝑝 𝑟𝑗 = 𝑒 𝑟𝑗 𝑒 𝑗 ̅ ,  
where  𝑒 𝑗 ̅ = ∑ 𝑒𝑟𝑗 𝑟,𝑗 ∑ 𝐶𝑟 𝑟  is the average number of edges of all chromosomes at the same network 
resolution. Subsequently, a predefined percentage value p (e.g.,  𝑝 = 0.02 ) can be set to find 
communities with at least p percentage (e.g., 2% in this work) of the average number of edges 
based on equation  𝑝 ∗ = 𝑎 𝑟 𝑔𝑚𝑖𝑛 𝑝𝑟𝑗 |𝑝 𝑟𝑗 − 𝑝|,  
where the closest value  𝑝 ∗  to  𝑝  is identified and the corresponding number of expected edges  𝑒 ∗  is calculated by  𝑒 ∗ = 𝑝 ∗ 𝑒 𝑗 ̅,  
which is the expected minimum size (or number of edges) of valid communities in a network 
at the chosen interaction resolution. Based on the predefined threshold value p, the expected 
minimum size of valid communities in a community interaction network can be determined by  𝑒 ∗  (e.g.,  the  number  of  edges  for  valid  communities  > 𝑒 ∗ )  and  estimated  automatically  for 
different Hi-C interaction networks at different resolutions (or window bin size). Please note 
that the edges of intercommunity interactions are not considered in the calculation. Thus, the 
number of identified total communities of chromosome r will be reduced from  𝐶 𝑟  to  𝑉 𝑟  after 
filtering communities that do not reach the minimum size (or the number of edges) of valid 
ones. 
 
Inferring communities and genomic feature matrices based on intra-chromosomal interactions 
First, an intra-chromosomal interaction matrix H (e.g., a normalized Hi-C contact count or Z-
score matrix) is converted to an adjacency matrix A based on a predefined resolution (e.g.,   12
window bin size = 500Kb). Then, multi-omics datasets  𝐻 (1) , ⋯ , 𝐻 (𝐺 )  from various genomic 
features are used to build genomic feature matrices  𝐹 (1) , ⋯ , 𝐹 (𝐺 )  follow the adjacency matrix 
A,  where  G  is  the  number  of  genomic  features  (e.g.,  gene  expressions,  nucleosome 
density/DNase-seq, histone modifications for enhancer “H3K4me1 and H3K27ac”, promoter 
“H3K4me3”, repressor “H3K9me3 and H3K27me3”, or insulator “CTCF”). The weights of 
genomic features are assigned to the edges of feature matrix F based on intra-chromosomal 
interactions in adjacency matrix A. Subsequently, SLM algorithm is applied on the unweighted 
adjacency  matrix  A  (e.g.,  with  a  window  bin  size=500kb)  to  predict  intra-chromosomal 
community interactions. After determining the minimum size (or number of edges) of valid 
communities, the valid communities V (or network clusters) are used in further data analysis 
such as the genomic feature enrichment and the differential network analysis.  
 
Evaluating genomic feature enrichment in intra-chromosomal interaction communities 
Based on the community network, an enrichment of various genomic features in valid intra-
chromosomal  communities  can  be  evaluated.  For  instance,  to  evaluate  the  difference  of 
genomic feature l between the inside community ( 𝐹 𝑣(𝑙) ) and the outside community ( 𝐹 𝑣′(𝑙) ), a 
two-sample T-test is applied on each valid community (or a network cluster), such as the T-
value of feature l is   𝑇 𝑖(𝑙) = 𝑇 𝑡 𝑒𝑠𝑡 (𝐹 𝑣(𝑙) , 𝐹 𝑣′(𝑙) ),  
where  𝐹 𝑣(𝑙)  represents the observed genomic feature l within the valid community v, while  𝐹 𝑣′(𝑙)  
is the expected genomic feature from randomly selected intra-chromosomal interactions (e.g., 
drawn from outside community v but with the same size as  𝐹 𝑣(𝑙) ). Meanwhile, the P-value (e.g.,  𝑃 𝑖(𝑙) ) of the corresponding  𝑇 𝑖(𝑙)  indicates whether there is a significant enrichment of genomic 
feature l at the community in sampling i or not (e.g., a significance level  𝛼 = 0.01 ). To obtain   13
an expected P-value ( 𝑃 (𝑙) ) for the enrichment of genomic feature l in a community, the random 
permutation test is repeated N times (e.g., resampling 1000 times) then calculate 
𝑃 (𝑙) = ∑ Ι {𝑃𝑖(𝑙)<𝛼} 𝑁 𝑖=1 𝑁 ,  
a small  𝑃 (𝑙)  indicates a significant feature enrichment in the community when it is compared 
against randomly drew samples outside the community. Subsequently, an expected T-value for 
the enrichment of genomic feature l in a community is  
𝑇 (𝑙) = ∑ 𝑇 𝑖(𝑙) 𝑁 𝑖=1 𝑁  
The permutation test is performed at each genomic feature l to access its overall enrichment in 
a  community  (or  network  cluster),  respectively.  Finally,  the  absolute  log10  transformed 
expected P-value with sign of expected T-value are visualized in a color coded heatmap, which 
illustrates both the significance and the direction (positive or negative) of genomic feature 
enrichment in intra-chromosomal communities. 
 
Differential network analysis - screening significant intra-chromosomal interaction changes 
Information  of  both  the  intra-chromosomal  interaction  network  from  Hi-C  data  and  the 
genomic feature enrichment in the communities are considered in differential network analysis. 
First,  the  network  structure  information  (e.g.,  node  centrality)  of  intra-chromosomal 
interactions are obtained in different conditions from the adjacency matrices A, respectively, 
which follows a specific resolution (e.g., window bin size =500kb) of Hi-C interactions. The 
degree  centrality of  a node  w  in  a network can tell  the importance of the node within the 
network based on its position[41], which is 
𝐶 𝑒 𝑛 𝑡 𝑟𝑎𝑙𝑖𝑡𝑦 (𝑤 ) = 𝑑 𝑒 𝑔𝑟𝑒𝑒 (𝑤 )
𝑛 − 1 ,    14
where  𝑤  is a given node in an intra-chromosomal interaction network,  𝑛  is the total number of 
nodes,  and  𝑑 𝑒 𝑔𝑟𝑒𝑒 (𝑤 )  is  the  number  of  interactions  connected  to  node  𝑤 .  A  significant 
change of a node centrality between the conditions, may indicate its importance in perturbating 
network structure follows condition change. Secondly, nodes located in the same community 
will share the same genomic feature l enrichment in the community, such as the expected T-
value of genomic feature l enrichment at a community is shared by nodes within the same 
community.  Subsequently,  both  the  node  centrality  and  the  expected  T-values  of  genomic 
feature  l  enrichment  at  the  node  𝑤  are  converted  to  Z-score  and  built  into  a  vector  𝑍 𝑤 (𝑧 𝐶 𝑒 𝑛 𝑡 𝑟𝑎𝑙𝑖𝑡𝑦 , ⋯ , 𝑧 𝑇(𝑙) , ⋯ )  for further differential network analysis. 
 
To identify significant intra-chromosomal interaction changes, the aforementioned procedure 
is repeated in two intra-chromosomal interaction networks derived from Hi-C data of MCF7 
cells and tamoxifen-resistant MCF7TR cells, respectively, where Z-score vectors (e.g.,  𝑍 𝑤  and  𝑍 𝑤 ′ ) of nodes are obtained. The Z-score vectors contain both the node centrality in a network 
and the genomic feature enrichments at the node, which can be used directly for differential 
network analysis, such as to compute a Euclidean distance between  𝑍 𝑤  and  𝑍 𝑤 ′    𝐷 𝑖𝑠𝑡 (𝑤 , 𝑤 ′ ) = ‖𝑍 𝑤 − 𝑍 𝑤 ′ ‖,  
where  ‖∙‖  represents  the  ℓ 2  norm.  Consequently,  nodes  with  significant  intra-chromosomal 
interaction changes between the conditions (e.g., MCF7 cells vs. tamoxifen-resistant MCF7TR 
cells) can be quantified by Euclidean distance. If the Euclidean distance  𝐷 𝑖𝑠𝑡 (𝑤 , 𝑤 ′ )  is mapped 
to a Gaussian distribution and applying a P-value threshold (e.g., P < 0.05), then a set of nodes 
with  significant  Euclidean  distance  changes  (or  network  interaction  changes)  between  the 
conditions will be identified. Subsequently, genes from these significantly changed nodes (e.g., 
P < 0.05) can be obtained for further biological interpretation such as to find differentially   15
expressed genes (DEG), to perform gene enrichment test in GO terms or biological pathways, 
and to visualize them in a network etc.  
  
Results 
Enrichment of genomic features in intra-chromosomal communities  
Generating genomic feature matrices based on intra-chromosomal interactions 
The proposed integrative data analysis pipeline (Figures 1A) is utilized to investigate Hi-C 
intra-chromosomal  interactions  in  untreated  and  one  hours  E2-treated  MCF7  cells, 
respectively.  First,  Z-scores  transformed  intra-chromosomal  interaction  matrices  are  built 
based on the normalized Hi-C contact count matrices[23] with previously published procedures 
[24].  Then,  significant  intra-chromosomal  interactions  were  identified  by  HOMER  [25]  in 
500kb resolution. The two results were converted to an adjacency matrices A with window bin 
size 500kb, respectively. In order to minimize the impact of random (or bias) interactions on 
the analysis[42], the first result was further filtered by removing the weakest 5, 10, 20, 40, and 
60  percentages  of  intra-chromosomal  interactions  before  building  adjacency  matrices, 
respectively. Based on these intra-chromosomal interaction matrices (or adjacency matrices), 
the contribution of genomic features in interactions (e.g., a pair of window bins or nodes with 
an edge in a graph) are estimated by a weighted edge approach [36]. Subsequently, genomic 
feature  matrices  F  based  on  the  adjacency  matrices  are  built  for  each  feature  in  every 
chromosome, respectively. For example, in Supplementary Figures (S1, 2 and 3), heatmaps of 
both  intra-chromosomal  interaction  matrices  (or  weighted  adjacency  matrix  A)  and  the 
corresponding genomic feature matrices F (e.g., gene expression, nucleosome density/DNase, 
insulator/CTCF,  and  histone  modifications  for  enhancer,  promoter  and  repressor)  before 
filtering weakest  interactions,  are shown for chromosomes 3, 17, and 20 in  zero hours E2 
treated MCF7 cells, respectively. Here, it is not easy to interpret the biological meaning of   16
these heatmaps due to the large size of networks, though all information come from the same 
cell line or condition. Thus, a network clustering or community prediction approach is needed. 
 
Identifying intra-chromosomal community interactions 
Here, a modularity function and a smart local moving (SLM) algorithm[38] were applied to 
predict intra-chromosomal community interactions, based on unweighted adjacency matrices 
A derived from either the significant intra-chromosomal interactions provided by HOMER or 
the  Hi-C  intra-chromosomal  interactions  after  filtering  the  weakest  5,  10,  20,  40,  and  60 
percentages of interactions, respectively. Figure 2 and Supplementary Figure S4 show that the 
median number of edges (or interactions) in A is between 8000 and 25000 for untreated and 
one hours E2-treated MCF7 cells, respectively. When filtering the weakest interactions, the 
modality  scores  of  identified  intra-chromosomal  community  interactions  increase  in  both 
conditions  as the percentage of removed weakest interactions rises. The score reaches 0.3 (e.g., 
a median of modularity score <0.3) when 60% of weakest intra-chromosomal interactions are 
filtered, which indicates there is a community structure in the network[39]. Nevertheless, the 
score is still smaller than that based on the significant interactions (e.g., a median of modularity 
scores>0.3) identified by HOMER. 
 
Figure 3 and Supplementary Figure S5 are violin plots of the number of valid communities (or 
network clusters) and the number of edges (or interactions) in valid communities for Hi-C data 
of untreated and one hours E2-treated MCF7 cells, respectively. The results indicate that there 
are  either  two  or  three  valid  communities  in  intra-chromosomal  interaction  networks, 
regardless of the filtering criteria for weak interactions. And the median number of edges in 
valid communities is between 5000 and 15000 (e.g., Figure 3 and SFigure 5), where more than 
half of intra-chromosomal interactions are assigned to valid communities (or clusters) in each   17
chromosome.  Therefore,  the  filtering  of  weak  (or  identification  of  significant)  intra-
chromosomal  interactions  is  essential  for  revealing  true  community  structure  of  intra-
chromosomal interactions. For the simplicity of current analysis, a manually selected cutoff 
value (e.g., minimum size of edges  𝑒 ∗ > 20 ) is used to define the valid intra-chromosomal 
communities  (or  network  clusters)  at  500kb  resolution.  More  information  of  the  intra-
chromosomal community interactions (e.g., the modularity scores and the number of valid 
clusters in a chromosome) based on the significant interactions from HOMER and the filtering 
of weakest interactions are shown in Supplementary Tables S1 to S12, respectively.  
  
Enrichment of genomic features in intra-chromosomal community interactions 
After intra-chromosomal interactions are summarized by community interactions or network 
clusters,  it  becomes  easy  to  evaluate  the  enrichment  of  genomic  features  in  communities. 
Figures 4 and 5 are heatmaps of identified intra-chromosomal community interactions (e.g., 
chromosomes 3, 17, and 20) in untreated and one hours E2-treated MCF7 cells, respectively, 
after  filtering the weakest  20% of intra-chromosomal interactions. The first  column of the 
heatmaps are intra-chromosomal interaction matrices (or weighted adjacency matrices A from 
Hi-C data), where valid communities (or network clusters) are color coded (e.g., the same 
community visualized by the same color). The plots reveal that the interactions from the same 
community (or network cluster) are physically close to each, which is similar to topologically 
associating domain (TAD). The second last column of Figures 4 and 5 are the heatmaps of 
absolute  log10  transformed  P-values  with  the  sign  of  the  T-values  from  genomic  features 
enrichment  test  (e.g.,  gene  expression,  nucleosome  density  –  DNase,  enhancer,  promoter, 
insulator, and repressor markers) in valid communities. These results indicate that there are 
two  major  communities  in  the  intra-chromosomal  interactions:  repressive  marks  (e.g., 
H3K27me3 and H3K9me3)  enriched  or active  marks  (e.g., enhancer  and promoter makers   18
H3K4me1, H3K4me3, and H3K27ac) enriched communities. This is analogous to the A/B 
compartments in the genome that are associated with either open or closed chromatin[43]. The 
same thing is observed in all results, regardless of the filtering criteria for intra-chromosomal 
interactions. Please refer to Supplementary Figures S6-11 for results of filtering the weakest 
40% and 60% interactions, as well as the significant interactions from HOMER, respectively. 
  
In Figures 4, 5 and Figures S6-11, the intra-chromosomal community interactions (e.g., the last 
column of the figures) are similar in both conditions, but the enrichment of genomic features 
in the communities (e.g., the second last column of the figures) are different from each other. 
For example, in chromosome 3 (Figure 4), the community C1 is enriched by repressor makers 
H3K27me3 and H3K9me3 (the second last column of Figure 4) at untreated MCF7 cell, but it 
is  disappeared  at  one  hours  E2  treated  condition  (the  second  last  column  of  Figure  5).  In 
chromosome  17,  a  community  C2  is  only  enriched  by  enhancer  markers  (H3K4me1  and 
H3K27ac)  at  untreated  MCF7  cell  (the  second  last  column  of  Figure  4),  but  a  similar 
community C3 is enriched by both the enhancer and the repressor markers (H3K27me3) at one 
hours  E2  treated  MCF7  cell  (Figure  5).  The  results  are  consistent  when  different  filtering 
conditions were applied on intra-chromosomal interactions (e.g., the weakest 20%, 40%, 60% 
of interactions were removed, or the significant interactions identified by HOMER).  
 
Optimizing the parameters for valid communities at different resolutions 
So far, the proposed analysis pipeline was only evaluated at a single resolution (e.g., window 
bin size = 500kb) with a fixed number of edges for valid communities (e.g.,  𝑒 ∗ > 20 ). The 
results indicate that the quality of inferred community structures is better by using significant 
chromosomal interactions identified by HOMER than that by an ad hoc filtering of the weakest 
interactions (Figures 2 and 3). For studying chromosomal interactions in different resolutions   19
(e.g., window bin size = 50kb, 100kb, and 500kb), new methods are needed to fine-tunning 
HOMER parameters and to estimate a minimum size of valid communities. Firstly, there are 
two essential  parameters (the window bin  resolution  and super resolution)  in  HOMER  for 
identifying  significant  chromosomal  interactions.  In  order  to  evaluate  the  impact  of  these 
parameters on HOMER, they are initially set to equal as recommended by the HOMER manual, 
then they are changed dynamically. For instance, in Figure 6A, the mean number of significant 
interactions found by HOMER (e.g., untreated and one hours E2 treated MCF7 cells) reaches 
a peak when the window bin size =50kb (e.g., the dashed lines in the left panel of Figure 6A), 
but it drops significantly after the bin size =10kb. At the same time, the median distance of 
identified significant Hi-C interactions drops sharply after the bin size reaches 100kb (e.g., the 
box plot in the left panel of Figure 6A), which indicates many long distance Hi-C interactions 
were excluded at here. These plots motivated us to find a way to estimate the optimal window 
bin size (or resolution) and the super resolution in HOMER: first, to fix the window resolution 
at  50kb  but  increase  the  super  resolution  from  50kb  to  500kb;  then,  the  higher  the  super 
resolution the larger the median distance for significant interactions (e.g., the box plot in right 
panel of Figure 6A), while the number of significant interactions is gradually decreasing such 
as a significant drop at the super resolution 200kb (e.g., the dash line plot in the right panel of 
Figure 6A). To have a tradeoff between the median distance of interactions and the number of 
interactions, the optimal window bin resolution and super resolution for HOMER at a specific 
interaction resolution can be estimated (e.g., window bin size = 50k and super resolution = 
200kb; Figure 6A). In other words, optimal HOMER parameters allow for the identification of 
as many significant interactions as possible, while remaining a reasonable median interaction 
distance. Based on Figure 6A and Supplementary Figure S12, if the window bin size is 50kb 
and 500kb, then the HOMER super resolution is 200kb and 500kb, respectively. Thus, this   20
strategy  can  be  used  by  HOMER  to  identify  significant  chromosomal  interactions,  while 
maintaining similar chromosomal interaction structures in different interaction resolutions. 
 
The  second  essential  parameter,  for  studying  intra-chromosomal  interactions  at  different 
resolutions  (or window bin  size), is the size of valid communities  -  the minimum size (or 
number of interactions) of valid communities. As described in the Method section,  𝑝 𝑟𝑗  is the 
normalized  number  of  edges  in  a  community  (or  network  cluster),  which  has  similar 
distributions and scales at different network resolutions (e.g., the light green bars in Figure 6B; 
bin size=500kb, 100kb, or 50kb). By setting a unified percentage cutoff value (e.g.,  𝑝 = 0.02  
; the black dashed line in bar plots of Figures 6B and 6C) for the minimum number of edges in 
a valid community (e.g., the light pink bars in Figure 6C and the bold row in Supplementary 
Table S13), the majority of small communities are excluded (e.g., have very small number of 
edges;  a  peak  of  light  pink  bars  at  the  left  most  X-axis  in  Figure  6C).  A  list  of  expected 
minimum  size  of  valid  communities  at  different  interaction  resolutions  is  shown  in 
Supplementary Table S13, which follow a real-world scale of window bin sizes. For example, 
if a window bin size is 500kb, then the minimum size of valid communities is ~20 edges (or 
interactions). Nevertheless, if the bin size is 100kb then the minimum size becomes ~100 edges, 
where both the resolution and the minimum size are increased five times. In other words, if the 
intra-chromosomal  interaction  resolution  reaches  50kb,  then  the  minimum  size  of  valid 
communities is at least 200 edges. In Supplementary Table S13, different percentage cutoff 
values p (e.g.,
 
𝑝 = 0.01 𝑜𝑟  0.05 ) were used to estimate the minimum size of valid communities 
at different resolutions, where  𝑝 = 0.02  gives the best results for identifying the minimum size 
of  valid  communities.  That  is  because  the  other  thresholds  remove  either  too  few 
(Supplementary Figure S13;  𝑝 = 0.01 ) or too many (Supplementary Figure S14;  𝑝 = 0.05 ) 
invalid communities.   21
 
Robustness of genomic feature enrichments at communities in different resolutions 
To identify significant interactions and to remove invalid communities in different interaction 
resolutions,  new  methods  have  been  developed  for  fine-tuning  HOMER  parameters  (e.g., 
window  resolution  and  super-resolution)  and  estimating  the  minimum  size  of  valid 
communities.  Thus,  the  proposed  analysis  pipeline  (Figure  1A)  was  repeated  in  different 
resolutions,  by  using  only  the  predicted  significant  intra-chromosomal  interactions  with 
HOMER  optimal  parameters  at  window  bin  size  =500kb,  100kb,  50kb,  respectively.  For 
instance, to identify significant interactions at window bin size = 50kb, the minimum size of 
valid communities is 480 (e.g., p=0.02; Supplementary STable 3) and the optimal HOMER 
window resolution and super resolution are 50kb and 200kb (Figure 6A), respectively. Figure 
7  and  Supplementary  Figures  S15  and  S16  are  results  of  chromosomes  17,  3,  and  20  in 
untreated  MCF7  cells,  respectively,  where  the  significant  intra-chromosomal  interactions 
obtained  from  HOMER  are  consistent  to  each  other  in  the  different  resolutions  such  as 
heatmaps of Figure 7A with bin size of 500kb, 100kb, and 50 kb. Subsequently, communities 
of these significant interactions are predicted by SLM algorithm, and a common percentage 
cutoff value (p in Supplementary Table S13) was applied to remove invalid communities. For 
instance if p=0.02 then the minimum size of valid communities is 23, 153, and 480 for 500kb, 
100kb, and 50kb resolutions, respectively.  
 
In Figure 7B, the black and the other color represent a pair of interactions (e.g., the two nodes 
with  one  edge)  are  clustered  in  the  different  communities  and  the  same  community, 
respectively. In Figure 7C, heatmaps of absolute log10 transformed P-values with sign of T-
values are results of genomic feature enrichment tests in communities, which were obtained by 
a random permutation t-test (e.g., 1000 times random sampling at a cutoff of P-value <0.01)   22
for evaluating the significance of genomic features inside a community versus that outside the 
community. In the heatmaps, the orange and blue colors represent the positive and negative 
enrichment of genomic features in the communities, respectively. These results are similar to 
each  other  in  different  interaction  resolutions  (e.g.,  bin  size  =  500kb,  100kb,  50kb).  For 
example, in chromosome 17, the intra-chromosomal interactions are mainly divided into two 
classes  of  communities  Figure  7:  one  with  low  gene  expression  but  high  enrichment  of 
repressive makers (H3K27me3 and H3K9me3), and the other with both high gene expression 
and high enrichment of activate markers (H3K4me1, H3K4me3, and H3K27ac. In the high 
resolutions (e.g., 100kb and 50kb), there are more than two communities but they are mainly 
divided into two clusters (e.g., the repressed and the activated genomic feature enrichment). 
This is consistent to known A/B compartments in the genome, which are associated with either 
euchromatin or heterochromatin[43]. Often, more detailed information of intra-chromosomal 
community interactions is revealed in the high resolution, such as an activated community C2 
in 500kb resolution (Figure 7C) was divided into two communities in 100kb resolution, and 
were  further  divided  into  three  communities  in 50kb  resolution.  Though  the  two  activated 
communities in 100kb resolution (e.g., C2 and C3)  are enriched differently by DNase and 
active markers (e.g., H3K4me3 and H3K27ac), they are similar in gene expression enrichment. 
The network graphs of valid community interactions are shown in Figure 7D, where the node 
size and the edge width indicate the community size and the number of community interactions, 
respectively.  In summary, both the heatmaps of genomic feature enrichments and the network 
graphs  of  valid  intra-chromosomal  community  interactions  are  consistent  in  different 
interaction resolutions (or window bin sizes).  
 
For the other chromosomes in untreated, one-hour E2-treated MCF7 cells, and MCF7TR cells, 
similar  results  are  obtained  (e.g.,  Supplementary  Figures  S15-22)  by  applying  the  newly   23
developed parameter tunning methods on Hi-C data in different resolutions. In this way, both 
the predicted intra-chromosomal community interactions and the identified genomic feature 
enrichments in valid communities are consistent in different window bin sizes (e.g., 500kb, 
100kb and 50kb; Figure 7). Thus, the results of high-resolution intra-chromosomal interactions 
(e.g., 50kb) can be used further for detailed investigation, such as to explore the association of 
dynamical changes of intra-chromosomal interactions with specific biological mechanism. 
 
Annotating genes involved in significant intra-chromosomal interaction changes 
Here,  a  differential  network  analysis  was  performed  to  identify  significant  Hi-C  intra-
chromosomal interaction changes between MCF7 cells and MCF7TR cells. All analyses were 
done  in  a  resolution  of  50kb  (window  bin  size  =  50kb),    ~2681  nodes  passed  filtering  of 
significant Euclidean distance changes (e.g., p-value<0.05) between the two conditions, and 
~1378 genes located in these selected nodes please refer to Supplementary Figure S23 and web 
Supplementary  Files  1  and  2  (e.g.,  all  web  supplementary  files  are  available  at  
https://dnaici.github.io/dnaici/).  After considering only differentially expressed genes (DEGs) 
(e.g., MCF7 vs. MCF7TR with absolute(relative ratio) > 0.66 [24]), the total number of selected 
genes were reduced to 515 that involve in both the differential intra-chromosomal interaction 
and the differential expression (or DIEGs; web Supplementary File 3). Subsequently, the gene 
ontology (GO) and biological pathway enrichment test were performed based on these 515 
DIEGs by using DAVID tool [44]. Enriched terms with the thresholds of Fold Enrichment >= 
2 and P-value < 0.05 were considered further (Figure 8, Supplementary Figure S24, and web 
Supplementary  File  4).  Interestingly,  we  observed  several  pathways  with  a  notable  fold 
enrichment, suggesting these pathways are potentially contributing to tamoxifen-resistance. 
For  example,  metabolism  pathways,  e.g.,  ascorbate  and  aldarate  metabolism,  porphyrin 
metabolism, as well as sphingolipid signaling pathway, are enriched, which indicated altered   24
cellular metabolism in resistant cells, aligning with findings from previous studies. [45, 46]. 
Besides, the PI3K-Akt signaling pathway and the EGFR signaling pathway, recognized for 
promoting breast cancer cell survival and could provide alternative growth signals that negate 
tamoxifen's  effectiveness,  exhibit  enrichment.  Furthermore,  the  mechanically  gated  and 
transmitter-gated ion channel activities were also found to be among the most significantly 
enriched GO terms, implicating that alteration of ion channels and transporters might play a 
role  in  contributing  to  endocrine-resistant  breast  cancer.  Not  surprisingly,  we  also  found 
endocrine resistance among the top pathways in KEGG, which further validated our method. 
Overall, the pathway enrichment test demonstrated looping-mediated mechanisms [11, 47] in 
regulating tamoxifen resistance breast cancer.  
 
Visualization of nodes involving in significant intra-chromosomal interaction changes 
For the selected 515 DIEGs, chromosome positions of the corresponding nodes (e.g., each node 
is  a  window  bin  with  size  50kb)  are  extracted  (web  Supplementary  File  5-8),  and  the 
interactions of nodes in each chromosome under MCF7 cells (Figure 9) and tamoxifen-resistant 
MCF7TR cells (web Supplementary File 9) were visualized by Cytoscape [48], respectively. 
In  Figure  9,  the  orange  and  blue  node  color  represent  active  (e.g.,  enriched  by  H3K4me1 
H3K4me3,  and  H3K27ac)  and  repressive  communities  (e.g.,  enriched  by  H3K27me3  and 
H3K9me3), respectively. If the nodes belong to the same community, then they will be marked 
by the same color at the outer ring of the network. Also, if there are edges between the nodes, 
then they interact to each other at the corresponding condition (e.g., significant Hi-C intra-
chromosomal interactions identified by HOMER). The differential network analysis between 
MCF7 and MCF7TR cells (e.g., Figure 9 and web Supplementary File 9) indicate that not only 
the genomic feature enrichments in nodes  are changed (e.g.,  Figure 9A;  from repressed to 
activated  community),  but  also  the  nodes  interactions  are  rewired  (or  change  of  intra-  25
chromosomal  interactions)  between  the  conditions.  Thus,  these  predicted  nodes  from 
differential network analysis are involved in significant intra-chromosomal interaction changes 
and are potentially important in breast cancer gene regulation.  
 
Especially,  in  chromosome  4,  the  genomic  feature  enrichment  in  selected  nodes  changed 
completely between MCF7 and MCF7TR cells: for example, two repressed communities (with 
green and yellow outer rings; left panel of Figure 9B) in MCF7 become active communities in 
MCF7TR (right panel of Figure 9B), and the interactions of nodes in repressed community 2 
at MCF7 (yellow outer ring in left panel of Figure 9B) are far less than that at MCF7TR (green 
outer ring in right panel of Figure 9B). It is also true for the other communities: for example, 
the only activated community in MCF7 was split to two repressed communities in MCF7TR, 
while the nodes interactions were also changed significantly (Figure 9C). Thus, genes in these 
selected  nodes  (e.g.,  chromosome  4;  Figure  9)  are  not  only  involved  in  significant  intra-
chromosomal interaction changes, but also might be responsible for the activation/repression 
of  the  region  under  different  conditions.  For  example,  CXXC4,  PITX2,  SLC16A10, 
TRAM1L1,  NDST3,  PRSS12,  FAM184A,which  normally  play  roles  in  maintaining  cell 
proliferation, estrogen response and metabolism were identified from the repressed community 
2 (left panel of Figure 9B). For instance, the downregulation of CXXC4 was reported as a 
contributing  factor  in  tamoxifen  resistance  through  the  activation  of  Wnt/β-catenin 
pathway[49].  Conversely,  25  genes  (e.g.,  CXCL2,  CXCL3,  CXCL8,  detailed  in  web 
Supplementary File 10) from the activated community (left panel of Figure 9C;) in the MCF7 
condition.  Specifically,  the  upregulation  of  CXCL2,  CXCL3,  and  CXCL8  is  positively 
associated  with  cancer  metastasis  and  chemo-resistance.  In  contrast,  downregulating  these 
chemokines can significantly suppress cancer cell motility[50]. Besides, several studies report 
that  CXCL3  and  CXCL8  play  multiple  roles  within  tumor  microenvironment,  including   26
recruiting  immunosuppressive  cells  to  the  tumor,  increasing  tumor  angiogenesis,  and 
promoting epithelial-to-mesenchymal transition (EMT), which are critical processes in tumor 
progression and metastasis[51, 52].  
 
Discussions and Conclusion 
In  this  work,  an  integrated  data  analysis  approach  (Figure  1)  is  developed  for  differential 
network analysis in intra-chromosomal interactions. First, the intra-chromosomal interaction 
networks  were  build  based  on  Hi-C  data  such  as  from  the  significant  intra-chromosomal 
interactions  identified  by  HOMER.  Usually,  the  higher  the  resolution  of  interactions,  the 
greater number of interactions is included in the network. Though there is neither a proper 
method to estimate optimal interaction resolution nor a solution to find optimal parameters (e.g., 
HOMER) for identifying significant interactions based on Hi-C dataset, Hi-C chromosomal 
interaction networks are frequently affected by these factors due to either the experimental 
limitation or the computational challenge, such as the sequencing depth of Hi-C experiments 
may affect the degree of Hi-C interaction resolution[53] that can be reached by a computer 
program. This motivated us to design a new approach (Figure 6) for fine-tuning parameters 
(e.g., window resolution and super-resolution) in HOMER at different interaction resolutions, 
while identifying significant  intra-chromosomal interactions.  It  guarantees  that  the inferred 
intra-chromosomal interaction networks are consistent in different resolutions (e.g., window 
bin size=500kb, 100kb, 50kb; Figures 6 and 7). For that reason, a proper high resolution intra-
chromosomal interaction network can be built (e.g., bin size =50kb) based on Hi-C data for 
further detailed network analysis.  
 
Secondly, intra-chromosomal interaction networks inferred from  Hi-C  data  often  contain  a 
huge number of interactions (e.g., tens thousand interactions with a window bin size 500kb;   27
Figures 2 and 3). This hinders the data analysis and biological interpretation of chromosomal 
interactions,  a  dimensionality  reduction  approach  (e.g.,  network  clustering  or  network 
community prediction by the SLM algorithm with local moving heuristic) to summarize the 
high  dimensional  intra-chromosomal  interactions  into  the  low  dimensional  community 
interactions (Figures 4 and 5) is definitely needed, where each predicated community (or super 
node) includes several nodes from the original interaction network. Application of SLM for 
network clustering has been successful in other similar works (e.g., Louvain algorithm) for 
identifying community structure in large physical, biological, and social networks. Usually, the 
SLM algorithm has more freedom than the other methods (e.g.,  Louvain algorithm[40]) in 
searching for high quality community structure because of an optimization of modularity score 
function. Particularly, if a modularity score reaches 0.3 then there is a significant community 
structure in the interaction network[39]. Nevertheless, the SLM algorithm frequently generates 
many small communities (or network clusters) when clustering a large network, due to either 
experimental error or computational bias. Thus, a new statistical method (Figure 6) is designed 
here (e.g., for finding optimal size of valid communities in a specific interaction resolution) to 
handle such problem in the network clustering. 
 
After  obtaining  valid  intra-chromosomal  community  interactions  in  different  resolutions 
(Figure 6), intra-chromosomal community network can be built for each chromosome (Figures 
4,  5  and  7)  in  breast  cancer  cells  (e.g.,  MCF7  and  MCF7TR)  based  on  Hi-C  data.  The 
community  network  is  a  dimensionally  reduced  super  network  for  intra-chromosomal 
interactions, where each node represents a community of intra-chromosomal interactions (or 
window bins) with similar network structure property. In the community (or super) network, 
the size of a node and the weight of an edge indicate the number of interactions within and 
between the communities, respectively. Heatmaps of inferred intra-chromosomal community   28
interactions  (Figures  4  and  5)  reveal  that  the  structure  of  intra-chromosomal  interactions 
resembles topologically associating domain (TAD) in the genome [54]. Subsequently, genomic 
feature enrichment (e.g., active/repressive histone markers, gene expression and nucleosome 
density etc.) was performed in predicted communities, by using genomic feature matrices F 
built  from  intra-chromosomal  interactions.  The  genomic  feature  matrices  can  be  either 
visualized  in  color  coded  heatmaps  or  used  to  assess  genomic  feature  enrichment  in  the 
communities, such as Figures 4, 5 and 7 for untreated/one-hours-E2-treated MCF7 cells and 
tamoxifen-resistant MCF7TR cells, respectively. The majority of identified communities are 
highly  enriched  by  either  the  active  markers  (H3K4me1  and  H3K27ac)  or  the  repressive 
markers (H3K27me3 and H3K9me3). This clearly resembles the nature of A/B compartments 
in the genome that associated with open/closed chromatin[43]. Thus, the network clustering 
(or  the  community  interaction  network)  not  only  significantly  reduces  the  size  of  intra-
chromosomal  interactions,  but  also  helps  tremulously  for  understanding  the  property  of 
network structures in a large intra-chromosomal interaction network[39, 55] after considering 
diverse genomic features. 
 
Finally,  differential  network  analysis  was  performed  by  considering  both  the  intra-
chromosomal interaction network and the genomic feature enrichment in the communities, and 
network nodes involve in significant intra-chromosomal interaction changes (e.g., MCF7 vs. 
MCF7TR) are identified. Genes associated with these significantly changed nodes are not only 
enriched  in  GO  and  pathways  consistent  with  existing  biological  perspectives  (e.g.,  the 
presence of looping-mediated mechanisms in endocrine-resistant breast cancer; Figure 8), but 
they also potentially relate to conditions that cause the physical and functional changes of intra-
chromosomal  interactions  (Figure  9).  These  results  indicate  that  the  proposed  differential 
network  analysis  in  intra-chromosomal  community  interactions  (Figure  1)  may  help  the   29
understanding  of  both  mechanism  and  structure  of  long-distance  gene  regulation  based  on 
chromatin loops[5].   
 
In conclusion, an integrated data analysis approach is developed for identifying significant 
intra-chromosomal interaction changes, by combining Hi-C data with other genomic features 
such as active/repressive histone modifications, gene expression, and nucleosome densities etc. 
And  novel  parameter  tunning  methods  are  designed  for  finding  either  a  proper  intra-
chromosomal interaction resolution or an optimal community size based on a specific Hi-C 
experiment. The proposed computational framework is evaluated in breast cancer cells (e.g., 
untreated or one-hours-E2-treated MCF7 cells and tamoxifen-resistant MCF7TR cells), where 
a  comprehensive  picture  of  chromatin  interactions  from  Mb-sized  intra-chromosomal 
communities  to  Kb-sized  foci  of  key  interactions  is  provided.  The  genes  located  in  nodes 
associated with significant intra-chromosomal interaction changes (e.g., MCF7 vs. MCF7TR), 
are not only enriched in GO and pathways consistent with existing biological perspectives, but 
also providing methodological and mechanistic insights for future research such as to explore 
significant genomic dynamics in cancer and beyond by integrating multi-omics datasets.  
 
 
 
 
 
 
 
 
   30
AUTHOR CONTRIBUTIONS 
ZY Visualization, Methodology, Validation, Formal analysis, Software, Writing-Reviewing 
and Manuscript preparation. KF Validation, Writing-Reviewing and Editing. JSD Validation, 
Web Documentation, Writing-Reviewing and Editing. MB Validation, Funding acquisition, 
Writing-Reviewing and Editing. VJ Conceptualization, Validation, and Writing-Reviewing.  
JW  Conceptualization,  Methodology,  Software,  Data  curation,  Writing-Reviewing  and 
Manuscript preparation, Visualization, Investigation, Formal analysis, Supervision, Validation, 
Project administration, Funding acquisition. 
 
ACKNOWLEDGMENTS  
JW  and YZ  are  supported by the South-Eastern Norway Regional  Health Authority  (HSØ 
2018107), Akershus University Hospital internal grant 2023 and 2024, and the Norwegian 
Research Council NOTUR and NOSTORE project (nn4605k and NS4605k). MB and YZ are 
supported by Norwegian Research Council (326101).  
 
DECLARATION OF INTERESTS  
The authors declare no competing interests.  
 
 
 
 
 
 
 
   31
Figures  
Figure  1.  A  graphical  illustration  of  integrative  data  analysis  with  both  the  Hi-C 
interactions and the multi-omics data. 
 
(A) For each chromosome, intra-chromosomal interactions under a particular condition/cell 
line  obtained  from  a  Hi-C  sequencing  experiment,  which  are  mapped  to  a  binary 
adjacency  matrix  A  with  a  predefined  resolution  (e.g.,  window  bin  size=500kb). 
Meanwhile,  multi-omics  data  such  as  histone  modifications  (e.g.,  enhancer  mark,   32
repressive  mark),  nucleosome  density,  and  gene  expressions  are  mapped  to  genomic 
feature matrices F, based on the Hi-C interactions (e.g., adjacency matrix A) with the same 
window bin size.  
(B)  For E2 treated MCF7 cells and tamoxifen-resistant MCF7TR cells, intra-chromosomal 
communities  (e.g.,  C 1 ,  C 2 ,  ….)  are  identified  by  applying  SLM  algorithm  on  the 
unweighted adjacency matrix A. Then, a random permutation T-test is used to evaluate 
the enrichment of genomic features in  the communities, where the log10  transformed 
expected P-value with the sign of expected T-value are converted to color coded heatmap, 
for assessing the significant positive or negative enrichment of genomic features in intra-
chromosomal communities. 
(C)  Differential network analysis is applied on the Hi-C intra-chromosomal interactions by 
combining both topological feature and genomic features. Euclidean distance is used to 
assessing  the  difference  between  nodes  from  E2  treated  MCF7  cells  and  tamoxifen-
resistant MCF7TR cells, respectively. Significant differential interacting expressed nodes, 
genes, and sub-networks are identified by setting thresholds of fitted Gaussian distribution 
and relative ratio (P-value < 0.05, Absolute(relative ratio) > 0.66). 
     33
Figure  2.  Violin  plot  of  modularity  scores  and  number  of  edges  based  on  different 
selection  of  intra-chromosomal  interactions  in  Hi-C  data  at  500kb  resolution  for 
untreated MCF7 cells. 
 
The upper panel of figure shows violin plots of modularity scores[39] for intra-chromosomal 
community interactions in 23 chromosomes that were calculated based on the filtering of 5, 10, 
20,  40,  60  percentages  of  the  weakest  intra-chromosomal  interactions  genome-widely, 
respectively. “HOMER” represents the significant intra-chromosomal interactions predicted 
by HOMER[25] program. The lower panel of figure displays violine plots of the number of 
edges (or intra-chromosomal interactions) in 23 chromosomes after filtering 5, 10, 20, 40, 60 
percentages  of  the  weakest  intra-chromosomal  interaction  genome-widely  or  selecting 
significant interactions by HOMER. A white dot in the violin plot represents the median value 
of the distribution.      34
Figure 3. Violin plot of number of valid clusters and number of edges in valid clusters 
based on different selection of intra-chromosomal  interactions in Hi-C data at 500kb 
resolution for untreated MCF7 cells. 
 
The upper panel of figure shows violin plots of the number of valid network clusters (e.g., 
communities with the number of edges > 20) in intra-chromosomal interactions  among 23 
chromosomes, which were predicted by filtering of 5, 10, 20, 40, 60 percentages of the weakest 
intra-chromosomal  interaction  genome-widely,  respectively.  “HOMER”  represents  the 
selection of significant intra-chromosomal interactions by HOMER[25] program. Here, a black 
line within each violin  plot represents every observation inside the  distribution. The lower 
panel  of  figure  displays  violine  plots  of  the  number  of  edges  (or  intra-chromosomal 
interactions) in valid clusters of intra-chromosomal interactions among 23 chromosomes after 
filtering of 5, 10, 20, 40, 60 percentages of the weakest intra-chromosomal interaction genome-
widely  or  selecting  significant  interactions  by  HOMER.  A  white  dot  in  the  violin  plot 
represents the median value of the distribution.      35
Figure 4. Heatmaps of intra-chromosomal community interactions for chromosomes 3, 
17 and 20 after removing the lowest 20 percentages of intra-chromosomal interactions at 
500kb resolution in untreated MCF7 cells. 
In  each  panel,  the  left  most  heatmaps  are  intra-chromosomal  interaction  matrices  where 
interactions are colored by their community (or network cluster) labels. For example, black 
indicates two nodes of an edge (or interaction) do not belong to the same community, but the 
other colors represent the community label of an interaction. The right most figures are network 
plots of predicted intra-chromosomal community interactions (only show communities with 
more than 20 edges) in a chromosome, where the size of nodes and the width of edges indicates 
the number of edges (interactions) in the community and the number communications between 
the communities, respectively. The color of nodes indicates the community label. In a panel, 
the second right most figures are heatmaps of log10 transformed P-values for enrichment of a 
genomic feature in a community (or network cluster) against that from randomly generated 
(e.g.,  the  number  of  P-value  <  0.01  in  1000  times  random  permutation  T-tests)  intra-
chromosomal interactions, where the orange and blue colors represent positive and negative   36
enrichment  of  the  feature,  respectively.  This  sign  of  the  enrichment  is  inferred  from  the 
expected T-values of random permutation T-test. The rest of red and green coded heatmaps are 
genomic features in intra-chromosomal community interactions: “Zscore” is the Z-score of 
intra-chromosomal  interactions;  “geneExp”,  “DNas”  and  “CTCF”  are  average  of  gene 
expressions, nucleosome densities, and insulator in the interactions, respectively; “H3K4me1” 
and “H3K27ac” are enhancer markers; “H3K4me3” is a promoter maker; “H3K27me3” and 
“H3K9me3” are repressor markers. Here, “cn1” and “cn2” are community number 1 and 2 in 
a chromosome, respectively, which are the same numbers as the node labels in the network 
plots.  
     37
Figure 5. Heatmaps of intra-chromosomal community interactions for chromosomes 3, 
17 and 20 after removing the lowest 20 percentages of intra-chromosomal interactions at 
500kb resolution in one-hours E2-treated MCF7 cells. 
 
In  each  panel,  the  left  most  heatmaps  are  intra-chromosomal  interaction  matrices  where 
interactions are colored by their community (or network cluster) labels. For example, black 
indicates two nodes of an edge (or interaction) do not belong to the same community, but the 
other colors represent the community label of an interaction. The right most figures are network 
plots of predicted intra-chromosomal community interactions (only show communities with 
more than 20 edges) in a chromosome, where the size of nodes and the width of edges indicates 
the number of edges (interactions) in the community and the number communications between 
the communities, respectively. The color of nodes indicates the community label. In a panel, 
the second right most figures are heatmaps of log10 transformed P-values for enrichment of a 
genomic feature in a community (or network cluster) against that from randomly generated 
(e.g.,  the  number  of  P-value  <  0.01  in  1000  times  random  permutation  T-tests)  intra-
chromosomal interactions, where the orange and blue colors represent positive and negative   38
enrichment  of  the  feature,  respectively.  This  sign  of  the  enrichment  is  inferred  from  the 
expected T-values of random permutation T-test. The rest of red and green coded heatmaps are 
genomic features in intra-chromosomal community interactions: “Zscore” is the Z-score of 
intra-chromosomal  interactions;  “geneExp”,  “DNas”  and  “CTCF”  are  average  of  gene 
expressions, nucleosome densities, and insulator in the interactions, respectively; “H3K4me1” 
and “H3K27ac” are enhancer markers; “H3K4me3” is a promoter maker; “H3K27me3” and 
“H3K9me3” are repressor markers. Here, “cn1”, “cn2”, and “cn3” are community number 1, 
2, and 3 in a chromosome, respectively, which are the same numbers as the node labels in the 
network plots.      39
Figure 6. Statistical results used to determine the optimal parameters for the proposed 
Hi-C integration pipeline. 
 
(A)  Box  plots  are  the  length  of  identified  significant  interactions  from  HOMER  in  all 
chromosomes.  The  dashed  lines  are  the  mean  number  of  significant  interactions  in  all 
chromosomes  with  different  resolution  combinations  (e.g.,  window  bin  size=super 
resolution=50kb). (B) Distribution of the normalized proportion of edges in each cluster and a 
unified cutoff p=0.02 (black dash lines) under different resolutions. (C) Distributions of the 
number of edges in each cluster and their cutoffs corresponding to p=0.02 (black dash lines) 
under different datasets (e.g., t0, t1, and tamr represent untreated, one-hours E2-treated MCF7 
cells,  and  tamoxifen-resistant  MCF7TR  cells,  respectively)  and  resolutions  (e.g.,  500kb, 
100kb, 50kb).   40
Figure 7. Heatmaps of intra-chromosomal community interactions for chromosomes 17 
in  three  different  resolutions  based  on  significant  intra-chromosomal  interactions 
identified by HOMER in untreated MCF7 cells. 
 
 
(A) The heatmaps of significant intra-chromosomal interactions identified by HOMER.  
(B)  The heatmaps  of color-coded community labels  for intra-chromosomal interactions 
predicted  by  SLM  algorithm,  where  the  black  and  the  other  colors  are  a  pair  of 
interactions (e.g., the two nodes with an edge) classified in the different and the same 
community, respectively.  
(C)  The  heatmaps  of  genomic  feature  enrichments  in  valid  communities,  where  the 
minimum size of valid communities is 23, 153, and 480 for window bin size equals 
500kb, 100kb, and 50kb, respectively. 
A B C
window resolution = 500kb, super resolution = 500kb
window resolution = 100kb, super resolution = 200kb
window resolution = 50kb, super resolution = 200kb
D  41
(D)  Interactions of valid community are represented by the super networks. 
     42
Figure 8. Top 15 enriched functional annotations of biological pathways on 515 DIEGs  
 
Enrichment analysis of KEGG pathway and REACTOME pathway on 515 DIEGs obtained by 
DAVID online tool. The pink bar represents the Fold Enrichment (top axis), and the green 
scatters  denotes  the  -log10  transformed  P-value  (bottom  axis).  The  top  15  functional 
annotations are selected under the thresholds of Fold Enrichment >= 2 and p-value =< 0.05, 
and  ranked  by  Fold  Enrichment.  The  515  DIEGs  were  filtered  through  both  differential 
network  analysis  and  Absolute(relative  ratio)>0.66  from  intra-chromosomal  interaction 
network between MCF7 cells and tamoxifen-resistant MCF7TR cells. 
 
   
A 
B
KEGG Pathway
REACTOME Pathway  43
Figure 9. Differential sub-networks in chromosome 4 between MCF7 cells and 
MCF7TR cells 
 
Sub-networks  in  chromosome  4  obtained  by  differential  network  analysis.  Nodes  in  the 
network were selected based on DIEGs between MCF7 cells and MCF7TR cells, with each 
node  representing  a  predefined  window  bin  in  chromosome  (i.e.,  resolution  =  50kb).  The 
orange and blue colors in the nodes represent whether the histone markers of the community 
are  active  or  not.  The  colors  of  the  outer  rings  represent  previously  identified  intra-
MCF7
A
B
MCF7TR
C
Repressed community 1
Repressed community 1
Repressed community 2
Repressed community 2
Activated community
Activated communitychromosomal  communities.  Gray  lines  between  nodes  represent  significant  Hi-C  intra-
chromosomal  interactions.  (A)  Overview  of  differential  network  in  chromosome  4.  (B) 
Subregion  1  of  the  chromosome  shows  different  genomic  features  and  intra-chromosomal 
interactions. (C) Subregion 2 of the chromosome shows different genomic features and intra-
chromosomal interactions.