求助：关于核苷酸和氨基酸筛选的问题（已解决，感谢！）

56836430

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

>lcl|Abi_c31169_g3_i1_m.75385 unnamed protein product
FRVSGAPSGKTKLGFRFKIAK-----PHKRSRETLLKEAGAEARGEDKAQVGRMTTRLKK
NRKKRGHVSAGHGRIGKHRKHPGGRGNAGGMHHHRILFDKYHPGYFGKVGMRHFHKLRNH
YHCPIVNIEKLWSLIPEDVKAKAT----AESAPLVDVTQFGYFKVLGKGELP-KQPIVVK
AKLISKLAERKIKEAGGAVMLTA
>lcl|Aco_Aquca_011_00529.1 unnamed protein product
-----------------------------------------------------MTTRFKK
NRKKRGHVSAGHGRIGKHRKHPGGRGNAGGMHHHRILFDKYHPG--------YFHKLRNK
FYCPTVNIDKLWSLVPQEIKDKANTT---GNAPLIDVTQLGYFKVLGKGMLPPSQPVVVK
AKLVSKNAEKKIKEAGGAVLLTA
>lcl|Ame_c52237_g2_i1_m.69793 unnamed protein product
-----------------------------------------------------MTTRLKK
HRKKRGHVSAGHGRIGKHRKHPGGRGNAGGLHHHRILFDKYHPGYFGKVGMRHFHKLKNH
YFCPIINIDKLWSLVPEDIKKKAT----ADNVPLLDVTHHGYFKVLGKGELP-PQPLVVK
AKLISKLAEKKIKEAGGAVLLTA

第二个是align好的核苷酸序列文件，fasta格式，如(b.fasta)：
>lcl|Abi_c31169_g3_i1_m.75385
TTTAGGGTTTCCGGCGCCCCGTCAGGGAAAACCAAATTAGGGTTTCGATTTAAAATCGCA
AAA---------------CCCCACAAGCGCTCGAGAGAAACCCTGTTGAAGGAGGCCGGG
GCAGAGGCGAGAGGGGAAGACAAGGCACAGGTTGGCAGGATGACGACCAGGTTGAAGAAG
AACAGGAAGAAGAGAGGGCACGTGAGCGCTGGCCACGGCCGTATCGGCAAGCACAGGAAG
CATCCGGGTGGTAGGGGTAACGCCGGAGGTATGCACCACCATCGCATTCTCTTCGATAAG
TACCATCCCGGTTACTTCGGTAAGGTTGGTATGAGGCATTTCCACAAGCTAAGGAATCAC
TATCACTGCCCTATTGTTAATATTGAAAAACTGTGGTCTCTGATTCCCGAGGATGTCAAG
GCCAAGGCCACT------------GCTGAGAGTGCCCCGCTCGTAGACGTTACTCAGTTT
GGGTATTTCAAGGTCTTGGGCAAGGGCGAGCTGCCC---AAGCAACCTATCGTAGTGAAG
GCCAAGCTTATCTCTAAGCTGGCCGAGAGGAAGATCAAGGAGGCCGGAGGGGCAGTCATG
TTGACTGCT
>lcl|Aco_Aquca_011_00529.1
------------------------------------------------------------
------------------------------------------------------------
---------------------------------------ATGACAACTCGATTCAAGAAG
AACAGGAAGAAGAGAGGACACGTGTCAGCTGGACATGGACGAATTGGAAAGCATAGAAAG
CATCCAGGAGGAAGAGGAAATGCAGGAGGAATGCATCATCATCGAATATTGTTTGATAAA
TACCATCCAGGG------------------------TATTTTCATAAGCTTCGTAATAAG
TTTTATTGCCCAACTGTTAATATTGATAAGTTATGGTCTTTGGTTCCTCAAGAAATCAAA
GATAAAGCTAATACTACC---------GGAAATGCCCCTTTAATTGATGTTACTCAATTA
GGGTATTTTAAGGTTTTGGGTAAGGGTATGTTACCACCTTCTCAACCTGTTGTTGTTAAG
GCTAAACTTGTTTCAAAGAATGCTGAGAAGAAGATTAAGGAAGCTGGTGGTGCTGTTCTA
CTTACTGCT
>lcl|Ame_c52237_g2_i1_m.69793
------------------------------------------------------------
------------------------------------------------------------
---------------------------------------ATGACGACTCGATTGAAGAAG
CACAGGAAGAAGAGAGGCCATGTGAGTGCAGGGCATGGTCGTATTGGAAAGCACAGGAAG
CATCCAGGAGGCAGGGGTAACGCCGGAGGTCTGCACCACCACCGCATTCTCTTTGATAAG
TACCATCCTGGATACTTCGGCAAGGTCGGTATGCGCCATTTCCACAAGCTCAAAAATCAC
TATTTCTGCCCCATTATTAATATCGACAAGCTCTGGTCTCTCGTCCCTGAAGATATCAAG
AAGAAAGCCACA------------GCTGATAATGTCCCCCTCTTAGATGTTACGCATCAC
GGCTATTTTAAGGTTTTGGGGAAGGGAGAATTGCCT---CCCCAGCCACTTGTTGTCAAG
GCCAAGCTTATTTCTAAGCTGGCCGAGAAGAAGATCAAGGAGGCAGGTGGGGCTGTTCTG
TTGACCGCT


另外还有一个文件是以上序列每个位点的score (c.txt):
0.000000
0.000000
0.000000
0.004008
0.003777
0.003657
0.003628
0.002800
0.000173
0.041396
0.041087
0.040902
0.040653
0.039747
0.038890
0.053296
0.051119
0.050429
0.042738
0.014122
0.011783
0.001929
0.007066
0.007922
0.008270
0.009295
0.020800
0.021852
0.022518
0.022070
0.022183
0.023023
0.022791
0.077558
0.099258
0.083209
0.034202
0.032664
0.032918
0.034097
0.035328
0.034721
0.135512
0.151489
0.163670
0.188523
0.245695
0.266098
0.285722
0.315185
0.363380
0.447316
0.541042
9.453623
9.482589
9.489742
9.493266
9.493580
9.493349
9.493314
9.493557
9.493586
9.493290
9.493337
9.493784
9.494013
9.494032
9.493960
9.493816
9.493815
9.493960
9.494031
9.494023
9.494006
9.494001
9.494003
9.494000
9.494005
9.493925
9.493927
9.494025
9.494033
9.493800
9.493784
9.493993
9.493990
9.493969
9.493882
9.493681
9.493720
9.996930
9.999200
9.999361
9.999575
9.999882
9.999759
9.999726
9.999896
9.999944
9.999930
9.999947
9.999947
9.999588
9.990940
9.417943
9.417428
9.417420
9.417448
9.417453
9.417495
9.417933
9.417899
9.526410
9.739041
9.951346
9.958738
9.981691
9.991584
9.995957
9.995432
9.995120
9.995554
9.999411
9.999811
9.998937
9.998946
9.999780
9.999766
9.999614
9.999510
9.999369
9.997099
9.947615
9.808741
9.075074
9.056618
8.150498
7.975688
6.985375
6.817541
6.363284
5.817629
5.532721
5.915577
5.878241
5.586190
5.512254
5.551737
5.988903
6.401943
7.864917
8.460415
8.868447
9.027906
9.056516
9.081413
9.082546
9.083229
9.083782
9.492370
9.572677
9.579134
9.579836
9.580672
9.820549
9.975557
9.994984
9.994225
9.974879
9.645265
8.934859
8.587425
2.708720
7.126993
9.018083
9.926556
9.970963
9.982491
9.990430
9.998740
9.999582
9.999713
9.999266
9.999191
9.999584
9.999491
9.999277
9.999483
9.999600
9.999212
9.998949
9.998173
9.991413
9.238468
9.238482
9.238677
9.238726
9.238749
9.237899
9.237364
9.238038
9.238764
9.238947


想根据score大小来筛选序列，当score大于5时，碱基得以保留，得到以下输出文件d.fasta和e.fasta：
d.fasta
>lcl|Abi_c31169_g3_i1_m.75385 unnamed protein product
MTTRLKKNRKKRGHVSAGHGRIGKHRKHPGGRGNAGGMHHHRILFDKYHP
GYFGKVGMRHFHKLRNHYHCPIVNIEKLWSLIPEDVKAKAT----AESAP
LVDVTQFGYFKVLGKGELP-KQPIVVKAKLISKLAERKIKEAGGAVMLTA

>lcl|Aco_Aquca_011_00529.1 unnamed protein product
MTTRFKKNRKKRGHVSAGHGRIGKHRKHPGGRGNAGGMHHHRILFDKYHP
G--------YFHKLRNKFYCPTVNIDKLWSLVPQEIKDKANTT---GNAP
LIDVTQLGYFKVLGKGMLPPSQPVVVKAKLVSKNAEKKIKEAGGAVLLTA

>lcl|Ame_c52237_g2_i1_m.69793 unnamed protein product
MTTRLKKHRKKRGHVSAGHGRIGKHRKHPGGRGNAGGLHHHRILFDKYHP
GYFGKVGMRHFHKLKNHYFCPIINIDKLWSLVPEDIKKKAT----ADNVP
LLDVTHHGYFKVLGKGELP-PQPLVVKAKLISKLAEKKIKEAGGAVLLTA

e.fasta:
>lcl|Abi_c31169_g3_i1_m.75385
ATGACGACCAGGTTGAAGAAGAACAGGAAGAAGAGAGGGCACGTGAGCGC
TGGCCACGGCCGTATCGGCAAGCACAGGAAGCATCCGGGTGGTAGGGGTA
ACGCCGGAGGTATGCACCACCATCGCATTCTCTTCGATAAGTACCATCCC
GGTTACTTCGGTAAGGTTGGTATGAGGCATTTCCACAAGCTAAGGAATCA
CTATCACTGCCCTATTGTTAATATTGAAAAACTGTGGTCTCTGATTCCCG
AGGATGTCAAGGCCAAGGCCACT------------GCTGAGAGTGCCCCG
CTCGTAGACGTTACTCAGTTTGGGTATTTCAAGGTCTTGGGCAAGGGCGA
GCTGCCC---AAGCAACCTATCGTAGTGAAGGCCAAGCTTATCTCTAAGC
TGGCCGAGAGGAAGATCAAGGAGGCCGGAGGGGCAGTCATGTTGACTGCT

>lcl|Aco_Aquca_011_00529.1
ATGACAACTCGATTCAAGAAGAACAGGAAGAAGAGAGGACACGTGTCAGC
TGGACATGGACGAATTGGAAAGCATAGAAAGCATCCAGGAGGAAGAGGAA
ATGCAGGAGGAATGCATCATCATCGAATATTGTTTGATAAATACCATCCA
GGG------------------------TATTTTCATAAGCTTCGTAATAA
GTTTTATTGCCCAACTGTTAATATTGATAAGTTATGGTCTTTGGTTCCTC
AAGAAATCAAAGATAAAGCTAATACTACC---------GGAAATGCCCCT
TTAATTGATGTTACTCAATTAGGGTATTTTAAGGTTTTGGGTAAGGGTAT
GTTACCACCTTCTCAACCTGTTGTTGTTAAGGCTAAACTTGTTTCAAAGA
ATGCTGAGAAGAAGATTAAGGAAGCTGGTGGTGCTGTTCTACTTACTGCT

>lcl|Ame_c52237_g2_i1_m.69793
ATGACGACTCGATTGAAGAAGCACAGGAAGAAGAGAGGCCATGTGAGTGC
AGGGCATGGTCGTATTGGAAAGCACAGGAAGCATCCAGGAGGCAGGGGTA
ACGCCGGAGGTCTGCACCACCACCGCATTCTCTTTGATAAGTACCATCCT
GGATACTTCGGCAAGGTCGGTATGCGCCATTTCCACAAGCTCAAAAATCA
CTATTTCTGCCCCATTATTAATATCGACAAGCTCTGGTCTCTCGTCCCTG
AAGATATCAAGAAGAAAGCCACA------------GCTGATAATGTCCCC
CTCTTAGATGTTACGCATCACGGCTATTTTAAGGTTTTGGGGAAGGGAGA
ATTGCCT---CCCCAGCCACTTGTTGTCAAGGCCAAGCTTATTTCTAAGC
TGGCCGAGAAGAAGATCAAGGAGGCAGGTGGGGCTGTTCTGTTGACCGCT



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

56836430

在网友的帮助下，现在能解决氨基酸的筛选问题，但是核苷酸的问题尚未解决。以下是解决氨基酸的文本：
#!/usr/bin/perl
use strict;
use Bio::SeqIO;

#限制程序参数
unless (@ARGV==2){
        print "Usage:\t$0\tScore_file\tFasta_file\n";
        sleep 3;
}
my $score_f=shift@ARGV;
my $fasta_f=shift@ARGV;

#读取Score数据
my $Scores={};
open (my $sc,"<$score_f") or die "Open score file failed!\n";
my $n=0;
while (<$sc>){
        chomp;
        if($_){
                $n++;
                $Scores->{$n}=$_;
        }
}
close $sc;

#读取Seq信息
my $Data={};
my $Fas=Bio::SeqIO->new(-file=>"$fasta_f",-format=>'fasta');
while(my $seq=$Fas->next_seq()){
        $Data->{$seq->id}=$seq;
}
$Fas->close();

#提取Score值大于5的位置对应的碱基
my $Seqs={};
foreach my $point(sort {$a<=>$b}keys %{$Scores}){
        if($Scores->{$point} >= 5){
                foreach my $seq_id(keys %{$Data}){
                        $Seqs->{$seq_id}.=$Data->{$seq_id}->subseq($point,$point) if $point <= $Data->{$seq_id}->length();
                }
        }
}

#输出数据
my $Out=Bio::SeqIO->new(-file=>">${fasta_f}_out.fasta",-format=>'fasta');
foreach my $seq_id(keys %{$Data}){
        print "Now is:\t$seq_id\n";
        $Data->{$seq_id}->seq($Seqs->{$seq_id});
        $Out->write_seq($Data->{$seq_id});
}
$Out->close();

print "Finished!\n";

求帮忙解决怎么筛选核苷酸

hztj2005

核苷酸序列文件中应该是三个字母算一个密码子？那么第二个文件中的点位是否是：
TTT-AGG-GTT-TCC-GGC-GCC-CCG-TCA-GGG-AAA-ACC-AAA-TTA-GGG-TTT-CGA-TTT-AAA-ATC-GCA-
AAA
如何与score.txt对应？

第一段中，氨基酸序列中如何算一个点位？又如何与score.txt对应？
FRVSGAPSGKTKLGFRFKIAK-----PHKRSRETLLKEAGAEARGEDKAQVGRMTTRLKK
NRKKRGHVSAGHGRIGKHRKHPGGRGNAGGMHHHRILFDKYHPGYFGKVGMRHFHKLRNH
YHCPIVNIEKLWSLIPEDVKAKAT----AESAPLVDVTQFGYFKVLGKGELP-KQPIVVK
AKLISKLAERKIKEAGGAVMLTA

hztj2005

score.txt只有203个值，而上面每个片段有609个字母或-，看来氨基酸序列也是三个字母一个点位了。
那么2楼的代码应该也可以处理氨基酸序列，为何还要代码呢?

56836430

回复 4# hztj2005


    氨基酸是一个位点一个值，所以能用。但是核苷酸是三个字母一个值，不知道怎么处理？万分感谢~~~

hztj2005

回复 5# 56836430

你是生物专业的，在表达问题时要让非生物专业的人能理解。

我关心的问题：“氨基酸是一个位点一个值”，你这里“位点”是指一个字母吗？那么上面的氨基酸片段有609个字母，而score.txt只有203个值，两者之间是如何对应的？你要把这个问题解释清楚。

至于“核苷酸是三个字母一个值”，那么609个字母构成203个“位点”，刚好对应score.txt只有203个值，反而容易处理。


   

56836430

回复 6# hztj2005


    每三个核苷酸组成一个氨基酸。score值（c文件）对应的是氨基酸的字母（a文件），也就是一个字母（203个字母）一个值（203个值 ）。对核苷酸， 如三楼所示，每三个字母对应一个值，即609个字母（b文件）对应203个值（c文件）。我现在需要处理的就是怎么对核苷酸进行处理。例如，当第一个score值小于5时，删除核苷酸前三个字母，以此类推。叩谢~

hztj2005

1. #!/usr/bin/perl
   
2. 
   
3. open($score_f, "score.txt");#打开文件
   
4. open($fasta_f, "Nucleotide.txt");#打开Nucleotide文件
   
5. open(Out, ">Nucleotidenew.txt");#输出结果文件
   
6. 
   
7. 
   
8. #读入score.txt
   
9. @Scores={};
   
10. while (<$score_f>){
    
11.          chomp;
    
12.          #print $_;
    
13.          if($_){
    
14.             $n++;
    
15.             push  @Scores,$_;
    
16.                
    
17.          }
    
18. }
    
19. 
    
20. @Scorestemp = splice (@Scores,1,203); #删除第一个原因不明的非法元素。
    
21. @Scores=@Scorestemp;
    
22. 
    
23. $Scorenum=@Scores;
    
24. print Out        ("$Scorenum\n");
    
25. 
    
26. #print Out        ("@Scores\n");
    
27. 
    
28. #读取Nucleotide.txt
    
29. 
    
30. my $fastaall="";
    
31. my $fastatemp="";
    
32. while (<$fasta_f>){
    
33.          chomp;
    
34.          $fastatemp = $fastaall.$_;
    
35.          $fastaall= $fastatemp ;         
    
36.          #print Out        $fastaall."\n";   
    
37.         
    
38. }
    
39. 
    
40. $fastaall =~ s/\s//g;#删除不当空白
    
41. $fastaall =~ s/(([-A-Z]){3})/$1 /g;#给每个位点添加一个空白
    
42. @Datafasta=split(/ /,$fastaall);#按空白切分
    
43. 
    
44. #检查读入数据分组是否正确
    
45. $Datafastanum=@Datafasta;
    
46. print Out        "\r\n";
    
47. print Out        ("Datafastanu:$Datafastanum\r\n");
    
48. $n=1;
    
49. foreach $item (@Datafasta) {
    
50.      print Out ("$n:$item\r\n");
    
51.      $n++;
    
52. }
    
53. 
    
54. 
    
55. #提取Score值大于5的位置对应的碱基
    
56. @Datanew ={};
    
57. for ($i = 0; $i <= @Datafasta; $i++) {
    
58. 
    
59.     if($Scores[$i] >= 5)
    
60.     {
    
61.             #print $Datafasta[$i]."\n";
    
62.             push  @Datanew, $Datafasta[$i];
    
63.            
    
64.     }
    
65. }
    
66. 
    
67. @Datanewtemp = splice (@Datanew,1,203); #删除第一个原因不明的非法元素。
    
68. @Datanew=@Datanewtemp;
    
69. 
    
70. #分组结果
    
71. print Out        ("@Datanew\n");
    
72. 
    
73. #最终结果
    
74. $result=join("",@Datanew);
    
75. print Out        "$result\n";
    
76. 
    
77. exit;
    
78. 
    
79. #打印最终结果个数：
    
80. $Datanewnum=@Datanew;
    
81. print Out        ("$Datanewnum\n");
    
82. 
    
83. exit;
    
84. 
    
85. ---------------------------------------------
    
86. 对perl不熟悉，上面代码写得不好，改一改：
    
87. #!/usr/bin/perl
    
88. 
    
89. open($score_f, "score.txt");#打开文件
    
90. open($fasta_f, "Nucleotide.txt");#打开Nucleotide文件
    
91. open(Out, ">Nucleotidenew.txt");#输出结果文件
    
92. 
    
93. #读入score.txt 并删除多余的空格
    
94. my @Scores = map {chomp;s/^ +//;$_} <$score_f>;
    
95. 
    
96. #读取Nucleotide.txt
    
97. my $fastaall="";
    
98. while (<$fasta_f>){
    
99.      chomp;
    
100.      $fastaall = $fastaall.$_;
     
101. }
     
102. 
     
103. $fastaall =~ s/\s//g;#删除不当空白
     
104. @Datafasta = $fastaall =~ m/\G[-A-Z]{3}/g;#按点位分割
     
105. 
     
106. #检查读入数据分组是否正确
     
107. $Datafastanum=@Datafasta;
     
108. print Out        ("Datafastanu:$Datafastanum\r\n");
     
109. $n=1;
     
110. foreach $item (@Datafasta) {
     
111.      print Out ("$n:$item\r\n");
     
112.      $n++;
     
113. }
     
114. 
     
115. #提取Score值大于5的位置对应的碱基
     
116. $i = -1;
     
117. @Datanew = grep {$i++;($Scores[$i] >= 5)}@Datafasta;
     
118. 
     
119. #分组结果
     
120. print Out        "@Datanew\n";
     
121. #最终结果
     
122. print Out        @Datanew;
     
123. exit;
     

复制代码

rubyish

回复 1# 56836430

c: total 203, [ score >= 5:  count = 149 ]
why? nide d length = 150 ?
why? nide e length = 450 ?

wode:
d: length 149
e: length 447

run: [ perl version >= 5.010 ]
perl abc.pl  a  b  c  d  e

1. #!/usr/bin/perl
   
2. use 5.010;
   
3. my $score = 5;
   
4. my @rw    = qw< < < < > > > >;
   
5. my ( $A, $B, $C, $D, $E ) = map { open my $f, shift @rw, $_; $f } @ARGV;
   
6. my @pick = map { state $i = -1; $i++; $_ >= $score ? $i : () } <$C>;
   
7. 
   
8. for ( [ $A, $D, 1 ], [ $B, $E, 3 ] ) {
   
9.     my ( $I, $O, $n ) = @$_;
   
10.     local $/ = '>'; <$I>;
    
11. 
    
12.     while (<$I>) {
    
13.         my ( $a, $b ) = split "\n", $_, 2;
    
14.         $b =~ tr/\n>//d;
    
15.         my @b = unpack "(A$n)*", $b;
    
16.         say $O '>', $a;
    
17.         say $O $_ for unpack '(A50)*', join '', @b[@pick];
    
18.         say $O '';
    
19.     }
    
20. }
    

复制代码

RE_HASH

open C, "c.txt";map {$T .= ($_ > 5)? "Y": "N";  } <C>; close C;#Read score

#make a RE template
$T =~ s/(N+)/{sprintf "(?:.{%d})", length $1} /ge;  $T =~ s/(Y+)/{$C++; $R .= "\$$C"; sprintf "(.{%d})", length $1; }/ge;

sub CutFile {
  my ($F, $T, $R) = @_;
  open A , $F;        local $_ = do {$/=undef; <A>}; s/(?<=\S{40})\n//g;  close $F; #Read File
  s/$T/qq(qq($R))/gee; s/(\S{50})/\1\n/g; #Cut string
  return $_;
}

printf "d.fasta\n%s", CutFile ("a.fasta", $T, $R) ;  
print "\n\n";
$T =~ s/(\d+)/{$1 * 3}/eg;  #Change Template X3
printf "e.fasta\n%s", CutFile ("b.fasta", $T, $R) ;