1．1开始

该示例数据集局限于第一个酵母染色体的前100kb。

1. 我们设置了一些必需的元变量

指定文件夹和实验设计路径wd <- system. library(GenoGAM)file("extdata/Set1"， package='GenoGAM')path(wd， "bam") expDesign <- file. path(wd， "bam")path(wd， "experimentDesign.txt") ## B. #寄存器并行后端(默认为"核数" - 2)BiocParallel::register(BiocParallel::MulticoreParam(worker = 2)) ## C. ##指定块和挂起大小。chunkSize <- 50000 overhangSize <- 1000 d# #实验设计，反映底层gam# # x =位置设计<- ~ s(x) + s(x, by =基因型)

2. 读入数据得到aGenoGAMDataSet对象。

##构建GenoGAMDataSet ggd <- GenoGAMDataSet(expDesign, directory = folder, chunkSize = chunkSize, overhangSize = overhangSize, design = design)

## INFO[2019-10-29 21:55:47]正在创建GenoGAMDataSet ## INFO[2019-10-29 21:55:48]读取数据## INFO[2019-10-29 21:55:49]读取wt_1 ## INFO[2019-10-29 21:55:49]读取wt_2 ## INFO[2019-10-29 21:55:50]读取mutant_1 ## INFO[2019-10-29 21:55:50]读取mutant_2 ## INFO[2019-10-29 21:55:51]已完成读取数据## INFO [2019-10-29 21:55:51] GenoGAMDataSet已创建

ggd

##类:GenoGAMDataSet ##维数:100000 4 ##样本(4):wt_1 wt_2 mutant_1 mutant_2 ##设计变量(1):基因型##瓦尺寸:52kbp ##瓦数量:2 ##染色体:chrI ##尺寸因素:## wt_1 wt_2 mutant_1 mutant_2 ## 0000 #公式:## ~ s(x) + s(x, by =基因型)

3. 计算尺寸因子

##计算大小因子ggd <- computeSizeFactors(ggd)

##信息[2019-10-29 21:55:52]计算大小因子完成

ggd

##类:GenoGAMDataSet ##维数:100000 4 ##样本(4):wt_1 wt_2 mutant_1 mutant_2 ##设计变量(1):基因型##瓦尺寸:52kbp ##瓦数量:2 ##染色体:chrI ##尺寸因素:## wt_1 wt_2 mutant_1 mutant_2 ## -0.04156987 0.04210819 -0.25589212 0.14328182 ##公式:## ~ s(x) + s(x, by =基因型)

##数据分析(ggd)

##数据帧100000行4列## wt_1 wt_2 mutant_1 mutant_2 ##     ## 10000 ## 2 00000 ## 3 00000 ## 4 00000 ## 5 00000 ## ... ... ... ... ...## 99996 00000 0 ## 99998 00000 ## 999999 0000 # 100000 1 1000

1．2模型拟合

1. 用固定超参数计算模型

##计算没有参数估计的模型，以节省时间在vignette结果<- genogam(ggd, lambda = 4601, theta = 4.51)

## INFO[2019-10-29 21:55:52]初始化模型## INFO[2019-10-29 21:55:52]已完成## INFO[2019-10-29 21:55:52]拟合模型## INFO[2019-10-29 21:56:00]已完成## INFO[2019-10-29 21:56:01]处理已完成## INFO[2019-10-29 21:56:01]已完成

结果

##家族:nb ##公式:~ s(x) + s(x, by =基因型)##类别:GenoGAM ##尺寸:100000 4 ##样本(4):wt_1 wt_2 mutant_1 mutant_2 ##设计变量(1):基因型##平滑函数(2):s(x) s(x):基因型##染色体:chrI ## ##大小因素:## wt_1 wt_2 mutant_1 mutant_2 ## -0.04156987 0.04210819 -0.25589212 0.14328182 ## ##交叉验证:未执行## ##样条参数:##结间距:20 ## b样条顺序:2 ##惩罚顺序:2 ## ## Tile设置:##块大小:50000 ## Tile大小:52000 ##悬挂大小:1000 ## tiles数量:2 ## #评估的基因组范围:## GRanges对象1范围和0元数据列:## seqnames ranges strand ##    ## [1] chrI 1-100000 * ## ------- ## seqinfo: 1个来自未指定基因组的序列

##匹配和标准错误匹配(结果)

##数据帧100000行2列## s(x) s(x):基因型## <数字> <数字> ## 1 -2.81190225503088 0.5320456038111 ## 2 -2.81128719366336 0.531848806972624 ## 3 -2.81067225244489 0.531651893033589 ## 4 -2.81005744734502 0.531454845516499 ## 5 - 2.809442774647943864 ## ... ... ...## 99996 -1.86772191646085 -0.110004776821936 ## 99997 -1.875276615882 -0.110350587546044 ## 99998 -1.88283132577617 -0.110696460554156 ## 99999 -1.89038604404279 -0.111042389677715 ## 100000 -1.89794076858129 -0.111388368748164

se(结果)

##数据帧100000行2列## s(x) s(x):基因型## <数字> <数字> ## 1 0.258647538432571 0.317642013209663 ## 2 0.257521587708836 0.316405090894515 ## 3 0.256401263130459 0.31517351476289 ## 4 0.255286594725364 0.31394731225843 ## 5 0.127609315401 0.312726508234327 ## ... ... ...0.209381149691971 ## 99997 0.157924553240842 0.210544652601457 ## 99998 0.158948126782311 0.211713323560403 ## 99999 0.159977531447117 0.212887133067987 ## 100000 0.161012724757534 0.214066047892396

2. 计算日志p值

result <- computessignificance (result)

## INFO[2019-10-29 21:56:01]计算位置p值## INFO[2019-10-29 21:56:01]完成

pvalue(结果)

##数据帧100000行2列## s(x) s(x):基因型## <数字> <数字> ## 1 0.0939371770243043 ## 2 0.095912533137095e -28 0.0927801973572222 ## 3 0.091647376361327e -28 0.0916313734527 ## 4 0.0904906101030138 ## 5 2.1184519790328816893581659938518 ## ... ... ...1.6057301625264e-32 0.600195329046113 ## 99998 0.601071576336622 ## 99999 3.2049678853473e-32 0.601947355046144 ## 100000 4.52673255409177e-32 0.602822579980591

中心10kb的绘图结果，其中两个瓷砖被连接在一起。

range = GRanges("chrI"， IRanges(45000, 55000)) plotGenoGAM(result, ranges = ranges)

3．1分析的目标

提供了一个小数据集来说明ChIP-Seq功能。这是Thornton等人发表的数据的一个子集(Thornton et al. 2014)他通过ChIP-Seq分析了组蛋白H3赖氨酸4三甲基化(H3K4me3)，比较了野生型酵母和具有截断形式的Set1(酵母H3赖氨酸4甲基化酶)的突变体。该分析的目标是鉴定突变株与野生型株相比甲基化显著差异的基因组位置。

为此，我们将建立一个GenoGAM对预期ChIP-seq片段计数的对数建模的模型\ (y \)作为基因组位置的平滑函数的和\ \ (x)．具体来说，我们写(用简化的符号):

其中基因型为1的数据来自突变型，和0的野生型。这里是函数\ \ (f (x))为参考水平，即野生型菌株的对数速率。这个函数f \ \(文本{突变/ wt} (x) \)是突变型与野生型的对数比。然后统计检验零假设\(f_\text{突变/wt}(x) = 0\)在每个位置\ \ (x)．在接下来的内容中，我们将展示如何构建数据集，执行模型拟合并执行测试。

3.2注册并行后端

方法注册并行后端BiocParallel包中。中的文档BiocParallel类，以便正确使用。还要注意，BiocParallel只是多个并行包的接口。例如，为了在集群上使用GenoGAM，BatchJobs软件包可能需要。并行后端可以随时注册，因为GenoGAM只会调用当前后端。

重要的是:根据这而且这发布在Bioconductor支持页面和R-devel邮件列表上的帖子，这是最重要的核心功能多核后端，共享内存，被Rs自己的垃圾收集器破坏，导致整个工作空间的复制跨所有核。考虑到内存中的数据量很大(没有HDF5后端)，它可能会使整个系统崩溃。我们强烈建议注册SnowParam后端为了避免这种情况如果研究的是较大的生物体．这种方式的开销略大，但只复制必要的数据到worker，保持内存消耗相对较低。我们从未经历过每个核心超过4GB的负载，通常人类基因组的负载在2GB左右，HDF5后端甚至更低。

在多核机器上，默认情况下可用核的数量- 2在默认后端' multicore ' BiocParallel::registered()[1]

## $MulticoreParam ##类:MulticoreParam ## bpisup: FALSE;bpnworkers: 2;bptasks: 0;BPJOB: BPJOB ## bplog: FALSE;bpthreshold:信息;bpstopOnError: TRUE ## bpRNGseed:bptimeout: 2592000;bpprogressbar: FALSE ## bpexportglobals: TRUE ## bplogdir: NA ## bpresultdir: NA ##集群类型:FORK

对于这个小例子，我们希望分配更少的工人。检查BiocParallel对于其他可能的后端和更多的选项SnowParam

BiocParallel::register(BiocParallel::SnowParam(workers = 3))

如果检查当前注册的后端，可以看到工人的数量和后端已经改变。

BiocParallel:注册()[1]

## $SnowParam ##类:SnowParam ## bpisup: FALSE;bpnworkers: 3;bptasks: 0;BPJOB: BPJOB ## bplog: FALSE;bpthreshold:信息;bpstopOnError: TRUE ## bpRNGseed:bptimeout: 2592000;bpprogressbar: FALSE ## bpexportglobals: TRUE ## bplogdir: NA ## bpresultdir: NA ##集群类型:SOCK

3．3构建GenoGAMDataSet

限制酵母染色体I的前100kb的BAM文件在本月/ extdata的文件夹GenoGAM包中。该文件夹还包含一个描述实验设计的平面文件。

我们首先从制表符分隔的文本文件加载实验设计experimentDesign.txt输入一个数据帧:

文件夹<- system。文件("extdata/Set1", package='GenoGAM') expDesign <- read.delim( file.path(folder, "experimentDesign.txt") ) expDesign

ID文件配对基因型## 1 wt_1 H3K4ME3_Full_length_Set1_Rep_1_chr1_100k。bam FALSE 0 ## 2 wt_2 H3K4ME3_Full_length_Set1_Rep_2_chr1_100k。bam FALSE 0 ## 3突变_1 H3K4ME3_aa762-1080_Set1_Rep_1_chr1_100k。bam FALSE 1 ## 4 mutant_2 H3K4ME3_aa762-1080_Set1_Rep_2_chr1_100k。bam FALSE 1

实验设计每行对应一个ChIP样本BAM格式的比对文件。对于多路复用排序，BAM文件必须已被解复用。实验设计有命名柱。三个列名具有固定的含义GenoGAMDataSet并且必须提供:ID，文件，配对．这个领域ID为每个对齐文件存储唯一标识符。建议使用简短和易于理解的标识符，因为它们随后用于标记数据和图表。这个领域文件存储没有目录的BAM文件名。这个领域配对布尔值真正的对于配对端测序数据，和假用于单端测序数据。所有其他柱均作为实验设计，且必须为二进制。这些列的名称必须与后面设计公式中提供的名称相同。它将允许我们以后建模不同的占用率或呼叫峰值。这里重要的是基因型列。

我们现在将计算每个基因组位置的测序片段中心，并将这些计数存储到GenoGAMDataSet类。GenoGAM将ChIP-Seq数据减少到片段中心计数，而不是全基覆盖，因此每个片段只计数一次。这减少了相邻核苷酸之间的人为相关性。对于单端库，通过将读端位置移动一个常数来估计片段中心。这个常数的估计利用了chipseq包(参见chipseq: estimate.mean.fraglen方法)。

创建对象所需的参数GenoGAMDataSet是:

• 的实验设计作为配置文件或data.frame的路径。
• 的目录的BAM文件
• 的设计公式

##构建GenoGAMDataSet wd_folder <- file。path(folder， "bam") ggd <- GenoGAMDataSet(expDesign, directory = wd_folder, design = ~ s(x) + s(x, by =基因型))

## INFO[2019-10-29 21:56:03]正在创建GenoGAMDataSet ## INFO[2019-10-29 21:56:04]读取数据## INFO[2019-10-29 21:56:04]读取wt_1 ## INFO[2019-10-29 21:56:05]读取wt_2 ## INFO[2019-10-29 21:56:05]读取mutant_1 ## INFO[2019-10-29 21:56:06]读取数据## INFO [2019-10-29 21:56:06] GenoGAMDataSet已创建

分析(ggd)

##数据帧100000行4列## wt_1 wt_2 mutant_1 mutant_2 ##     ## 10000 ## 2 00000 ## 3 00000 ## 4 00000 ## 5 00000 ## ... ... ... ... ...## 99996 00000 0 ## 99998 00000 ## 999999 0000 # 100000 1 1000

getOverhangSize(ggd)

## [1] 1000

getChunkSize(ggd)

## [1] 98000

的GenoGAMDataSet类将此计数数据存储到索引基因组位置的结构中瓷砖，定义为chunkSize而且overhangSize．了解这些参数需要一些背景知识。GenoGAM中的平滑是基于样条(图1)，这是分段多项式。的结是多项式连接的位置。根据我们的经验，在典型应用中，每20到50 bp就需要一个结，以获得足够的平滑配合分辨率。广义加性模型的拟合涉及的步骤要求与结数的平方成比例的操作，防止沿整个染色体拟合。为了对GAMs进行全基因组拟合，GenoGAM在重叠区间进行拟合(瓷砖)，并在连续贴图重叠的中点连接拟合。的参数chunkSize而且overhangSize定义瓦片，其中块是不重叠其他瓦片的瓦片的核心部分，悬垂是两个重叠部分。

这两个变量都影响了计算速度和内存使用，因为它们定义了整个基因组被分割用于并行计算的块的数量(以及大小)。此外，过小的悬垂尺寸可能会导致接合点的跳汰配合，其中块被粘在一起。因此建议将悬垂大小保持在1kb左右(默认情况下是这样设置的)。块大小在默认情况下大致估计，以适应整体建模问题，而不使用太多内存。然而，这种估计是一种非常简单的启发式方法，因此如果遇到问题，建议自己设置块大小。我们对50 - 100kb之间的数字有很好的经验。

图1:样条示例。显示的是7个三次b样条基函数，它们共同构成了完整的样条(上面的深蓝色)。结点在每个基函数的底部中心被描绘成深灰色的圆点。所有基函数都乘以了它们各自的系数，因此不同。

3.4GenoGAMSettings和HDF5后端提供更多选项

对于酵母基因组以外的大多数基因组大小，将所有数据保存在内存中是不切实际的，甚至是不可能的。因此，使用HDF5后端将所有相关数据存储在硬盘上。这可以通过设置参数简单地完成hdf5为true。这还将触发将底层数据按染色体拆分为数据帧列表。理论上，拆分数据是不必要的。然而，由于在Bioconductor基础设施中表示长整数的问题，大于2^31(约2.14千兆字节对)的数据帧不能轻松存储，或者只能通过需要更多存储空间的替代方法存储。此外，由于HDF5数据存储在硬盘驱动器上，它可以很容易地移动到不同的设备，但必须再次读入。为了简化，参数fromHDF5是否可以设置为true，以触发GenoGAMDataSet读取已经存在的HDF5文件。

ggd <- GenoGAMDataSet(expDesign, directory = wd_folder, design = ~ s(x) + s(x, by =基因型)，hdf5 = TRUE)

## INFO[2019-10-29 21:56:06]正在创建GenoGAMDataSet ## INFO[2019-10-29 21:56:07]读取wt_1 ## INFO[2019-10-29 21:56:07]读取wt_2 ## INFO[2019-10-29 21:56:07]读取mutant_1 ## INFO[2019-10-29 21:56:08]读取mutant_2 ## INFO[2019-10-29 21:56:08]已完成读取data ## INFO[2019-10-29 21:56:09]将chrI写入HDF5 ## INFO [2019-10-29 21:56:09] chrI已写入## INFO [2019-10-29 21:56:10] GenoGAMDataSet已创建

##注意与上述试验的DataFrame相比数据类型的差异(ggd)

## $chrI ## # <100000 x 4>矩阵类DelayedMatrix和类型“integer”:## wt_1 wt_2 mutant_1 mutant_2 ## [1，] 00000 ## [2，] 00000 ## [3，] 00000 ## [4，] 0000 ## [5，] 00000 00 ## # ... ... .## [99996，] 0000 ## [99997，] 0000 ## [99998，] 0000 ## [99999，] 0000 ## [100000，] 1 100000

默认情况下，HDF5文件存储在临时文件夹中，如果必须保存和共享数据，这是不切实际的。因此HDF5文件夹可以在GenoGAMSettings对象。该对象包含许多元参数，这些元参数负责的功能GenoGAM．它特别保护以下领域:

• 读入和处理数据
• 模型拟合和交叉验证过程中的参数估计
• 数据和存储控制

所有设置都可以通过打印GenoGAMSettings对象。除了元参数之外，它还将显示并行后端信息和记录器阈值，尽管它们都不是GenoGAMSettings对象。

GenoGAMSettings ()

## -------------------- 记录参数  ----------------- ## 中心:真正的染色体# #:# #定制过程功能:禁用  ## ## -------------------- BAM参数  --------------------- ## 类:ScanBamParam # # bamFlag (NA除非指定):# # bamSimpleCigar:假# # bamReverseComplement:假# # bamTag: # # bamTagFilter: # # bamWhich: 0 # # bamWhat范围:pos, qwidth # # bamMapqFilter: NA  ## ## -------------------- 并行端  ------------------- ## 类:SnowParam # # bpisup:虚假的;bpnworkers: 3;bptasks: 0;BPJOB: BPJOB ## bplog: FALSE;bpthreshold:信息;bpstopOnError: TRUE ## bpRNGseed:bptimeout: 2592000;bpprogressbar:假# # bpexportglobals:真正的# # bplogdir: NA # # bpresultdir: NA # #集群类型:袜子  ## ## -------------------- 日志设置  -------------------- ## 记录器阈值:信息  ## ## -------------------- HDF5设置  ---------------------- ## HDF5目录:/ tmp / RtmpJ73Hov / HDF5Array_dump # # HDF5压缩级别(0 - 9):6  ## ## -------------------- 数据设置  ---------------------- ## 区域大小:4000 # # # #完全结:20 ## ## --------------------优化参数------------ ##优化方法:Nelder-Mead ##优化控制:## maxit: 60 ## fnscale: -1 ## trace: 1 ##参数估计控制:## eps: 1e-06 ## maxiter: 1000

range <- GRanges("chrI"， IRanges(30000, 50000)) params <- Rsamtools::ScanBamParam(which = range) settings <- GenoGAMSettings(center = FALSE, bamParams = params) ggd <- GenoGAMDataSet(expDesign, directory = wd_folder, design = ~ s(x) + s(x, by = genotype)， settings = settings)

## INFO[2019-10-29 21:56:11]正在创建GenoGAMDataSet ## INFO[2019-10-29 21:56:11]读取数据## INFO[2019-10-29 21:56:11]读取wt_1 ## INFO[2019-10-29 21:56:12]读取wt_2 ## INFO[2019-10-29 21:56:12]读取mutant_1 ## INFO[2019-10-29 21:56:12]读取mutant_2 ## INFO[2019-10-29 21:56:12]已完成读取数据## INFO [2019-10-29 21:56:12] GenoGAMDataSet已创建

##注意高计数测定(ggd)

##数据帧与20001行和4列## wt_1 wt_2 mutant_1 mutant_2 ##     ## 1 20 6 3 5 ## 2 20 6 3 5 ## 4 20 6 3 5 ## 5 20 6 3 5 ## 5 20 6 3 5 ## # ... ... ... ... ...## 19997 14 5 6 7 ## 19998 14 5 6 7 ## 19999 14 5 6 7 ## 20000 13 5 6 7 ## 20001 13 5 6 7

rowRanges (ggd)

## seqnames pos strand ##    [1] chrI 30000 * ## [2] chrI 30001 * ## [3] chrI 30002 * ## [4] chrI 30003 * ## [5] chrI 30004 * ## ... ... ... ...## [19997] chrI 49996 * ## [19999] chrI 49998 * ## [20000] chrI 49999 * ## [20001] chrI 50000 * ## ------- ## seqinfo:来自未知基因组的1个序列

现在是染色体说明。由于在示例文件中只有chrI存在##，读入将记录一个错误并生成一个空的GenoGAMDataSet对象settings <- GenoGAMSettings(chromosomeList = c("chrII"， "chrIII")) ggd <- GenoGAMDataSet(expDesign, directory = wd_folder, design = ~ s(x) + s(x, by = genotype)， settings = settings)

## INFO[2019-10-29 21:56:12]正在创建GenoGAMDataSet ## ERROR[2019-10-29 21:56:12]没有染色体来读取。检查指定的设置或BAM文件头## INFO [2019-10-29 21:56:12] GenoGAMDataSet已创建

##空的GenoGAMDataSet，因为在例子BAM文件长度中缺少数据(ggd)

## [1] 0

ggd

##类:GenoGAMDataSet ##维度:0 0 ##样本(0):##设计变量(0):##尺寸因素:##数值(0)##公式:## ~ s(x)

注意，optimControl和estimControl是一个包含更多参数的列表。然而，除了所显示的以外，没有一个是重要的设置<- GenoGAMSettings(optimMethod = "L-BFGS-B"， optimControl = list(maxit = 100L)， estimControl = list(eps = 1e-10, maxiter = 500L))设置

## -------------------- 记录参数  ----------------- ## 中心:真正的染色体# #:# #定制过程功能:禁用  ## ## -------------------- BAM参数  --------------------- ## 类:ScanBamParam # # bamFlag (NA除非指定):# # bamSimpleCigar:假# # bamReverseComplement:假# # bamTag: # # bamTagFilter: # # bamWhich: 0 # # bamWhat范围:pos, qwidth # # bamMapqFilter: NA  ## ## -------------------- 并行端  ------------------- ## 类:SnowParam # # bpisup:虚假的;bpnworkers: 3;bptasks: 0;BPJOB: BPJOB ## bplog: FALSE;bpthreshold:信息;bpstopOnError: TRUE ## bpRNGseed:bptimeout: 2592000;bpprogressbar:假# # bpexportglobals:真正的# # bplogdir: NA # # bpresultdir: NA # #集群类型:袜子  ## ## -------------------- 日志设置  -------------------- ## 记录器阈值:信息  ## ## -------------------- HDF5设置  ---------------------- ## HDF5目录:/ tmp / RtmpJ73Hov / HDF5Array_dump # # HDF5压缩级别(0 - 9):6  ## ## -------------------- 数据设置  ---------------------- ## 区域大小:4000 # # # #完全结:20 ## ## --------------------优化参数------------ ##优化方法:L-BFGS-B ##优化控制:## maxit: 100 ## fnscale: -1 ## trace: 1 ##参数估计控制:## eps: 1e-10 ## maxiter: 500

##设置HDF5文件夹为“myFolder”tmp <- tempdir() settings <- GenoGAMSettings(hdf5Control = list(dir = file. conf))path(tmp， "myFolder")， level = 0))

"/tmp/RtmpJ73Hov/myFolder"

HDF5Array: getHDF5DumpCompressionLevel ()

## [1] 0

##另一种设置方法是通过HDF5Array包HDF5Array::setHDF5DumpDir(文件。path(tmp， "myOtherFolder"))设置

## -------------------- 记录参数  ----------------- ## 中心:真正的染色体# #:# #定制过程功能:禁用  ## ## -------------------- BAM参数  --------------------- ## 类:ScanBamParam # # bamFlag (NA除非指定):# # bamSimpleCigar:假# # bamReverseComplement:假# # bamTag: # # bamTagFilter: # # bamWhich: 0 # # bamWhat范围:pos, qwidth # # bamMapqFilter: NA  ## ## -------------------- 并行端  ------------------- ## 类:SnowParam # # bpisup:虚假的;bpnworkers: 3;bptasks: 0;BPJOB: BPJOB ## bplog: FALSE;bpthreshold:信息;bpstopOnError: TRUE ## bpRNGseed:bptimeout: 2592000;bpprogressbar:假# # bpexportglobals:真正的# # bplogdir: NA # # bpresultdir: NA # #集群类型:袜子  ## ## -------------------- 日志设置  -------------------- ## 记录器阈值:信息  ## ## -------------------- HDF5设置  ---------------------- ## HDF5目录:/ tmp / RtmpJ73Hov / myFolder # # HDF5压缩级别(0 - 9):0  ## ## -------------------- 数据设置  ---------------------- ## 区域大小:4000 # # # #完全结:20 ## ## --------------------优化参数------------ ##优化方法:Nelder-Mead ##优化控制:## maxit: 60 ## fnscale: -1 ## trace: 1 ##参数估计控制:## eps: 1e-06 ## maxiter: 1000

例如，我们选择了两倍长的区域，但也减少了两倍的结(两倍的结间距)，这导致每个贴图的结数与之前相同。- GenoGAMSettings(dat控件= list(regionSize = 8000, bpknots = 40))设置

## -------------------- 记录参数  ----------------- ## 中心:真正的染色体# #:# #定制过程功能:禁用  ## ## -------------------- BAM参数  --------------------- ## 类:ScanBamParam # # bamFlag (NA除非指定):# # bamSimpleCigar:假# # bamReverseComplement:假# # bamTag: # # bamTagFilter: # # bamWhich: 0 # # bamWhat范围:pos, qwidth # # bamMapqFilter: NA  ## ## -------------------- 并行端  ------------------- ## 类:SnowParam # # bpisup:虚假的;bpnworkers: 3;bptasks: 0;BPJOB: BPJOB ## bplog: FALSE;bpthreshold:信息;bpstopOnError: TRUE ## bpRNGseed:bptimeout: 2592000;bpprogressbar:假# # bpexportglobals:真正的# # bplogdir: NA # # bpresultdir: NA # #集群类型:袜子  ## ## -------------------- 日志设置  -------------------- ## 记录器阈值:信息  ## ## -------------------- HDF5设置  ---------------------- ## HDF5目录:/ tmp / RtmpJ73Hov / myOtherFolder # # HDF5压缩级别(0 - 9):0  ## ## -------------------- 数据设置  ---------------------- ## 区域大小:8000 # # # #完全结:40 ## ## --------------------优化参数------------ ##优化方法:Nelder-Mead ##优化控制:## maxit: 60 ## fnscale: -1 ## trace: 1 ##参数估计控制:## eps: 1e-06 ## maxiter: 1000

3．5尺寸因子估计

测序库通常在测序深度上有所不同。在GenoGAM中，这种变化是通过在方程的右边项中添加一个样本特定常数来控制的(1）．这些常数的估计由函数来完成computeSizeFactor ()如下:

##再次读入数据，因为最后一次读入是一个空的GenoGAM ggd <- GenoGAMDataSet(expDesign, directory = wd_folder, design = ~ s(x) + s(x, by =基因型))

## INFO[2019-10-29 21:56:12]正在创建GenoGAMDataSet ## INFO[2019-10-29 21:56:13]读取数据## INFO[2019-10-29 21:56:13]读取wt_1 ## INFO[2019-10-29 21:56:14]读取wt_2 ## INFO[2019-10-29 21:56:14]读取mutant_1 ## INFO[2019-10-29 21:56:14]读取mutant_2 ## INFO[2019-10-29 21:56:14]已完成读取数据## INFO [2019-10-29 21:56:14] GenoGAMDataSet已创建

ggd <- computeSizeFactors(ggd)

##信息[2019-10-29 21:56:15]计算大小因子完成

sizeFactors (ggd)

## wt_1 wt_2 mutant_1 mutant_2 ## -0.04156987 0.04210819 -0.25589212 0.14328182

注意:GenoGAM中的大小因子为自然对数刻度。此外，从设计中自动识别因素组。考虑到更复杂的设计，这可能会失败。或者，可能需要不同的因子组。在这种情况下，组可以通过factorGroups论点。

3.6模型拟合

GenoGAM模型需要进一步拟合两个参数:正则化参数，\λ(\ \)，为色散参数\θ(\ \)．正则化参数\λ(\ \)控制平滑的数量。更大的\λ(\ \)就是，平滑函数越平滑。色散参数\θ(\ \)控制观察到的计数偏离方程所模拟的期望值的程度(1）．离散度捕捉了通常在重复样本中看到的生物和技术变化，但也捕捉了模型的错误。在GenoGAM中，离散度是通过假设计数遵循负二项分布和均值来建模的\ \(μ= \ operatorname {E} (y) \)它的对数是由方程(1）还有方差\(v = \mu + \mu^2/\theta\)．

如果没有提供，则提供参数\λ(\ \)而且\θ(\ \)都是通过交叉验证得到的。这一步对于小插图来说太耗时了。为了快速浏览这个例子，我们手动提供这些值:

没有参数估计的拟合模型拟合<- genogam(ggd, lambda = 4601, theta = 4.51)

## INFO[2019-10-29 21:56:15]初始化模型## INFO[2019-10-29 21:56:15]已完成## INFO[2019-10-29 21:56:25]拟合模型## INFO[2019-10-29 21:56:25]已完成## INFO[2019-10-29 21:56:25]处理拟合## INFO[2019-10-29 21:56:25]处理已完成## INFO[2019-10-29 21:56:25]已完成

适合

##家族:nb ##公式:~ s(x) + s(x, by =基因型)##类别:GenoGAM ##尺寸:100000 4 ##样本(4):wt_1 wt_2 mutant_1 mutant_2 ##设计变量(1):基因型##平滑函数(2):s(x) s(x):基因型##染色体:chrI ## ##大小因素:## wt_1 wt_2 mutant_1 mutant_2 ## -0.04156987 0.04210819 -0.25589212 0.14328182 ## ##交叉验证:未执行## ##样条参数:##结间距:20 ## b样条顺序:2 ##惩罚顺序:2 ## ## Tile设置:##块大小:98000 ## Tile大小:1e+05 ##悬挂大小:1000 ## tiles数量:1 ##评估的基因组范围:## GRanges对象1范围和0元数据列:## seqnames ranges strand ##    ## [1] chrI 1-100000 * ## ------- ## seqinfo: 1个来自未指定基因组的序列

参数估计说明:估计参数\λ(\ \)而且\θ(\ \)通过交叉验证，调用genogam ()没有设定他们的价值。这将对每个具有初始参数值的瓷砖执行k-fold交叉验证，并迭代直到收敛，通常约为30次迭代。可通过参数设置折叠次数kfolds．我们建议对20到40个不同的区域执行此操作(数量可以通过参数设置)地区)．在此过程中，底层优化算法(默认为Nelder-Mead)将探索可能的参数空间。在低计数区域的情况下，这可能导致不稳定的矩阵运算，特别是在线性系统的分辨率。为了稳定计算，可以通过参数添加一个小的一阶正则化每股收益．这将改变处罚条款从\(\beta^T \textbf{S} \beta\)来\(β\ ^ T \ textbf{年代}\β+ \ε\ textbf{我}\)，其中单位矩阵\ (\ textbf{我}\)与罚矩阵的维数相同吗\ (\ textbf{年代}\)．的参数每股收益代表\ε(\ \)在上面的方程中。因此，a太大了\ε(\ \)也会对一般的适应度有显著的影响，而小的呢\ε(\ \)将主要影响接近于零的小值。我们建议\ε(\ \)取值范围为1e-4 ~ 1e-3。

在GenoGAM中交叉验证的一个重要区别是，它没有遗漏数据点，而是忽略了间隔。这通过考虑局部相关性(省略单个点可能导致过拟合)来改进参数估计。如果数据中存在重复，交叉验证可以通过在ChIP-Seq片段大小附近使用更大的间隔(默认为20bp)，即200bp来进一步改进。因此，复制中的数据将捕获遗漏的大间隔。

##不计算fit_CV <- genogam(ggd, intervalSize = 200)

关于并行计算的说明:如果使用SnowParams后台，启动一个worker和复制相关数据可能需要几秒钟的时间。特别是如果这是连续完成的。在交叉验证期间，计算是在小块上完成的，因此计算可能只是总开销的一小部分。因此，在交叉验证和之后重置期间，工作人员的数量将自动设置为最优的较低数量。

3．7策划的结果

按照包装的惯例mgcv而且gam类定义的平滑函数的拟合名称通过设计公式中的变量(与配置中的列名相同)遵循此模式s (x):变．这里是我们感兴趣的平滑函数f \ \(文本{突变/ wt} (x) \)因此命名为s (x):基因型．

range = GRanges("chrI"， IRanges(45000, 55000)) plotGenoGAM(fit, ranges = ranges)

因为数据通常太大，无法绘制整体。绘图函数需要一个特定的范围，最好由农庄对象提供给范围论点。此外，还可以指定平滑是否应该在相同的刻度上绘制(规模)，则数值应以日志格式显示(日志)或是否应包括原始数据计数数据(提供GenoGAMDataSet来ggd不幸的是，在小插图中很难表现出来)。

plotGenoGAM(适合，范围=范围，比例= FALSE)

3.8统计测试

我们在每个平滑和每个位置进行测试\ \ (x)零假设，它的值为0computeSignificance ()．这就给出了逐点的p值。

fit <- computessignificance (fit)

## INFO[2019-10-29 21:56:27]计算位置p值

pvalue(适合)

##数据帧100000行2列## s(x) s(x):基因型## <数字> <数字> ## 11 0.094066901907606 ## 2 0.094066901907606 ## 2 0.0929472532979708e -28 0.0929572815942151 ## 3 0.09935352252415942151 ## 4 2.164084441135602e -28 0.0894804406588641 ## ... ... ...2.31805863114621e-32 0.60198351364854 ## 99999 3.27429977098754e-32 0.602859394685401 ## 100000 4.62426303076389e-32 0.60373476562017

3.9微分绑定

如果需要区域显著性，如在差分绑定的情况下，则调用computeRegionSignificance ()．这将返回提供的农庄由p-value和FRD列更新的对象。如果光滑的如果它是为所有平滑计算的，则不为。

gr <- GRanges("chrI"， IRanges(c(1000, 7000)， c(4000, 9000))) db <- computeRegionSignificance(fit, regions = gr, smooth = "s(x):基因型")

## INFO[2019-10-29 21:56:27]估计区域p值和FDR ## INFO[2019-10-29 21:56:27]计算位置日志p值## INFO[2019-10-29 21:56:28]完成## INFO[2019-10-29 21:56:28]执行多重修正## INFO[2019-10-29 21:56:28]完成

db

## seqnames ranges strand | pvalue FDR ##    | <数字> <数字> ## [1]chrI 1000-4000 * | 0.249909462695555 0.24990946462695555 ## [2] chrI 7000-9000 * | 0.249535616079423 0.249909462695555 ## ------- # seqinfo: 1个来自未指定基因组的序列;没有seqlengths

3.10峰打电话

所计算的拟合允许通过函数轻松调用窄峰和宽峰callPeaks ()．由此产生的data.table对象提供了类似于ENCODE窄峰和宽峰格式的结构。请注意,分数是p值的负自然对数。此外，没有提供峰值边界(就像在GenoGAM 1.0中实现的那样，一个计算它们的函数)。由于这不是峰值呼叫数据，我们使用高阈值来演示功能。

<- callPeaks(fit, threshold = 1, smooth = "s(x):基因型")

##信息[2019-10-29 21:56:28]完成

山峰

## $ ' s(x):基因型' ## seqnames pos zscore score fdr summit ## 1: chrI 77623 3.953526539 10.1647521 0.0000000 0.8218586 ## 2: chrI 59362 2.299714606 4.53455556 0.5555556 1.4055236 ## 3: chrI 62226 2.605805754 5.3854186 0.5714286 1.7176303 ## 4: chrI 59883 2.632129470 5.4625846 0.6000000 1.3510637 ## 5: chrI 81446 2.233273750 4.3610131 0.6000000 1.3675995 ## 6: chrI 83493 2.307658693 4.5555672 0.6250000 1.55556027 ## 7: chrI 90239 2.611568704 5.40225708: 0.6666667 0.9923707 ## 8:11: chrI 22332 0.274332932 0.9367119 0.8367347 4.2479926 ## 12: chrI 2357 0.010564995 0.7016124 0.8474576 3.0423462 ## 13: chrI 57733 0.184091053 0.8510392 0.8490566 1.4130322 ## 15: chrI 44866 0.130774568 0.8030138 0.8518519 2.0274223 ## 16:19: chrI 38415 0.381625912 1.0459171 0.8604690 1.5213175 ## 20: chrI 23137 0.015640543 0.7057045 0.8620690 4.6110303 ## 21: chrI 21262 0.200954529 0.8666270 0.8653846 1.7848703 ## 22: chrI 14576 0.338250353 1.0007945 0.8695652 3.6128455 ## 23: chrI 27850 0.257966317 0.9207596 0.8800000 1.9101958 ## 24:chrI 15951 0.387984743 1.0526442 0.8809524 1.5720949 ## 25: chrI 61001 0.204287531 0.8697307 0.8823529 0.883636 1.2718125 ## 27: chrI 31237 0.343267025 1.0059453 0.8888889 3.1725112 ## 28: chrI 13524 0.395600685 1.0607392 0.9024390 1.6710578 ## 30: chrI 23903 0.623616735 1.3226075 0.9354839 1.3912920 ## 31: chrI 49699 1.245915796 2.2405717 0.9375000 1.3353167 ## 32:chrI 4797 0.455874666 1.1262711 0.9459459 1.4687162 ## 33: chrI 41399 0.422117971 1.0892479 0.9487179 2.0248993 ## 35: chrI 76493 -0.246390343 0.5153191 0.1196190 0.9571429 2.2251183 ## 36: chrI 39240 -1.212076272 0.1196190 0.9571429 2.2251183 ## 37: chrI 57191 0.624069757 1.3231660 0.9666667 0.7640392 ## 37: chrI 76594 -0.229948090 0.5260510 0.9672131 2.2555678 ## 38: chrI 28556 -0.338869459 0.4578443 0.9682540 1.3313444 ## 40: chrI 9026 -0.388672158 0.4288760 0.9692308 4.1543744 ## 40:chrI 51485 0.557091881 0.9696970 2.0802139 ## 41: chrI 936 -0.814219812 0.2328903 0.9705882 1.9832077 ## 43: chrI 37176 -1.054367667 0.1576570 0.9710145 1.9079031 ## 44: chrI 47888 0.469061169 1.1409575 0.9722222 1.2486934 ## 44: chrI 73443 0.434946481 1.1032210 0.9736842 1.4434126 ## 45: chrI 68581 -0.375111593 0.4366288 0.9843750 3.6005765 ## 46: chrI 34811 -0.795611070 0.2396910 0.9850746 1.8922657 ## 48:chrI 86031 1.509538353 2.7244749 1.0000000 1.5380176 ## 49: chrI 98729 1.479515619 2.6664101 1.0000000 1.2249742 ## 50: chrI 75279 1.2946257689506 1.3612746 1.8438616 ## 52: chrI 12615 0.627263028 1.3271073 1.0000000 0.8576204 ## 53: chrI 17336 0.604812635 1.2995599 1.0000000 0.7174821 ## 54: chrI 46915 0.506105669 1.1828924 1.0000000 1.2082224 ## 56: chrI 87700 0.470411916 1.1424690 1.0000000 1.2082224 ## #* * * * * * * * * * * *

<- callPeaks(fit, threshold = 1, peakType = "broad"， cutoff = 0.75, smooth = "s(x):基因型")

##信息[2019-10-29 21:56:29]完成

山峰

# # $ ' s (x):基因型的# # seqnames开始结束宽度链得分meanSignal罗斯福# # 1:chrI 1 578 578 * 0.2881752 1.3225292 0.8687306 # # 2: chrI 1937 5520 3584 * 0.2880604 1.0858642 0.8687306 # # 3: chrI 6557 11060 4504 * 0.2460828 1.1853691 0.8687306 # # 4: chrI 11900 26025 14126 * 0.2880586 1.2344399 0.8687306 # # 5: chrI 26974 29868 2895 * 0.2883320 1.1872678 0.8687306 # # 6: chrI 30637 33619 2983 * 0.2880309 1.3303773 0.8687306 # # 7: chrI 37978 38914 937 * 0.2883072 1.1866105 0.8687306 # # 8:12: chrI 46322 48244 1923 * 0.2884650 1.3329063 0.8687306 ## 13: chrI 68366 68770 405 * 0.2883727 0.9572166 0.8687306 ## 14: chrI 69998 71392 1395 * 0.2880892 1.4932166 0.8687306 ## 15: chrI 72905 79495 6591 * 1.3726813 2.2360777 0.8687306 ## 16:17: chrI 84745 86766 2022 * 0.2884386 1.6236602 0.8687306 ## 18: chrI 87351 92044 4694 * 0.2886375 2.3349403 0.8687306 ## 19: chrI 52845 66840 13996 * 0.1385248 1.7610177 0.8706416 ## 20: chrI 92840 100000 7161 * 0.1497634 1.5247567 0.8706416

这个函数writeToBEDFile ()提供了一种简单的方法来写入峰值数据。表到narrowPeaks或broadPeaks文件。后缀将自动确定

writeToBEDFile(peaks, file = 'myPeaks')

4．1犰狳错误

问题:在模型拟合过程中会出现以下错误:

错误:矩阵乘法:不兼容矩阵尺寸:22333830147200x5360120267008000和4294972496x1

解决方案:首先，确保正确安装了所有Armadillo依赖项。看到在这里

其次，错误很可能与Armadillo使用32位矩阵有关，从而导致GenoGAM使用的大型矩阵出现问题。解决方案是启用ARMA_64BIT_WORD，在RcppArmadillo中默认不启用。这应该在编译期间完成，但如果由于某种原因失败了，您可以手动执行#定义ARMA_64BIT_WORD在my_R_Directory / lib / R /图书馆/ RcppArmadillo / include / RcppArmadilloConfig.h．看到在这里．