1概述

本文档描述了如何使用DEScan2来检测表观基因组测序数据中差异富集的区域。DEScan2是一个基于R/Bioconductor的工具，为Bioconductor v3.6和R v3.4.3开发。

2一个典型的差分富集分析工作流程

首先，让我们加载包并设置串行计算。

库("DEScan2")库("RUVSeq")库("edgeR") BiocParallel::register(BiocParallel::SerialParam())

# #示例数据

该分析将使用小鼠19号染色体的小鼠Sono-seq数据。该数据是通过情境恐惧条件反射学习后从小鼠海马中获得的数据的一小部分。试验设4个恐惧条件组(FC)和4个家养对照(HC)。数据使用bowtie2与小鼠基因组(mm9)对齐，允许多映射读取。如果从相对于基因组位置的重复reads中分离出来，则删除重复。本例中的对齐数据是作为床文件提供的，但是DEScan接受bam格式的对齐。

bam。files <- list.files(system.file(file.path("extdata"，"bam")， package="DEScan2")， pattern="bam$"， full.names=TRUE)

调用每个样本的峰值首先，对样本的峰值调用将使用findPeaks函数。在这种情况下，Bam文件被用作输入，但bed文件也可以(在任何情况下，我们建议指定一个基因组代码作为数据的参考)。为了节省时间，我们将在单个样本上执行这一步;但是，请注意，可以提供文件名的向量。如果保存flag设置为TRUE时，该函数将为文件夹中的每个输入文件生成一个排序的床文件outputFolder(“peaks”目录是默认目录)。findPeaks实现自适应窗口大小扫描来查找峰值，并需要2个参数来定义要测试的重叠窗口的大小。每个窗口的浓缩计算相对于minCompWinWidth5 kb,maxCompWinWidth10kb(局部)或整个染色体，三者中最大的被报道。

参数描述:

文件

包含待分析文件路径的字符向量。

文件类型

字符，“bam”或“bed”表示输入文件的格式。

binSize

扫描最小窗口的整数大小(以碱基对为单位)，默认为50。

minWin

整数，表示以碱基对为单位的最小窗口大小，默认为50。(binSize的倍数)

maxWin

整数，表示允许的最大窗口大小，单位为50 bp，默认为1000。(binSize的倍数)

minCompWinWidth

Z-score比较背景窗口的最小基线宽度(默认为5000)。

maxCompWinWidth

Z-score比较背景窗口的最大宽度(默认值是10000)。

zthresh

z分数的截止值，默认为10。只有z分数更高的窗口才会被保留。

minCount

一个小常数(通常不大于1)要在日志转换之前添加到计数中，以避免log(0)的问题。

sigwin

一个整数值，用于计算峰值信号的长度(默认值为10)。

保存

标志，表示是否需要保存输出

outputFolder

字符串，保存山峰的文件夹名称(可选)，如果该目录不存在，将创建该目录。默认为“峰值”。

力

一个布尔标志，指示是否强制输出覆盖。

genomeName

用来作为输入文件参考的基因组代码。我们强烈建议(特别是床上文件)指定基因组代码，但是，如果没有互联网连接或与服务器的联系存在问题，则不能使用此参数。

onlyStdChrs

只与标准染色体工作的标志。函数参数)

空空的如果不是NULL，像“chr#”这样的字符表示要使用的染色体

peaksGRL <- DEScan2::findPeaks(files=bam。文件[1], filetype="bam", genomeName="mm9", binSize=50, minWin=50, maxWin=1000, zthresh=10, minCount=0.1, sigwin=10, minCompWinWidth=5000, maxCompWinWidth=10000, save=FALSE, outputFolder="peaks", force=TRUE, onlyStdChrs=TRUE, chr=NULL, verbose=FALSE)

##调整复制之间的峰值以生成区域findPeaks已经在每个样本上进行了测试finalRegions函数可以用来对齐在多个样本中发现的重叠峰。所有对齐文件的峰值文件都可以在“extdata/peaks”文件夹中找到。finalRegions将产生一个GRanges对象，其中包含对齐的峰值的位置。如果saveFlag设置为TRUE时，该函数将在outputFolder．

参数描述:

peakSamplesGRangesList

命名为GRangesList，其中每个元素都是一个被称为峰值的样本。每个grange都需要一个z分数mcols值。(通常由findPeaks函数返回)。

zThreshold

考虑对齐峰值的z分数，它需要等于或高于findPeaks中使用的ztresh，默认值为20。

minCarriers

对准峰值后，需要重叠的最小生物重复数。

saveFlag

布尔值，指示是否将结果保存在以制表符分隔的文件中

outputFolder

存储输出文件的目录名。

峰值。文件<- system。file(file.path("extdata"，"peaks"，"RData"，"peaksGRL_all_files.rds")， package="DEScan2") peaksGRL <- readRDS(peaks.file) regionsGR <- DEScan2::finalRegions(peakSamplesGRangesList=peaksGRL, zThreshold=10, minCarriers=3, saveFlag=FALSE, outputFolder=NULL, verbose=FALSE) head(regionsGR)

## seqnames ranges strand | z-score n-peaks ##    |   ## chr19。* | 16.3522 3 ## chr19;6 ## chr19;P023_np06_k3 chr19 57226800-57228199 * | 18.3281 6 ## chr19。* | 25.2732 6 ## chr19;P028_np04_k4 chr19 57235750-57236399 * | 13.3579 4 ## chr19。4 ## k-carrier ## <数字> ## chr19;3 ## r;# 8226;3 ## r;# 8226;4 ## ch19;P032_np04_k3 3 ## ------- ## seqinfo: 1个来自mm9基因组的序列

该函数的输出是一个类似床的文件，其中的列表示基因组坐标，以及其他列:AvgZ，组合形成一个公共区域的峰值的平均z分数，NumCarriers，一个区域存在的样本数量。

结果区域可用于生成一个计数矩阵countFinalRegions函数。该函数获取要计数的区域(可以是任何类似于床的数据结构)，以及包含要计数的读取的文件的路径。所有对齐文件的Bam文件都可以在“extdata/ Bam /”文件夹中找到。为了加快进程，还可以指定最小载波数。在这种情况下，我们将不指定最小载波数，并在计数后进行过滤。的包装器summarizeOverlaps．

bam。path <- system.file(file.path("extdata"，"bam")， package="DEScan2") finalRegions <- DEScan2::countFinalRegions(regionsGRanges=regionsGR, readsFilePath=bam。path, fileType="bam"， minCarriers=1, genomeName="mm9"， onlyStdChrs=TRUE, saveFlag=FALSE, verbose=FALSE) counts <- summarizeexperimental::assay(finalRegions) regions <- summarizeexperimental::rowRanges(finalRegions)

它返回一个SummarizedExperiment对象。为了得到计数矩阵，有必要使用分析函数。得到的计数矩阵为每个区域包含一行，为每个样本包含一列。这种结构类似于常见的RNA-seq数据，可以用类似的工具进行标准化和分析。使用函数rowRanges地区的画眉会被归还。在区域名称上使用计数矩阵的行名，可以访问峰值的坐标。此外，区域GRanges对象可用于过滤出计数矩阵的行。实际上，为了可读性，我们将重命名并重新排序列，并过滤至少4个载体，以便仅测试相关区域。

counts <- counts[regions$ ' k-carriers ' >= 4，] counts <- counts[rowsum (counts) > 0，] colnames(counts) <- c("FC1"， "FC4"， "HC1"， "HC4"， "FC6"， "FC9"， "HC6"， "HC9") counts <- counts[，order(colnames(counts)))] head(counts)

## FC1 FC4 FC6 FC9 HC1 HC4 HC6 HC9 ## chr19。P021_np06_k6 129 9 50 26 51 46 27 20 ## chr19。15 4 21 341 53 56 11 6 ## chr19。P028_np04_k4 28 2 33 11 16 25 27 12 ## chr19。P036_np04_k4 25 5 18 9 5 43 10 2 ## chr19。36 20 42 18 28 48 9 12 ## chr19。P048_np04_k4 21 2 4 3 54 36 28

##使用RUV的规范化

为了控制“不需要的变化”，例如，批处理、库准备和其他讨厌的影响，可以使用RUVSeq包中的样本间归一化方法ruv。任何基于总文库计数的归一化方法都不适用于表观遗传测序实验，因为计数矩阵中总计数的差异可能是由于兴趣的处理造成的。

library("RColorBrewer") colors <- brewer。pal(3， "Set2") set <- EDASeq::betweenLaneNormalization(counts, which = "upper") groups <- matrix(c(1:8)， nrow=2, byrow=TRUE) trt <- factor(c(rep("FC"， 4)， rep("HC"， 4))))

相对日志表达式的箱形图(RLE =样本中读计数与中位数读计数的对数比值)和主成分图(PC)显示了样本间归一化的明显需要。

EDASeq::plotRLE(set, outline=FALSE, ylim=c(- 4,4)， col=colors[trt]， main="No Normalization RLE") (set, col=colors[trt]， main="No Normalization PCA"， labels=FALSE, pch=19)

的参数k指定需要移除的不需要变量的数量，在本例中我们使用4，但这由用户决定。我们可以在PCA图中看到，在RUVs归一化后，前2个主成分将两组分开，这表明治疗是变异的主要来源。

k <- 4 s <- RUVSeq::RUVs(set, cIdx=rownames(set)， scIdx=groups, k=k) EDASeq::plotRLE(s$normalizedCounts, outline=FALSE, ylim=c(- 4,4)， col=colors[trt]， main="Normalized RLE") EDASeq::plotPCA(s$normalizedCounts, col=colors[trt]， main="Normalized PCA"， labels=FALSE, pch=19)

##区域差异富集测试

现在，我们准备使用在edgeR中实现的负二项式拟似然GLM方法来寻找差异富集区域(详细信息请参阅edgeR包小图)。这是通过考虑一个设计矩阵来实现的，该设计矩阵既包括感兴趣的协变量(这里是治疗状态)，也包括不需要的变化因素。最后通过对之前计算的区域进行子集化，得到了差异富集区域的坐标。

设计<-模型。matrix(~0 + trt + s$W) colnames(design) <- c(levels(trt)， paste0("W"， 1:k)) y <- edgeR::DGEList(counts=counts, group=trt) y <- edgeR::estimateDisp(y, design) fit <- edgeR::glmQLFit(y, design, robust=TRUE) con <- limma::makeContrasts(FC - HC, levels=design) qlf <- edgeR::glmQLFTest(fit，对比度=con) res <- edgeR::topTags(qlf, n=Inf, p.value=0.05) head(res$table)

## logFC logCPM F PValue FDR ## chr19。P362_np04_k4 2.587 11.83 17.49 0.00003929 0.005186 ## chr19。P506_np07_k7 2.005 12.59 11.13 0.00097202 0.047489 ## chr19。P254_np07_k5 2.067 13.05 10.93 0.00107929 0.047489

暗(res美元表)

## [1] 3 5

地区(rownames (res表美元))

## seqnames ranges strand | z-score n-peaks ##    |   ## chr19。P362_np04_k4 chr19 57993950-57994699 * | 18.6357 4 ## chr19。* | 19.7085 7 ## chr19;7 ## k-carrier ## <数字> ## chr19;4 ## c;7 # c;P254_np07_k5 5 ## ------- ## seqinfo: 1个来自mm9基因组的序列

3.关于性能和串行计算的说明

DEScan2包使用BiocParallel包来支持并行计算。在这里，我们使用register命令来确保小插图使用串行计算运行。

然而，在真实的数据集中，并行计算可以大大加快计算速度，在许多基因和/或许多细胞的存在。

在DEScan2中允许并行计算有两种方法，第一种是register()一个并行后端(参见?BiocParallel::register了解详细信息)。或者，可以将BPPARAM对象传递给findPeaks和finalRegions函数，例如:

library("BiocParallel") peaksGRL <- DEScan2::findPeaks(files=bam。文件[1], filetype="bam", genomeName="mm9", binSize=50, minWin=50, maxWin=1000, zthresh=10, minCount=0.1, sigwin=10, minCompWinWidth=5000, maxCompWinWidth=10000, save=FALSE, outputFolder="peaks", force=TRUE, onlyStdChrs=TRUE, chr=NULL, verbose=FALSE, BPPARAM=BiocParallel::MulticoreParam(2))

我们发现MulticoreParam()在Mac上可能有一些性能问题;因此，我们建议在Mac上工作时使用DoparParam()。