abr*_*ham 2 bash awk r fasta str-replace
我有一个基因组的 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
Run Code Online (Sandbox Code Playgroud)
我想添加一个基于 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
Run Code Online (Sandbox Code Playgroud)
例如,对于 fastaStrain-01.faa文件,其IMEHDJCA_03186id 对应于Strain-01,因此我想将ID_01列 ID 的编号(来自file.txt)添加到序列的标题中,如下所示:
ID_01对应于IMEHDJCA_03186ID_02对应于IMEHDJCA_03187ID_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
Run Code Online (Sandbox Code Playgroud)
我只想将 ID 代码添加file.txt到 fasta 文件的标头中。
有想法这样做吗?或bash,R或任何其他方式?
非常感谢
更新#2:通过以下方式消除特定菌株(名称)处理awk:
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
Run Code Online (Sandbox Code Playgroud)
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
Run Code Online (Sandbox Code Playgroud)
这会生成:
>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
Run Code Online (Sandbox Code Playgroud)
注意:此最新更新将为文件中的所有行执行 ID 插入Strain-{01,02}.faa(请参阅下面的更新 #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
Run Code Online (Sandbox Code Playgroud)
对原始答案进行一些修改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
Run Code Online (Sandbox Code Playgroud)
这会生成:
############### 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.
Run Code Online (Sandbox Code Playgroud)
原答案
一个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
Run Code Online (Sandbox Code Playgroud)
对于给定的样本数据,这会生成:
>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
Run Code Online (Sandbox Code Playgroud)
| 归档时间: |
|
| 查看次数: |
854 次 |
| 最近记录: |