求帮忙，写的脚本不知错在哪里

cxrfly

1、

我有一个文件a，包括91个蛋白序列ID；另外一个b包括2万多蛋白序列。

File a:
jgi|Physo3|320426|gm1.1_g
jgi|Physo3|551829|estExt_Genewise1Plus.C_1_t10001


File b:
>jgi|Physo3|320426|gm1.1_g
MAILAYVNFKREGLVSSYASDPVLALGAMKVCAGSGEVGEMVARILLLLAMDKCVMGDKDFSQCYHTIGQ
FVPVQQFLDVLQGNKERPIYRKKMDSEITKGPEFKVRRSKWEGWYVGFTHFVQLHLEPNEDTLWFLLGRR
AAGIFPRNQSGADLLIPMFKKRSNSRPVGGETTRGEQEGHEVLMMLIQNQLSFTGVHDTIRIYMTARGQE
ADGNCIYTTTMKDSLQKSPHQRGEATAKEDGMEAVGGNTVGASEMRAESLAKNETNPGVTIISEETQEGK
EKEPRRILVPERSMLANMVKSQYDPLGLVVDDRGDDEKQSVSSYAKLSVTANSAKLLKSARKALACAPGC
SSKNFGVNAVDDGSKGTSHCHNKGGDTIKRPREKKLEQSQGRKRAKCER*
>jgi|Physo3|551829|estExt_Genewise1Plus.C_1_t10001
MTTPARSPLSAPSSLLRRVASPVPAGPGSGGGARSQTLLSRLFRRKDRSPAAPRPLSRSLYGEEEETKTH
VDSPMRALLRSLSLQRKSRPQEVRTPPADAAQVLLEGVLTELAQTRYWRSKPCYFFTFTLRQFHFLPETR
RPGRLVAIATLSPDDVVSDLRTHRQVKHALELQLAAAHQRVEIAGARDADDDTREESASQLLRLAATDER
TLVKWLQMLNSCISRLVLQAEERQYRHGDAAGWNPVTEPAAKIVPMFGDEGRRSSMLVMTRDRDHNLSAI
SFTSSECSNNNHYEEDIPAEEIRAARKDAVGEEDEDETKPPADERYDRVTPALEDEELPARRAIEYCCTP
EYRSSGSGSSNSQNLVDTDNNNAADEKEGNNQEDNQGDDDSEEDSLAGDLATQSSVAPYKSAVLAGRGVA
AYQSVGLEDEYTKAGSRSAVANGRGRLPRAAFGATTSSSGLFPTLRPQTHSFNANGQVFTLDTRYQLVKP
IGNGAYGAVIAVKDVVNGGENLAVKKITNIFEDLVDAKRILREVRLLGHFNHKNITRLLDLSPPPSRKQF
DDMYIITELMETDLHQVIYSMQPMSDDHVKYFLYQMLCALHHIHSAGVLHRDMKPSNILLNANCDLKVCD
FGLARGGVGGQEELLQPGELTEYVVTRWYRAPEIMLNCLHYTTAIDVWAVGCIFAEMLLREPLFPGNDYL
HQLKLIIKFLGTPKQEDIGFVKNTKAMRFLTKLAISKPKKWRDVFAGSGTENAVSSEAIDLLSKMLFFNP
EKRISVDAALRHPYLATFFDENDLVVSTAFDFSFDLPDDQLSKDALINLLCEDIEQFHPPVPTSVTPMPL
SAAASRFFRMGMTASAS*


2、

我想利用91个ID，从另外的文件中抽取这些ID对应的氨基酸序列。
脚本是这样写的，不对，没有结果文件产生。烦您指点下！


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

#welcome words
print "\n######Dear,";

print " Welcome you to use this code!\n";

print " Please give me your feedback after using.######\n\n\n";



my $proteinfile='file b';
unless (open(PROTEINFILE, $proteinfile)){
print "Could not open the file $proteinfile!\n";
}
#reset line seperator '$/' to '>' read in each set of genes one by one
$/= ">";

while (<PROTEINFILE>){
my $protein_ID='';
my $protein_sequence='';

if ($_=~/^>/){
next;}

if ($_=~/(.*\n)/){
$protein_ID=$1;
}
$_=~s/^[^\n]+\n//;

$protein_sequence=$_;
$protein_sequence=~s/>$//;


my $inputfile="file a";
my @ID_file=get_file ($inputfile);

my $ID='';

foreach $ID (@ID_file){
chomp $ID;
if ($protein_ID eq $ID){
my $data=">$protein_sequence\n";

my $outputfile = 'test.output';
open (OUTPUT,">>$outputfile");

print OUTPUT  "$data\n";

close OUTPUT;

   }
else {
exit;
       }

                              }
}
close PROTEINFILE;
exit;

#####################################################

#read in data from one file you input
sub get_file {
my ($filename)=@_;
unless (open (GET_FILE_DATA,$filename)){
print STDERR "Cann\'t open \"$filename\"\n";
exit;
}
my @filedata = ();
@filedata=<GET_FILE_DATA>;
close GET_FILE_DATA;
return @filedata;
}

huang6894

不考虑顺序且文件a ID唯一的话：


否则：

http://bbs.chinaunix.net/forum.p ... mp;fromuid=29097174

欧阳西风

回复 1# cxrfly


    已在原shell版给你回复了，这儿就不再复制一遍了。

cxrfly

回复 2# huang6894


    您写的好简洁，好多我还没有学，呵呵，需要消化下。

cxrfly

回复 3# 欧阳西风


    谢谢你的提醒。我现在发现确实是我由于修改了行符，导致ID文件读入后不对。
我已经修改了，程序正常运行了，谢谢你。

cxrfly

谢谢huang6894和欧阳西风两位的帮忙。我努力看了点书，自己把我之前写的修改了下，现在可以准确运行了，也得到了希望的结果。可能我的程序不是最简洁最好的，但是毕竟它可以运行了。
贴出来，希望有同样困惑的人可以借鉴下，也请指导我如何提高代码写作水平。

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

#welcome words
print "\n######Dear, ";

print " Welcome you to use this code!\n";

print " Please give me your feedback after using.######\n\n\n";

my $inputfile="file a";

unless (open (GET_ID,$inputfile)){
print STDERR "Cann\'t open \"$inputfile\"\n";
exit;
}

my @id1 = ();
my @id2 = ();

@id1=<GET_ID>;
foreach my $fileid (@id1){
$fileid=~ s/^\s*|\s*$//g;
push @id2, $fileid;
}
close GET_ID;


my $proteinfile="file b";
unless (open(GET_PROTEIN, $proteinfile)){
print "Could not open the file $proteinfile!\n";
}


#reset line seperator '$/' to '>' read in each set of genes one by one
$/= ">";

while (<GET_PROTEIN>){
my $protein_ID='';
my $protein_sequence='';

if ($_=~/^>/){
next;}

if ($_=~/(.*)\n/){
$protein_ID=$1;
}

$_=~s/^[^\n]+\n//;
$protein_sequence=$_;
$protein_sequence=~s/>$//;
$protein_sequence=~ s/\s*//g;

foreach my $ID (@id2){

if ($ID eq $protein_ID){
my $data=">$protein_ID\n$protein_sequence\n";

my $outputfile = 'test.output';
open (OUTPUT,">>$outputfile");

print OUTPUT  "$data\n";

close OUTPUT;
print "$ID\n";
   }
}
}

close GET_PROTEIN;
exit;

youzhengchuan

什么叫蛋白序列id，什么是氨基酸序列呢