☰

我想获取另一个 fasta 文件中存在的整个标头的名称

tag-icon • tag-icon text-processing awk sed bioinformatics
我想获取另一个 fasta 文件中存在的整个标头的名称

我有一个这样的 fasta 文件：

>TRINITY_DN100_c0_g1_i1 len=242 path=[0:0-241]
AGCAATTCAAACTGCTGCAATCTGGGCTCGTGAAAACGATATGGTATTACACTTACATCGTGCAAGACAAACAACAAGTGCGTTACGTCACGTTACGTTAATTCGTTCGTTCGTCTCTCTATTGCGTTGCGTTACGCTCTCGCCGCGGACCCAGCACCGCATACCCGCCCATACAAGTCATACAGTACACAACACAAAAACACAACCAACTATTTCCTTGGAAGAGAAAGCAAGCCCAAAAC
>TRINITY_DN105_c0_g1_i1 len=260 path=[0:0-259]
TGAAAATATTAATTCTCAACCTTTTATGCGTTGGAGAGAAAGATATTTATTTTCAATAGAGGGAGTTAATCGTGCAGCGGCAGCAAGTGGTGAAATTAAAGGACATTATTTAAATGTTACTGCAGGTACAATGGAAGACATGTACGAACGCGCCAAGATTGATGTGCCTGAGAACCACATGAATAACGAGGAGCAATACACACTTCACTACCAAGAGTACCTTGTGGGTAGCTCGGCTGGTGTGCCCAAGGATATGAAGG
>TRINITY_DN103_c0_g1_i1 len=260 path=[0:0-259]
GTTCTCTTCGGTGGCAGCCTTACGGCCGACCACCTGGTATTGTCGCATAATTCCCGCAGCAGTCATGATGTCTATTGTTTGTCGTGAAAAGAAATGAATTAAGAGAGTCATAGTTACTCCCGCCGTTTACCCGCGCTTGGTTGAATTCCTTCACTTTGACATTCAGAGCACTGGGCAGAAATCACATTGCGTCAACACCATCTCTGTTTCAACGAAATCAGCAGTATCTGTAGAAGTGTAGTTAAAACTAATATCTTTCC

我想使用标题的开头将物种名称添加到包含在其他 fasta 文件中的标题中，并且不更改第一个 fasta 的序列。

>TRINITY_DN100_c0_g1_i1  HQ912515@Guinardia_delicatula
AGCAATCCAAACTGCTGCAATCTGGGCTCGTGAAAACGATATGGTATTACACTTACACCGTGC
>TRINITY_DN105_c0_g1_i1 KR048205@Mougeotia_transeaui
AGCAATTCAAAGTGCTGCAATCTGGGCTCGTGAAAACGATATGTTATTACACTTACACCGTGCA
>TRINITY_DN103_c0_g1_i1 RP957897@Luticola_sparsipunctata
AGCAATTCAAAGTGCTGCAATCTGGGCTCGTGAAAACGATATGATTTTACACTTACACCGTGCA

我想要这样的东西：

>TRINITY_DN100_c0_g1_i1 len=242 path=[0:0-241]@Guinardia_delicatula
AGCAATTCAAACTGCTGCAATCTGGGCTCGTGAAAACGATATGGTATTACACTTACATCGTGCAAGACAAACAACAAGTGCGTTACGTCACGTTACGTTAATTCGTTCGTTCGTCTCTCTATTGCGTTGCGTTACGCTCTCGCCGCGGACCCAGCACCGCATACCCGCCCATACAAGTCATACAGTACACAACACAAAAACACAACCAACTATTTCCTTGGAAGAGAAAGCAAGCCCAAAAC
>TRINITY_DN100_c0_g1_i1 len=260 path=[0:0-259]@Mougeotia_transeaui 
TGAAAATATTAATTCTCAACCTTTTATGCGTTGGAGAGAAAGATATTTATTTTCAATAGAGGGAGTTAATCGTGCAGCGGCAGCAAGTGGTGAAATTAAAGGACATTATTTAAATGTTACTGCAGGTACAATGGAAGACATGTACGAACGCGCCAAGATTGATGTGCCTGAGAACCACATGAATAACGAGGAGCAATACACACTTCACTACCAAGAGTACCTTGTGGGTAGCTCGGCTGGTGTGCCCAAGGATATGAAGG
>TRINITY_DN103_c0_g1_i1 len=260 path=[0:0-259]@Luticola_sparsipunctata
GTTCTCTTCGGTGGCAGCCTTACGGCCGACCACCTGGTATTGTCGCATAATTCCCGCAGCAGTCATGATGTCTATTGTTTGTCGTGAAAAGAAATGAATTAAGAGAGTCATAGTTACTCCCGCCGTTTACCCGCGCTTGGTTGAATTCCTTCACTTTGACATTCAGAGCACTGGGCAGAAATCACATTGCGTCAACACCATCTCTGTTTCAACGAAATCAGCAGTATCTGTAGAAGTGTAGTTAAAACTAATATCTTTCC

你知道我该怎么做吗？

答案1

实现您所描述的目标的一种方法是：

  1. 使用带有标头的第二个文件创建映射
  2. 读取第一个文件并用映射中的匹配项替换匹配行

使用 awk，例如：

awk -F@ 'NR == FNR { if ($0 ~ /^>/) { a[$1] = $0 }; next } { if ($0 ~ /^>/ && $0 in a) $0 = a[$0]; print }' second_file_with_headers first_file

例子：

$ cat second_file_with_headers
>TRINITY_DN100_c0_g1_i1 len=242 path=[0:0-241]@Guinardia_delicatula
AGCAATCCAAACTGCTGCAATCTGGGCTCGTGAAAACGATATGGTATTACACTTACACCGTGC
>TRINITY_DN100_c0_g1_i1 KR048205>TRINITY_DN105_c0_g1_i1 len=260 path=[0:0-259]AGCAATTCAAAGTGCTGCAATCTGGGCTCGTGAAAACGATATGTTATTACACTTACACCGTGCA
>TRINITY_DN103_c0_g1_i1 len=260 path=[0:0-259]@Luticola_sparsipunctata
AGCAATTCAAAGTGCTGCAATCTGGGCTCGTGAAAACGATATGATTTTACACTTACACCGTGCA

$ cat first_file
>TRINITY_DN100_c0_g1_i1 len=242 path=[0:0-241]
AGCAATTCAAACTGCTGCAATCTGGGCTCGTGAAAACGATATGGTATTACACTTACATCGTGCAAGACAAACAACAAGTGCGTTACGTCACGTTACGTTAATTCGTTCGTTCGTCTCTCTATTGCGTTGCGTTACGCTCTCGCCGCGGACCCAGCACCGCATACCCGCCCATACAAGTCATACAGTACACAACACAAAAACACAACCAACTATTTCCTTGGAAGAGAAAGCAAGCCCAAAAC
>TRINITY_DN105_c0_g1_i1 len=260 path=[0:0-259]
TGAAAATATTAATTCTCAACCTTTTATGCGTTGGAGAGAAAGATATTTATTTTCAATAGAGGGAGTTAATCGTGCAGCGGCAGCAAGTGGTGAAATTAAAGGACATTATTTAAATGTTACTGCAGGTACAATGGAAGACATGTACGAACGCGCCAAGATTGATGTGCCTGAGAACCACATGAATAACGAGGAGCAATACACACTTCACTACCAAGAGTACCTTGTGGGTAGCTCGGCTGGTGTGCCCAAGGATATGAAGG
>TRINITY_DN103_c0_g1_i1 len=260 path=[0:0-259]
GTTCTCTTCGGTGGCAGCCTTACGGCCGACCACCTGGTATTGTCGCATAATTCCCGCAGCAGTCATGATGTCTATTGTTTGTCGTGAAAAGAAATGAATTAAGAGAGTCATAGTTACTCCCGCCGTTTACCCGCGCTTGGTTGAATTCCTTCACTTTGACATTCAGAGCACTGGGCAGAAATCACATTGCGTCAACACCATCTCTGTTTCAACGAAATCAGCAGTATCTGTAGAAGTGTAGTTAAAACTAATATCTTTCC

$ awk -F@ 'NR == FNR { if ($0 ~ /^>/) { a[$1] = $0 }; next } { if ($0 ~ /^>/ && $0 in a) $0 = a[$0]; print }' second_file_with_headers first_file
>TRINITY_DN100_c0_g1_i1 len=242 path=[0:0-241]@Guinardia_delicatula
AGCAATTCAAACTGCTGCAATCTGGGCTCGTGAAAACGATATGGTATTACACTTACATCGTGCAAGACAAACAACAAGTGCGTTACGTCACGTTACGTTAATTCGTTCGTTCGTCTCTCTATTGCGTTGCGTTACGCTCTCGCCGCGGACCCAGCACCGCATACCCGCCCATACAAGTCATACAGTACACAACACAAAAACACAACCAACTATTTCCTTGGAAGAGAAAGCAAGCCCAAAAC
>TRINITY_DN105_c0_g1_i1 len=260 path=[0:0-259]
TGAAAATATTAATTCTCAACCTTTTATGCGTTGGAGAGAAAGATATTTATTTTCAATAGAGGGAGTTAATCGTGCAGCGGCAGCAAGTGGTGAAATTAAAGGACATTATTTAAATGTTACTGCAGGTACAATGGAAGACATGTACGAACGCGCCAAGATTGATGTGCCTGAGAACCACATGAATAACGAGGAGCAATACACACTTCACTACCAAGAGTACCTTGTGGGTAGCTCGGCTGGTGTGCCCAAGGATATGAAGG
>TRINITY_DN103_c0_g1_i1 len=260 path=[0:0-259]@Luticola_sparsipunctata
GTTCTCTTCGGTGGCAGCCTTACGGCCGACCACCTGGTATTGTCGCATAATTCCCGCAGCAGTCATGATGTCTATTGTTTGTCGTGAAAAGAAATGAATTAAGAGAGTCATAGTTACTCCCGCCGTTTACCCGCGCTTGGTTGAATTCCTTCACTTTGACATTCAGAGCACTGGGCAGAAATCACATTGCGTCAACACCATCTCTGTTTCAACGAAATCAGCAGTATCTGTAGAAGTGTAGTTAAAACTAATATCTTTCC

$

awk 脚本可以很容易地解释：

awk         # the awk program, can be gawk, mawk or any other implementation
-F@         # Sets FS="@" . This tells awk to split on the @ character.
NR == FNR { if ($0 ~ /^>/) { a[$1] = $0 }; next }
   # While handling the first argument (here, second_file_with_headers)
   # For lines starting with >, store the line into an array with the part before @ as the key
   # Discard all lines in this file with 'next'
{ if ($0 ~ /^>/ && $0 in a) $0 = a[$0]; print }
   # For all other lines, meaning that we are now working on the second argument (first_file)
   # For lines starting with >, change the line into the stored version from the array if it exists
   # then print the line
second_file_with_headers
   # Use the file with the headers as the first argument
first_file
   # Use the file which needs to be updated as the second argument

答案2

awk -F'[@[:blank:]]+' '
NR==FNR { if(/^>/) seen[$1]=$NF; next }
        { if($1 in seen) $0=$0 "@" seen[$1]; print }' second.fasta first.fasta

我们-F可以指定F我们设置为@字符的字段分隔符以及[:blank:]可以出现一次或多次的空格（制表符/空格）[...]+；因此字段将用这些字符分隔。

NF==FNR是一个始终为真的条件awk对于第一个输入文件总是如此（第二.fasta文件在此）；NR并FNR针对每行输入递增，但FNR针对下一个输入文件重置回 1。

这if(/^>/)是一个条件，检查该行是否以字符开头（^在正则表达式中指向该行的开头）>，如果是，则将最后一个字段保存$NF到关联的数组中（我们将其命名seen）并归档#1$1将用作键。

该next语句跳过处理其余代码并再次跳转到代码的开头；如果NR==FNR条件仍然为真，它将重复以下代码NR==FNR { ... }，否则将执行下一个块，并且正在第二个输入文件上执行第一.fasta我们检查是否可以在数组中找到其中的第一个字段seen，如果是，则在该行的末尾附加找到的键的值，然后执行print整行中的任何操作$0（无论是否更新） ）。

或者，如果两个文件中的所有标头均已表示并且也位于同一行号中，则您可以执行以下操作而不是缓冲第二.fasta文件存入内存：

paste -d@ first.fasta second.fasta |awk -F@ '/^>/{ $1=$1 $NF } { print $1 }'

相关内容