染色体状态分类及注释(ChromHMM)¶
作者：中国科学院分子细胞科学卓越创新中心
基因组中非编码区域在生理、进化、疾病等过程中起着重要的调节作用。表观基因组标记的映射，如组蛋白修饰、组蛋白变异、开放染色质区域和相关标记，已经成为一种注释基因组的强有力的手段，可用来识别潜在的调控元件以及研究它们在不同细胞类型和人类疾病中的活性变化。本文描述了如何使ChromHMM软件来对整个染色体区域(包含编码区和非编码区)做状态分类、富集用以寻找具有生物学意义的基因组区间，并使用HOMER来对不同染色体状态进行基因注释，然后比较不同状态下基因的表达谱。使用的测试数据集为IMR90细胞系的19种表观数据和表达量数据。
预处理¶
根据基因组染色体大小(chrlen.txt)、组蛋白标记的样本表格(design_sheet.txt)和样本比对生成的bam文件，使用ChromHMM的子命令BinarizeBam将bam文件转化为二值化的信号文件：
chrlen.txt格式类似如下：
chr1 248956422 chr2 242193529 chr3 198295559 chr4 190214555 chr5 181538259 chr6 170805979 chr7 159345973 chr8 145138636 chr9 138394717 chr10 133797422design_sheet.txt格式类似如下：
IMR90 H3K4me1 GSM521895_sort_uniq.bam IMR90 H3K4me3 GSM521901_sort_uniq.bam IMR90 H3K9me3 GSM469974_sort_uniq.bam IMR90 H3K27ac GSM469966_sort_uniq.bam IMR90 H3K27me3 GSM469968_sort_uniq.bam IMR90 H3K36me3 GSM521890_sort_uniq.bam IMR90 H3K9ac GSM469973_sort_uniq.bam IMR90 H3K4ac GSM521893_sort_uniq.bam IMR90 H3K14ac GSM521881_sort_uniq.bam IMR90 H3K18ac GSM521884_sort_uniq.bam IMR90 H3K23ac GSM521885_sort_uniq.bam IMR90 H3K56ac GSM521902_sort_uniq.bam IMR90 H4K5ac GSM469975_sort_uniq.bam IMR90 H2AK5ac GSM521866_sort_uniq.bam IMR90 H2BK120ac GSM521869_sort_uniq.bam IMR90 H2BK12ac GSM521871_sort_uniq.bam IMR90 H2BK20ac GSM521879_sort_uniq.bam IMR90 H3K4me2 GSM521899_sort_uniq.bam IMR90 H3K79me1 GSM521904_sort_uniq.bam所有输入样本的bam文件置于一个文件夹下，然后将该路径作为输入，目录格式如下：
system('ls /sibcb2/bioinformatics2/chenyulong/testdir/19.test_chromhmm/bamfile',intern=T) [1] "GSM469966_sort_uniq.bam" "GSM469968_sort_uniq.bam" [3] "GSM469973_sort_uniq.bam" "GSM469974_sort_uniq.bam" [5] "GSM469975_sort_uniq.bam" "GSM521866_sort_uniq.bam" [7] "GSM521869_sort_uniq.bam" "GSM521871_sort_uniq.bam" [9] "GSM521879_sort_uniq.bam" "GSM521881_sort_uniq.bam" [11] "GSM521884_sort_uniq.bam" "GSM521885_sort_uniq.bam" [13] "GSM521890_sort_uniq.bam" "GSM521893_sort_uniq.bam" [15] "GSM521895_sort_uniq.bam" "GSM521899_sort_uniq.bam" [17] "GSM521901_sort_uniq.bam" "GSM521902_sort_uniq.bam" [19] "GSM521904_sort_uniq.bam"将bam文件转化为信号的二值化文件，命令行代码如下：
java -mx4000M -jar ChromHMM.jar BinarizeBam -gzip -b 200 -f 0 -g 0 -p 0.0001 chrlen.txt design_sheet.txt bamfile binarization
结果如下：binarization/ ├── IMR90_chr10_binary.txt.gz ├── IMR90_chr11_binary.txt.gz ├── IMR90_chr12_binary.txt.gz ├── IMR90_chr13_binary.txt.gz ├── IMR90_chr14_binary.txt.gz ├── IMR90_chr15_binary.txt.gz ├── IMR90_chr16_binary.txt.gz ├── IMR90_chr17_binary.txt.gz ├── IMR90_chr18_binary.txt.gz ├── IMR90_chr19_binary.txt.gz ├── IMR90_chr1_binary.txt.gz ├── IMR90_chr20_binary.txt.gz ├── IMR90_chr21_binary.txt.gz ├── IMR90_chr22_binary.txt.gz ├── IMR90_chr2_binary.txt.gz ├── IMR90_chr3_binary.txt.gz ├── IMR90_chr4_binary.txt.gz ├── IMR90_chr5_binary.txt.gz ├── IMR90_chr6_binary.txt.gz ├── IMR90_chr7_binary.txt.gz ├── IMR90_chr8_binary.txt.gz ├── IMR90_chr9_binary.txt.gz ├── IMR90_chrX_binary.txt.gz └── IMR90_chrY_binary.txt.gz
染色体状态的分类及富集¶
以上一步得到的二值化的信号文件作为输入，使用LearnModel子命令来完成染色体的分类及富集分析，命令行代码如下：
java -mx4000M -jar ChromHMM.jar LearnModel -gzip -d 0.001 -color 0,0,255 -p 5 -i chrhmm binarization learnmodel
结果如下：
state18/learnmodel
├── emissions_18_chrhmm.png #染色体状态的分类热图，png格式
├── emissions_18_chrhmm.svg #染色体状态的分类热图，svg格式
├── emissions_18_chrhmm.txt #染色体状态分类热图的作图数据
├── IMR90_18_chrhmm_dense.bed.gz #可在IGV浏览器中展示分类结果的bed文件，所有状态在同一个track里面
├── IMR90_18_chrhmm_expanded.bed.gz #可在IGV浏览器中展示分类结果的bed文件，不同状态在不同的track里面
├── IMR90_18_chrhmm_overlap.png #染色体状态的富集热图，png格式
├── IMR90_18_chrhmm_overlap.svg #染色体状态的富集热图，svg格式
├── IMR90_18_chrhmm_overlap.txt #染色体状态富集热图的作图数据
├── IMR90_18_chrhmm_RefSeqTES_neighborhood.png #染色体状态在TES区域的富集热图，png格式
├── IMR90_18_chrhmm_RefSeqTES_neighborhood.svg #染色体状态在TES区域的富集热图，svg格式
├── IMR90_18_chrhmm_RefSeqTES_neighborhood.txt #染色体状态在TES区域的富集热图的作图数据
├── IMR90_18_chrhmm_RefSeqTSS_neighborhood.png #染色体状态在TSS区域的富集热图，png格式
├── IMR90_18_chrhmm_RefSeqTSS_neighborhood.svg #染色体状态在TSS区域的富集热图，svg格式
├── IMR90_18_chrhmm_RefSeqTSS_neighborhood.txt #染色体状态在TSS区域的富集热图的作图数据
├── IMR90_18_chrhmm_segments.bed.gz #染色体状体的分类结果的bed文件
├── model_18_chrhmm.txt #包含learnmodel子命令使用参数的文件
├── transitions_18_chrhmm.png #状态转化的热图，png格式
├── transitions_18_chrhmm.svg #状态转化的热图，svg格式
├── transitions_18_chrhmm.txt #状态转化热图的作图数据
└── webpage_18_chrhmm.html #染色体状态分类的网页报告染色体分类热图：
染色体状态的分类统计¶
染色体状态的注释
使用homer软件的子命令annotatePeaks.pl来对注释染色体状态，示例代码如下：
annotatePeaks.pl peakbed fasta -gtf gtf -cpu 5 -gid >prefix_annostate.txt
结果如下：PeakID Chr Start End Strand Annotation Nearest PromoterID Gene Name E15-10799 chr14 103473601 103475600 + intron ENSG00000166165.13 CKB E8-66258 chr7 79987001 79987200 + intron ENSG00000206818.1 RNU6-849P E1-20958 chr2 80861801 80864400 + intron ENSG00000230975.1 AC012355.1 E8-60630 chr6 37210601 37210800 + Intergenic ENSG00000137193.14 PIM1 E2-55960 chr2 65287601 65291400 + Intergenic ENSG00000138071.13 ACTR2 E13-57227 chr22 44019601 44021200 + intron ENSG00000188677.14 PARVB E18-5662 chr15 58932201 58932400 + intron ENSG00000137776.17 SLTM E11-41293 chr2 200610201 200610400 + intron ENSG00000138356.14 AOX1 E9-60679 chr3 149461801 149462200 + TTS ENSG00000282502.1 FKBP1AP4染色体各状态的统计
使用脚本chromhmm_box_hist.R来统计不同染色体状态下基因的表达量，以及不同染色体状态下功能结构分布，示例代码如下：
Rscript chromhmm_box_hist.R -a annostate.txt -e expression.txt [-o outdir]
染色体不同状态下基因的表达量分布：
染色体不同状态下功能结构的分布：
结果说明¶
从分析目的来说，ChromHMM利用不同组蛋白标记将染色体分成潜在的有生物学意义的状态，例如上面的示例，通过19种组蛋白标记将每一条染色体分成18种不同的状态(具体见上面的染色体分类热图)，然后结合每一种状态在基因组功能区域的富集情况来预测染色体状态可能存在的生物学意义。如已知H3k4me3修饰主要发生在基因组的启动子区，而H3k27me3修饰除了发生在启动区外也主要集中在基因富集区域，H3K4me3、H3k27me3修饰的主要生物学功能分别是活化和阻遏基因的表达。从上面的染色体分类热图可以看出状态17、18中的H3k4me3修饰信号值很强，同时富集结果也显示这两种状态在启动子区域有较强的信号，故可以预测这两种状态可能与基因的高表达有关，而状态4中H3k27me3修饰有很强的信号，同时也在启动子区域和基因body区有较强的富集信号，故可以预测该状态可能与基因的沉默有关。同时，从表达量的分布箱线图上也可以看出状态17、18中基因的表达量也比较高而状态4的表达量则比较低，结合功能分布柱状图来看，状态17、18包含大量的启动子区域而状态4主要包含大量的启动子区及基因富集区域，这些结果与前面的分类富集结果也比较吻合，也进一步说明了染色体分类结果的有效性。