问题标签 [samtools]

我在 Linux 上使用 IntelliJ IDEA 2017.1.5，无法使用 IDE 编译 htsjdk。也许有人已经想出办法了？

需要明确的是，htsjdk 从命令行编译得很好，而不是在 IDE 中。

要编译，我单击 Gradle 选项卡中的小回收图标（“刷新外部项目”）。经过一些工作，它给了我：“编译完成，出现 8 个错误”。

构建选项卡显示此错误：

错误:(3, 28) java: 包 org.scalatest.testng 不存在

我之前没有将 Scala 与 IntelliJ 一起使用，所以我不确定如何继续。

我有一个 RNA-seq bam 文件，但很少有读数让我感到困惑。

根据bam头，这个bam文件是按坐标排序的，使用tophat创建的，markduplicate步骤没有做。但是一些读取在 samflag 中被标记为重复。更糟糕的是，当我运行 picard markduplicate 时，这些读取的 pcr 重复标志被切换，将它们标记为不重复。此外，我手动找到了此读取的副本（具有相同起始位置和配对起始位置的相同读取），因此初始标记看起来是正确的。

所以我的问题是：
知道为什么会发生这种情况吗？
Tophat 是否标记为重复？（我不这么认为）
如果读取已经被标记为重复，那么 picard markduplicate 是否会切换？

以下是标记重复步骤之前和之后读取的外观。
Before:
C0RTF 1187 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

After Markduplicate
C0RTF 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

谢谢

我希望在 awk 的每一行上打印 $fname 的以下输出：

电流输出：

期望的输出：

这样输出的平均深度是可追溯的。我已经尝试了 FILENAME 和各种各样的东西，但没有成功。谢谢！

这里的关键是这是一个巨大的文件。我的目标是避免一次将整个文件读入内存，并避免解析循环中的每一行以到达我需要的行（因为它需要很长时间。该文件实际上有 1500 万行长）。

我目前正在做的是将文件打开为...

...将指针直接移动到所需行的位置（使用许多恶作剧，但它有效）并在单独的行中读取。

类似于下面的行（尽管这些示例来自文件的开头，为了方便和信息起见）......

有些人可能会注意到这是一个BAM文件，所以如果有更好的方法可以做到这一点，欢迎提出建议……尽管samtools过滤器无法满足我的需要。我必须按行搜索，而不是按数据搜索。

我正在处理一个双端全基因组测序的 bamfile，并且想要过滤掉来自特定基因组区域的未映射到正确配对中的读数（这些有时表示结构变异）。我正在使用 samtools，并尝试使用“标志”选项过滤读取，以选择未映射成正确对的读取。如果我是正确的，这些读数的标志值不应该有 2。（https://broadinstitute.github.io/picard/explain-flags.html）

但是，根据 samtools，我的所有读取都没有映射成正确的对。当我计算（-c）我指定的区域中的所有读取时，没有过滤器，它给我的总数为 179：

当我过滤正确配对的读取时，即标志包含“2”（-f 2），计数为零：

我检查了读取是否被识别为配对（-f 1），以及配对是否被映射（-F 8），它们都是：

我还尝试了基因组中的其他区域，并且到处遇到同样的问题。我使用相同的 BAM 文件检查了 IGV 中的区域，IGV 告诉我大多数读取都正确配对，只有少数不是。有谁知道这里发生了什么？BAM 文件是否以不考虑正确配对映射的方式标记？

欢迎任何帮助！非常感谢。

我正在尝试使用samtools两个文件 File1 和 File2 来制作堆积文件。

我已按染色体将 File1 和 File2 分开，导致有 44 个文件以下列格式命名：

其中 ${c} 是 1 到 22 之间的一个数字，$TISSUE 是结肠或肌肉——22 条染色体代表结肠，22 条染色体代表肌肉。IE; chr1.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY chr2.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY

这些文件由两列组成，第一列只显示染色体编号，第二列是该染色体上的位置。IE;

对于文件中的每一行（例如"chr2.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY"），我需要从第 2 列中获取位置，将其称为“x”，并使用它来获取a-b、 wherea=x-5和的范围b=x+5。然后，我会将这些值插入以下脚本：

例如，假设我正在查看 2 号染色体，位置 103977（上面的第 1 行）。那么我的脚本将是

所以基本上，它是一个循环中的一个循环。就像是，

提前致谢。我是使用 Linux 的新手，而且我基本上没有受过计算机科学培训。

我正在尝试限制并行化脚本。该脚本的目的是获取 10 个样本/文件夹内的列表，并使用列表的记录执行 samtools 命令，这是最苛刻的部分。

这是简单的版本：

为了使用我们的本地服务器，该脚本包含一个可以工作的分叉命令。但它会分叉，直到服务器的所有资源都用完并且没有其他人可以处理它。

因此，我想parallel -j 50用 gnu 并行实现类似的东西。我在待分叉的samtools命令前面试了一下，比如

这不起作用（也尝试使用反引号），我得到了

或者以某种方式甚至vim被调用。但我也不确定这是否是parallel脚本中命令的正确位置。您是否知道如何解决此问题，从而限制分叉的进程数量？

我还想过使用基于 FIFO 的信号量实现这里提到的https://unix.stackexchange.com/questions/103920/parallelize-a-bash-for-loop/103921，但我希望gnu parallel能做我正在寻找的为了？我检查了更多页面，例如https://zvfak.blogspot.de/2012/02/samtools-in-parallel.html和https://davetang.org/muse/2013/11/18/using-gnu-parallel/但通常不是这种组合问题。

这是脚本的更详细版本，以防其中的任何命令可能相关（我听说 awk 反引号和新行通常可能是一个问题？）

我想过滤一个文件，以便我可以获得在第 1 列中匹配但在第 2 列中不匹配的行。在以下示例中：

我想在第 1 列中找到恰好出现两次的整体，并且在第 2 列中不匹配，即应该忽略组合 column1、column2 中的重复项。所以预期的输出是：

第 3、4、5 等列中的内容对于过滤并不重要，但我确实需要保留这些信息。

我还需要从另一个输出中将其输入，这是读取文件并保留标题所必需的。所以我需要一些格式：

我尝试了几个版本的awk，但都无济于事。任何帮助将非常感激。

我正在尝试使用已下载到 main.cpp 文件文件夹中的samtools C API ( https://github.com/samtools/samtools ) 编译（在 Linux 上，使用 G++）一个简单的 main.cpp 程序。我想要一个非常简单的makefile来编译main.cpp（并最终编译samtools代码）。但是，由于我对 makefile 知之甚少，我可能做错了什么。

这是我的生成文件：

这是我的 cpp 主要内容：

当我运行“make”时，它编译时没有警告，但在运行时，它告诉我：“加载共享库时出错：libhts.so.2 无法打开共享对象文件”

任何帮助都非常受欢迎！提前致谢。

我有一个包含起点和终点的区域列表。

我使用了samtools faidx ref.fa <region>命令。这个命令给了我那个区域的正链序列。

在 samtools 手册中有一个提取反向链的选项，但我不知道如何使用它。

有人知道如何在 samtools 中为反向链运行此命令吗？

我的地区是这样的：