一個處理生物數(shù)據(jù)的perl程序的問題，研究了好長時間沒有找出來，煩高手�。�！

aids260

#!/usr/local/bin/perl
#done by zmyin aids260@163.com
($input1,$input2,$output)=@ARGV;
open I1,"$input1"||die"$!";
open I2,"$input2"||die"$!";
open O,">>$output"||die"$!";
chomp(@name=<I1>);
$n;
$/=">";
@a;$b;$c;$d;

@cds_seq=<I2>;
$i=0; #cds計數(shù)
foreach $n(@name){
           $d=$n;            
            foreach $b(@cds_seq){
                        $b=~/(\w+-TA)/;
                        $c=$1;
                                    if($d eq $1){
                         $b=~s/>$//;
              $b=">".$b;
              print O $i;
                          $i++;
                     last;
                }
                }
    }
print O "$i full cds\n";       
close I1;
close I2;
close O;



上邊是我的源程序，使用命令是   exact-cds-sequence.pl   cds-name.txt   silkcds.fa   cds-sqe.txt
其中exact-cds-sequence.pl是程序名稱  后邊三個參數(shù)，前兩個是已有的文件名，最后一個是生成的結(jié)果
我把圖片和 cds-name.txt   silkcds.fa   兩個文件都傳上去了，我的目的是（如上傳的圖片所示）把cds-name.txt   的一些基因的名字
比如“BGIBMGA000062-TA”
在silkcds.fa   文件中找到 BGIBMGA000062-TA的名字以及它的全部序列，比如是：

>BGIBMGA000062-TA  cds:novel  sequence:nscaf1071:1239819:1245122:+  gene:BGIBMGA000062  protein:BGIBMGA000062-PA
ATGAAGTTAGTTCAGTTTTCATACAAAGATAGTCCAAAAAATATACGTGT
GGGCTACCTGGAAGGAGATGATATTGTAGATATTAATAAGGCGGACTCCA
GTTTGCCGACCACTTTGCTCCAAATACTCAGGAATGGAGACTTAGAAAAA
GTGAAGAAGTTGAAATCAACAAAACCAGCAACTATACCACTATCATCAGT
CACTCTAACTGCACCCATACATGGTGTAGATAAAATCCTCTGTATCGGCT
TGAACTACAAGGATCACTGCCAAGAGCAGAATCTCACCCCACCTCCTGTG
CCGATGGTGTTCAGTAAATTTTCAAGCACCATCATTGGACCTGATCAGCC
TGTTAGGATCAGAACTGATGTTACTAAGAAGGTGGACTGGGAGGTGGAGC
TGTGCGTGGTGGTGGGGCGCGAGGCCAGCTGCGTGCGCGAGGAGGACGCG
CTGCAGCACGTGGCCGGGTACACCGTCGCGCAGGACATCAGCGCCAGGGA
CTGGCAGAAAGAGAAGAACATGGGGCAGTTCCTGCTAGGGAAGTCCATGG
ACACGTTCTGTCCGCTGGGCCCGTGCGTGCTGACGAGCGACGAGGTGGGC
GCGGCCGTGGAGCTGCGCGTGTCCTGCTCGCTCAACGGGGTCCTCAAGCA
GAGCAGCAGCACGGCGCAGCTAGTGCACTCCATCCCGAGCCTGCTGCACA
GGATCTCCTCCGTGATGACCCTGGTCCCCGGCGACCTGATCCTGACGGGC
ACCCCGGGGGGCGTGGGCATGTACCGGCAGCCCCCGGAGTACCTGCAGCC
CGGGGACGTGCTCACCAGCGAGATCGAGAAGATCGGCGCCTTCGACGTTC
GCATCGAGAAGTTTTAG


把這些輸入到文件cds-sqe.txt中

longbow0

沒看到圖，你是不是想從一個 fasta 文件 中找到含有指定字符串的序列，然后輸出到文件中？

1. 
   
2. # 部分代碼
   
3. 
   
4. # 讀fasta文件
   
5. my $o_seqi = Bio::SeqIO->new(
   
6.     -file => 'silkcds.fa',
   
7.     -format => 'fasta',
   
8. );
   
9. 
   
10. # 輸出fasta文件
    
11. my $o_seqo = Bio::SeqIO->new(
    
12.     -file => ">cds-seq.txt",
    
13.     -format => 'fasta',
    
14. );
    
15. 
    
16. # 在 description 中查找字符串 $cds_name
    
17. while (my $o_seq = $o_seqi->next_seq) {
    
18.     if ($o_seq->description =~ /$cds_name/ ) {
    
19.         $o_seqo->write_seq($o_seq);
    
20.     }
    
21. }
    

復(fù)制代碼

另外，最好用

1. 
   
2. use strict;
   
3. use warnings;
   

復(fù)制代碼

toniz

回復(fù) 1# aids260


    不知道你想說些啥。。。最好是能表達(dá)清楚點，把文件的數(shù)據(jù)格式也貼出來，然后遇到什么問題 。

aids260

回復(fù) 2# longbow0


    我把圖片上傳了，我的目的是想從指定fasta文件中找出指定名字的基因序列，把它輸出到指定文件中��！  謝謝了啊，再幫我看一下啊

aids260

回復(fù) 3# toniz


    麻煩幫我看一下

longbow0

看看 Bio::SeqIO， Bio::Seq

1. 
   
2. while (my $o_seq = $o_seqi->next_seq) {
   
3.      if ($o_seq->id =~ /$cds_name/ ) {
   
4.          $o_seqo->write_seq($o_seq);
   
5.      }
   
6. }
   

復(fù)制代碼

x9x9

學(xué)習(xí)著寫了一個，不知道對不對；

1. use strict;
   
2. use warnings;
   
3. use 5.010;
   
4. 
   
5. (my $input1,my $input2,my $output)=@ARGV;
   
6. open I1,"$input1"||die"$!";
   
7. open I2,"$input2"||die"$!";
   
8. open O,">$output"||die"$!";
   
9. 
   
10. my @cds=<I1>;
    
11. 
    
12. $/="";
    
13. my $lines=<I2>.">";
    
14. 
    
15. foreach (@cds) {
    
16. chomp;
    
17. print O $1,"\n" if $lines~~/($_.*?)>/s;
    
18. }
    
19. close I1;
    
20. close I2;
    
21. close O;

復(fù)制代碼

toniz

回復(fù) 1# aids260


    其實樓主你的代碼可以實現(xiàn)這個功能的。不過有點小毛病而已。改改就Ok了。

1. foreach $n(@name){
   
2. $d=$n;   
   
3. $d =~s/\s+//g;

復(fù)制代碼

這里加個去空格的。我用perl，是很不相信chomp這東東的。

然后

1.         if($d eq $1){
   
2.           $b=~s/>$//;
   
3.                        $b=">".$b;
   
4.                        print O $b;

復(fù)制代碼

要輸出的是 $b 不是 $i

改這兩個地方就可以輸出到數(shù)據(jù)了。

aids260

回復(fù) 8# toniz


    謝謝你，你說的就是問題之所在��！