用雪茄推断序列的长度

Question

给你一点上下文：我正在尝试将一个sam文件转换成bam

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

出现以下错误

^{pr2}$

冒犯的底线是这样的：

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文件时可能有一个bug，报告了雪茄中的101个字符。我可以给自己一个奢侈的机会来丢失几次读取，而且我目前无法访问生成sam文件的源代码，因此没有机会查找bug并重新运行对齐。换句话说，我需要一个务实的解决方案来继续（目前）。因此，我设计了一个python脚本来计算我的核苷酸字符串的长度，将其与雪茄中注册的内容进行比较，并将“正常”的行保存在一个新文件中。在

#!/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。老实说，我对它的行为有点谨慎。在非常明显的情况下，这个模块似乎计算错了长度。是否有其他模块或更明确的策略来将雪茄转化为序列的长度？在

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

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

Answer 1

SAM spec为我们提供了这个雪茄操作表，指明哪些操作“使用”查询或引用，以及如何从雪茄字符串计算序列长度的明确说明：

                                                             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” and “consumes reference” indicate whether the CIGAR operation causes the alignment to step along the query sequence and the reference sequence respectively.

...

• Sum of lengths of the M/I/S/=/X operations shall equal the length of SEQ.

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

如果我们从（已经非常短的）雪茄模块中提取出相关的代码，那么我们只剩下这样一个简短的Python实现，它是上面引述的求和操作的一个简短的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

Answer 2

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

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

这就是说，如果它碰巧是一个琐碎的正确处理（也许你的破坏的对齐过程总是给第一个或最后一个雪茄操作添加额外的基），那么你需要知道的是根据雪茄字符串计算查询长度的正确方法，这样你就知道从中减去什么。在

samtools使用htslib函数bam_cigar2qlen计算这个值。在

bam_cigar2qlen调用的其他函数在sam.h中定义，包括一个有用的注释，该注释显示了操作使用查询序列和引用序列的真值表。在

简而言之，要像samtools（实际上是htslib）那样计算雪茄字符串的查询长度，您应该为雪茄操作M、I、S、=、或{{}添加给定长度，并忽略任何其他操作的雪茄操作长度。在

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

看起来您应该能够使用雪茄python模块使用mask_left或{}方法和mask="H"从左端或右端进行硬剪辑。在