如何使用不同文件中的序列 ID 从文件中提取 FASTA 序列？

Question

我有两个文件：

sequence.fasta - 包含多个 FASTA 序列的大文件

ids.txt - 由制表符分隔格式的序列 ID 组成。

我想将这些序列提取到 sequence.fasta 的另一个文件中，其 ID 在 ids.txt 中匹配。

sequence.fasta

的样本

>AUP4056.1
MFKSLIQFFKSKSNTSNIKKENAVQRQERQDIEGWITPYSGQELLNTELRQHHLGLLWQQVSMTREMFEH
LYQKPIERYAEMVQLLPASESHHHSHLGGMLDHGLEVISFAAKLRQNYVLPLNAAPEDQAKQKDAWTAAV
IYLALVHDIGKSIVDIEIQLQDGKRWLAWHGIPTLPYKFRYIKQRDYELHPVLGGFIANQLIAKETFDWL
ATYPEVFSALMYAMAGHYDKANVLAEIVQKADQNSVALALGGDITKLVQKPVISFAKQLI`

>XIM5213.2
FKISSKGPGDGWLTEDGLWLMSKTTADQIRAYLMGQGISVPSDNRKLFDEMQAHRVIESTSEGNAIWYCQ
LSADAGWKPKDKFSLLRIKPEVIWDNIDDRPELFAGTICVVEKENEAEEKISNTVNEVQDTVPINKKENI
ELTSNLQEENTALQSLNPSQNPEVVVENCDNNSVDFLLNMFSDNNEQQVMNIPSADAEAGTTMILKSEPE
NLNTHIEVEANAIPKLPTNDDTHLKSEGQKFVDWLKD

ids.txt

的样本

AUP4056.1 GUP5213.2 ARD5364.5 HAE6893.7
JIK6023.5 YUP7086.9

我需要输出如下

>AUP4056.1
MFKSLIQFFKSKSNTSNIKKENAVQRQERQDIEGWITPYSGQELLNTELRQHHLGLLWQQVSMTREMFEH
LYQKPIERYAEMVQLLPASESHHHSHLGGMLDHGLEVISFAAKLRQNYVLPLNAAPEDQAKQKDAWTAAV
IYLALVHDIGKSIVDIEIQLQDGKRWLAWHGIPTLPYKFRYIKQRDYELHPVLGGFIANQLIAKETFDWL
ATYPEVFSALMYAMAGHYDKANVLAEIVQKADQNSVALALGGDITKLVQKPVISFAKQLI

>GUP5213.2
ELTSNLQEENTALQSLNPSQNPEVVVENCDNNSVDFLLNMFSDNNEQQVMNIPSADAEAGTTMILKSEPE
NLNTHIEVEANAIPKLPTNDDTHLKSEGQKFVDWLKDKLFKKQLTFNDRTAKVHIVNDCLFIVSPSSFEL
YLQEKGESYDEECINNLQYEFQALGLHRKRIIKNDTINFWRCKVIGPKKESFLVGYLVPNTRLFFGDKIL
INNRHLLLEE

我已经尝试过 Perl 单行代码，但这行不通。既不给出任何错误也不给出任何输出。

perl -ne 'if(/^>(\S+)/){$c=$i{}}$c?print:chomp;$i{$_}=1 if @ARGV' ids.txt sequence.fasta

谁能帮我更正这段代码，或者是否还有其他 Perl 脚本？

Answer 1

这里的问题是 one-liners 很难遵循、理解和理清。

所以写出来'long hand':

#!/usr/bin/env perl

use strict;
use warnings;

open ( my $id_file, '<', 'ids.txt' ) or die $!;
#use split here, to split any lines on whitespace. 
chomp ( my @ids = map { split } <$id_file> );
close ( $id_file );

my %sequences;

open ( my $input, '<', 'sequence.fasta' ) or die $!;
{
   local $/ = '';    #paragraph mode; Read until blank line

   while ( <$input> ) {
      my ( $id, $sequence ) = m/>\s*(\S+)\n(.*)/ms;
      $sequences{$id} = $sequence;
   }
}

foreach my $id (@ids) {
   if ( $sequences{$id} ) {
      print ">$id\n";
      print "$sequences{$id}\n";
   }
}

如果你想从@ARGV读取文件名：

my ( $ids_file, $sequence_file ) = @ARGV;

我不会尝试将它压缩回一个衬里 - 你可能可以，但当你回来时会很难理解它。

Answer 2

如果你想要一个衬里 - 你的 post 事实上建议 - 这就是你可以做的：

perl -pe '$i=if/^>(\S+)/;map$i{$_}++,split;$i{$i}or$_=""' ids.txt seq.fasta

如何使用不同文件中的序列 ID 从文件中提取 FASTA 序列？

How to extract FASTA sequences from a file using sequence IDs in adifferent file?

perl

fasta