怎么用python计算序列长度_python – 使用CIGAR推断序列的长度

为了给你一些上下文：我试图将一个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个字符,但在创建在CIGAR中报告101的sam文件时可能存在错误.我可以给自己丢失一些读取的奢侈品,而且我目前无法访问产生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

最佳答案 我怀疑没有工具来解决这个问题的原因是因为没有通用的解决方案,除了使用没有出现此问题的软件再次执行对齐.在您的示例中,查询序列与引用完全对齐,因此在这种情况下,CIGAR字符串不是很有趣(只是以整个查询长度为前缀的单个Match操作).在这种情况下,修复只需要将101M更改为98M.

但是,对于更复杂的CIGAR字符串(例如包含插入,删除或任何其他操作的字符串),您将无法知道CIGAR字符串的哪个部分太长.如果从CIGAR字符串的错误部分中减去,则会留下未对齐的读数,这对于您的下游分析而言可能比仅留下整个读数更糟糕.

也就是说,如果它恰好是正确的(可能你的破坏的对齐程序总是为第一个或最后一个CIGAR操作添加额外的基础),那么你需要知道的是根据这个计算查询长度的正确方法. CIGAR字符串,以便您知道从中减去什么.

samtools使用htslib函数bam_cigar2qlen计算它.

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

简而言之,要按照samtools(真正的htslib)的方式计算CIGAR字符串的查询长度,您应该为CIGAR操作M,I,S,=或X添加给定长度,并忽略CIGAR操作的长度.任何其他操作.

python雪茄模块的当前版本似乎使用same set of operations,并且用于计算查询长度的算法(这是len(雪茄)将返回的)看起来对我来说是正确的.是什么让你认为它没有给出正确的结果？

看起来你应该能够使用带有mask =“H”的mask_left或mask_right方法,使用雪茄python模块从左端或右端硬夹.