自2009年Lieberman及其合作者发明Hi-C技术至今，Hi-C互作分析中Compartmant的主流检测方法一直是基于PCA降维后，用第一维度(通常用第一维度)来区分A/B compartment。直到2018年有了一个新的方法，基于Cscore，它也是本文主要讨论的对象。

按照套路，Cscore文章[1]中开篇提到了基于PCA检测Compartment时存在的问题：

1). 因为采用PCA降维后，第一维度的生物学意义难以解释(elusive)，因此不同的数据集之间无法进行比较

2). 当分析的互作矩阵很大时，基于PCA降维的办法非常消耗内存，使整个计算过程变得非常慢

基于Cscore开发的CscoreTools可以很好地解决以上两个问题。

1. 软件准备

有关HiCPro的安装、使用及统计结果解读，可参见我的另几篇博文：

1) Hi-C数据比对软件HiCPro的安装与使用

2) HiCPro分析流程详解

3) HiCPro统计结果解读

CscoreTools的安装：

按照github官网上的说明，CscoreTools下载后修改权限即可用。如果出现以下问题

./CscoreTool1.1: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.20' not found (required by ./CscoreTool1.1)

./CscoreTool1.1: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.21' not found (required by ./CscoreTool1.1)

按照github官网上的说明，编译一下也可以运行

g++ CscoreTool1.1.cpp twister.cpp -fopenmp -O3 -o CscoreTool1.1

注意：

CscoreTools自带generateEqualLengthBed.cpp程序，这个程序在后面会有用，这里需先编译一下：

g++ generateEqualLengthBed.cpp -o generateEqualLengthBed

2. CscoreTools的使用

按照CscoreTools使用说明，需先准备5个信息(文件)

CscoreTool1.1 <windows.bed> <input.summary> <outputPrefix> <session> <minDis> [chrName]

2.1 windows.bed

这个可以通过CscoreTools自带的generateEqualLengthBed生成

usage: generateEqualLengthBed <chrSizes.txt> <out.bed> <windowsize>

chrSizes.txt是染色体长度记录文件可以通过samtools faidx生成，如下：

samtools faidx /path/to/reference.fasta

#reference.fasta是提供的参考基因组，做互作分析需要染色体水平的基因组。这一步完成后会在/path/to/reference.fasta目录下生成一个/path/to/reference.fasta.fai

再运行generateEqualLengthBed

generateEqualLengthBed  /path/to/reference.fasta.fai  windows500k.bed  500000

完成后生成windows.bed 文件

2.2 input.summary

这个文件记录着Hi-C互作信息，是最主要的输入文件。格式借用HOMER软件规定的Hi-C Summary格式，具体如下
1). Read Name (can be blank)
2). chromosome for read 1
3). positions for read 1 (5' end of read, one-indexed)
4). strand of read 1 (+ or -)
5). chromosome for read 2
6). positions for read 2 (5' end of read, one-indexed)
7). strand of read 2 (+ or -)

我们可以利用HiCPro生成的allValidPairs文件来构建这个Hi-C Summary文件，这也是这篇博文最主要的目的！

HiCPro生成的allValidPairs文件在hicpro_output/hic_results/data/XX/XX.allValidPairs，具体如何运行HiCPro可以参见我的博文Hi-C数据比对软件HiCPro的安装与使用

allValidPairs文件格式如下：

1). read name

2). chr_reads1

3). pos_reads1

4). strand_reads1

5). chr_reads2

6). pos_reads2

7). strand_reads2

8). fragment_size

[9). allele_specific_tag] #如果提供Allele specific相关输入信息，则会生成这列

因此去掉最后几列即可：

cut -f 1-7 allValidPairs > input.summary

2.3 outputPrefix

输出文件名前缀

2.4 session

是一个整数，程序会将input.summary文件分成这么多份，同时来跑。相当于线程数

2.5 minDis

因为近程互作很可能会受到Loops、TAD的影响，这对于检测Compartment来说是干扰信息，CscoreTools会屏蔽掉这些近程互作。考虑到TAD通常在Mb水平，CscoreTools作者建议这个值设为1Mb，即1000000

运行完成后会生成四个文件：

1) XXX_bias.txt它是评估出的每个window的bias因子，即公式中的B_i（后面第3节会有详解），这里需要说明的是，如果运行时没有指定染色体，程序会按染色体顺序依次划分window，不会分染色体

2) XXX_hh.txt它是评估出的距离曲线，即后续公式中的H(d_ij)。如果设文件中第一列为n，则表示对应的距离为10^0.04*n

3) XXX_cscore.txt是基因组上每个window评估出的Cscore值，即后续公式中的C_i

4) XXX_cscore.bedgraph，这个bedgraph文件用于后续可视化分析，低比对率的windows(bias<0.2)或者太高拷贝数的windows(bias>5)会被过滤掉。因此我们后续来做的时候，有一些windows为空

通过检查发现按XXX_bias.txt文件去掉绝对值大于5和绝对值小于0.2的bias后，剩下的bin可以和XXX_cscore.bedgraph对应。表明XXX_bias.txt、XXX_hh.txt和XXX_cscore.txt是原始文件，XXX_cscore.bedgraph是过滤之后最终的结果。

我们可以通XXX_cscore.bedgraph文件中Cscore值，以0为分界，区分为A/B compartment。当然，还需要结合基因密度等其它信息综合评估正值为A compartment还是B compartment。

3. CscoreTools的算法

核心原理：

如果两个bin（在CscoreTools引文中作者称为window，为了尊重原文，我们下面都称为window）来自同一种类型的compartment (例如window i和window j都属于compartment A)，它们之间的互作比来自不同类型compartment (例如window m和window n分别属于compartment A和compartment B)的互作要强一些(i和j之间的互作强于m和n之间的互作)。

3.1 F score & C score

作者先定义了一个F_ij Score，它表示window i和window j有多少可能属于同一个compartment，假设window i和window j都在compartment A中，则其概率为P_iP_j，而window i和window j都来自compartment B的概率为(1-P_i)*(1-P_j)，那么

F_ij = P_ij + (1-P_i)(1-P_j) = [1+(2P_i-1)(2P_j-1)]/2

如何理解这个F Score？

假设window i来自compartment A的概率是1，那么如果window j来自compartment A的概率为1，则F_ij = 1；如果window j来自compartment A的概率为0，则F_ij=0。同时如果window i来自compartment B，结果相同

F score不仅做到了两个window属于相同compartment (无论是A&A还是B&B)时概率最高，两个window属于不同compartment概率越低。而且，将这一值限定在0-1之间！

从F score中可以看出i和j是对称的，即无论i属于compartment A还是compartment B，只要i和j一致，结果不变。因此，作者定义了C score：

C_i = 2P_i -1

3.2 E score

F score虽然很好，但是还是会受到一些干扰条件的影响，还不是最优。干扰因素如

=> 基因组距离越远window之间的互作越少。基于此，作者又提出了一个H(d_ij)，它是互作随距离变化的换算函数

=> 互作会受到PCR偏好性和基因组mappability的影响。针对这点，作者也提出了一个指标B_i，它表示互作bias

综上，对F score进行改进，得E score：

E_ij = B_iB_jH(d_ij)(1+C_iC_j)

3.3 Poisson loglikelihood function

在第3节开头我们提到，作者观测到来自相同类型compartment之间的互作要比来自不同类型compartment之间的互作强，那么连接数和来自相同类型compartment的概率是怎样的拟合关系呢？

作者认为符合泊松分布。

注意：这里有一个前提条件——两个window要相距很远。作者认为至少要在1M以上，这也是考虑到TAD会对compartment分析造成干扰。

泊松分布的概率函数如下：

泊松分布图像如下，有关泊松分布更多内容，可以参见泊松分布的wiki页面：

(https://en.wikipedia.org/wiki/Poisson_distribution)

回到compartment分析中，n_ij对应k，它是自变量；因变量符合泊松分布，其中参数λ= E_ij。

那么泊松分布的极大似然函数如下：

这里面参数λ是未知的，为了找到参数λ，使泊松分布的似然函数最大化，通常会对似然函数等号两边都取对数，即

回到compartment分析中来：

即求E_ij使得泊松分布的似然函数最大化，而E_ij有三个未知参数B_i，H(d_ij)和C_i。即求泊松分布的似然函数最大时三个参数B_i，H(d_ij)和C_i的值。

综上，函数中有三个未知数，现求函数值最大时这三个未知数的值。

如果之前有看过《EM算法(期望最大化算法)简介》一文，可能这时可能很快想到分别求偏导，再逐步迭代的办法。

CscoreTools也确实是这样做的：

1) 先设定初始值

2) 固定B_i和C_i，求H(d_ij)的偏导；

3) 固定C_i和上步求得的H(d_ij)值，求B_i

4) 固定上述求得的B_i值和H(d_ij)值，求C_i

5) 再将C_i值也更新，这样就完成一轮迭代

6) 不断重复上述过程，直到收敛为止

注意：文章中在这里提到迭代的顺序是H-B-C，关于H-B的顺序文章中没有详说

结果收敛后就可以得到Cscore，及两个window i和j来自相同类型Compartment的概率。

参考材料：

[1] Xiaobin Zheng and Yixian Zheng. CscoreTool: fast Hi-C compartment analysis at high resolution. Bioinformatics. 2018

[2] http://homer.ucsd.edu/homer/interactions/HiCtagDirectory.html

[3] http://nservant.github.io/HiC-Pro/MANUAL.html

[4] Nicolas Servant, Nelle Varoquaux, Bryan R. Lajoie, et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biology. 2015

[5] Lieberman-Aiden, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science. 2009