luria的个人博客分享 http://blog.sciencenet.cn/u/luria

博文

HiCPro接CscoreTools检测Compartment

已有 498 次阅读 2019-9-11 13:22 |个人分类:Hi-C|系统分类:科研笔记| HiCPro, Compartment, Cscore, HiC

2009Lieberman及其合作者发明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

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

 

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

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

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

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

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


通过检查发现按XXX_bias.txt文件去掉绝对值大于5和绝对值小于0.2bias后,剩下的bin可以和XXX_cscore.bedgraph对应。表明XXX_bias.txtXXX_hh.txtXXX_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 iwindow j都属于compartment A),它们之间的互作比来自不同类型compartment (例如window mwindow n分别属于compartment Acompartment B)的互作要强一些(ij之间的互作强于mn之间的互作)

 

3.1 F score & C score

作者先定义了一个Fij Score,它表示window iwindow j有多少可能属于同一个compartment,假设window iwindow j都在compartment A中,则其概率为PiPj,而window iwindow j都来自compartment B的概率为(1-Pi)*(1-Pj),那么

Fij = Pij + (1-Pi)(1-Pj) = [1+(2Pi-1)(2Pj-1)]/2

 

如何理解这个F Score

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

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

 

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

Ci = 2Pi -1

 

3.2 E score

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

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

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

综上,对F score进行改进,得E score

Eij = BiBjH(dij)(1+CiCj)

 

3.3 Poisson loglikelihood function

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

作者认为符合泊松分布。

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

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


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

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

 

回到compartment分析中,nij对应k,它是自变量;因变量符合泊松分布,其中参数λ= Eij

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


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


回到compartment分析中来:

即求Eij使得泊松分布的似然函数最大化,而Eij有三个未知参数BiH(dij)Ci。即求泊松分布的似然函数最大时三个参数BiH(dij)Ci的值。

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

 

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

CscoreTools也确实是这样做的:

1) 先设定初始值

2) 固定BiCi,求H(dij)的偏导;

3) 固定Ci和上步求得的H(dij)值,求Bi

4) 固定上述求得的Bi值和H(dij)值,求Ci

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

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

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

 

结果收敛后就可以得到Cscore,及两个window ij来自相同类型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

 




http://blog.sciencenet.cn/blog-2970729-1197554.html

上一篇:HiCPro统计结果解读
下一篇:宏基因组Binning软件MetaBat2的安装与使用

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-11-13 02:39

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部