bamcalib 工具

bamcalib 是一个小型 Rust 程序，使用基于 Hu 等人原始方法的归一化过程对校准的 Chip-Seq 数据进行归一化。

安装

您需要已安装 cargo

• curl --proto '=https' --tlsv1.2 -sSfhttps://sh.rustup.rs | sh
• sudoapt install cargo
• sudopacman -Srust
• brewinstall rust

然后运行以下命令

cargo install bamcalib --version 0.1.5

您还可以使用 lbmc/bamcalib:0.1.5 Docker 镜像

docker run -it lbmc/bamcalib:0.1.5 bamcalib

用法

示例

bamcalib \
  -i IP_11_sorted.bam \
  -w T_11_sorted.bam \
  -bigwig-ip IP_11.bigwig \
  -bigwig-WCE T_11.bigwig \
  -t 8

bamcalib 预期以下参数

• -i, --ip-bam <bam_ip> 用于 IP 数据的排序 bam 文件
• -w, --wce-bam <bam_wce> 用于 WCE 数据的排序 bam 文件
• -bigwig-ip, --bigwig-ip <bigwig> 输出的归一化 bigwig 文件名
• -bigwig-wce, --bigwig-wce <bigwig> 输出的归一化 bigwig 文件名

可选参数

• 为校准基因组（默认：calib_）的染色体名称添加校准前缀的代码-c，--cal-prefix <cal_prefix>
• 打印帮助信息的命令-h，--help
• 将OR值乘以此值（必须在所有样本中相同）的命令-s，--scaling <SCALING> [默认：1000]
• 移除与该位代码匹配的读取的命令-n，--no_bits <no_bits>。您可以使用此网站获取适当的位代码。默认值是1540，用于
  • 段未映射
  • 未通过质量控制
  • PCR或光学重复
• 解析bam文件和写入bigwig文件的线程数的命令-t，--thread <thread>
• 不进行片段估计的命令--no-fragment-estimation从读取信息计算覆盖率而不是片段估计
• 输出normIP和normWCE bigwig文件的命令--intermediate-bigwig（而不仅仅是两个比率）
• 打印版本信息的命令-V，--version

方法

归一化

Hu等人提出了以下归一化程序，以从

• 在参考基因组上的IP样本位置$t$的覆盖率$IP_{x}\left(t\right)$
• 在IP样本上的校准基因组位置$t$的覆盖率$IP_{c}\left(t\right)$
• 在WCE样本上的参考基因组位置$t$的覆盖率$WCE_{x}\left(t\right)$
• 在校准基因组上的WCE样本位置$t$的覆盖率$WCE_{c}\left(t\right)$

大小为$T_x$的参考基因组（忽略染色体分段）和大小为$T_c$的校准基因组中，他们提出计算以下$\text{OR}$因子

\text{OR} = \frac{10^6 \sum_{t^{'}=1}^{T_c}WCE_{c}\left(t^{'}\right)}{\sum_{t^{'}=1}^{T_c}IP_{c}\left(t^{'}\right)\sum_{t^{'}=1}^{T_x}WCE_{x}\left(t^{'}\right)}

相反，我们提出以下公式，其中$\beta$（默认为$10^3$）是一个任意缩放因子。

\text{norm}IP\left(t\right) = IP_{x}\left(t\right) \frac{\beta \frac{1}{T_c}\sum_{t^{'}=1}^{T_c}WCE_{c}\left(t^{'}\right)}{\frac{1}{T_c}\sum_{t^{'}=1}^{T_c}IP_{c}\left(t^{'}\right)\frac{1}{T_x}\sum_{t^{'}=1}^{T_x}WCE_{x}\left(t^{'}\right)}

我们希望通过缩放校准基因组覆盖率来消除Chip-Seq实验对$IP_{x}\left(t\right)$的技术影响

\frac{1}{\frac{1}{T_c}\sum_{t^{'}=1}^{T_c}IP_{c}\left(t^{'}\right)}

我们希望通过缩放缩放后的WCE覆盖率来消除Chip-Seq实验对$IP_{x}\left(t\right)$的细胞比例差异影响

\frac{\frac{1}{T_c}\sum_{t^{'}=1}^{T_c}WCE_{c}\left(t^{'}\right)}{\frac{1}{T_x}\sum_{t^{'}=1}^{T_x}WCE_{x}\left(t^{'}\right)}

为了更好地分析基因组重复区域的覆盖率信息，我们建议将以前的归一化扩展到按核苷酸归一化信号，并使用以下以下OR比率

\text{ratio}IP\left(t\right) = \frac{\text{norm}IP\left(t\right)}{\text{norm}WCE\left(t\right)}

与

\text{norm}WCE\left(t\right) = WCE_{x}\left(t\right) \frac{1}{\frac{1}{T_c}\sum_{t^{'}=1}^{T_c}WCE_{c}\left(t^{'}\right)}  \alpha

我们找到$\alpha$，使得

E\left(\text{norm}IP\left(t\right)\right)= E\left(\frac{\text{norm}IP\left(t\right)}{\text{norm}WCE\left(t\right)}\right)

这给出

\text{norm}WCE\left(t\right) = WCE_{x}\left(t\right) \frac{1}{\sum_{t^{'}=1}^{T_x}IP_{x}\left(t^{'}\right)} \sum_{t^{'}=1}^{T_x}\frac{IP_{x}\left(t^{'}\right)}{WCE_{x}\left(t^{'}\right)} 

注意，所有计算都按基因组大小缩放，而不是按读取数量进行，如Hu等人中所述，这也是为什么我们添加缩放因子$\beta$（默认为$10^3$）的原因。

使用这种方法，我们保留了Hu等人归一化在样本间平均读取密度上的有趣特性（即，我们可以以定量方式比较两个不同的样本）并且我们考虑了WCE样本中观察到的局部读取密度偏差（差异染色质可及性、重复、低可映射区域等）。

基因组覆盖率密度

为了计算覆盖密度 $X_y(t)$，其中 $X \in \left[IP, WCE\right]$ 且 $y \in \left[c, x\right]$，我们统计与位置 $t$ 重叠的读取数 $r(t)$。对于正确配对的读取（具有在同一染色体上的配对读取，并且起始位置在读取末尾之后结束），我们还计算第一个读取的末尾和其配对读取的起始之间密度为 1。$X_y(t) = r(t) + g(t)$。

一些片段可能被人为地拉长，因此，我们通过去除片段长度 0.1 的上下值来计算两个配对读取之间片段大小稳健均值 $\mu$。大小超过 $\phi^{-1}(0.95, \mu, 1.0)$ 的片段将被设置为在 $\phi^{-1}(0.95, \mu, 1.0)$ 值处结束，其中 $\phi()$ 是正态累积分布函数。

一些片段可能比读取长度短，在这种情况下，我们不将重叠读取区域计为 2 个片段的覆盖，而是计为 1 个片段的覆盖。