-1

我有一个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];
}

我的文件如下所示：

注释.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
4

2 回答 2

0

您的代码中有很多警告，您的方法效率低下。让我首先向您展示一个有效的 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;
}

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

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

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

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

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

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

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;
于 2013-03-23T10:38:40.457 回答
0

考虑使用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 记录。

希望这可以帮助！

于 2013-03-23T17:32:12.317 回答