问题标签 [samtools]

在SAM格式中，每条对齐线代表一个段的线性对齐，每条线有11个必填字段，即QNAME、FLAG、RNAME、POS、MAPQ等。

假设我想要一个给定 BAM 文件中所有“QNAMES”的 NumPy 数组。或者，可以采用几列并将它们导入 Pandas Dataframe。

pysam 可以实现此功能吗？

人们可以很自然地用 来打开给定的 BAM 文件，pysam.AlignmentFile()然后用 来访问各个段pysam.AlignmentSegment()，例如

但是，您可以将所有 QNAMES 保存到 NumPy 数组中吗？

是linux，使用conda安装pysam，pip install pysam一直失败。成功安装 pysam 后，pysam 显示conda list并出现，anaconda2/pkgs/ 但是import pysam在 python 2.7.12 中，它失败了Traceback (most recent call last): File "<stdin>", line 1, in <module> ImportError: No module named pysam

请帮忙。

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

退出并出现以下错误

违规行如下所示：

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

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

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

最初的问题

我正在用python（3.5）编写一个生物信息学脚本，它解析一个表示在基因组上对齐的测序读数的大型（排序和索引） bam文件，将基因组信息（“注释”）与这些读数相关联，并计算遇到的注释类型. 我正在测量我的脚本处理对齐读取的速度（超过 1000 次读取），我获得了以下速度变化：

什么可以解释这种模式？

我的直觉会让我押注一些数据结构，随着它变得越来越密集，它会逐渐变慢，但它会不时扩展。

不过，内存使用似乎并不显着（运行近 2 小时后，我的脚本仍然只使用了我计算机内存的 0.1%，根据htop）。

我的代码是如何工作的（见最后的实际代码）

我正在使用该pysam模块进行 bam 文件解析。AlignmentFile.fetch方法为我提供了一个迭代器，以AlignedSegment对象的形式提供有关连续对齐读取的信息。

我根据它们的对齐坐标和gtf格式的注释文件（使用 bgzip 压缩并使用 tabix 索引）将注释与读取相关联。我使用TabixFile.fetch方法（仍然 from pysam）来获取这些注释，我过滤它们并以frozenset字符串的形式生成它们的摘要（process_annotations，下面未显示，返回这样的 a frozenset），在内部循环的生成器函数中在 AlignedSegment 迭代器上。

我将生成的frozensets 提供给一个Counter对象。计数器能否对观察到的速度行为负责？

我怎样才能知道如何避免这些定期减速？

附加测试

根据评论中的建议，我分析了我的整个分析cProfile，发现在访问注释数据时花费了最多的运行时间（pysam/ctabix.pyx:579(__cnext__)请参阅稍后的调用图），如果我理解正确的话，这是一些与 samtools C 库交互的 Cython 代码. 观察到的减速的原因似乎很难理解。

为了加快我的脚本速度，我尝试了另一种基于pybedtoolspython 接口和 bedtools 的解决方案，它还可以从 gtf 文件 ( https://daler.github.io/pybedtools/index.html ) 中检索注释。

速度

速度提升非常重要。以下是实际的命令和计时结果（两者实际上是并行运行的）：

（需要注意：两种方法的注释结果略有不同。也许我应该检查一下我如何处理坐标。）

速度曲线没有先前观察到的长期下降：

我的速度问题已经解决，但我仍然对为什么基于 tabix 的方法有这些速度下降感兴趣。为此，我添加了“生物信息学”和“samtools”标签。

调用图

作为记录，我在分析结果上使用 gprof2dot 生成了调用图：

以下是基于 tabix 方法的调用图：

对于基于 pybedtools 的方法：

编码

这是我当前代码的主要部分：

（我使用了上下文管理器和其他高阶函数（make_annotation_processor以及这个函数调用的函数），以便更容易在同一个脚本中使用各种注释检索方法。）

我有以下更改染色体名称的 sed 命令：

我的问题是如何将结果插入新路径，同时将变量 $filename 作为每个新文件名的一部分？它总是将结果插入 /myoldpath/ 或 /mynewpath/ 中的字面意思“filename.chr.bam”。我在该部分的语法中遗漏了什么$file > /mynewpath/${filename}_chr.bam吗？

任何帮助，将不胜感激

我是序列分析的新手，我正在做一些练习来帮助我学习使用 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)。读取部分重叠删除的雪茄串应该是什么样子？三）。读取完全重叠的雪茄串（也就是说，其末端映射在删除的任一侧）删除看起来像什么？

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

我目前正在尝试将 .fasta 文件导入 bwa 以使用参考基因​​组将我的读取映射到。但是，我目前收到此错误：

有什么帮助吗？这是我的代码：

关于如何考虑可以使用需要（稍微）不同的 shell 调用的软件版本的管道的想法？

有时，在带有conda的版本之间切换，shell 调用是不同的，例如 Samtools 0.1.18 和 Samtools 1.3.1。前缀需要不同的格式。

我已经想到了三种方法，并正在寻找其他建议：

1. 我将在我的 YAML 中添加一个配置变量，用它来设置版本号。有一个条件，使得只包含匹配的版本规则。与我在这里使用软件选择标志的方式类似。

2. 我将在规则中进行版本检测，并相应地更改 shell 调用。在上面的示例中，唯一的区别是添加了“-T”参数。这很简单，但我担心最终我会遇到一种情况，它不仅仅是一个额外的论点。

3. 编写完全独立的规则，规则名称中包含版本，让用户有责任选择正确的版本。这将不可避免地导致模棱两可的冲突，但我可以处理这些。

反思： 我不相信嵌套条件是最好的方法，因为它使维护更加困难，而且真的不是那么优雅。在规则中进行检测并没有那么糟糕，但我不喜欢将控制流推送到规则中的想法（现在它已经从我的管道规则中删除了）。我想避免产生歧义冲突。

我缺少支持功能吗？起初我以为这是Snakemake 的 'Version' 指令，但实际上并不是我想要的。要么，要么我错过了如何使用它。

想法？

支持文档

山姆工具 1.3.1

1.3.1 ---> config["samtools"] + ' sort -no -m ' + str(config["sortMem"]) + ' - -T' + str(log.stdErr)

山姆工具 0.1.18

0.1.18 --> config["samtools"] + ' sort -no -m ' + str(config["sortMem"]) + ' - ' + str(log.stdErr)

不同之处在于 str(log.stdErr) 之前的“-T”。

I'm trying to install samtools with conda (Ubuntu 16.04):

But i get this index error:

Somebody knows about this? I have tried unsuccessfully to update it..

我正在编写一个类的函数来将 sam 文件转换为 bam 文件，然后对其进行排序并在 python 中对其进行索引（我知道在 bash 中很容易实现，但了解更多信息并没有什么坏处）。

类的另一个函数已经在定义的目录 (cwd) 中生成了一个 sam 文件。

这是代码：

但是什么都没有出来

我试过了：

但也没有。尝试将其保存为变量并执行 .communicate()，什么也没有。

所以我不知道我做错了什么？

谢谢，

Xp