我是Perl的新手,我对如何完成这项任务感到困惑。我有两个文件:
它包含许多序列(database)
这里我给出我的两个文件的一小部分:
Seq.txt包含:
'>' SCO0700, probable ABC transporter protein, ATP-binding component.
MASSMEKPLDHRYRGEHPIRTLVYLFRADRRRLAGAVAVFTVKHSPIWLLPLVTAAIVDT
VVQHGPITDLWTSTGLIMFILVVNYPLHLLYVRLLYGSVRRMGTALRSALCTRMQQLSIG
'>' SCO0755, putative ABC transporter 797720:799942 forward MW:79858
VSTAQETRGRRRAAPPRRSVPKSRARTVRTPTVLQMEAVECGAASLAMVLGHYGRHVPLE
ELRIACGVSRDGSRASNLLKAARSYGFTAKGMQMDLAALAEVTAPAILFWEFNHYVVYDG
'>' SCO2305,putative ABC transporter ATP-binding subunit 2474063:2474989 forward MW:32345
MRPTEGTTPAVAFTGAAKAYGDVRAVDGVDLRIGCGETVALLGRNGAGKSTTIALLLGLC
PPDAGTVELFGGPAERAVRAGRVGAMLQEARAVPRVTVGELVAFVAGRYPAPMPVGQALE
'>' SCO1144, putative ABC transporter ATP-binding protein 1202480:1204282 reverse MW:64637
MHPDRESAWTAPADAVEQPRQVRRILKLFRPYRGRLAVVGLLVGAASLVSVATPFLLKEI
LDVAIPEGRTGLLSLLALGMIFGAVLTSVFGVLQTLISTTVGQRVMHDLRTAVYGRLQQM
'>' SCO1148, putative ABC transporter 1207772:1209553 forward MW:63721
MIGVAPPSYDPAAPTTANTLPVGARPTVRAYVGELLRRHRRAFLFLVTVNTVAVIASMAG
PYLLGGLVERVSDDARELRLGLTATLFVLALVVQAVFVREVRLRGAVLGERMLADLREDF
PID.txt包含:
SCO0755
SCO1144
我写的代码如下:
open (PID, 'PID.txt');
my @PID = '<'PID'>';
close(PID);
open (MSD, 'Seq.txt');
my @MSD = '<'MSD'>';
close(MSD);
chomp(@MSD);
my $MSD=join (' ', @MSD);
print "$MSD \n";
for ($i = 0; $i<=2; $i++) {
my $a=$PID[$i];
if ($MSD =~ m/$a(.*?)>/) # ">" end of the string
{
print "$1 \n";
$output= ">".$a.$1;
print $output;
open (MYFILE, '>>data.txt');
print MYFILE "$output\n";
close (MYFILE);
}
}
为什么它不能识别$a
?如果我放入$a,那么绑定操作符会识别$a
,但不会返回我想要的序列( ID存储在$a
中),而是返回第一个序列。
我期望的结果是:
'>' SCO0755, putative ABC transporter 797720:799942 forward MW:79858
VSTAQETRGRRRAAPPRRSVPKSRARTVRTPTVLQMEAVECGAASLAMVLGHYGRHVPLE
ELRIACGVSRDGSRASNLLKAARSYGFTAKGMQMDLAALAEVTAPAILFWEFNHYVVYDG
'>' SCO1144, putative ABC transporter ATP-binding protein 1202480:1204282 reverse MW:64637
MHPDRESAWTAPADAVEQPRQVRRILKLFRPYRGRLAVVGLLVGAASLVSVATPFLLKEI
LDVAIPEGRTGLLSLLALGMIFGAVLTSVFGVLQTLISTTVGQRVMHDLRTAVYGRLQQM
发布于 2011-02-22 18:02:45
我不能告诉您为什么会得到这样的输出,因为您发布的代码不是有效的Perl,不能编译或运行- my @PID = '<'PID'>';
在语法上是无效的。(应为my @PID = <PID>;
,不带任何引号。)因此,这显然不是您运行来产生这些结果的代码。
您没有得到任何匹配的原因是,尽管发布的代码执行了chomp(@MSD)
,但它也没有执行chomp(@PID)
,所以PID只有在后面跟一个换行符时才会匹配。在发布的数据中,它们不是。(即使它们是,chomp(@MSD)
也会删除它们。)
解决这个问题会让您更进一步,但仍然不会产生您想要的结果,因为您的正则表达式是错误的。试试这个(将$a
重命名为$target
,因为a:它是一个更有意义的名称,b:$a
和$b
很神奇,所以您不应该使用它们):m/'>' $target([^']*)/
最后,您的for ($i...)
循环是不正确的,这是在C风格的for
中很容易犯的错误。使用for (list)
要好得多。
修复所有这些东西,以及切换到词法文件句柄和open
的三参数形式(正如David Precious已经提到的),并进行一些通用的代码清理,给了我们:
#!/usr/bin/env perl
use strict;
use warnings;
open my $pid_fh, '<', 'PID.txt';
my @PID = <$pid_fh>;
close $pid_fh;
chomp(@PID);
open my $msd_fh, '<', 'Seq.txt';
my @MSD = <$msd_fh>;
close $msd_fh;
chomp(@MSD);
my $msd = join(' ', @MSD);
my $output;
open my $outfile, '>>', 'data.txt';
for my $target (@PID) {
if ($msd =~ m/'>' $target([^']*)/) {
$output = ">" . $target . $1;
print $output, "\n";
print $outfile "$output\n";
}
}
...which生成以下输出:
>SCO0755, putative ABC transporter 797720:799942 forward MW:79858 VSTAQETRGRRRAAPPRRSVPKSRARTVRTPTVLQMEAVECGAASLAMVLGHYGRHVPLE ELRIACGVSRDGSRASNLLKAARSYGFTAKGMQMDLAALAEVTAPAILFWEFNHYVVYDG
>SCO1144, putative ABC transporter ATP-binding protein 1202480:1204282 reverse MW:64637 MHPDRESAWTAPADAVEQPRQVRRILKLFRPYRGRLAVVGLLVGAASLVSVATPFLLKEI LDVAIPEGRTGLLSLLALGMIFGAVLTSVFGVLQTLISTTVGQRVMHDLRTAVYGRLQQM
选择了正确的序列;我将把它们完全按照您的要求格式化,作为读者的练习。
发布于 2011-02-22 17:37:08
首先,不要在代码中使用$a
和$b
。它们是特殊的变量,只有在sort
块中才有意义;在其他地方避免使用它们,而使用有意义的变量名。
第二,
my @PID = '<'PID'>';
假设您正在尝试将文件句柄PID的内容读入一个数组,您的意思是:
my @PID = <PID>;
第三,目前常见的最佳实践是使用3-arg open和lexical文件句柄,例如:
open(my $pidfh, '<', 'PID.txt') or die "...";
my @PID = <$pidfh>;
close $pidfh;
你的脚本的顶部有没有use strict;
?
为了方便查找,我会将您感兴趣的PID读入一个哈希表中,然后循环遍历Seq.txt;记住您正在查看并存储其内容的条目;每次看到新条目时,查看您构建的前一个条目是否是您想要的条目,如果是,则将其打印出来。这样你就不需要在内存中保存文件的内容,这在文件非常大的情况下很有用。
大致如下所示:
#!/usr/bin/perl
use strict;
# Read in a list of PIDs we're interested in
my %want_pid;
open(my $pidfh, '<', 'PID.txt') or die "Failed to open PID.txt - $!";
while (my($pid) = <$pidfh> =~ m{([A-Z0-9]+)}) {
$want_pid{$pid}++ if $pid;
}
# Now process the file and print entries we want
open(my $seqfh, '<', 'Seq.txt') or die "Failed to open Seq.txt - $!";
my $current_pid;
my $current_text;
while (my $line = <$seqfh>) {
if (my ($new_pid) = $line =~ m{^ '>' \s+ ([A-Z0-9]+) , }x) {
# We're at the start of a new entry; if the last one is one we want,
# print it.
if ($want_pid{$current_pid}) {
print $current_text;
}
$current_pid = $new_pid;
$current_text = $line;
} else {
# It's a continuation of an entry
$current_text .= $line;
}
}
close $seqfh;
(还有改进的空间,但这应该会让你走上正轨。)
发布于 2011-02-22 17:59:09
测试一下这对你是否有效:
use warnings;
use strict;
die "Usage: $0 <pid file> <seq file>\n" unless @ARGV == 2;
open my $pid, "<", $ARGV[0] or die "Error: Cannot open file $ARGV[0]: $!\n";
open my $seq, "<", $ARGV[1] or die "Error: Cannot open file $ARGV[1]: $!\n";
my %pid = ();
while ( <$pid> ) {
chomp;
s/^\s*(\S*)\s*$/$1/;
++$pid{$_};
}
$/ = "\'>\'";
foreach ( <$seq> ) {
$_ = substr $_, 0, -3;
my ($p) = split /\,/;
$p =~ /(\S+)/;
print "'>'", $_ if exists $pid{$1};
}
致以敬意,
https://stackoverflow.com/questions/5075889
复制相似问题