包:csaw
作者:Aaron Lun (alun@wehi.edu.au）
编辑日期:2015-07-17
我假设这里的每个人都知道ChIP-seq是什么，但如果不知道，那么你可以看看这张可爱的图片。
为什么要在生物条件之间寻找差异结合(DB) ?
为什么要寻找新创网站?
为什么使用windows?
这是一个相当模块化的过程:
但在我们开始之前，我们必须建立一些数据。我们将使用Tiwari等人2012年的研究。这比较了NF-YA转录因子在胚胎干细胞(ESCs)和末端神经元(TNs)中的结合。
bam。文件<-文件。path("/home/ubuntu/data/aaronlun", c("es_1.bam", "es_2.bam", "tn_1.bam", "tn_2.bam"))这些是BAM文件，因此我们已经将读取映射到mm10小鼠基因组的构建。对于那些感兴趣的人，这是用subread在SAMtools和Picard的帮助下。完整的处理管道可以在下面的位置找到。
系统。File ("doc"， "sra2bam.sh"， package="csaw")## [1] "/home/ubuntu/R-libs/csaw/doc/sra2bam.sh"第一步是计算读取窗口的次数。Reads被定向扩展到平均片段长度，以表示免疫沉淀过程中存在的原始片段。大小恒定(10 bp)的窗口以固定的间隔(50 bp)平铺在整个基因组上。
碎片弹。Len <- 110窗口。宽度<- 10间距<- 50然后在每个库中计算重叠每个窗口的估算片段的数量。我们使用windowCounts函数来进行实际计数。我们过滤掉低质量的比对(比对质量分数小于50)，除了标准染色体组，我们忽略了其他的东西。
param <- readParam(minq=50, restrict=paste0('chr'， c(1:19， 'X'， 'Y'))) data <- windowCounts(bam。文件,ext =碎片弹。len宽度=窗口。Width, spacing=spacing, param=param)数据##类:rangedsummarizeexperiment ## dim: 3950319 4 ##元数据(4):间距宽度shift final。ext ## assays(1):计数## rownames: NULL ## rowRanges元数据列名(0):## colnames: NULL ## colData names(4): bam。文件总数ext参数这给了我们一个RangedSummarizedExperiment我们可以玩的对象。看看rowRanges(数据)，为窗口的基因组坐标;分析(数据)，表示每个窗口的计数;而且colData(数据)，以查阅图书馆的个别资料。
Max.delay <- 500 dedup。on <- readParam(dedup=TRUE, minq=50) x <- correlateReads(bam。(0:max.delay, x, type="l"， ylab="CCF"， xlab="Delay (bp)")
类中可以指定许多有用的读提取参数readParam对象。这些选项允许我们:
strandedCounts。这个想法是定义一个单一的readParam对象，以便在整个分析过程中使用，以确保每次提取相同的读取(有一些例外，例如，correlateReads)．例如，看看我们之前构造的那个。
参数##在单端模式下提取读取##重复删除被关闭##允许的最小映射分数为50 ##从两条链中提取读取##读取提取限制为21个序列##没有指定区域来丢弃读取有几个不同的标准化选项csaw，但我们只关注消除组合偏差。当DB区域在一个库中占据了更大的测序蛋糕时，就会引入这种方法，从而抑制了其他所有内容的覆盖(例如，非DB站点，背景区域)。我们不想检测由于这种抑制而产生的差异，因此我们基于大多数基因组是非db背景的假设进行标准化。首先，我们将读取数计算到10 kbp的大箱子中，以量化背景区域的覆盖范围。
binned <- windowCounts(bam。files, bin=TRUE, width=10000，参数=参数)这些后台库之间的任何系统差异都必须是人为的，应该通过缩放规范化来消除。的正常化函数是一个S4方法normalizeCounts函数，该函数反过来包装calcNormFactors函数从刨边机包中。方法执行实际的规范化m值修剪均值(TMM)法。
Normfacs <- normalize(binned) Normfacs## [1] 1.006 0.973 1.015 1.007我们可以在MA图上看看这些因素。下面的代码比较了库1和库2的分箱计数，即两个ES复制。注意这个大斑点，它主要对应于基因组的背景区域。红线的位置表示对数缩放因子，用于将一个库与另一个库规范化。
adj.counts <- cpm(asDGEList(binned)， log=TRUE) curt .x <- adj.counts[，1] curt .y <- adj.counts[，2] smoothScatter(x=(curt .x+ curt .y)/2+6*log2(10)， y= curt .x-cur. counts)， y=y, xlab="A"， ylab="M"， main="1 vs 2")区域<- diff(log2(normfacs[c(2,1)]))) abline(h=所有。dist,坳=“红色”)
另一个偏见来源是不同文库之间免疫沉淀效率的差异。如果存在这样的差异，我们不会想要检测它们，因为它们在生物学上不有趣。因此，我们可以通过消除高丰度窗口(即绑定位点)中丰度的任何系统差异来规范化。这涉及到对剩余部分进行过滤(请参阅下一节)和规范化。
keep <- aveLogCPM(asDGEList(data)) > -1 normalize(data[keep，])## [1] 0.610 0.799 1.533 1.338注意:我们只能选择一组归一化因子——要么去除成分偏差，要么去除效率偏差。选择取决于实验的生物学背景。具体来说，你是否期望条件之间的绑定发生大规模变化?
在这种情况下，我们会考虑成分偏差，因为在细胞类型之间，NF-YA结合不应该是恒定的。
趋势偏差也可以作为MA图中的趋势观察到。这些问题不能通过像TMM这样的缩放方法来消除，而是需要使用非线性方法。在csaw，这可以通过设置实现类型=“黄土”在对正常化。这将执行基于黄土的归一化以去除趋势，并产生一个偏移矩阵以用于以后的模型拟合。
过滤是使用每个窗口的平均丰度来完成的，如使用aveLogCPM在刨边机。这可以解释为每个窗口的log-(平均计数)/ million。在我们的统计模型中，平均丰度大致独立于每个窗口的DB状态，所以当我们选择高丰度的窗口时，我们避免了数据窥探。
- asDGEList(data))然后，根据某个阈值，我们保留被认为是高丰度的窗口。的RangedSummarizedExperiment对象可以很容易地进行子集来实现这一点(我们将保留原始的，以防万一)。
Keep <- abundance > -1 original <- data data <- data[Keep，] nrow(data)## [1] 13823大量的低丰度窗口是有问题的，因为它们会弄乱一些刨边机的统计方法和假设。过滤消除了大多数这些窗口及其相关问题，同时还减少了计算工作和多次测试校正的严重性。
扔进垃圾箱。2 <- windowCounts(bam。filter.stat.bg <- filterWindows(original, background=binned. files, bin=TRUE, width=2000, param=param)2, type="global")背景。Keep <- filter.stat。Bg $filter > log2(3) summary(background.keep)##模式FALSE TRUE NA的##逻辑3930588 19731 0我们可以模仿呼叫高峰苹果电脑，其中每个窗口的背景在本地估计。在本例中，我们将本地背景定义为每个窗口周围的2 kbp间隔。
surround <- 2000 neighbour <- suppressWarnings(resize(rowRanges(original)， surround, fix="center")) wider <- regionCounts(bam。文件，区域=邻居，ext=碎片。len参数=参数)然后我们使用filterWindows函数将这些处理为丰富值，用于在其本地背景上的每个窗口。我们保留了那些比它的邻居多3倍或更多的窗户。
filter.stat.loc <- filterWindows(original, wider, type="local")Keep <- filter.stat。Loc $filter > log2(3) summary(local.keep)##模式FALSE TRUE NA的##逻辑3940821 9498 0我们首先定义实验设计。这对于我们的设置来说是相当容易的，它包含两组两个生物重复。更复杂的设计也可以容纳。
分组< -因子(c(“西文”、“西文”,tn, tn))设计< -模型。矩阵(~0 +分组)colnames(设计)<- levels(分组)我们使用负二项式(NB)框架在刨边机包中。这解释了在生物重复之间表现出过度分散的低离散计数。具体地说，重复之间的变异性是使用NB分散参数建模的，该参数通过estimateDisp函数。
y <- asDGEList(data, norm.factors=normfacs) y <- estimateDisp(y, design)我们可以用准似然(QL)方法。QL分散解释了窗口特定的变异性，而NB分散则模拟了重复之间的生物变异性。我们为每个窗口的计数拟合一个广义线性模型(GLM)，其中该窗口的QL离散度由GLM偏差估计。
fit <- glmQLFit(y, design, robust=TRUE)最后，我们使用QL F-test来测试两种单元格类型之间的DB。
对比<- makeContrasts(es - tn, levels=design) results <- glmQLFTest(fit，对比度=对比度)norep。fit <- glmFit(y, design, dispersion=0.05)结果<- glmLRT(norep。健康,对比=对比)健壮的= TRUE？箱子adc <- cpm(asDGEList(bin .2)， log=TRUE) plotMDS(bin. 2)时会有,标签=分组)
这不是很高级，但是我们可以画出NB的分散plotBCV命令。在QL框架中，我们只关注趋势，因为只有趋势NB分散将在以后使用。我们通常会看到丰度下降的趋势，可能会在非常大的丰度处稳定下来。
plotBCV (y)
类似地，对于QL分散，我们可以使用plotQLDisp函数。这表明，当原始值被压缩到趋势值时，经验贝叶斯收缩对估计的影响。我们对所有下游加工都使用缩水的估算。其思想是减少存在有限重复的不确定性，以提高DB测试的检测能力。
plotQLDisp(适合)
如果你有更复杂的实验(例如，多组，阻塞因素)，你可能需要更仔细地考虑设计矩阵和对比。请查看刨边机如何处理它们的用户指南。
为了纠正跨越数千个窗口的多次测试，我们可以考虑控制错误发现率(FDR)。然而，在所有窗口中控制FDR并不是那么有用，因为很少有人会从窗口的角度来解释ChIP-seq结果。相反，我们想要控制所有基因组区域的FDR，即区域层次上罗斯福。
为了做到这一点，我们可以使用mergeWindows函数。小于托尔分开的被扔进同一个集群。合并相邻或重叠的窗口可以减少冗余并简化对结果的解释。每个簇现在代表一个不同的基因组区域。
cluster <- mergeWindows(rowRanges(data)， tol=1000) cluster $region##具有4238个范围和0个元数据列的GRanges对象:## seqnames ranges strand ## ## [1] chr1 [6466701, 6466760] * ## [2] chr1 [7088951, 7088960] * ## [3] chr1 [7397901, 7398110] * ## [5] chr1 [9545251, 9545360] * ## ... ... ... ...## [4234] chrY [259151, 259360] * ## [4235] chrY [90761101, 90761110] * ## [4236] chrY [90805151, 90805160] * ## [4237] chrY [90808801, 90808910] * ## [4238] chrY [90812101, 90813810] * ## ------- ## seqinfo:来自未指定基因组的21条序列 我们可以结合p-values为每个集群中的所有窗口使用Simes的方法。这是计算单个组合p-value为每个集群/区域，表示针对全局空值的证据。我们还得到了一些其他的统计数据，比如每个集群中的窗口总数，以及在每个方向上发生重大变化的数量。
tabcom <- combineTests(群集$id，结果$表)头(tabcom)## nWindows logFC。logFC。down PValue FDR ## 12 00 0.7041 0.7698 ## 21 01 0.0179 0.0527 ## 3 5 05 0.0800 0.1462 ## 4 3 0 3 0.0442 0.0961 ## 5 3 0 3 0.0589 0.1181 ## 62 0 2 0.0121 0.0409输出表的每一行对应一个集群/区域。控制跨区域的FDR很简单，因为我们可以对组合应用Benjamini-Hochberg校正p值。
托尔？我们如何避免连锁效应?clustered2 <- mergeWindows(rowRanges(data)， tol=100, max.width=5000)选项卡。getBestTest(clustered$id, results$table)头(tab.best)## best logFC logCPM F PValue FDR ## 12 -0.295 -0.9971 0.155 1.0000 1.0000 ## 2 3 -1.858 -0.8812 6.091 0.0179 0.0546 ## 3 6 -1.211 1.1286 5.999 0.0934 0.1716 ## 4 10 -1.454 -0.0262 5.291 0.0798 0.1524 ## 5 13 -1.514 -0.7097 4.550 0.1169 0.2025 ## 6 16 -2.204 -0.8219 8.383 0.0121 0.0432我们可以使用许多不同的集群策略，只要它们对窗口的DB状态是不可见的。一种特别有用的策略是根据窗口所在的启动子对窗口进行聚类。这允许与差分表达式分析进行直接协调。
基因。体<-基因(TxDb.Mmusculus.UCSC.mm10.knownGene) prom <- promoter(基因。船体，上游=3000，下游=1000##具有6个范围和1个元数据列的GRanges对象:## seqnames ranges strand | gene_id ## | ## 100009600 chr9 [21074497, 21078496] - | 100009600 ## 100009609 chr7 [84963010, 84967009] - | 100009609 ## 100009614 chr10 [77708446,77712445] + | 100009614 ## 100009664 chr11 [45805083,45809082] + | 100009664 ## 100012 chr4 [144161652,134771004] - | 100017 ## ------- ## seqinfo:mm10基因组66个序列(1个圆形) 然后我们识别出重叠每个启动子的窗口。每个基因启动子的窗口集被定义为一个簇。
olap <- finoverllaps (prom, rowRanges(data)##命中12282次，0元数据列的对象:## queryHits subjectHits ## ## [1] 7 8461 ## [2] 7 8462 ## [3] 18 10101 ## [4] 18 10102 ## [5] 18 10103 ## ... ... ...## [12278] 23643 7903 ## [12279] 23643 7904 ## [12280] 23646 8274 ## [12281] 23646 8275 ## [12282] 23646 8276 ## ------- ## queryLength: 23653 ## subjectLength: 13823 这是一个简单的问题，结合p-values，就像我们之前做的那样。这是使用包装器函数完成的，该包装器函数直接接受支安打对象。我们最终得到一个表，其中的每个条目都有一行毕业舞会。
tabbroad <- combineOverlaps(olap, results$table) head(tabbroad[!is.na(tabbroad$PValue)，])## nWindows logFC。logFC。down PValue FDR ## 7 20 2 0.000120 0.00261 ## 18 5 05 0.000924 0.00744 ## 22 20 2 0.179911 0.24410 ## 25 40 0 0.935199 0.95201 ## 28 3 0 3 0.005547 0.02282 ## 36 40 4 0.016083 0.04409这意味着我们可以对启动子区域的DB做出直接的陈述。注意NA的项表示毕业舞会它们没有重叠的窗口，这里没有显示出来，因为它们不是特别令人兴奋。
属性可以向区域添加简单的基于基因的注释detailRanges函数。这可以识别每个区域重叠的注释基因组特征(外显子、启动子、内含子)，以及每个区域侧翼的基因组特征(加上每个侧翼特征与区域之间的间隙距离)。
anno <- detailRanges(clustered$region, orgdb=org. mm . example .db, txdb= txdb . mmusculus . ucc .mm10. knowngene) head(anno$overlap)# #[1]”“Pcmtd1 | 0 - 1 | +”# # (3 ] "" "" ## [ 5]“Rrs1 | 0 | + Adhfe1 | 0 | +”“Vcpip1 | 0 | - 1700034 p13rik | 0 - 1 | +”这里可能需要一些解释。每个字符串表示一个逗号分隔的列表，其中每个条目的格式为符号|外显子|链，即相关基因的符号，重叠的外显子(我对于内含子，和0对于启动子)和基因链。对于侧翼区域，另加(距离)字段表示注释的特征与每个区域之间的距离。
头(伊斯兰教纪元左)美元## [1] "" "" "" "" ## [5] "" Vcpip1|1|-[19]"现在我们可以把一个表保存到磁盘上，供其他人阅读等等。当然，如果你愿意，你可以保存为其他格式，使用包rtracklayer。
所有人。结果<- data.frame(as.data.frame(clustered$region)[，1:3]， tabcom, anno)所有。results <- all.results[order(all.results$PValue)，]结果,文件= "保存。tsv"， row.names=FALSE, quote=FALSE, sep="\t")控件可以执行简单的可视化Gviz包中。方法从BAM文件中提取读取extractReads函数，我们以每百万阅读数为单位绘制覆盖率。分别用蓝色和红色标记每条线。
current .region <- GRanges("chr18"， IRanges(77806807,77807165)) collected <- list() for (i in 1:length(bam.files)) {reads <- extractReads(bam. files)files[i]， curm .region) pcov <- as(coverage(reads[strand(reads)=="+"])/data$totals[i]*1e6， "GRanges") ncov <- as(coverage(reads[strand(reads)=="-"])/data$totals[i]*1e6， "GRanges") ptrack <- DataTrack(pcov, type="直方图"，lwd=0, fill=rgb(0,0,1，.4)， ylim=c(0,1)， name=bam。files[i]， col.axis="black"， col.title="black") ntrack <- DataTrack(ncov, type="histogram"， lwd=0, fill=rgb(1,0,0，.4)， ylim=c(0,1)) collected[[i]] <- OverlayTrack(trackList=list(ptrack,ntrack))} gax <- GenomeAxisTrack(col="black") plotTracks(c(gax, collected)， from=start(curt .region)， to=end(curt .region)))
anno2 <- detailRanges(clustered$region, orgdb=org. mm . exe .db, txdb= txdb . mmusculus . ucsc .mm10. txt)knownGene, dist=3000，启动子=c(0,2000))将结果存储为农庄为了进一步的操作。这可以通过将各种统计信息填充到元数据字段中来实现。
所有人。regions <- clustered$region elementMetadata(all.regions) <- tabcom . tabcom . txt您刚刚完成了一个基于窗口的操作新创微分结合分析。有关更多信息，您可以查看csaw用户指南，或每个功能的文档。
csawUsersGuide ()作为奖励，你可以猜一猜含有蛋白质的食物，如图所示。