ChIP-seq 基本原理 染色体免疫共沉淀（ChIP ,Chromatin Immunoprecipitation）结合高通量测序技术（NGS）用于分析转录因子，以及其他染色质相关蛋白影响生物体表型变化的机制。研究DNA与蛋白质互作如何调节基因表达对于理解许多生物学过程以及疾病的发生机制具有重要意义。ChIP-seq当前主要应用于两个方面：一方面，转录因子与DNA序列的相互作用识别位点（启动子，增强子等各种顺式作用元件）的研究；另一方面主要应用在表观遗传学领域，组蛋白修饰与核小体定位等。
对于组蛋白修饰语ChIP-seq的研究 。组蛋白修饰是基因激活或者抑制的重要指示器，组蛋白修饰存在于基因表达的调控区域如启动子，增子等。基因的激活表达一般与转录起始位点组蛋白H3K4me3和H3K27ac修饰密切相关，且活跃的增强子能够被H3K4me1和H3K27ac富集鉴定。活跃转录基因的基因体与H3K36me3相关。基因抑制能够通过包含H3K9me3和H3K27me3的两种不同机制介导。基因的表达与否以及表达水平是通过多个组蛋白修饰相互拮抗共同作用的结果。利用染色体免疫共沉淀对已知具有激活或抑制功能的组蛋白修饰marker进行靶定，根据相应组蛋白修饰表达水平的上调或下调来研究预测基因的激活或者抑制。
ChIP-seq 建库测序 染色体免疫共沉淀建库测序流程
ChIP-seq 建库流程 a. 甲醛交联整个细胞系使得目标蛋白与染色质紧密连接
b. 使用超声波或者MNase酶降解等方法将基因组DNA打成小片段
c. 添加与目标蛋白特异的抗体进行孵育，与目标蛋白形成免疫沉淀结合复合物结合在beads上
d. 沉淀洗脱，富集与目标蛋白结合的DNA序列
e. 蛋白酶K解交联，纯化DNA
f-g. PCR扩增准备测序
CHIP-seq数据分析 基本流程
ChIP-seq 分析流程 （1） 质控。与RNA-seq数据分析流程类似，ChIP-seq获得的原始数据首先需要使用fastQC进行质量检测，并对不合格的数据进行进一步处理
（2） Mapping reads。Bowtie将测序获得的reads序列比对到参考基因组上，获得sam或者bam格式的比对文件，获得reads的位置信息。使用samtools软件对比对结果进行排序，以及去除PCR重复，获得最终可以进行下游分析的数据
（3） Peak calling。使用MACS2获得ChIP-seq富集信息,生成bed文件中详细列出了每一个reads的位置信息与可信度
（4） IGV可视化基因位置信息.为了 减少bam文件的大小从而在IGV可视化过程中降低计算机资源的消耗，需要将bam文件转化为BigWig文件或者tdf格式文件
（5） 下游信息挖掘。使用R-packages Chipseeker对处理好的数据进行进一步的信息挖掘与可视化；Metascape进行GO Annotation
分析代码 1、数据获取 fastq- dump -- split - 3 -gzip SRR*
2、质控检测 fastqc -t 20 -o $CHIPseq /fastqc $CHIPseq /rawdata/MCF *.fastq.gz 3、匹配基因组序列 bowtie-build /picb/epigenome/data/genome/hg19/hg19.fa index #构建比对模板 #保存文件$CHIPseq/Align/bowtie.sh for i in `cat id`; do gzip -d $CHIPseq /rawdata/ $i .fastq.gz ; done #将下载的fastq.gz文件解压（-d） for i in `cat id`; do bowtie -S -q -p 20 /picb/epigenome/data/genome/bowtie_indexes/hg19 $CHIPseq /rawdata/ $i .fastq > $i .sam ; done #-S 输出sam文件 #-q 输入文件 for i in `cat id`; do samtools sort $i .sam > $i .bam; done #将sam文件进行排序并生成bam文件输出 #-U 表示为single单端测序 #-x 后面输入index文件 #-S 表示输出的sam文件 for i in `cat id`; do samtools rmdup $i .bam rm $i .bam; done #去重复 for i in `cat id`; do samtools index rm $i .bam ; done #生成后缀为.bai的index文件用于后面macs2进行peak calling
4、Call Peaks 本实验数据主要是在研究组蛋白H3K27ac修饰对于调节基因的表达情况，组蛋白修饰作为表观遗传学研究的重要标签对于基因表达的激活或抑制具有重要的指示作用。H3K27ac修饰是非常重要的基因激活表达标签，对数据进行peak-calling就是要找出H3K27ac修饰增加的位点从而找到相对于正常的乳腺细胞，在乳腺癌细胞中表达上调的基因。
对于热点的定义有着严格的统计学意义，这些热点位置多次被测得的reads所覆盖且可近似认为其服从泊松分布.
macs2 callpeak -t rmMCF7H3K27ac_N1 .bam rmMCF7H3K27ac_N2 .bam -c rmMCF7_IN_N .bam -g hs -n H3K27 macs2 callpeak -t rmMCF7_CTCF_N1 .bam rmMCF7_CTCF_N2 .bam -c rmMCF7_IN_N .bam -g hs -n CTCF 5、Peak Annotation library (ChIPseeker) #读入数据 CTCF<-readPeakFile( 'CTCF_peaks.narrowPeak' ,head= F ) CTCF<-CTCF[,c( 1 , 6 )] H3K27<-readPeakFile( 'H3K27_peaks.narrowPeak' ,head= F ) H3K27<-H3K27[,c( 1 , 6 )] peaklist<-list( 'CTCF' =CTCF, 'H3K27' =H3K27) #对promoter进行注释 library (TxDb.Hsapiens.UCSC.hg19.knownGene) promoter<-getPromoters(TxDb =TxDb.Hsapiens.UCSC.hg19.knownGene ) #把peak Align到窗口上生成矩阵；单独作图 CTCFMatrix<-getTagMatrix(CTCF,windows = promoter) H3K27Matrix<-getTagMatrix(H3K27,windows = promoter) #distribution of promoter tagHeatmap(CTCFMatrix,xlim = c(- 1000 , 1000 ),color = '#e79686' ,title = 'CTCF' ) tagHeatmap(H3K27Matrix,xlim = c(- 1000 , 1000 ),color = '#E7475E' ,title = 'H3K27' ) #合并比较作图 pdf( 'Promoter Distribution.pdf' ) peakHeatmap(peaklist,weightCol = 'V9' ,TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene,upstream = 4000 ,downstream = 4000 ,color = '#4e709d' ) dev.off() #degree of combination展示结合强度,peak结合的数目越多说明结合强度越大 plotAvgProf(CTCFMatrix,xlim = c(- 1000 , 1000 ),xlab= 'Genomic Region (5`->3`)' ,ylab= 'CTCF Read Count frequency' ) plotAvgProf(H3K27Matrix,xlim = c(- 1000 , 1000 ),xlab= 'Genomic Region (5`->3`)' ,ylab= 'h3k27 Read Count frequency' ) #合并作图 tagMatrixlist<-lapply(peaklist,getTagMatrix,windows=promoter) pdf() plotAvgProf(tagMatrixlist,xlim = c(- 1000 , 1000 ),facet = 'col' ) dev.off() #读入peak数据 library (ChIPseeker) CTCF<-readPeakFile( 'CTCF_peaks.narrowPeak' ,head= F ) H3K27<-readPeakFile( 'H3K27_peaks.narrowPeak' ,head= F ) #生成注释文件 CTCFAnno <-annotatePeak(CTCF,tssRegion = c(- 2000 , 2000 ),TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene) H3K27Anno <-annotatePeak(H3K27,tssRegion = c(- 2000 , 2000 ),TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene) Annolist<-list( 'CTCF' =CTCF, 'H3K27' =H3K27) #绘制饼图可视化peakannotation分布 plotAnnoPie(CTCFAnno) plotAnnoPie(H3K27Anno) #multiply Annotation/peak vennpie(CTCFAnno) vennpie(H3K27Anno) pdf( 'CTCFAnno.pdf' ,width = 15 ,height = 15 ) upsetplot(CTCFAnno,vennpie= TRUE ) dev.off() pdf( 'H3K27Anno.pdf' ) upsetplot(H3K27Anno,vennpie= TRUE ) dev.off() #the nearist gene annotation plotDistToTSS(CTCFAnno,title = 'Distribution of transtcription factor-binding loci\nrelative to TSS $CTCF' ) plotDistToTSS(H3K27Anno,title = 'Distribution of transtcription factor-binding loci\nrelative to TSS $H3K27' ) #一张图展示数据 peakAnnoList<-lapply(peaklist,annotatePeak,TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene,tssRegion=c(- 1000 , 1000 )) plotDistToTSS(peakAnnoList) ###富集分析绘制气泡图 library (clusterProfiler) genes=lapply(peaklist, function (i) seq2gene(i,c(- 1000 , 1000 ), 1000 ,TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene)) cc=compareCluster(geneClusters = genes,fun= 'enrichKEGG' ,organism= 'hsa' ) dotplot(cc,showCategory= 10 )
结果与分析 Promoter Distribution 组蛋白H3K27ac作为激活基因表达的组蛋白Marker一搬位于启动子区域用于调控激活基因的表达，CTCF是一种转录因子其结合位点一般位于绝缘子区域。故对序列比对之后的reads在启动子区域的分布情况进行可视化展示，如Figure11。对图A展示了CTCF转录因子与H3K27ac修饰在启动子周围的分布。他们有着不同的分布模式，H3K27ac在TSS转录起始位点分布比较集中且覆盖范围较广；CTC转录因子在TSS区域也存在结合但是结合位点分布范围较窄；图B展示了CTCF转录因子与H3K27ac修饰在TSS附近的结合强度。CTCF转录因子在TSS附近启动子结合强度较大且集中，在TSS附近结合的范围较窄；H3K27ac在TSS附近一个相对比较宽的区域都有较强的结合.图C与图B传递了相同的信息，CTCF转录因子与H3K27ac修饰相对于TSS不同距离的分布情况，相比较于CTCF转录因子，H3K27ac修饰在TSS区域有更广的结合范围与更强的结合强度.
Pie plot 对CTCF转录因子结合位点与H3K27ac修饰位点在基因结构上分布情况进行注释；A图为CTCF转录因子结合位点的分布，大部分注释在远端基因间区，启动子区；B图为H3K27ac修饰结果与CTCF转录因子相似，但是在启动子区域注释更加显著。
UpSet plot 考虑同一个基因可能存在多个reads注释到基因结构的不同位点的情况，只是进行饼图的绘制不能够详细的体现这种情况的发生。进一步做了分布位点注释图Figure13。图A表示CTCF转录因子结合位点 的注释情况，纵轴表示注释到的基因数目，横轴相应的基因结构位点 。比如对于A图CTCF转录因子结合位点有3511个基因同时在5`UTR，启动子区，外显子和内含子区域都有相应的reads进行注释。
Bubble plot for GO enrichment
对获得的数据进行功能的注释，该热图展示了CTCF转录因子与H3K27ac修饰所调节基因的功能富集情况。气泡的大小代表所占基因的比例，气泡越大说明所富集的基因数目越多；气泡的颜色表示可信度，红色代表可信度较高，蓝色代表具有较低的可信度。
CTCF转录因子所调控的基因与H3K27ac修饰调控的基因表达在人类乳头瘤病毒感染，MAPK信号通路中都有较高的富集。
