1

我是序列分析的新手，我正在做一些练习来帮助我学习使用 pysam 和 samtools 进行 WGS 数据分析。我想做的一件事是从二维牛津纳米孔数据（大读数）中检测（相当大的）缺失。为此，我从大肠杆菌基因组中提取了前 10kb 以及覆盖该区域的测序读数。调用原始基因组 A。然后我通过在 A 的中间插入 1kb 序列来创建基因组 A'，并使用 A' 作为参考来对齐 A 的读取以模仿序列中的删除。我现在想编写一个程序来检测我的“删除”的位置。我的问题是我读取的 CIGAR 字符串不符合我的期望，我认为这一定是错误的。

假设我有一个序列 ....GTTGCA ---1kb 删除--- GAACGT... 并且读取与该序列对齐。我做出以下假设：

案例 1. 删除左侧且不与删除重叠的读取可以以 aHbS（a 和 b 为常数，a,b >=0）开始，后跟一系列 Ms、Is、Ds，然后以 cSdH 结束。我不希望在这些读取中找到大段 Is 和 Ds。

案例 2. 从左侧部分与删除重叠的读取应与 (1) 中的读取相同，但应以 rS 结尾，常数 r 的大小取决于读取与删除重叠的程度。

案例 3. 读取与删除完全重叠（请记住，我有很长的读取，所以存在这样的读取）应该与 (1) 中的读取相同，但我希望在我的 CIGAR 字符串中看到 1000D 或类似的东西，然后读取应与 (1) 中的读取相同。这是我在数据中没有观察到的。我的“删除”从 5kb 开始，但具有 4500 < POS < 5000 且长度超过 2kb 的读取实际上似乎包含相同的 Ms、Is 和 Ds 序列，就好像它们与参考对齐一样。

我的问题，我希望不是离题，因为我宁愿询问数据格式而不是实际编程，是 i)。我上面的哪个假设是错误的 ii)。读取部分重叠删除的雪茄串应该是什么样子？三）。读取完全重叠的雪茄串（也就是说，其末端映射在删除的任一侧）删除看起来像什么？

我附上了一个图，希望能帮助说明我的三个案例。

在此处输入图像描述

4

2 回答 2

0

如果我正确理解您的问题，我认为您的假设听起来不错。您所描述的内容（CIGAR 字符串中未表示删除）让我认为，虽然您更改了参考基因组（可能是 *.fasta 文件？），但您可能没有重新运行您的读取对齐到那个参考。

您可能从中获取 CIGAR 字符串的 BAM/SAM 文件是比对的结果，当您仅更改参考基因组时，它本身不会改变。更改参考基因组后，您现在需要重新进行比对以获取新的 BAM/SAM 文件，然后查看其中的 CIGAR 字符串，现在应该反映模拟删除。请让我知道这是否是一个正确的评估。

（我知道这将是一个更值得评论的帖子，但我还没有写评论的权限。）

于 2019-11-21T19:01:36.123 回答
0

我不是这个主题的专家，但我敢打赌，案例 1 和 2 就像你说的那样。案例 3 可能因使用的矫正器而异。一种表示这种对齐方式的方法如您所料：xHxSxM1000DxMxSxH，但是拆分映射器可以提供另一种方式：xHxSxM1000SxH 加上 xH1000SxMxSxH 在两个不同的条目中，其中一个标记为主要对齐，另一个标记为补充对齐（一些旧的对齐器可以将其标记为次要，因为补充在标准后面出现）。对准器在其表示方式中起着至关重要的作用。

您是否检查过您的完全重叠读取在您的 sam/bam 文件中只表示一次？

于 2018-06-29T10:10:19.433 回答