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

Question

我有一个带有起始点和结束点的区域列表。

我使用了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)

Answer 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

Answer 2

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