使用 R 根据 FASTA 文件中名称的初始段使它们唯一

Question

我的目标是根据 FASTA 文件中名称初始部分的唯一性来排除序列。但是，我在创建并使用下面的代码时遇到了错误。

FASTA 文件通常包含 4,333 个序列，我的目标是每个生物体仅保留一个序列。例如，如果智人有三个序列，我打算只选择一个。我利用下面的代码根据初始部分使名称变得唯一，但它未能实现所需的唯一性。您能指导我如何增强代码吗？非常感谢您的帮助

>Homo_sapiens_NP_004698.3_protein ATG12 isoform 1
MAEEPQSVLQLPTSIAAGGEGLTDVSPETTTPEPPSSAAVSPGTEEPAGDTKKKIDILLKAVGDTPIMKTKKWAVERTRT
IQGLIDFIKKFLKLVASEQLFIYVNQSFAPSPDQEVGTLYECFGSDGKLVLHYCKSQAWG
>Homo_sapiens_AAH12266.2_autophagy related 12 homolog (S. cerevisiae)
MTSREHQVSLCNCVPLLRRLLCDAPWRKARPLHALSRYFRSRVSPSKMAEEPQSVLQLPTSIAAGGEGLTDVSPETTTPE
PPSSAAVSPGTEEPAGDTKKKIDILLKAVGDTPIMKTKKWAVERTRTIQGLIDFIKKFLKLVASEQLFIYVNQSFAPSPD
QEVGTLYECFGSDGKLVLHYCKSQAWG
>Homo_sapiens_BAG36125.1_protein product
MAEEPQSVLQLPTSIAAGGEGLTDVSPETTTPEPPSSAAVSPGTEEPAGDTKKKIDILLKAVGVTPIMKTKKWAVERTRT
IQGLIDFIKKFLKLVASEQLFIYVNQSFAPSPDQEVGTLYECFGSDGKLVLHYCKSQAWG
>Pan_troglodytes_PNI35014.1_protein CK820_G0038523
MAEEPQSVLQLPPSSAAGGEGLTDVSPETTTPEPPSSAAVSPGTEEPAGDTKKKIDILLKAVGDTPIMKTKKWAVERTRT
IQGLIDFIKKFLKLVASEQLFIYVNQSFAPSPDQEVGTLYECFGSDGKLVLHYCKSQAWG
>Pan_paniscus_XP_034816260.1_protein ATG12 isoform X1
MTSREHQVSLCNCVPLLRRLLCDAPWRKARPLHALSRYFRSRVFPSKMAEEPQSVLQLPPSSAAGGEGLTDVSPETTTPE
PPSSAAVSPGTEEPAGDTKKKIDILLKAVGDTPIMKTKKWAVERTRTIQGLIDFIKKFLKLVASEQLFIYVNQSFAPSPD
QEVGTLYECFGSDGKLVLHYCKSQAWG

>Pan_paniscus_XP_000000H_protein ATG12 isoform X1
MTSREHQVSLCNCVPLLRRLLCDAPWRKARPLHALSRYFRSRVFPSKMAEEPQSVLQLPPSSAAGGEGLTDVSPETTTPE
PPSSAAVSPGTEEPAGDTKKKIDILLKAVGDTPIMKTKKWAVERTRTIQGLIDFIKKFLKLVASEQLFIYVNQSFAPSPD
QEVGTLYECFGSDGKLVLHYCKSQAWG

#Code used
library(Biostrings)
library(stringr)

modified_sequences <- readAAStringSet( "path/modified_sequences.fasta")
modified_sequences

first_part <- sub("^>([^_]+)_.*$", "\\1", names(modified_sequences))
unique_first_part <- make.unique(first_part)
name_mapping <- setNames(unique_first_part, first_part)
keep_indices <- !duplicated(unique_first_part)
unique_sequences <- modified_sequences[keep_indices]
writeXStringSet(unique_sequences, file = "path/unique_sequences.fasta")
unique_sequences

Answer 1

使用您提供的示例数据，我无法重现实际错误，但代码确实无法返回唯一序列，即每个属一个序列（id 的第一部分）。

我注意到 R 中 AAStringSet 对象中的序列 id 不是以“">”开头，因此

sub

函数调用没有效果，因为它与指定的模式不匹配。为了确保这不是使用不同软件包版本的影响，这些是我使用的：

library(Biostrings)
library(stringr)

> packageVersion("Biostrings")
[1] ‘2.68.1’
> packageVersion("stringr")
[1] ‘1.5.0’

不过，您的代码不会生成唯一的名称，因为您正在使用

make.unique

以及随后的

!duplicated

。由于您首先使名称唯一，因此不会有任何重复的名称，因此选择所有序列。

这是代码的修改版本，它选择与属（“Homo”和“Pan”）相关的第一个序列。

library(Biostrings)
library(stringr)

modified_sequences <- readAAStringSet( "path/modified_sequences.fasta")
modified_sequences

first_part <- sub("^([^_]+)_.*$", "\\1", names(modified_sequences))

keep_indices <- !duplicated(first_part)

unique_sequences <- modified_sequences[keep_indices]
unique_sequences

writeXStringSet(unique_sequences, file = "path/unique_sequences.fasta")

保留的独特序列是：

AAStringSet object of length 2:
    width seq                                                  names               
[1]   140 MAEEPQSVLQLPTSIAAGGEGLTDV...GTLYECFGSDGKLVLHYCKSQAWG Homo_sapiens_NP_0...
[2]   140 MAEEPQSVLQLPPSSAAGGEGLTDV...GTLYECFGSDGKLVLHYCKSQAWG Pan_troglodytes_P...