求幫忙，寫的腳本不知錯(cuò)在哪里

cxrfly

1、

我有一個(gè)文件a，包括91個(gè)蛋白序列ID；另外一個(gè)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個(gè)ID，從另外的文件中抽取這些ID對應(yīng)的氨基酸序列。
腳本是這樣寫的，不對，沒有結(jié)果文件產(chǎn)生。煩您指點(diǎn)下！


#!/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://www.72891.cn/forum.p ... mp;fromuid=29097174

歐陽西風(fēng)

回復(fù) 1# cxrfly


    已在原shell版給你回復(fù)了，這兒就不再復(fù)制一遍了。

cxrfly

回復(fù) 2# huang6894


    您寫的好簡潔，好多我還沒有學(xué)，呵呵，需要消化下。

cxrfly

回復(fù) 3# 歐陽西風(fēng)


    謝謝你的提醒。我現(xiàn)在發(fā)現(xiàn)確實(shí)是我由于修改了行符，導(dǎo)致ID文件讀入后不對。
我已經(jīng)修改了，程序正常運(yùn)行了，謝謝你。

cxrfly

謝謝huang6894和歐陽西風(fēng)兩位的幫忙。我努力看了點(diǎn)書，自己把我之前寫的修改了下，現(xiàn)在可以準(zhǔn)確運(yùn)行了，也得到了希望的結(jié)果�？赡芪业某绦虿皇亲詈啙嵶詈玫�，但是畢竟它可以運(yùn)行了。
貼出來，希望有同樣困惑的人可以借鑒下，也請指導(dǎo)我如何提高代碼寫作水平。

#!/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，什么是氨基酸序列呢