HiC数据分析实战之通过文章来了解流程

通过第一讲：三维基因组学习笔记，我们了解了3D基因组研究范围，然后根据我在生信技能树发布的生信工程师标准提炼出基础技能，也就是第二讲：生信基础技能 最后提炼出了数据分析流程，并且安装好了对应的软件，也就是第3讲：流程及软件

本来准备直接实战了，但是在看一些新的paper 时候发现我漏掉了hic技术应用的文章解读，我还是需要带领大家看看那些已经发表的好文章到底是如何处理hic数据的。

癌细胞的HIC文章

文章是 ：3D genome of multiple myeloma reveals spatial genome disorganization associated with copy number variations 数据公布在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE87585

北大李程课题组的研究人员比较了骨髓瘤细胞与正常B细胞之间的TAD的差异，在GM12878，RPMI8226与U266三个细胞系中，其分别得到了2756，3457，3342个TAD，其中有1281个TAD在三个细胞系中保守存在，740个TAD特异的存在于两种骨髓瘤细胞系中，这些数据表明在癌症细胞中TAD的结构会发生相当比例的改变，同时TAD的长度变小，数目增多。

作者进一步对骨髓瘤细胞与正常B细胞相比发生compartment改变区域内的基因进行了信号通路富集分析，结果表明，富集到的信号通路与骨髓瘤都密切相关，包括了MAPK，TNF，cytokine-cytokine受体相互作用等信号通路。

虽然本次我们讲解HiC，但事实上这个文章利用的各种数据比较多，包括：

我们关心的HiC数据

主要是4个HiC样本，如下：

GSM2334835: Hi-C U266 MboI; Homo sapiens; OTHER

GSM2334834: Hi-C U266 HindIII; Homo sapiens; OTHER

GSM2334833: Hi-C RPMI-8226 MboI; Homo sapiens; OTHER

GSM2334832: Hi-C RPMI-8226 HindIII; Homo sapiens; OTHER

查看其中一个数据：

数据量不小，想下载全部的4个hic样本来完全重复出来该文章的分析过程及结果对服务器计算资源的考验很大，

其分析结果包括：

然后再看其文章描述的数据处理步骤，作者使用了2013任兵教授的nature文章的数据分析方法：GSE43070. 简单点说，就是：

• all Hi-C sequencing reads were mapped to the human reference genome (hg19) using Bowtie2

• The two ends of paired-end reads were mapped independently using the first 36 bases of each read.

• We filtered out redundant and non-uniquely mapped reads, and kept the reads within 500 bp upstream of enzyme cutting sites (HindIII or Mbol) due to the size selection.

• We utilized the iterative correction and eigenvector decomposition (ICE) method and HiCNorm to normalize raw interaction matrices

好奇怪，里面没有用的hiclib也没有用hicpro软件，而且也没有走完我们第三讲总结好的那些流程。

看看数据处理的中间文件

我尝试下载了 HindIII_HiC_TAD_40kb.tar.gz 文件和HindIII_HiC_ice_matrix_500kb文件并且简单查看，如下:

mkdir -p ~/project/hic/data/myelom
cd ~/project/hic/data/data/myelom  
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2334nnn/GSM2334834/suppl/GSM2334834_U266_HindIII_HiC_TAD_40kb.tar.gz
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2334nnn/GSM2334834/suppl/GSM2334834_U266_HindIII_HiC_ice_matrix_500kb.tar.gz
tar zxvf GSM2334834_U266_HindIII_HiC_TAD_40kb.tar.gz
tar zxvf GSM2334834_U266_HindIII_HiC_ice_matrix_500kb.tar.gz

作者选择了40kb的分辨率来看这个CNVs and TADs 的关系。

使用R包HiTC来看500kb的分辨率下的compartments A/B switches and gene expression的关系。作者给出的分析结果文件是；

resolution_500k/cis/ice_normalization/
|-- [1007K]  chr10_500k_normalized_matrix.txt
|-- [1.0M]  chr11_500k_normalized_matrix.txt
|-- [1.0M]  chr12_500k_normalized_matrix.txt
|-- [621K]  chr13_500k_normalized_matrix.txt
|-- [508K]  chr14_500k_normalized_matrix.txt
|-- [3.1M]  chr1_500k_normalized_matrix.txt
|-- [469K]  chr15_500k_normalized_matrix.txt
|-- [393K]  chr16_500k_normalized_matrix.txt
|-- [373K]  chr17_500k_normalized_matrix.txt
|-- [366K]  chr18_500k_normalized_matrix.txt
|-- [198K]  chr19_500k_normalized_matrix.txt
|-- [229K]  chr20_500k_normalized_matrix.txt
|-- [ 93K]  chr21_500k_normalized_matrix.txt
|-- [ 90K]  chr22_500k_normalized_matrix.txt
|-- [1.4M]  chr23_500k_normalized_matrix.txt
|-- [ 68K]  chr24_500k_normalized_matrix.txt
|-- [3.5M]  chr2_500k_normalized_matrix.txt
|-- [2.3M]  chr3_500k_normalized_matrix.txt
|-- [1.9M]  chr4_500k_normalized_matrix.txt
|-- [1.9M]  chr5_500k_normalized_matrix.txt
|-- [1.7M]  chr6_500k_normalized_matrix.txt
|-- [1.5M]  chr7_500k_normalized_matrix.txt
|-- [1.2M]  chr8_500k_normalized_matrix.txt
`-- [940K]  chr9_500k_normalized_matrix.txt

这些txt文件总共是6206行，乘以500Kb的分辨率，也就对应着人类的3Gb的基因组大小。

其中的任何一个文件，都是可以拿出去画热图的，本身就是一个矩阵，我画21号染色体如下：

代码很简单，但是不知道图对不对以及图后面蕴含的生物学意义。

rm(list=ls())
options(stringsAsFactors = F)
a=read.table('~/GitBook/qc/chr21_500k_normalized_matrix.txt')
library(pheatmap)
pheatmap(a,cluster_rows = F,cluster_cols = F,labels_row = '',labels_col ='')

所以还需要慢慢学。

当然，作者还比较不同细胞系找到的TADs区别。

不过我们后面的实战演练，暂时不使用这个数据集。