从fasta文件中删除多个序列

Question

我有一个由两行组成的字符序列的文本文件：一个标题，以及下一行中的序列本身。该文件的结构如下：

>header1
aaaaaaaaa
>header2
bbbbbbbbbbb
>header3
aaabbbaaaa
[...]
>headerN
aaabbaabaa

在另一个文件中，我有一个我想删除的序列标题列表，如下所示：

>header1
>header5
>header12
[...]
>header145

我们的想法是从第一个文件中删除这些序列，因此所有这些标题+以下行。我使用sed做了以下，

while read line; do sed -i "/$line/,+1d" first_file.txt; done < second_file.txt

它工作但需要很长时间，因为我用sed多次加载整个文件，而且它非常大。关于如何加快这个过程的任何想法？

Answer 1

使用第二个文件中的delete命令创建脚本：

sed 's#\(.*\)#/\1/,+1d#' secondFile.txt > commands.sed

然后将该文件应用于第一个文件

sed -f commands.sed firstFile.txt 

Answer 2

$ awk 'NR==FNR{a[$0];next} $0 in a{c=2} !(c&&c--)' list file
>header2
bbbbbbbbbbb
>header3
aaabbbaaaa
[...]
>headerN
aaabbaabaa

c是你想跳过多少行，从刚刚匹配的那一行开始。见https://stackoverflow.com/a/17914105/1745001。

或者：

$ awk 'NR==FNR{a[$0];next} /^>/{f=($0 in a ? 1 : 0)} !f' list file
>header2
bbbbbbbbbbb
>header3
aaabbbaaaa
[...]
>headerN
aaabbaabaa

f是否在目标数组>...中找到最近读取的a[]线。 f=($0 in a ? 1 : 0)可以缩写为f=($0 in a)，但为了清晰起见，我更喜欢三元表达式。

第一个脚本依赖于您知道每个记录长多少行，而第二个依赖于每个记录以>开头的记录。如果您知道哪一个，那么您使用的是样式选择。

Answer 3

你可以使用这个awk：

awk 'NR == FNR{seen[$0]; next} /^>/{p = !($0 in seen)} p' hdr.txt details.txt

Answer 4

您遇到的问题很容易回答，但在处理通用fasta文件时无法帮助您。 Fasta文件有一个序列标题，后跟一行或多行，可以连接起来表示序列。 Fasta文件格式大致遵循以下规则：

• 描述行（defline）或标题/标识符行以<greater-then>字符（>）开头，为序列提供名称和/或唯一标识符，还可能包含其他信息。
• 在描述行之后是标准单字母字符串中的实际序列。除了有效字符之外的任何内容都将被忽略（包括空格，制表符，星号等）。
• 序列可以跨越多行。
• 通过在公共文件中连接几个单序列FASTA文件，通常通过在两个后续序列之间留下空行来获得多序列FASTA格式。

大多数提出的方法将在具有多线序列的多快速方法上失败

以下将始终有效：

awk '(NR==FNR) { toRemove[$1]; next }
     /^>/ { p=1; for(h in toRemove) if ( h ~ $0) p=0 }
    p' headers.txt file.fasta

这与EdMorton和Anubahuva的答案非常相似，但这里的区别是文件headers.txt只能包含标题的一部分。

Answer 5

这个awk可能适合你：

awk 'FNR==NR{a[$0]=1;next}a[$0]{getline;next}1' input2 input1

Answer 6

一种选择是创建一个long sed表达式：

sedcmd=
while read line; do sedcmd+="/^$line\$/,+1d;"; done < second_file.txt
echo "sedcmd:$sedcmd"
sed $sedcmd first_file.txt

这只会读取一次文件。请注意，我将^和$添加到sed模式（所以>header1与>header123不匹配......）

---

如果您有数千个文件，使用文件（如@daniu所建议的）可能会更好，因为您可能会使用此方法命中命令行最大计数。

Answer 7

但尝试GNU;

sed -E ':s $!N;s/\n/\|/;ts ;s~.*~/&/\{N;d\}~' second_file.txt| sed -E -f -  first_file.txt

将time命令添加到两个脚本以比较速度， 看看time while read line;do...和time sed -....导致我的测试这是在不到OP的一半时间内完成的