如何消除fasta文件中的重复序列

Question

我试图用所有发布的序列来构建细菌数据库类型，使用bowtie2来计算我的读取数据的覆盖率，并利用fasta_library进行映射。为此，我将从ncbi下载的所有基因组序列合并到一个fasta_library中(我在fasta文件中合并了74个文件)，问题是在这个fasta文件(我创建的库)中，我有很多重复的序列，这对覆盖范围有很大的影响，所以我想问是否有任何方法来消除我的Library_File中的重复，或者是否有任何方法来合并没有重复的序列，或者还有其他方法来计算我对引用序列的读取覆盖率

我希望我说的够清楚了，如果有什么不清楚的地方请告诉我。

Answer 1

如果您控制了您的设置，那么您可以安装塞克吉并在FASTA文件上运行以下命令：

$ seqkit rmdup -s < in.fa > out.fa

如果您有多个文件，您可以将它们连接起来并作为标准输入输入它们：

$ seqkit rmdup -s < <(cat inA.fa ... inN.fa) > out.fa

rmdup选项移除重复项，-s选项根据序列调用重复项，忽略标头中的差异。我不确定输出中保留了哪个标头，但这可能是要考虑的问题。

为了避免第三方依赖，并了解如何删除dups，可以使用awk。

这样做的目的是将所有FASTA记录逐个读入一个关联数组(或者哈希表，在Python中也称为“字典”)，前提是序列还没有在数组中。

例如，从如下所示的单行FASTA文件in.fa开始：

>test1
ATAT
>test2
CGCG
>test3
ATAT
>test4
GCCT

我们可以删除重复项，保留第一个标头，如下所示：

$ awk 'BEGIN {i = 1;} { if ($1 ~ /^>/) { tmp = h[i]; h[i] = $1; } else if (!a[$1]) { s[i] = $1; a[$1] = "1"; i++; } else { h[i] = tmp; } } END { for (j = 1; j < i; j++) { print h[j]; print s[j]; } }' < in.fa > out.fa
$ cat out.fa
>test1
ATAT
>test2
CGCG
>test4
GCCT

如果您需要修改，它需要一些关于awk的知识。这种方法还取决于您的FASTA文件是如何构造的(包含一行或多行序列的记录，等等)，尽管将FASTA文件修改为上述结构非常容易(头和序列各有一行)。

任何哈希表方法都会使用相当大的内存(我猜想seqkit可能在这个特定的任务上做出了同样的妥协，但我还没有查看源代码)。对于非常大的FASTA文件来说，这可能是一个问题。

如果您有一个可以安装软件的本地环境，那么最好使用seqkit。如果您有一个it锁定设置，那么awk也会为这个任务工作，因为它是随大多数Unixes一起出现的。