使用 CIGAR 推断序列的长度答案

【问题标题】：Infer the length of a sequence using the CIGAR使用 CIGAR 推断序列的长度
【发布时间】：2017-02-04 06:24:27
【问题描述】：

为您提供一些背景信息：我正在尝试将 sam 文件转换为 bam

samtools view -bT reference.fasta sequences.sam > sequences.bam

退出并出现以下错误

[E::sam_parse1] CIGAR and query sequence are of different length
[W::sam_read1] parse error at line 102
[main_samview] truncated file

而有问题的行看起来像这样：

SRR808297.2571281       99      gi|309056|gb|L20934.1|MSQMTCG   747     80      101M    =       790     142     TTGGTATAAAATTTAATAATCCCTTATTAATTAATAAACTTCGGCTTCCTATTCGTTCATAAGAAATATTAGCTAAACAAAATAAACCAGAAGAACAT      @@CFDDFD?HFDHIGEGGIEEJIIJJIIJIGIDGIGDCHJJCHIGIJIJIIJJGIGHIGICHIICGAHDGEGGGGACGHHGEEEFDC@=?CACC>CCC      NM:i:2  MD:Z:98A1A

我的序列长度为 98 个字符，但在创建 sam 文件时可能存在错误，在 CIGAR 中报告了 101。我可以让自己奢侈地丢失几次读取，并且我目前无法访问生成 sam 文件的源代码，因此没有机会寻找错误并重新运行对齐。换句话说，我需要一个务实的解决方案来继续前进（目前）。因此，我设计了一个 Python 脚本，它计算我的核苷酸字符串的长度，将其与 CIGAR 中注册的内容进行比较，并将“正常”的行保存在一个新文件中。

#!/usr/bin/python
import itertools
import cigar

with open('myfile.sam', 'r') as f:
    for line in itertools.islice(f,3,None): #Loop through the file and skip the first three lines
            cigar=line.split("\t")[5]
            cigarlength = len(Cigar(cigar)) #Use module Cigar to obtain the length reported in the CIGAR string
            seqlength = len(line.split("\t")[9])

            if (cigarlength == seqlength):
                    ...Preserve the line in a new file...

如您所见，为了将 CIGAR 转换为显示长度的整数，我使用了模块 CIGAR。老实说，我对它的行为有点警惕。在非常明显的情况下，这个模块似乎错误地计算了长度。是否有其他模块或更明确的策略将 CIGAR 转换为序列的长度？

旁注：有趣的是，至少可以说，这个问题已被广泛报道，但在互联网上找不到实用的解决方案。请参阅以下链接：

https://github.com/COMBINE-lab/RapMap/issues/9
http://seqanswers.com/forums/showthread.php?t=67253
http://seqanswers.com/forums/showthread.php?t=21120
https://groups.google.com/forum/#!msg/snap-user/FoDsGeNBDE0/nRFq-GhlAQAJ

【问题讨论】：

“说实话，我对它的行为有点警惕。这个模块在非常明显的情况下似乎错误地计算了长度。” - 比如？ underlying code 短小精悍；很难想象它会以何种方式失败。

标签： python module bioinformatics samtools bam

【解决方案1】：

SAM spec 为我们提供了这张 CIGAR 操作表，它指示哪些“消耗”查询或引用，并附有关于如何从 CIGAR 字符串计算序列长度的明确说明：

                                                             Consumes  Consumes
Op  BAM Description                                             query  reference
M   0   alignment match (can be a sequence match or mismatch)   yes   yes
I   1   insertion to the reference                              yes   no
D   2   deletion from the reference                             no    yes
N   3   skipped region from the reference                       no    yes
S   4   soft clipping (clipped sequences present in SEQ)        yes   no
H   5   hard clipping (clipped sequences NOT present in SEQ)    no    no
P   6   padding (silent deletion from padded reference)         no    no
=   7   sequence match                                          yes   yes
X   8   sequence mismatch                                       yes   yes

“Consumes query”和“consumes reference”表示 CIGAR 操作是否导致对齐分别沿着查询序列和参考序列进行。

...

M/I/S/=/X 操作的长度之和应等于 SEQ 的长度。

这使我们可以通过将 CIGAR 中所有“消费查询”操作的长度相加来简单地从其 CIGAR 计算序列的长度。这正是雪茄模块中发生的事情（参见https://github.com/brentp/cigar/blob/754cfed348364d390ec1aa40c951362ca1041f7a/cigar.py#L88-L93），所以我不知道为什么这里的 OP 认为该模块的实现是错误的。

如果我们从（已经很短的）雪茄模块中提取相关代码，我们会得到类似这样的东西，作为上面引用中描述的求和操作的简短 Python 实现：

from itertools import groupby

def query_len(cigar_string):
    """
    Given a CIGAR string, return the number of bases consumed from the
    query sequence.
    """
    read_consuming_ops = ("M", "I", "S", "=", "X")
    result = 0
    cig_iter = groupby(cigar_string, lambda chr: chr.isdigit())
    for _, length_digits in cig_iter:
        length = int(''.join(length_digits))
        op = next(next(cig_iter)[1])
        if op in read_consuming_ops:
            result += length
    return result

【讨论】：

【解决方案2】：

我怀疑没有解决此问题的工具的原因是因为没有通用解决方案，除了使用不会出现此问题的软件再次执行对齐之外。在您的示例中，查询序列与参考完全对齐，因此在这种情况下，CIGAR 字符串不是很有趣（只是一个以整体查询长度为前缀的 Match 操作）。在这种情况下，修复只需将101M 更改为98M。

但是，对于更复杂的 CIGAR 字符串（例如，包含 Insertions、Deletions 或任何其他操作的字符串），您将无法知道 CIGAR 字符串的哪一部分太长。如果您从 CIGAR 字符串的错误部分中减去，您将得到一个未对齐的读数，这对于您的下游分析可能比只留下整个读数更糟糕。

也就是说，如果它恰好是微不足道的（也许您损坏的对齐过程总是在第一个或最后一个 CIGAR 操作中添加额外的碱基），那么您需要知道的是计算查询长度的正确方法根据 CIGAR 字符串，以便您知道要从中减去什么。

samtools 使用htslib 函数bam_cigar2qlen 进行计算。

bam_cigar2qlen 调用的其他函数在sam.h 中定义，包括一个有用的注释，显示操作消耗查询序列与参考序列的真值表。

简而言之，要以 samtools（实际上是 htslib）的方式计算 CIGAR 字符串的查询长度，您应该为 CIGAR 操作 M、I、S、= 添加给定的长度, 或 X 并忽略任何其他操作的 CIGAR 操作的长度。

python cigar 模块的当前版本似乎正在使用same set of operations，并且计算查询长度的算法（这是len(Cigar(cigar)) 将返回的）对我来说看起来很合适。是什么让您认为它没有给出正确的结果？

看起来您应该能够使用 cigar python 模块使用 mask_left 或 mask_right 方法和 mask="H" 从左端或右端硬剪辑。

【讨论】：