bam转fq

常用步骤¶

踩坑日记|bam 文件转换为 fastq 文件 - 知乎 (zhihu.com)

1.排序¶

根据 Read Name 排序

为什么要排序，在后面的常见问题里进行了说明

samtools sort -n SAMPLE.bam -o SAMPLE_sorted.bam

补充 1：使用 Sambamba 替代 samtools 排序

更快的处理 bam 数据—Sambamba-腾讯云开发者社区-腾讯云 (tencent.com)

Sambamba 能多线程。就测试来看，对同一个文件排序，由原来的近 20m，缩短为了 2m25s

sambamba sort -n -t 线程数 SAMPLE.bam -o SAMPLE_sorted.bam

补充 2：如何判断是否已经排序

生信分析常见文件——BAM 文件-CSDN 博客

1.查看 bam 文件 samtools 直接查看 bam 文件_samtools 查看 bam 文件-CSDN 博客

samtools view -h 文件名 |less -S

2.判断：

方法一：从头部判断

方法二：如果头部没有找到 HD 标签，手动核对任一 Read Name 是否相邻
文件的第一列是 Read Name，双端测序时，相同的 Read Name 会出现两次
所以如果发现两个相同的 Read Name 没有挨在一块，就说明没有按 Read Name 排序过

2.转换¶

Bam 中提取出序列（即 Bam 转 fq）

samtools fastq -@ 8 SAMPLE_sorted.bam \
-1 SAMPLE_R1.fastq.gz -2 SAMPLE_R2.fastq.gz \
-0 /dev/null -s /dev/null -n

1. samtools fastq:

这是 samtools 工具集中的一个命令，用于将 SAM 或 BAM 格式的比对文件转换为 Fastq 格式的序列文件。

2. -@ 8:

-@: 指定线程数。
8: 使用 8 个线程进行转换。这可以加快转换速度，尤其是在处理大型 BAM 文件时。根据你的 CPU 核心数调整此值。建议设置为 CPU 核心数或略小于 CPU 核心数。

3. SAMPLE_sorted.bam:

这是输入的 BAM 文件名。
SAMPLE_sorted.bam: 假设你的 BAM 文件名为 SAMPLE_sorted.bam。
非常重要: 这个 BAM 文件应该是已经按 Read Name 排序的 (使用 samtools sort -n 命令)。否则，输出的 R1 和 R2 Fastq 文件中的 reads 可能无法正确配对。

4. -1 SAMPLE_R1.fastq.gz:

-1: 指定输出的 R1 (正向 reads) Fastq 文件的文件名。
SAMPLE_R1.fastq.gz: 输出的 R1 Fastq 文件名为 SAMPLE_R1.fastq.gz。 .gz 后缀表示输出的文件将被 gzip 压缩。

5. -2 SAMPLE_R2.fastq.gz:

-2: 指定输出的 R2 (反向 reads) Fastq 文件的文件名。
SAMPLE_R2.fastq.gz: 输出的 R2 Fastq 文件名为 SAMPLE_R2.fastq.gz。 .gz 后缀表示输出的文件将被 gzip 压缩。

6. -0 /dev/null:

-0: 指定输出的未比对上的 reads (unmapped reads) 的 Fastq 文件的文件名。
/dev/null: 这是一个特殊的设备文件，表示 "空设备"。任何写入 /dev/null 的数据都会被丢弃。
作用： 这条命令会将所有未比对上的 reads 丢弃，不保存到任何文件中。这可以减少输出文件的大小。
Hard Clipping 策略的产物，在“为什么可以 BAM 转 Fastq-潜在问题”中有具体说

7. -s /dev/null:

-s: 指定输出的单端 reads (single-end reads) 的 Fastq 文件的文件名。
/dev/null: 同样是 "空设备"。
作用： 这条命令会将所有单端 reads 丢弃，不保存到任何文件中。

8. -n:

-n: 根据 Read Name 输出 Fastq 文件的 Read ID。也就是保留 Read Name 的完整信息。
作用： 确保输出的 Fastq 文件中的 Read ID 与 BAM 文件中的 Read Name 保持一致，这对于后续的分析非常重要。

总结：

这条 samtools fastq 命令的作用是：

使用 8 个线程将按 Read Name 排序的 BAM 文件 SAMPLE_sorted.bam 转换为 Fastq 文件。
将正向 reads (R1) 保存到压缩的 Fastq 文件 SAMPLE_R1.fastq.gz 中。
将反向 reads (R2) 保存到压缩的 Fastq 文件 SAMPLE_R2.fastq.gz 中。
丢弃所有未比对上的 reads 和单端 reads。
确保输出的 Fastq 文件中的 Read ID 与 BAM 文件中的 Read Name 保持一致。

注意事项：

排序： 再次强调，输入的 BAM 文件必须是按 Read Name 排序的。
线程数： 根据你的 CPU 核心数调整 -@ 选项的值。
丢弃 reads： 确保你知道你在丢弃未比对上的 reads 和单端 reads。如果你后续的分析需要用到这些 reads，就不要使用 -0 /dev/null 和 -s /dev/null 选项。
输出文件名： 根据你的需求修改输出文件名。

snakemake 流程化¶

详情见snakemake 流程化 - bam 转 fq

常见问题¶

1.为什么要排序¶

理论解释¶

为什么将 BAM 文件转换为 Fastq 需要先按 Read Name 排序？

根本原因在于，双端测序的 reads (R1 和 R2) 在 BAM 文件中可能不是成对相邻存放的，尤其是当 BAM 文件是按染色体坐标排序的时候（详见后面的例子）
所以转换为 Fastq 时，我们需要确保每个 read 的 R1 和 R2 正确配对，以便后续的分析。

以下是更详细的解释：

1.BAM 文件的排序方式：

按坐标排序 (Coordinate-sorted): 这是 BAM 文件的常见排序方式。 BAM 文件中的 reads 按照它们在参考基因组上的位置 (染色体和坐标) 进行排序。这种排序方式适合于基因组比对、变异检测等需要考虑 reads 在基因组上位置的分析。

按 Read Name 排序 (Name-sorted): 在这种排序方式中，BAM 文件中的 reads 按照它们的 Read Name 进行排序。对于双端测序数据，具有相同 Read Name 的 R1 和 R2 reads 会相邻存放。

2.双端测序与 Fastq 文件:

双端测序生成两个文件，一个包含正向 reads (R1)，另一个包含反向 reads (R2)。 Fastq 文件的格式是 reads 成对出现，通常按照原始测序的顺序排列。每个 read 对应四行： Read Name, 序列, 占位符 (+), 质量值。

对于 paired-end reads, R1 和 R2 必须在各自的 Fastq 文件中保持顺序一致，并且 read name 要对应。

3.转换过程中的配对问题:

当 BAM 文件按坐标排序时，属于同一对的 R1 和 R2 reads 可能在 BAM 文件中相隔很远，甚至分布在不同的染色体上（如果发生了嵌合）。

如果直接将按坐标排序的 BAM 文件转换为 Fastq，无法保证 R1 和 R2 reads 在 Fastq 文件中的顺序是配对的，这会严重影响下游分析的准确性。

4.Read Name 排序的作用:

通过 samtools sort -n 将 BAM 文件按 Read Name 排序后，确保了具有相同 Read Name 的 R1 和 R2 reads 在 BAM 文件中是相邻的。

这样，在后续使用 bedtools bamtofastq 或 samtools fastq 转换时，就可正确地将成对的 reads 提取到对应的 R1 和 R2 Fastq 文件中，保证了 paired-end 信息的完整性。

总结：

先按 Read Name 排序是为了保证在 BAM 文件转换为 Fastq 文件时，双端测序的 reads (R1 和 R2) 能够正确配对，从而避免下游分析出现错误。
如果 BAM 文件已经是按 Read Name 排序的，则可以跳过 samtools sort -n 这一步。但是，如果不能确定 BAM 文件的排序方式，最好还是先进行排序，以确保结果的正确性。

补充说明：

bedtools bamtofastq 和 samtools fastq 都是用于将 BAM 文件转换为 Fastq 文件的工具。 samtools fastq 功能更强大，可以直接生成 gzip 压缩的 Fastq 文件，并且可以丢弃未配对的 reads。

例子辅助理解¶

简化假设：

我们只有 3 条 reads。
这些 reads 来自一个简单的基因组，只有一条染色体 (chr1)，长度为 100bp。
这些 reads 是双端测序的，因此每条 read 都有 R1 和 R2。

BAM 文件（按坐标排序）：

samtools 看一下是否排序了

一个按坐标排序的 BAM 文件大致会按照 染色体 和 比对位置 进行排序。

Read ID | Read Name | 染色体 | 比对位置
------- | -------- | -------- | --------
1       | read_1/1    | chr1     | 10
3       | read_2/1    | chr1     | 20
5       | read_3/1    | chr1     | 30
2       | read_1/2    | chr1     | 60
4       | read_2/2    | chr1     | 70
6       | read_3/2    | chr1     | 80

这里：

Read Name 是 read 的名称（例如 read_1）。双端测序 reads 的 Read Name 相同，用于标识它们是一对。
那么怎么区分哪个是正向 R，哪个是反向 R？
Read ID 是一个唯一的标识符，用于区分不同的 reads (方便我们理解)。
染色体 是 read 比对到的染色体。
比对位置 是 read 比对到染色体上的起始位置。

可以看到：

Reads 按照比对位置升序排列。
同一对 reads (例如 read_1 的 R1 和 R2) 在文件中不是相邻的，它们被其他 reads 隔开。
BAM 文件还包含序列、质量值等信息，这里为了简化而省略。

BAM 文件（按 Read Name 排序）：

一个按 Read Name 排序的 BAM 文件会按照 Read Name 进行排序。在这个例子中，BAM 文件中的 reads 的顺序可能是这样的：

Read ID | Read Name | 染色体 | 比对位置
------- | -------- | -------- | --------
1       | read_1/1    | chr1     | 10
2       | read_1/2    | chr1     | 60
3       | read_2/1    | chr1     | 20
4       | read_2/2    | chr1     | 70
5       | read_3/1    | chr1     | 30
6       | read_3/2    | chr1     | 80

可以看到：

Reads 按照 Read Name 排序。
同一对 reads (例如 read_1 的 R1 和 R2) 在文件中是相邻的。

为什么按 Read Name 排序后更容易转换到 Fastq？

现在假设我们要将这些 BAM 文件转换为两个 Fastq 文件 (R1.fastq 和 R2.fastq)。

按坐标排序的 BAM： 如果直接转换，工具需要记住每个 read 的 Read Name，并找到与之配对的 read。这需要大量的内存和计算资源，而且容易出错。
按 Read Name 排序的 BAM： 工具可以简单地按顺序读取 BAM 文件，将每对相邻的 reads 分别写入 R1.fastq 和 R2.fastq。由于同一对 reads 是相邻的，所以配对非常简单高效。

实际 BAM 文件：

实际的 BAM 文件远比这个例子复杂。它包含大量的 reads，以及各种 metadata 信息 (例如比对质量、mapping flags 等)。 BAM 文件通常是二进制格式，不能直接用文本编辑器打开。你需要使用 samtools 等工具来查看 BAM 文件的内容。例如：

samtools view -H input.bam  ## 查看 BAM 文件的 header
samtools view input.bam | head -n 20 ## 查看 BAM 文件的前 20 行 (文本格式)

补充：

虽然双端测序 reads 的 Read Name 相同，但通常会在 Read Name 的末尾添加一个后缀来区分 R1 和 R2。常用的后缀包括：

1./1 和 /2 后缀：这是最常见的区分方式。

R1 的 Read Name 结尾是 /1。 R2 的 Read Name 结尾是 /2。

例如：read_1/1 和 read_1/2。

2..1 和 .2 后缀：有些测序平台或数据处理流程会使用 .1 和 .2 作为后缀。

R1 的 Read Name 结尾是 .1。 R2 的 Read Name 结尾是 .2。

例如：read_1.1 和 read_1.2。

3._1 和 _2 后缀：也有些情况使用 _1 和 _2 作为后缀。

R1 的 Read Name 结尾是 _1。 R2 的 Read Name 结尾是 _2。

例如：read_1_1 和 read_1_2。

4.Read Flag (BAM 文件特有):

在 BAM 文件中，有一个叫做 "flag" 的字段，其中包含了 reads 的各种信息，包括是否是 paired-end read、是 R1 还是 R2 等。这个 flag 信息是区分 R1 和 R2 的最可靠的方式。

工具如何利用这些信息：

bedtools bamtofastq 和 samtools fastq 在转换 BAM 文件时，会读取 Read Name 或 Read Flag，根据后缀或 flag 信息来判断哪个 read 是 R1，哪个 read 是 R2，并将它们分别写入对应的 Fastq 文件。

重要：检查你的数据！

你需要检查你的 BAM 文件，确定 Read Name 的后缀是什么 (是 /1 和 /2，还是 .1 和 .2，还是其他？)。
如果 Read Name 没有后缀，那么你需要依赖 Read Flag 来区分 R1 和 R2 (这需要更高级的 BAM 文件处理技巧)。

2.为什么可以 Bam 转 Fastq¶

Fastq (通常): 通常是从 SRA 文件转换而来，包含原始 reads 序列和质量信息，未经比对和排序。
BAM: 是 Fastq 文件经过比对和处理后的结果文件，包含了 reads 的比对信息。

既然 BAM 已经是比对后的数据，为什么要转回 Fastq 呢？

答案是：可以这样做，而且在某些情况下是合理的。

为什么可以将 BAM 转回 Fastq？¶

1.数据可移植性和互操作性:

Fastq 是一种非常通用的格式，几乎所有的基因组学工具都支持 Fastq 作为输入。

虽然 BAM 文件包含了比对信息，但有时候你可能只需要 reads 的序列信息，而不需要比对信息。将 BAM 转回 Fastq 可以方便地将数据提供给那些只接受 Fastq 作为输入的工具。

2.修改比对或重新比对:

比对是一个计算密集型的过程。如果你对比对结果不满意 (例如，使用了错误的参考基因组、使用了不合适的比对参数)，你可以将 BAM 转回 Fastq，然后使用不同的参数或工具重新进行比对。

例如，你可能需要使用更新的参考基因组版本，或者需要使用更敏感的比对参数来检测更多的变异。

3.数据子集提取:

有时候你只需要 BAM 文件中一部分 reads 的序列信息。例如，你可能只对未比对上的 reads 感兴趣 (unmapped reads)。

将 BAM 转回 Fastq，然后只保留未比对上的 reads，可以减小数据量，提高后续分析的效率。

4.数据共享:

虽然 BAM 文件包含了比对信息，但共享 BAM 文件可能涉及一些隐私问题，因为 BAM 文件包含了 reads 在基因组上的精确位置。

将 BAM 转回 Fastq 可以去除比对信息，从而在一定程度上保护了数据的隐私。当然，这种方法并不能完全匿名化数据，因为 reads 的序列信息仍然可以用于推断样本的来源。

5.存储成本考虑:

虽然 BAM 通常是压缩格式，但如果原始数据只需要序列信息，fastq 可能更节省空间。尤其对于需要长期存档的数据，可能需要权衡考虑。但是，这种考虑现在越来越少，因为存储成本持续下降。

BAM 转回 Fastq 的潜在问题¶

丢失比对信息： 将 BAM 转回 Fastq 会丢失比对信息【这是必然的，因为我们只需要序列信息】，包括比对位置、比对质量、CIGAR 字符串等。如果你后续的分析需要用到这些比对信息，那么这种转换就不可取。
Read Name 信息： BAM 转回 Fastq 时，需要确保 Read Name 信息能够正确地保留，并且 R1 和 R2 能够正确配对。不同的转换工具可能会有不同的处理方式。
计算开销： 虽然比对是一个计算密集型的过程，但将 BAM 转回 Fastq 也需要一定的计算开销。

关键点：这种转换会不会丢失序列信息？？

答案是： 理想情况下，将 BAM 转回 Fastq **不应该**丢失序列信息。 samtools fastq 等工具的设计目标就是从 BAM 文件中提取出原始的 reads 序列和质量值，并将其恢复为 Fastq 格式。
但是，在实际操作中，可能会发生以下情况导致序列信息丢失或改变：

1.Hard Clipping：

Hard clipping 是一种比对策略，它会将 reads 的部分序列从比对结果中 "剪掉" (移除)，通常是因为这些序列无法比对到参考基因组上。 BAM 文件会记录 hard clipping 的信息。

问题： 如果在 BAM 文件中存在 hard clipping，那么在转回 Fastq 时，被 hard clipping 掉的序列信息可能无法完全恢复。不同的工具对 hard clipping 的处理方式可能不同。有些工具可能会直接丢弃被 hard clipping 掉的序列，导致序列长度变短。

解决方案： 在进行 BAM 转 Fastq 之前，最好检查 BAM 文件中是否存在大量的 hard clipping。如果存在，需要考虑是否去除 hard clipping 信息，或者选择能够正确处理 hard clipping 的转换工具。

2.质量值编码：

不同的测序平台可能使用不同的质量值编码方式 (例如，Phred33, Phred64)。BAM 文件会记录质量值编码方式。

问题： 如果转换工具无法正确识别 BAM 文件中的质量值编码方式，那么在转回 Fastq 时，质量值可能会被错误地转换，导致质量信息失真。

解决方案： 确保转换工具能够正确识别 BAM 文件中的质量值编码方式。如果需要，可以手动指定质量值编码方式。

3.软件 Bug：

虽然 samtools fastq 等工具是经过广泛测试的，但仍然可能存在 bug，导致序列信息丢失或改变。

解决方案： 使用最新版本的软件，并仔细阅读软件文档，了解其已知的问题和限制。

4.数据损坏：

罕见情况下，BAM 文件本身可能已经损坏，导致无法正确提取序列信息。

总结：

虽然理论上将 BAM 转回 Fastq 不应该丢失序列信息，但在实际操作中，需要注意 hard clipping、质量值编码等问题，并选择合适的工具和参数，以确保转换的准确性。

最好的方法，是在转换前后对数据进行验证，例如：

比较序列总数： 比较原始 Fastq 文件和从 BAM 文件转回的 Fastq 文件中的 reads 总数，看是否一致。
比较序列长度： 检查从 BAM 文件转回的 Fastq 文件中的 reads 长度是否与原始 Fastq 文件一致。
随机抽样比对： 从原始 Fastq 文件和从 BAM 文件转回的 Fastq 文件中随机抽取一部分 reads，然后将它们比对到参考基因组上，比较比对结果是否一致。

通过这些验证步骤，可以确保将 BAM 转回 Fastq 的过程没有引入错误。记住，数据验证是基因组学分析中非常重要的一环！永远不要盲目相信工具的结果，一定要进行验证。

如何将 BAM 转回 Fastq？¶

可以使用 samtools fastq 命令将 BAM 文件转换为 Fastq 文件：

samtools fastq input.bam -1 output_R1.fastq -2 output_R2.fastq

input.bam: 输入的 BAM 文件。
output_R1.fastq: 输出的 R1 Fastq 文件。
output_R2.fastq: 输出的 R2 Fastq 文件。

总结：

将 BAM 转回 Fastq 是一个可行的操作，并且在某些情况下是合理的。但是，需要权衡利弊，并确保转换过程中不会丢失重要信息。

【 只有当你确定后续的分析不需要比对信息时，才可以考虑将 BAM 转回 Fastq】

还需要注意 Read Name 信息的正确保留和 R1/R2 的正确配对。考虑到存储成本，这种做法已经比较少见了。

本文阅读量次
本站总访问量次

Authors: wangshangjian