这一章主要讲的是Chip-seq的基本分析方法
Chip-seq
（引用原文）
ChIP-seq是使用染色质免疫共沉淀技术将蛋白质结合的DNA样本测序后分析全基因组蛋白质结合位点的分析方法。研究对象包括转录因子，组蛋白修饰分析
ChIP-seq实验主要分为如下几步：
1.Sample fragmentation, 将基因组打断
2.Immunoprecipitation, 免疫共沉淀
3.DNA purification, DNA提纯
4.sequence, 测序
通常打断的基因组有150-300bp长，那么我们测的是蛋白或者TF结合的那个片段，当这些reads被mapping回到基因组上时，就会形成一个peak，那么peak的位置即为结合的部分
常用的软件有：
1.mapping： bowtie1/2, Rsubread
2.peak calling：CCAT, SICER, MACS, ZINBA, BayesPeak, chipseq, ChIPseqR, CSAR, csaw, GenoGAM, iSeq, PICS
3.Peak annotation and analysis: ChIPpeakAnno
4.Gene network building: GeneNetworkBuilder, ChIPXpress
5.Motif enrichment analysis: The MEME Suite, Homer
质控
当我们完成上游分析，我们可以在R里面进行QC
library(ChIPQC)
experiment <- ChIPQC(samples)
ChIPQCreport(experiment)
峰注释
我们call完peak以后，我们也许知道peak的位置，但是却不知道是哪些位置的（TSS）或者功能等，不过我们要下载好注释好的包
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
annoData <- toGRanges(TxDb.Hsapiens.UCSC.hg19.knownGene, feature="gene")
#导入peak文件
library(ChIPpeakAnno)
packagePath <- system.file("extdata", package = "ChIPpeakAnno")
dir(packagePath, "gff|bed|narrowPeak|broadPeak|gz")
#导入峰的GFF文件
toGRanges(file.path(packagePath, "GFF_peaks.gff"), format = "GFF")
#导入bed文件
toGRanges(file.path(packagePath, "WS220.bed"), format = "BED")
#导入narrowpeak
toGRanges(file.path(packagePath, "peaks.narrowPeak"), format = "narrowPeak")
#导入broadpeak
toGRanges(file.path(packagePath, "TAF.broadPeak"), format = "broadPeak")
上述介绍了如何导入我们需要的文件
接下来就可以进行注释了
#以narrowpeak为例
packagePath <- system.file("extdata", "macs", package = "ChIPseqStepByStep")
macs2.files <- dir(packagePath, pattern="*.q1_peaks.narrowPeak$")
peaks <- sapply(macs2.files, function(.ele)
toGRanges(file.path(packagePath, .ele), format="narrowPeak"))
# 寻找最近的TSS
anno <- annotatePeakInBatch(gr1,
AnnotationData=annoData)
head(anno, n=2)
差异peak分析
类似于RNA-seq的差异分析一样，我们可以寻找差异binding site
DiffBind
文件是根据实验设计和数据存储地址等信息创建的一个csv格式文件，包含的表头信息有"SampleID"、 "Tissue"、 "Factor"、 "Condition" 、"Treatment"、"Replicate" 、"bamReads" 、"ControlID"、 "bamControl" "Peaks"、 "PeakCaller"（bam,peak文件分别在比对和call peak的步骤产生）
//www.greatytc.com/p/f849bd55ac27
library(DiffBind)
## 准备好样品文件
samples <- read.csv(file.path(system.file("extra", package="DiffBind"),
"tamoxifen.csv"))
samples[1:2, ]
pf <- system.file("extra", package="DiffBind")
samples$bamReads <- file.path(pf, samples$bamReads)
samples$bamControl <- file.path(pf, samples$bamControl)
samples$Peaks <- file.path(pf, samples$Peaks)
tmpfile <- tempfile()
write.csv(samples, tmpfile)
##读取文件
tamoxifen <- dba(sampleSheet=tmpfile)
##查看样品间的关系
plot(tamoxifen)
##计数
tamoxifen <- dba.count(tamoxifen, summits=250) ##只对峰上下游250bp(总计500bp)进行计数
##差异分析
tamoxifen <- dba.contrast(tamoxifen, categories=DBA_CONDITION) ##比对的条件为samples中的Condition列
##因为没有bam文件，所以这个例子其实并不能跑起来
tamoxifen <- dba.analyze(tamoxifen)
##提取结果
tamoxifen.DB <- dba.report(tamoxifen)
这个包的详细教程：http://www.bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf
其中输入文件为：
csaw
library(csaw)
## 1. Loading in data from BAM files.
library(csaw)
param <- readParam(minq=50)
data <- windowCounts(bam.files, ext=110, width=10, param=param)
## 2. Filtering out uninteresting regions.
library(edgeR)
keep <- aveLogCPM(asDGEList(data)) >= -1
data <- data[keep,]
## 3. Calculating normalization factors.
binned <- windowCounts(bam.files, bin=TRUE, width=10000, param=param)
data <- normOffsets(binned, se.out=data)
## 4. Identifying DB windows.
y <- asDGEList(data)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design, robust=TRUE)
results <- glmQLFTest(fit)
## 5. Correcting for multiple testing.
merged <- mergeWindows(rowRanges(data), tol=1000L)
tabcom <- combineTests(merged$id, results$table)
