我想用另一个名称替换 fasta 文件中的序列名称

Question

我有一个fasta文件和一个文本文件fasta文件包含fasta格式的序列，文本文件包含基因名称现在我想用文本文件中的基因名称替换fasta文件中'>'符号后的序列名称 尽管我已经编写了一个脚本，但我是 perl 的新手，但我不知道为什么它不起作用，任何人都可以帮我解决这个问题 以下是我的脚本：

print"Enter annotated file...";
$f1=<STDIN>;
print"Enter sequence file...";
$f2=<STDIN>;
open(FILE1,$f1) || die"Can't open $f1";
@annotfile=<FILE1>;
open(FILE2,$f2) || die"Can't open $f2";
@seqfile=<FILE2>;
@d=split('\t',@annotfile[0]);

for($i=0;$i<scalar(@annotfile);$i++)
{
@curr_all=split('\t',@annotfile[$i]);
@curr_id[$i]=@curr_all[0];
@gene_nm[$i]=@curr_all[1];
}
for($j=0;$j<scalar(@seqfile);$j++)
{   
$id=@curr_id[$j];
$gene=@gene_nm[$j];


@seqfile[$j]=~s/$id[$j]/$gene[$j]/g;
print @seqfile[$j];
}   

我的文件如下所示：

annot.txt

pool75_contig_389 泛素连接酶 e3a
pool75_contig_704 肿瘤易感性
pool75_contig_1977 丝氨酸苏氨酸蛋白磷酸酶 4 催化亚基
pool75_contig_3064 bardet-biedl 综合征 2 蛋白 P
pool75_contig_2499 琥珀酰连接酶

goat300.fasta

goat300.fasta


>pool75_contig_704
CCCTTTCTCCCTTCCCAACATTCAGAGATACTGAATCGAAACTCTTACTGTCTGTTAGAT
GACAAAGAGTTATCCATCCTACATACTCCAATTTCCTTCCGCAACTTGTGATTTCGCCGC
TTGAATCTTGACGCCGTGCGTCCACAGTTTGTTGTGTTTTATCAATCAAGGTCATTATCA
ACCGAAGACGCTATCTATTTTCTTGGCGAAGCTCTCGGAAAGGAGCCATCGAAATGGAAG
TATTTCTCAAGAAAGTCCGCGAGTTATCCCGGAAGCAGTTC
>pool75_contig_389
GACCTATACCGGACCGTCACTGAAAGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
ACGATCCAGGCATGGAGTTGTGGTGACGAGTAGGAGGGTCACCGTGGTGAGCGGGAAGCC
TCGGGCGTGAGCCTGGGTGGAGCCGCCACGGGTGCAGATCTTGGTGGTAGTAGCAAATAT
TCAAGTGAGAACCTTGAAGGCCGAGGTGGAGAAGGNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTCATTTGTAT
CGCCCGGAAAACGTCACAAGAACGGGAGTTGCGTACAGAA
>pool75_contig_1977
AAGGGACACCGTTGGGTGAGGCGAGCTGCGTTCCTCGAACCATGGCTTCAAAAAGCGACT
TAGACCGTCAGATTGAACAGCTCAGGGCCTGCAAGCTCATTACAGAGGATGAGGTTAAGG
CACTCTGCGCTAAGGCGCGTGAGATTTTAATTGAAGAGAGTAATGTCCAGTGCGTGGACT
CACCTGTCACGGTTTGTGGCGATATCCACGGCCAGTTTTACGACTTGATTGAACTGTTTA
AAGTGGGCGGAGATGTTC
>pool75_contig_3064
TTACTATTTCTGGGCCTTAAGACTGGCTTAGTCGCTTACGACCCTTATAACAATGTAGAT
GTATATTATAAGGATCTTCCTGATGGTGCTAACGCTATGTTAATTTATTCAAACTCACCG
ACAAAGGAACAGAATATGCTTTGGCAGGTGGAAACTGTTCGATAATTGGATTGAACGACG
GCGGATGCGAGGTATTTTGGACAGTCACTGGCGACTCCGTTTGCTCTCTTTGCTCGATTA
AATCCGACAGCGATAAGTCAAGAGATTTTGTGGTTGGCTCTGAAGATTTTGACATCCGAA
TCTTCCATGGGGATGCCATAATATATGAAATCACGGAGTCTGATG
>pool75_contig_2499
AAGAGAAGAGGTGAGTTTGAGTATTGTTTGTGTGTGTGTGGTTGGGTGAGTGTGTGGTAT
GTGGTGTATGTGTGTGATGAATGTATGTGAAAGAGAGTGATGAATCTCATGGATATGTTC
GAGTTCGTGGTTTCCATTGATCGGTTATAGCCGAGATGATGGATGTGTTCCATGTGTCTG
ATTTCAGTTTAGGATTGTGTTGATGATGTTGATGATGAAAATTGTTGATGGTGATGACGA
TAGTGATGATGATGACGATGTTTCGGATAATGGTGATGATGATGATGGTTCCGACGATGA
TGTTTCGCTTGATGATGGTGATAATGATGACTCCGAAAATAACGTTGACTCGGATGAG

Answer 1

考虑使用 Bio::SeqIO 解析您的 Fasta 数据集，而不是自己进行。 Bio::SeqIO 为这项任务而生，并且为此而开发得很好。此外，如果您从事生物信息学，了解 Bio::SeqIO 对您很有帮助。鉴于此，请考虑以下事项：

use strict;
use warnings;
use Bio::SeqIO;

open my $fh, '<', 'annot.txt' or die $!;
my %annot = map { /(\S+)\s+(.+)/; $1 => $2 } <$fh>;
close $fh;

my $in = Bio::SeqIO->new( -file => 'goat300.fasta', -format => 'Fasta' );

while ( my $seq = $in->next_seq() ) {
    my $seqID = $annot{ $seq->id } // $seq->id;
    print "$seqID\n" . $seq->seq . "\n";
}

数据集的输出：

tumor susceptibility
CCCTTTCTCCCTTCCCAACATTCAGAGATACTGAATCGAAACTCTTACTGTCTGTTAGATGACAAAGAGTTATCCATCCTACATACTCCAATTTCCTTCCGCAACTTGTGATTTCGCCGCTTGAATCTTGACGCCGTGCGTCCACAGTTTGTTGTGTTTTATCAATCAAGGTCATTATCAACCGAAGACGCTATCTATTTTCTTGGCGAAGCTCTCGGAAAGGAGCCATCGAAATGGAAGTATTTCTCAAGAAAGTCCGCGAGTTATCCCGGAAGCAGTTC
ubiquitin ligase e3a
GACCTATACCGGACCGTCACTGAAAGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNACGATCCAGGCATGGAGTTGTGGTGACGAGTAGGAGGGTCACCGTGGTGAGCGGGAAGCCTCGGGCGTGAGCCTGGGTGGAGCCGCCACGGGTGCAGATCTTGGTGGTAGTAGCAAATATTCAAGTGAGAACCTTGAAGGCCGAGGTGGAGAAGGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTCATTTGTATCGCCCGGAAAACGTCACAAGAACGGGAGTTGCGTACAGAA
serine threonine-protein phosphatase 4 catalytic subunit
AAGGGACACCGTTGGGTGAGGCGAGCTGCGTTCCTCGAACCATGGCTTCAAAAAGCGACTTAGACCGTCAGATTGAACAGCTCAGGGCCTGCAAGCTCATTACAGAGGATGAGGTTAAGGCACTCTGCGCTAAGGCGCGTGAGATTTTAATTGAAGAGAGTAATGTCCAGTGCGTGGACTCACCTGTCACGGTTTGTGGCGATATCCACGGCCAGTTTTACGACTTGATTGAACTGTTTAAAGTGGGCGGAGATGTTC
bardet-biedl syndrome 2 protein P
TTACTATTTCTGGGCCTTAAGACTGGCTTAGTCGCTTACGACCCTTATAACAATGTAGATGTATATTATAAGGATCTTCCTGATGGTGCTAACGCTATGTTAATTTATTCAAACTCACCGACAAAGGAACAGAATATGCTTTGGCAGGTGGAAACTGTTCGATAATTGGATTGAACGACGGCGGATGCGAGGTATTTTGGACAGTCACTGGCGACTCCGTTTGCTCTCTTTGCTCGATTAAATCCGACAGCGATAAGTCAAGAGATTTTGTGGTTGGCTCTGAAGATTTTGACATCCGAATCTTCCATGGGGATGCCATAATATATGAAATCACGGAGTCTGATG
succinyl- ligase
AAGAGAAGAGGTGAGTTTGAGTATTGTTTGTGTGTGTGTGGTTGGGTGAGTGTGTGGTATGTGGTGTATGTGTGTGATGAATGTATGTGAAAGAGAGTGATGAATCTCATGGATATGTTCGAGTTCGTGGTTTCCATTGATCGGTTATAGCCGAGATGATGGATGTGTTCCATGTGTCTGATTTCAGTTTAGGATTGTGTTGATGATGTTGATGATGAAAATTGTTGATGGTGATGACGATAGTGATGATGATGACGATGTTTCGGATAATGGTGATGATGATGATGGTTCCGACGATGATGTTTCGCTTGATGATGGTGATAATGATGACTCCGAAAATAACGTTGACTCGGATGAG

哈希%annot 是通过读取和捕获annot.txt 数据的内容来初始化的。 Bio::SeqIO 对象是使用您的goat300.fasta 文件数据创建的。 while 循环遍历您的 fasta 序列。变量$seqID 要么采用%annot 散列中键的关联值，要么保持当前序列ID（// 表示法表示已定义或，因此确保 $seqID 将是定义）。最后打印出 Fasta 记录。

希望这会有所帮助！

Answer 2

您的代码中有很多警告，您的方法效率低下。让我首先向您展示一个有效的 Perl 程序。后面我会解释的。

#!/usr/bin/perl
use strict;
use warnings;

# Read the annotations file
print"Enter annotated file...\n";
# my $f1 = <STDIN>;
my $f1 = 'annot.txt';
open(my $fh_annotations, '<', $f1) or die "Can't open $f1";
my @annotfile = <$fh_annotations>;
close $fh_annotations;

# Read the sequence file
print"Enter sequence file...\n";
# my $f2 = <STDIN>;
my $f2 = 'goat300.fasta';
open(my $fh_genes, '<', $f2) or die "Can't open $f2";
my @seqfile = <$fh_genes>;
close $fh_genes;

# Process the annotations data
my %names; # this hash is going to hold the names
foreach my $line (@annotfile) {
  chomp $line;                      # remove newline
  my @fields = split /\t/, $line;   # split into array
  $names{$fields[0]} = $fields[1];  # save in the hash as key->value pair
}

# Process the sequence data
foreach my $line (@seqfile) {
  # Look at each line
  if ($line =~ m/>(.+)$/) {
    # If there is a heading there, remember it...
    if (exists $names{$1}) {
      # ... check if we know a name for it and replace it in the line
      $line =~ s/($1)/$names{$1}/;
    }
  }
  # output the line (this would be done to another filehandle)
  print $line;
}

这会读取这两个文件并将它们保存在内存中，就像您的一样。但我没有尝试为名称构建两个数组，而是选择了a hash，这是一个键/值对。可以把它想象成一个有名称而不是数字并且没有特定排序的数组。

一旦设置了这些名称，我就可以处理序列文件了。我simply look at each line 并通过寻找> 标志来检查那里是否有标题。如果它在那里（它进入$1，因为括号），我看看我们的%names哈希中是否有一个哈希条目（带有exists）。如果这样做，我们可以用正确的名称替换标题。

之后，我们可以将其写入一个新文件。我只是打印它。

我使用了一些其他技术。不幸的是，人们在 BioPerl 环境中获得的文献已经过时了。请采纳这个建议，它会让您的生活更轻松。

• 始终使用strict 和warnings。他们会告诉您代码存在的问题。
• 始终使用my 声明您的变量。这与其他语言不同，您需要在问题的顶部设置一个变量。您可以在需要的地方声明它。 var 仅存在于特定范围内，这意味着在最近的封闭 { 和 } 括号或块之间。
• 使用三参数打开和词法文件句柄来保证安全。阅读更多here。
• Perl 提供 foreach 作为 C for 循环的替代方案。在这种情况下，它让事情变得容易多了。

关于这个程序的另一件事：虽然这个示例数据很短，但我相信您的实际数据可能要大得多。考虑在读取序列文件时对其进行处理，以免内存不足。没有必要保存所有的行，除非你想对它们做其他事情。

open my $fh_out, '>', $filename_out or die $!;
open my $fh_in, '<', $filename_in or die $!;
while (my $line = <$fh_in>) {
  # do stuff with the line, like your regex
  print $fh_out $line;
}
close $fh_in;
close $fh_out;