问题标签 [samtools]

我正在尝试通过使用此处的 seqff 脚本来研究胎儿分数。该指令说我需要一个无头 sam 文件才能工作，因此在与 对齐后bwa，我使用以下方法处理我的 sam 文件samtools：

由于没有参考文件，此过程还会删除 @SG 标头。我还使用grep删除对齐后得到的 target.sam 的标题以进行比较。然而，当我比较时，我注意到对处理的.sam 的 ENet 计算返回了 NA 结果，因此该脚本中计算胎儿分数的方法的一半已经消失。

我已经在多种情况下对其进行了测试，似乎没有使用samtools，脚本运行良好，但samtoolsEnet calculaton 不起作用。有人可以告诉我脚本有什么问题吗？

我正在使用 R 4.0，我对脚本所做的所有修改都在 args 上，因为它无法识别我输入的 args。

我想使用pysam库将我的 SAM 文件转换为 BAM。我正在尝试从 samtools 编写此命令：

作为这样的python代码：

根据本文档。但不幸的是，它不起作用。我该怎么做？我应该先创建一个空的 BAM 文件还是自动创建它？最后我有这个错误：

SamtoolsError: 'samtools 返回错误 1: stdout=b'\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff\x06\x00BC\x02\x004\x01\x9dQMO\x83@\x10 \xa5G\xfa#\xccz\xb3\x87\xfd\xa2\xa0\x86\x83\xa9\xad\xd86i\t\x16\xc3u\xb3P\xc0&,PvcC\xff\xa6\x7f\xc8\xc5\ x88&\xe8\xc1xx\x99\xc9\xbcyy\x937\xf3\xfb\xed\xe8md\x18\xb3\xf0\xc9\x0c}\xd7"\x90\x90)a\xa1\xc5\xe8.\xcdX \xd2&\x8cPv\xb9\xf2"\xb8\x8av\xac\xac\xd6\xa5L\x1ben\xf4\xee\xed\xf5x\x16,\xcd\xf5\x83\x1b\x9f\xb8\x19\xf8 \x1f%\xf2]\x82n\x10u`C\xa9M\xcc\xc5\xa6\x1b\x03\x91\n\x00\x95E\x00\x9c[\x00>;\x00\xcb\xa4\xe1* y\xc1\xbcP\x82\x97e\x8c\xfde\xc8\x1e\x97\xd4\xc1\xbd\x99\x97W\r\x1e\\\x84\x89\xcd\x9a\x94\xef%\x13 \xbc\xae\xd3=\xde\xf2z\xd1&\x84~5l\xa0@\x19\x97\x8a\xff\xdfo\xfa\xe9\x97\x1dJ^\x0ci\x16z(;\xa2\xfc\ xdc\'!\xb9PUU\xc8.\x8e\xef>\xe8\xa3\xa1\x88\x12pU\xb7\x9a\x9at\xd1\xf4+\xe0\xf5\x90\x9e\x00\x0c\x01\x8c\x01\xac\x00V\xa2\xee\x10g\xb6 <\xe6\x9a\xfei\x1bxH\x8b\xc1\xdd\xefT\xcc\xc5X\x7f\xd5\x98h\xfc\xf5\xa7\xc6\x85V\xbc\x03\xb7nn4\x10\x02\x00\ x00\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff\x06\x00BC\x02\x00\x1b\x00\x03\x00\x00\x00\x00\x00\x00\x00\ x00\x00', stderr=[main_samview] 随机对齐检索仅适用于索引 BAM 或 CRAM 文件。\n'

先感谢您！

我正在尝试将我的工作流程从 bwa 升级到 bwa-mem2。该命令通常有效。

但是，一旦我开始使用 bwa-mem2，管道就会中断。

我得到以下输出：

我也将 fasta 文件从常规 bwa 索引更改为 bwa-mem2 索引。

我想提取读取名称列表，其中一个是唯一映射的，另一个是多映射的。

但是，经过一些尝试和错误，我发现

第一个生的

总是这样坏：

或者像这样：

（即部分 sam 文件行）

我怀疑bam文件本身坏了，然后我做

但它正常返回

是否有任何理由插入“句末”或其他工作人员？

先感谢您，

[修改问题]

感谢您的友好指导和解释。

好的，我重新解释一下情况。由于原始文件很重，我无法掌握许多工具所做的整个过程，所以让我解释一下。

“samtools view”命令采用 .bam 或 .sam 文件并返回超过 11 个以 \t 作为分隔符的字段。

我想用“XS:i:”提取行，所以“grep XS:i:”和 awk '/XS:i:/' 都可以。

尽管预计每个完整的行都将输出，但已经出现了具有完整字段集的行。

从另一个试验中，我可以理解该行没有损坏，我想知道这是因为 samtools 或 grep 和/或 awk。

是否可以通过 grep 或 awk '/' 部分提取太长的行？? 谢谢大家，我知道 grep -q 不适合我的目的。

我正在尝试使用 for 循环从包含可能字符串列表的文件中提取包含字符串的文件的所有行。我还想将 grep 的结果导出到文件名中包含变量的新文件。

这是我所拥有的：

这段代码所做的只是为每个变量创建一个空白文件。为什么 grep 不提取包含该变量的行并将其放入这些文件中？

作为参考，variables.txt文件如下所示：

这是samtools view输出的样子：

对于那些可能不熟悉的人，samtools view只需读出.bam文件即可。你可以把它想象成cat.

提前致谢！

我正在尝试编写一个循环遍历目录中的 SAM 文件并用于samtools view将这些 SAM 文件转换为 BAM 文件的 python 脚本。这就是我所拥有的，但我正在努力解决如何将 SAM 文件 (i) 输入到最后一行的 bash 命令中。

感谢您的帮助！

我通过 velocyto.py 将 10X bam 文件传输到 loom 文件中，并使用 Scanpy 进行细胞聚类。但是，当我用 Scanpy 进行数据处理时，我发现 loom 文件中没有 mito 基因，见下文。 在此处输入图像描述

运行 velocyto.py 时，它显示：“警告 - .bam 文件引用注释 (.gtf) 文件中不存在的染色体 'MT+'” “警告 - .bam 文件引用染色体 'MT-' 不存在于注释 (.gtf) 文件中”

BAM 文件来自 Cell Ranger v3.1 生成的 10X。参考基因组是 10X 的 refdata-gex-mm10-2020-A.tar.gz。重复注释文件是来自 UCSC 的 mm10_rmsk.gtf。

是不是因为 cell ranger v3.1 太旧，没有与 refdata-gex-mm10-2020-A.tar.gz 和 mm10_rmsk.gtf 相同的染色体名称？但是当我检查BAM中的染色体名称时，我发现它是ChrM，与refdata-gex-mm10-2020-A.tar.gz和mm10_rmsk.gtf一样？那么为什么提示显示 .bam 文件指的是染色体“MT-”？

这是检查 bam 中染色体名称的代码。$ samtools 查看/home/hyjforesight/mybam.bam | 切-f3 | perl -ne '打印除非 $seen{$_}++' - chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrM chrX chrY

我想可视化具有多个基因组的 ba​​m 文件。我使用 samtools 对其进行了排序，并得到了一个索引文件。

我首先尝试了 Qualimap，但是当生成跨参考的覆盖范围时，结果并不好，因为值在轴末端重叠。

我需要获得更好的可视化。我使用了各种工具，IGV、GenomeWorkbench 甚至 ngscat，但其中很多工具不接受我的 bam 文件，因为它有几种/多种类型的基因组（病毒和人类）。我需要一个可以接受它并给我一个更清晰的图表的工具

谢谢！

我正在使用“samtools quietd”将 MD 标签添加回 BAM 文件。原始 BAM 的大小约为 50Gb（使用 pacbio HIFI 读取的全基因组序列）。我遇到的问题是“冷静”的速度非常慢！这些作业已经运行了 12 个小时，并且只生成了 600MB 带有 MD 标签的 BAM。这样，50GB BAM 需要 30 天才能完成！

这是我用来添加MD标签的代码（很正常）：

我使用的samtools版本是v1.10。

顺便说一句，我使用 16 个核心来运行平静，但是，看起来 samtools 仍在使用 1 个核心来运行它：

我可以知道如何让平静更快吗？或者有没有其他工具可以更有效地完成同样的工作？

非常感谢

我正在尝试创建用于分析生物数据的自动命令链。
为此，我在 Slurm 集群中使用 Samtools。下面这一行是我为分析运行的命令之一：
samtools view -h file.sam | awk '$6 ~ /N/ || $1 ~ /^@/' | samtools view -h > spliced.file.sam
使用它，我得到了预期的输出（简单）。
但是，当我想将命令插入作业时，--wrap会出现语法错误。
如图所示：

srun在命令的开头和结尾使用&，在提交时非常有用，但是当我想创建命令管道时可以使用它吗？我可以为这个命令添加依赖项吗？有没有可能使用
--wrap这个命令的方法？

我的目标是创建一个自动命令管道，如下面的链接所示 - https://gencore.bio.nyu.edu/building-an-analysis-pipeline-for-hpc-using-python/

提前致谢。