我有一个包含起点和终点的区域列表。
我使用了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的，可以选择--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
只想提一下，如果遇到问题，您可能需要将 samtools 更新到最新版本。就我而言，samtools V1.2 不起作用，而 V1.10 起作用。