【发布时间】：2021-06-23 22:14:06
【问题描述】：
我正在尝试将FASTA 文件读入哈希。
问题
my $primer_seq_14bp = q(); // --> works
my $primer_seq_14bp; // --> error
Use of uninitialized value $primer_seq_14bp in hash element at ./script.pl line 24, <FWDPRIMER> line 1.
为什么这个变量需要声明为空字符串？
这是我第一个“我自己”脚本的一部分，因此欢迎就格式、效率等提供反馈！
代码
#!/usr/bin/perl
use strict; use warnings;
my %primer_name_for;
my $primer_name;
my $primer_seq_14bp;
# Store fwd primer sequences in a sequence{name} hash.
open FWDPRIMER, '<', $ARGV[0] or die "error reading $ARGV[0]\n";
while (<FWDPRIMER>) {
chomp;
if ( /^(>.+)/ ) {
$primer_name = $1;
}
else {
$primer_seq = $_;
$primer_seq_14bp = substr $primer_seq, 0, 14;
}
$primer_name_for{$primer_seq_14bp} = $primer_name;
}
close FWDPRIMER;
输入
Fasta 文件，已格式化：
>Sequence_name
DNA_sequence
>Sequence_name
DNA_sequence
等等
整个文件都是这种格式，所以我想我可以使用$_，因为如果该行不是名称，它必须是一个序列。
目标
一个哈希，其中DNA_sequence 是Sequence_name 的键。最终的操作是使用 read 序列搜索 hash 以获得 read 的生物学名称（reads 带有默认名称）。
示例
输入引物
>snp_fwd_primer
AAGCTCCTGCAGGTCATCTC
输入读取
>read_name
AAGCTCCTGCAGGTCATCTCTAGTTGACACCTTTGCTGACAATTATTGTG
期望的输出
请注意，read 的前 20 bp 与引物序列匹配。该脚本将读取序列切割为 14 bp，因为这是我最短的引物序列的长度。
>snp_fwd_primer
AAGCTCCTGCAGGTCATCTCTAGTTGACACCTTTGCTGACAATTATTGTG
更新
我有一个解决方案，但是，下面的答案要好得多。我正在为其他迷路的生物学家分享它（我的解决方案基于生物信息学论坛）。
我更改了记录分隔符 $/ 以一次读取输入引物 fasta 一个条目（我总是将多行 fasta 转换为单行；对于 multi-fasta，请参阅@ikegami 的答案）。然后使用正则表达式将标题与序列分开。主要的缺点是我必须从输入引物 fasta 中手动删除第一个 >（了解到您不能在 perl 脚本中关闭并重新打开同一个文件）。
这不适用于输入 read fasta，因为我只想要序列行，所以我跳过了任何以 > 开头的行。
#!/usr/bin/perl
# Usage: thisScript.pl fwdPrimers.fasta reads.fasta > reads_with_primer_names.fasta
# Fasta files need to be in single-line format.
# Before running, remove the first instance of '>' from the fwdPrimerFastaFile!
# Purpose: for reads that have marker forward primer sequences, rename the read with the marker name.
use strict; use warnings;
open PRIMERFILE, '<', $ARGV[0] or die "error reading $ARGV[0]\n";
$/ = '>';
# Store fwd primer sequences in a sequence{name} hash.
my %primer_name_for;
while (<PRIMERFILE>) {
/(.+)\n(.+)\n/ and my $primer_name = $1 and my $primer_seq = $2;
my $primer_seq_14bp = substr $primer_seq, 0, 14;
$primer_name_for{$primer_seq_14bp} = $primer_name;
}
close PRIMERFILE;
# Store reads in a full_length_read{trimmed_read} hash.
# Use trimmed read to search the fwd primer hash.
# Combine primer name and full length read.
open READFILE, '<', $ARGV[1] or die "error reading $ARGV[1]\n";
$/ = "\n";
my %read_fullSeq_for;
my $read_seq;
while (<READFILE>) {
chomp;
if($_ =~ /^>/){
next;
}
else {$read_seq = $_};
my $read_seq_14bp = substr $read_seq, 0, 14;
$read_fullSeq_for{$read_seq_14bp} = $read_seq;
if (exists $primer_name_for{$read_seq_14bp} and $read_fullSeq_for{$read_seq_14bp}) {
print ">$primer_name_for{$read_seq_14bp}\n$read_fullSeq_for{$read_seq_14bp}\n";
}
}
close READFILE;
__END__
【问题讨论】：
-
您正在使用
undef作为哈希键。这就是警告的意思。undef将被转换为空字符串。不是一个非常有意义的哈希键。解决方法是在未定义变量时不使用该变量作为哈希键。 -
从某种意义上说，当你说它“有效”时，它实际上并不有效。你只是没有收到关于它不起作用的警告。
-
欢迎来到 SO！根本问题是这些变量无缘无故地被前向删除。为什么不在使用它们的地方声明/限定它们？更严格的范围 = 更少的错误，更容易理解代码。
-
似乎最好不要重新发明轮子来阅读 FASTA：metacpan.org/dist/TIGR/view/lib/TIGR/FASTA/Reader.pm。使用嵌套循环可以避免作用域，外循环读取键，内循环读取值，将外循环留在下一个键。维基百科似乎建议 FASTA 值可以是多行的，所以我仍然不确定 OP 的基于伪代码的文件的实际规范是什么，尽管现在有更多的上下文（谢谢，OP！）。
-
@Kip 您可能应该使用一个旨在读取 fasta 文件的模块。比如可能 bioperl，但我没有这样的经验。
标签： perl