根据score筛选序列，求帮助！ (已解决，谢谢b114213903!)

56836430

我有两个文件，一个是align好的序列文件，fasta格式，如(a.fasta)：

>lcl|Abi_c29393_g1_i1_m.53755 unnamed protein product
MIHFVLLISRQGKVRLTKWYSPYSQKERSKVIREISGVILTRGPKLCNFI
>lcl|Aco_Aquca_038_00013.1 unnamed protein product
MIHFVLLISRQGKVRLTKWYSPYTQKERTKVIRELSGVILTRGPKLCNFV
>lcl|Aco_Aquca_038_00013.2 unnamed protein product
MFL----------------------QVIRELSGVILTRGPKLCNFVEWRG

另外一个文件是以上序列每个位点的score (b.txt):
0.033658
0.034218
0.038550
0.042490
0.150432
0.214897
0.466222
8.644571
8.904594
8.971139
9.004972
9.013395
9.014203
9.014603
9.012165
9.012979
9.013093
9.013061
9.013095
9.298905
9.300411
9.300800
9.300818
9.300811
9.300889
9.300984
9.300952
9.300918
9.300968
9.300942
9.300397
9.300424
9.300866
9.301543
9.302279
9.303268
9.934049
9.978881
9.993884
9.999056
9.999716
9.999755
9.999635
9.999471
9.998696
9.998483
9.998806
9.999327
9.999733
9.999860
9.999929
9.999871
9.999862

想根据score大小来筛选序列，当score大于5时，碱基得以保留，得到以下输出文件c.fasta：
>lcl|Abi_c29393_g1_i1_m.53755 unnamed protein product
ISRQGKVRLTKWYSPYSQKERSKVIREISGVILTRGPKLCNFI
>lcl|Aco_Aquca_038_00013.1 unnamed protein product
ISRQGKVRLTKWYSPYTQKERTKVIRELSGVILTRGPKLCNFV
>lcl|Aco_Aquca_038_00013.2 unnamed protein product
------------------QVIRELSGVILTRGPKLCNFVEWRG

请大家帮个忙，应该怎么解决?

b114213903

1. #!/usr/bin/perl
   
2. use strict;
   
3. use Bio::SeqIO;
   
4. 
   
5. #限制程序参数
   
6. unless (@ARGV==2){
   
7.         print "Usage:\t$0\tScore_file\tFasta_file\n";
   
8.         sleep 3;
   
9. }
   
10. my $score_f=shift@ARGV;
    
11. my $fasta_f=shift@ARGV;
    
12. 
    
13. #读取Score数据
    
14. my $Scores={};
    
15. open (my $sc,"<$score_f") or die "Open score file failed!\n";
    
16. my $n=0;
    
17. while (<$sc>){
    
18.         chomp;
    
19.         if($_){
    
20.                 $n++;
    
21.                 $Scores->{$n}=$_;
    
22.         }
    
23. }
    
24. close $sc;
    
25. 
    
26. #读取Seq信息
    
27. my $Data={};
    
28. my $Fas=Bio::SeqIO->new(-file=>"$fasta_f",-format=>'fasta');
    
29. while(my $seq=$Fas->next_seq()){
    
30.         $Data->{$seq->id}=$seq;
    
31. }
    
32. $Fas->close();
    
33. 
    
34. #提取Score值大于5的位置对应的碱基
    
35. my $Seqs={};
    
36. foreach my $point(sort {$a<=>$b}keys %{$Scores}){
    
37.         if($Scores->{$point} > 5){
    
38.                 foreach my $seq_id(keys %{$Data}){
    
39.                         $Seqs->{$seq_id}.=$Data->{$seq_id}->subseq($point,$point) if $point <= $Data->{$seq_id}->length();
    
40.                 }
    
41.         }
    
42. }
    
43. 
    
44. #输出数据
    
45. my $Out=Bio::SeqIO->new(-file=>">out.fasta",-format=>'fasta');
    
46. foreach my $seq_id(keys %{$Data}){
    
47.         print "Now is:\t$seq_id\n";
    
48.         $Data->{$seq_id}->seq($Seqs->{$seq_id});
    
49.         $Out->write_seq($Data->{$seq_id});
    
50. }
    
51. $Out->close();
    
52. 
    
53. print "Finished!\n";

复制代码

Data.rar (694 Bytes, 下载次数: 13)

2015-10-08 19:28 上传

点击文件名下载附件
数据文件

56836430

万分感谢！！！我再看看