从 CIGAR 字符串中检测大量删除答案

【问题标题】：Detect large deletions from CIGAR string从 CIGAR 字符串中检测大量删除
【发布时间】：2017-09-10 04:17:05
【问题描述】：

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

假设我有一个序列 ....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

我的问题是一世）。我上面的哪个假设是错误的 ii)。读取部分重叠删除的雪茄串应该是什么样子？三）。读取完全重叠的雪茄串（也就是说，其末端映射在删除的任一侧）删除看起来像什么？

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

【问题讨论】：

你应该问biostars.org

标签： bioinformatics samtools pysam

【解决方案1】：

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

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

【讨论】：

【解决方案2】：

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

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

（我知道这篇文章更值得发表评论，但我还没有编写 cmets 的权限。）

【讨论】：