许可

在艺术许可2.0下，您可以自由使用和重新发布此软件。但是，如果您在自己的出版物中使用此软件包，我们会要求您引用以下论文:

\ start {itemize} \item[]张晓霞，张志强，张志强，张志强。(2009)。PICS: Chip-seq的概率推断。生物识别67,151 -163。结束\{逐条列记}

简介

ChIP-Seq将染色质免疫沉淀与大规模并行短读测序相结合，能够以比ChIP-chip更高的灵敏度、特异性和空间分辨率分析体内全基因组转录因子- dna关联。虽然ChIP-Seq为研究提供了新的机会，但它也为统计分析提出了新的挑战，这些挑战来自生物系统的复杂性及其数字序列数据中的可变性和偏差。我们提出了一种称为PICS (ChIP-Seq概率推断)的方法，用于从ChIP-Seq对齐读取数据中提取信息，以识别由转录因子结合的区域。PICS通过模拟定向读的局部浓度来识别富集区域，并使用DNA片段长度先验信息通过贝叶斯分层t混合模型来区分紧密相邻的结合事件。它的每个事件片段长度估计还允许它从具有非典型片段长度的分析区域中移除。PICS使用预先计算的全基因组读拟合性剖面和截断的t分布来调整绑定事件模型，以补偿由于局部基因组重复而丢失的读。PICS估计模型参数中的不确定性，这些不确定性可用于定义绑定事件位置上的置信区域并过滤估计。

一个典型的PICS分析包括以下步骤:\begin{enumerate} \item将数据转换为农庄用于高效处理\项目基因组分割通过segmentPICS估计绑定位点位置和其他PICS参数通过图片通过交换IP和控制样本重复1-2，使用1-3估计FDR viapicsFDR\item输出丰富的地区和图片带床和假发文件的配置文件\end{枚举}

数据输入和格式化

与任何包一样，首先需要使用以下命令加载它

库(图片)

PICS管道的第一步包括将对齐的读取(来自IP和对照样本)转换为一种格式，允许将基因组有效分割为一组具有足够正向和反向读取的候选区域。数据格式可以是床类型的数据帧以及AlignedReads函数返回的对象readAligned'可以读取各种文件格式，包括Eland, MAQ, Bowtie等。有关如何将数据读入R(而不是从BED文件)的更多细节，请参阅ShortRead的小插图。在下面列出的示例中，我们使用床类型文件，它由制表符分隔的文件组成，包含以下列:空格、开始、结束、链，其中每一行表示一次读取。

除了IP和Control数据集之外，我们还使用了另一个参数，该参数由所询问的基因组的可映射性配置文件组成。对于每条染色体，特定读长度的可映射性配置文件由0和1组成的向量组成，它给出了染色体中每个碱基对的估计读可映射性得分。在某个位置的分数为1意味着我们应该能够在该位置唯一对齐该长度的读取，而分数为0表示该长度的读取不应该在该位置唯一对齐。如上所述，不能映射到唯一基因组位置的reads通常被丢弃。为了方便起见，并且由于可映射区域和不可映射区域之间的转换通常比这些区域短得多，我们将每个染色体的可映射配置文件紧凑地总结为仅指定零值配置文件的不可映射间隔的不相交并集。我们将此信息存储为床文件，其中每一行对应一个不可映射的间隔。PICS使用这种可映射性概要来估计由于无法对齐而丢失的读取。一系列读取长度和种类的可映射性配置文件可从:

http: / wiki.rglab.org/index.php ? title =公共:Mappability_Profile

要创建自己的可映射性配置文件，请使用以下脚本:

https://github.com/SRenan/proMap/

一旦数据被读取，我们就可以创建GRanges对象，它存储基因组位置序列和相关注释。列表中的每个元素都包含IP和控制读以及对应染色体的可映射性配置文件。在读取数据时，我们还可以对重复的读取进行排序和删除，建议在PICS分析中这样做(排序实际上是必需的)。

本文档中数据的存放路径为:…/PICS/inst/extdata文件夹。

路径<- system. Path。file("extdata"， package = "PICS") dataIP <- read.table(file. data, package = "PICS")path(path， "Treatment_tags_chr21_sort.bed")， header=TRUE, colClasses=c("factor"， "integer"， "integer"， "factor")) dataIP <- as(dataIP， "GRanges") dataCont <- read.table(文件。path(path， "Input_tags_chr21_sort.bed")， header = TRUE, colClasses = c("factor"， "integer"， "integer"， "factor")) dataCont <- as(dataCont， "GRanges")

Map <- read.table(file. table)path(path， "mapProfileShort")， header=TRUE, colClasses=c("factor"， "integer"， "integer"， "NULL")) map <- as(map， "GRanges")

图片分析

基因组分割

我们通过预处理来自单个ChIP-Seq实验的双向对齐读取数据，将基因组分割为候选区域，以检测具有最少数量的正向和反向读取的候选区域。然后，PICS将处理这些候选区域。的minReads参数将严重影响返回的候选区域数量。原则上，用小的值为好minReads为了不错过任何真正的结合位点。但是，过小的值可能会返回许多假阳性区域，从而导致\texttt{PICS}函数的计算时间变长。如果有一个控制样本，你也可以选择一个值minReads这将给你一个控制中候选区域数量与IP中候选区域数量的比例，约为25% \%。当使用\texttt{picsFDR}时，大于25 %的比率很可能会导致较大的FDR，因此没有必要使用过低的阈值。请看下面的说明。\换行符

seg <- segmentPICS(dataIP, dataC = dataCont, map = map, minReads = 1)

参数估计

在将基因组分割成候选区域后，我们使用\texttt{PICS}函数以概率的方式检测结合区域，并返回结合位点估计值、置信区间等。在下面的代码中，我们假设您正在处理转录因子数据，如所指定的数据类型选项设置为特遣部队转录因子。我们还计划在未来的PICS版本中支持组蛋白修饰数据，该版本将通过指定打开数据类型=“H”．

为了提高PICS包的计算效率，我们建议使用\texttt{parallel}包，它允许简单的并行计算。在接下来的内容中，我们假设您的计算机上安装了\texttt{parallel}，并且您至少有两个内核。如果不是，可以省略第一行和参数nCores，计算将在单个CPU上进行。默认情况下，函数将只使用一个核心。

库(并行)

在我们的例子中，假设我们已经使用\texttt{segmentPICS}分割了基因组，我们可以继续执行以下命令:

pics <- pics (seg, nCores = 2)

罗斯福估计

为了估计FDR，我们需要在交换IP和控制样本后重新运行我们的分析。注意，这需要一个对照样本的存在。在去除带有滤波器指定的噪声参数的绑定站点估计后，我们继续执行以下命令来计算FDR。

segC <- segmentPICS(dataCont, dataC = dataIP, map = map, minReads = 1) picsC <- PICS(segC) fdr <- picsFDR(PICS, picsC, filter = list(delta = c(50, Inf)， se = c(0,50)， sigmaSqF = c(0,22500)， sigmaSqR = c(0,22500)))

然后，一旦估算出FDR，就可以将FDR视为分数的函数，以选择一个阈值，从而获得所需的FDR。

plot(pics, picsC, xlim = c(2,8)， ylim = c(0， .2)， filter = list(delta = c(50,300)， se = c(0,50)， sigmaSqF = c(0,22500)， sigmaSqR = c(0,22500))， type = "l")

块plotFDR的图

您还可以将FDR可视化为区域数量的函数。

图(fdr[， c(3,1)])

块plotFDR2的图

结果输出

为了便于数据处理和数据输出，我们使用\texttt{IRanges}包将我们的结果总结为农庄对象。由此产生的农庄对象可以使用\texttt{IRanges}包进行分析和/或使用\texttt{rtracklayer}包导出到bed/wig文件。在每一种情况下，我们可以使用滤波器噪声区域和或得分过低的区域。特别是\texttt{picsFDR}函数可以用来设置得分阈值，从而达到所需的分数罗斯福．

将丰富的区域写入BED文件

myFilter = list(delta = c(50,300)， se = c(0,50)， sigmaSqF = c(0, 22500)， sigmaSqR = c(0, 22500)) rdBed <- makeGRangesOutput(pics, type="床"，filter = c(myFilter, list(score = c(1, Inf))))

下面的命令将创建适当的床文件，因为它需要\texttt{rtracklayer}包，所以我们在这里不运行它，只是为了您的目的而显示它。

library(rtracklayer)导出(rdBed， "myfile.bed")

将密度分数写入WIG文件

下面的命令将创建适当的农庄和相应的wig文件，因为它需要\texttt{rtracklayer}包，我们在这里不运行它，只是为了您的目的而显示它。

rdWig <- makeGRangesOutput(pics, type="wig"， filter=c(myFilter, list(score=c(1,Inf)))) export(rdWig， "myfile.wig"))