使用 samtools 从 FASTA 文件的反向链中提取用户指定的序列答案

【问题标题】：Extract user specified sequence from reverse strand of from FASTA file Using samtools使用 samtools 从 FASTA 文件的反向链中提取用户指定的序列
【发布时间】：2020-03-03 08:19:41
【问题描述】：

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

我使用了samtools faidx ref.fa <region> 命令。该命令为我提供了该区域的正向链序列。

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

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

我的地区是这样的：

 LG2:124522-124572 (Forward)
 LG3:250022-250072 (Reverse)
 LG29:4822278-4822318 (Reverse)
 LG12:2,595,915-2,596,240 (Forward)
 LG16:5,405,500-5,405,828 (Reverse)

【问题讨论】：

标签： samtools

【解决方案1】：

如您所见，samtools 可以选择--reverse-complement（或-i）从反向链输出序列。

据我所知，samtools 不支持允许指定链的区域表示法。

一个快速的解决方案是将您的区域文件分成正向和反向位置并运行samtools 两次。

下面的步骤相当冗长，只是为了让步骤清晰。例如，使用 bash 中的进程替换来清理它是相当简单的。

# Separate the strand regions.

# Use grep and sed twice, or awk (below).
grep -F '(Forward)' regions.txt | sed 's/ (Forward)//' > forward-regions.txt
grep -F '(Reverse)' regions.txt | sed 's/ (Reverse)//' > reverse-regions.txt

# Above as an awk one-liner.
awk '{ strand=($2 == "(Forward)") ? "forward" : "reverse"; print $1 > strand"-regions.txt" }' regions.txt

# Run samtools, marking the strand as +/- in the FASTA output.
samtools faidx ref.fa -r forward-regions.txt --mark-strand sign -o forward-sequences.fa 
samtools faidx ref.fa -r reverse-regions.txt --mark-strand sign -o reverse-sequences.fa --reverse-complement

# Combine the FASTA output to a single file.
cat forward-sequences.fa reverse-sequences.fa > sequences.fa
rm forward-sequences.fa reverse-sequences.fa

【讨论】：

【解决方案2】：

只想提一下，如果遇到问题，您可能需要将 samtools 更新到最新版本。就我而言，samtools V1.2 不起作用，而 V1.10 起作用。

【讨论】：