将文本添加到 fasta 文件的标题答案

【问题标题】：add text to a header of a fasta file将文本添加到 fasta 文件的标题
【发布时间】：2021-07-18 22:04:08
【问题描述】：

我有一个基因组的 fasta 文件 (txt)，类似于：

$ cat Strain-01.faa
>IMEHDJCA_03186 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>IMEHDJCA_03189 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

我想根据 file.txt 中的列表添加一个额外的 ID。

$ cat file.txt
ID      Gene        Strain-01       Strain-02       Strain-03
ID_01   pphB        IMEHDJCA_03186  DIBHEKPI_01648  LLMDBGDK_00598
ID_02   group_1001  IMEHDJCA_03187  DIBHEKPI_01635  LLMDBGDK_00611
ID_03   group_1002  IMEHDJCA_03189  DIBHEKPI_01628  LLMDBGDK_00616

例如对于fasta Strain-01.faa文件有IMEHDJCA_03186 id对应Strain-01，所以我想将列ID的ID_01编号（来自file.txt）添加到标题序列，类似于：

ID_01 对应于IMEHDJCA_03186
ID_02 对应于IMEHDJCA_03187
ID_03 对应于IMEHDJCA_03189

结果会是这样的：

$cat Strain-01_edited.faa
>ID_01 IMEHDJCA_03186 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>ID_02 IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>ID_03 IMEHDJCA_03189 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

我只想在fasta文件的头部添加一个file.txt的ID码。

有什么想法吗？在bash 或R，或任何其他方式？

非常感谢

【问题讨论】：

bash 或任何其他 shell 对于此或任何其他文本操作来说都是一个糟糕的选择，请参阅why-is-using-a-shell-loop-to-process-text-considered-bad-practice。我希望您的意思是 Unix tool（例如 awk、sed 等）而不是 bash。

标签： r bash awk str-replace fasta

【解决方案1】：

更新 #2： 消除awk 处理的特定菌株（名称）：

我们会将所有可能的应变/ID 映射加载到awk
这将允许处理任何*.faa 文件，而无需知道菌株名称
这将允许处理带有混合菌株的 *.faa 文件（不知道这是否是 OP 必须解决的问题）
降低了 awk 代码的复杂性（与 UPDATE #1 相比），但会为更多 id[] 数组条目消耗额外的内存

样本数据（第一个字段中的菌株混合）：

# for this (nonsensical?) file the first 3 blocks include a strain
# from each of the 3 columns (of strain names) from file.txt; the
# 4th block contains a nonsensical strain that doesn't exist in
# file.txt (ie, 4th line should not see an insertion of a ID value)

$ cat Strain-mixed.faa
>IMEHDJCA_03186 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>DIBHEKPI_01635 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>LLMDBGDK_00616 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD
>NO_MATCH hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

新的awk 代码将所有菌株加载到id[] 数组中：

awk '
NR==1   { next }                           # skip 1st line of 1st file
FNR==NR { for (i=3; i<=NF; i++)            # for rest of 1st file load id[] with ...
              id[$i]=$1                    # all strain/ID combos
          next
        }
/^>/    {                                   # for 2nd file, if 1st column is ">"
          ndx=substr($1,2)                  # strip off ">"
          if ( ndx in id )                  # if 1st field (sans ">") is an index in id[] then ...
             ( $1=">" id[ndx] " " ndx )     # rewrite 1st field to include our id[] value
        }
1                                           # print current line (of 2nd file)
' file.txt Strain-mixed.faa

这会生成：

>ID_01 IMEHDJCA_03186 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>ID_02 DIBHEKPI_01635 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>ID_03 LLMDBGDK_00616 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD
>NO_MATCH hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

注意：此最新更新将为 Strain-{01,02}.faa 文件中的所有行执行 ID 插入（请参阅下面的 更新 #1）。

更新 #1： 扩展原始答案以解决（我认为）Paul Hodges 关于概括答案以支持可变菌株名称的问题：

从file.txt动态确定使用哪一列菌株
动态处理匹配的<strain>.faa文件

样本数据：

$ cat Strain-01.faa
>IMEHDJCA_03186 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>IMEHDJCA_03189 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

# for this next file I simply copied data from OP's Strain-01.faa and
# modified the initial field for blocks 1 & 3; net result is we should
# see 2 of the blocks receive insertions of ID values

$ cat Strain-02.faa
>DIBHEKPI_01635 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>DIBHEKPI_01648 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

$ cat Strain-XX.faa
cat: Strain-XX.faa: No such file or directory

对原始的awk 答案进行一些修改，并包装在 (bash) for 循环中以处理不同的菌株：

for strain in Strain-01 Strain-02 Strain-XX
do
    printf "\n############### ${strain} / ${strain}.faa\n\n"

    awk -v strain="${strain}" '                 # pass bash variable in as awk variable (same name)

    NR==1   { for (i=3; i<=NF; i++)             # 1st row of 1st file, look for matching strain name
                  { if ( $i == strain )         # if we find a match then ...
                       { strain_ndx=i           # make note of the column and ...
                         next                   # skip to next line from 1st file
                       }
                  }

              # if we got here we did not find a matching strain name so 
              # print a message and exit from our awk script 

              print "Unable to locate entry for "strain" in "FILENAME". Aborting."
              exit
            }

    FNR==NR { id[$(strain_ndx)]=$1              # for rest of 1st file build array of ids
              next
            }

    /^>/    {                                   # for 2nd file, if 1st column is ">"
              ndx=substr($1,2)                  # strip off ">"
              if ( ndx in id )                  # if 1st field (sans ">") is an index in id[] then ...
                 ( $1=">" id[ndx] " " ndx )     # rewrite 1st field to include our id[] value
            }
    1                                           # print current line (of 2nd file)
    ' file.txt "${strain}.faa"
done

这会生成：

############### Strain-01 / Strain-01.faa

>ID_01 IMEHDJCA_03186 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>ID_02 IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>ID_03 IMEHDJCA_03189 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

############### Strain-02 / Strain-02.faa

>ID_02 DIBHEKPI_01635 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>ID_01 DIBHEKPI_01648 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

############### Strain-XX / Strain-XX.faa

Unable to locate entry for Strain-XX in file.txt. Aborting.

原始答案

一个awk想法：

awk '
FNR==NR { id[$3]=$1 ; next }                # for 1st file build array of ids
/^>/    {                                   # for 2nd file, if 1st column is ">"
          ndx=substr($1,2)                  # strip off ">"
          if ( ndx in id )                  # if 1st field (sans ">") is an index in id[] then ...
             ( $1=">" id[ndx] " " ndx )     # rewrite 1st field to include our id[] value
        }
1                                           # print current line (of 2nd file)
' file.txt fasta.dat

对于给定的样本数据，这会生成：

>ID_01 IMEHDJCA_03186 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>ID_02 IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>ID_03 IMEHDJCA_03189 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

【讨论】：

不错。您可以一次扩展它以按列和文件名匹配菌株吗？
@PaulHodges ：我认为这些更新解决了您的问题/建议......？
这应该对 OP 有所帮助。我希望在awk 脚本中看到这一切，但是当问题解决后不要浪费时间来招待我，哈哈；）

【解决方案2】：

R 方法要长一些，但如果你想看看你在 REPL 中做了什么，这可能是一个不错的选择。

# read in file and split into groups, each starting with a line like "> ..."
strain <- readLines(con = 'Strain-01.faa')
strain <- split(strain, cumsum(grepl('^>', strain)))

# extract ids from line 1 of each group
ids <- sapply(strain, function(x) gsub('^>(\\w+).*', '\\1', x[1]))

# read in ID lookup table and match to the extracted IDs
id_lkp <- read.table('file.txt', header = TRUE)
id_num <- with(id_lkp, ID[match(ids, Strain.01)])

# for each group, append the id after > to the first line
for(i in seq_along(strain)){
  strain[[i]][1] <- sub('^>(\\w+)', paste0('>', id_num[i], ' \\1'), strain[[i]][1])
}

# write to output file
writeLines(unlist(strain), file('output.txt'))

^{由reprex package (v2.0.0) 于 2021-07-19 创建}

现在 output.txt 看起来像：

>ID_01 IMEHDJCA_03186 Serine/threonine-protein phosphatase 2
MEFKHRFIDGSRYQRIFVIGDIHGKLALLQDTLKRVDFHGERDLLISVGDLIDRGPDSVG
VLDYYQTHDWFEAVMGNHEWMMVNALDAQNKLERSEKEAYFIKIWHRNGCEWSQNL
>ID_02 IMEHDJCA_03187 Serine transporter
MKESRETLNFSDTLPTETWTKHDTHWVLSLFGTAVGAGILFLPINLGIGGFWPLVLLALL
AFPMTFWGHRALARFVLSSKQADADFTDVVEEHFGAKAGRLISLLYFLSIFPILLIYGVG
>ID_03 IMEHDJCA_03189 hypothetical protein
MNNQRHGITFGIERIGSQTILVFKATGTLTHQDYQAIAPVLEAALAGINRQQMNMLADIS
EFSGWEPRAAWDDFQLGLKIGFSVNKVAVYGDKNWQELAAKVGSWFISGEMKSFGD

【讨论】：