学渣求学霸帮忙改改作业。。。

pitonas

{:2_172:} 大神! 这个, 内存占用小了。

1. #!/usr/bin/perl
   
2. 
   
3. my %code = (
   
4.     'TCA' => 'S',    # Serine
   
5.     'TCC' => 'S',    # Serine
   
6.     'TCG' => 'S',    # Serine
   
7.     'TCT' => 'S',    # Serine
   
8.     'TTC' => 'F',    # Phenylalanine
   
9.     'TTT' => 'F',    # Phenylalanine
   
10.     'TTA' => 'L',    # Leucine
    
11.     'TTG' => 'L',    # Leucine
    
12.     'TAC' => 'Y',    # Tyrosine
    
13.     'TAT' => 'Y',    # Tyrosine
    
14.     'TAA' => '_',    # Stop
    
15.     'TAG' => '_',    # Stop
    
16.     'TGC' => 'C',    # Cysteine
    
17.     'TGT' => 'C',    # Cysteine
    
18.     'TGA' => '_',    # Stop
    
19.     'TGG' => 'W',    # Tryptophan
    
20.     'CTA' => 'L',    # Leucine
    
21.     'CTC' => 'L',    # Leucine
    
22.     'CTG' => 'L',    # Leucine
    
23.     'CTT' => 'L',    # Leucine
    
24.     'CCA' => 'P',    # Proline
    
25.     'CCC' => 'P',    # Proline
    
26.     'CCG' => 'P',    # Proline
    
27.     'CCT' => 'P',    # Proline
    
28.     'CAC' => 'H',    # Histidine
    
29.     'CAT' => 'H',    # Histidine
    
30.     'CAA' => 'Q',    # Glutamine
    
31.     'CAG' => 'Q',    # Glutamine
    
32.     'CGA' => 'R',    # Arginine
    
33.     'CGC' => 'R',    # Arginine
    
34.     'CGG' => 'R',    # Arginine
    
35.     'CGT' => 'R',    # Arginine
    
36.     'ATA' => 'I',    # Isoleucine
    
37.     'ATC' => 'I',    # Isoleucine
    
38.     'ATT' => 'I',    # Isoleucine
    
39.     'ATG' => 'M',    # Methionine
    
40.     'ACA' => 'T',    # Threonine
    
41.     'ACC' => 'T',    # Threonine
    
42.     'ACG' => 'T',    # Threonine
    
43.     'ACT' => 'T',    # Threonine
    
44.     'AAC' => 'N',    # Asparagine
    
45.     'AAT' => 'N',    # Asparagine
    
46.     'AAA' => 'K',    # Lysine
    
47.     'AAG' => 'K',    # Lysine
    
48.     'AGC' => 'S',    # Serine
    
49.     'AGT' => 'S',    # Serine
    
50.     'AGA' => 'R',    # Arginine
    
51.     'AGG' => 'R',    # Arginine
    
52.     'GTA' => 'V',    # Valine
    
53.     'GTC' => 'V',    # Valine
    
54.     'GTG' => 'V',    # Valine
    
55.     'GTT' => 'V',    # Valine
    
56.     'GCA' => 'A',    # Alanine
    
57.     'GCC' => 'A',    # Alanine
    
58.     'GCG' => 'A',    # Alanine
    
59.     'GCT' => 'A',    # Alanine
    
60.     'GAC' => 'D',    # Aspartic Acid
    
61.     'GAT' => 'D',    # Aspartic Acid
    
62.     'GAA' => 'E',    # Glutamic Acid
    
63.     'GAG' => 'E',    # Glutamic Acid
    
64.     'GGA' => 'G',    # Glycine
    
65.     'GGC' => 'G',    # Glycine
    
66.     'GGG' => 'G',    # Glycine
    
67.     'GGT' => 'G',    # Glycine
    
68. );
    
69. # @ARGV = qw/sca cd/;
    
70. my ( $one, $two ) = @ARGV;
    
71. open my ($TWO), $two;
    
72. my %cds;
    
73. 
    
74. while (<$TWO>) {
    
75.     my ( $chr, $s, $e, $p ) = ( split /\s+/, $_, 6 )[ 0, 3, 4, 2 ];
    
76.     next if $p ne 'CDS';
    
77.     my $l = $e - $s + 1;
    
78.     push @{ $cds{$chr} }, [ $s, $l, $_ ];
    
79. }
    
80. 
    
81. open OUT, '>', 'Sequence_extraction_part.fa' or die "don't worry";
    
82. open AA,  '>', 'aa_extraction_part.fa'       or die "be happy";
    
83. open my ($ONE), $one;
    
84. 
    
85. $/ = '>';
    
86. <$ONE>;
    
87. 
    
88. while (<$ONE>) {
    
89.     my ( $key, @seq ) = split /\n/;
    
90.     ($key) = $key =~ /^(\w+)/;
    
91.     next unless exists $cds{$key};
    
92.     my $seq = join '', @seq;
    
93.     chop $seq if $seq =~ />$/;
    
94.     for my $k ( @{ $cds{$key} } ) {
    
95.         extract( $seq, @$k );
    
96.     }
    
97. 
    
98. }
    
99. 
    
100. close OUT;
     
101. close AA;
     
102. 
     
103. sub extract {
     
104.     my ( $SEQ, $start, $len, $m ) = @_;
     
105.     my $Seq = substr $SEQ, $start - 1, $len;
     
106.     my $aa;
     
107.     while ( $Seq =~ /(.{3})/g ) {
     
108.         $aa .= $code{$1} || bad_codon($1);
     
109.     }
     
110.     my ( @OUT, @AA );
     
111.     my ( $lenseq, $lenaa ) = ( length($Seq), length($aa) );
     
112.     for ( my $j = 0 ; $j * 60 < $lenseq ; $j++ ) {
     
113.         push @OUT, substr( $Seq, $j * 60, 60 ), "\n";
     
114.     }
     
115.     for ( my $i = 0 ; $i * 60 < $lenaa ; $i++ ) {
     
116.         push @AA, substr( $aa, $i * 60, 60 ), "\n";
     
117.     }
     
118.     print OUT ">$m", @OUT;
     
119.     print AA ">$m",  @AA;
     
120. }
     
121. 
     
122. sub bad_codon {
     
123.     my $bad = shift;
     
124.     print qq[Bad codon "$bad"!!\n];
     
125.     exit;
     
126. }
     

复制代码

回复 10# huang6894


   

huang6894

回复 11# pitonas


    吐司大神。。。。 谢谢，谢谢

huang6894

回复 11# pitonas


   吐司大神。。。我做测试的时候。。。发现一个问题：
不知道我的操作哪错了。。。

我的输入文件（第二个文本）：

1. 
   
2. 》head -3 Prunus_mume_scaffold.gff
   
3. scaffold101     GLEAN   mRNA    1862    8809    0.997087        -       .       ID=Pm020658;
   
4. scaffold101     GLEAN   CDS     1862    2248    .       -       0       Parent=Pm020658;
   
5. scaffold101     GLEAN   CDS     2709    2861    .       -       0       Parent=Pm020658;
   
6. 》cut -f1 Prunus_mume_scaffold.gff|sort|uniq|wc -l
   
7. 41
   
8. 》 grep 'scaffold' Prunus_mume_scaffold.gff|grep CDS|wc -l
   
9. 47536
   

复制代码

我的输入文件（第一个文本）：

1. 
   
2. 格式：
   
3. >scaffold16 36.5
   
4. AATTTTCGCTGGGTTGGGTGTGTCAATTTTTGGCCCAAACACAAACCCACCGTGAGGGTACCCCCGAGTTGCCTTGAGACAGTCAGTTATTTTTCCTCACGAGTTTCGGGGACGCCTGAACCGTGTAGTGCCAGGTTGCGGCTCGGATTGAGATGTG
   
5. >scaffold35 33.2
   
6. GTGTCTTTATATACATATAT
   
7. 》grep scaffold plum_0630.scafSeq.FG|wc -l
   
8. 41
   

复制代码

结果文件：

1. >未优化程序：
   
2. 》grep '>scaffold' Sequence_extraction_part.fa|wc -l
   
3. 47536
   
4. >大神优化后的程序一：
   
5. 》grep '>' Sequence_extraction_part.fa|wc -l
   
6. 1723
   
7. 大神优化的程序二“
   
8. 》grep '>scaffold' Sequence_extraction_part.fa|wc -l
   
9. 427

复制代码

按道理来说最终结果应该是生成47536行的结果的呀，看了一下结果文件，大神的两个脚本第一个脚本只存在输入的第二个文本第一行建立的hash，第二个程序只存在输入的第一个文本的前两段建立的hash。。。大神能不能指导我一下呢？

afukada

雖然黑吐司大神有解了

還是請問一下

你想做的事情是不是把序列(sequence)取出來之後

把CDS的部份給截選出來

再將他轉錄成胺基酸(amino-acid)

huang6894

回复 14# afukada


    是的。。。就这么简单

pitonas

{:2_172:} 大神，测试下

1. #!/usr/bin/perl -w
   
2. 
   
3. my $two = 'Prunus_mume_scaffold.gff';
   
4. open my ($TWO), $two;
   
5. my %cds;
   
6. 
   
7. while (<$TWO>) {
   
8.     my ( $chr, $s, $e, $p ) = ( split /\s+/, $_, 6 )[ 0, 3, 4, 2 ];
   
9.     next if $p ne 'CDS';
   
10.     my $l = $e - $s + 1;
    
11.     push @{ $cds{$chr} }, [ $s, $l, $_ ];
    
12. }
    
13. 
    
14. my $cdssum;
    
15. $cdssum += $_ for map { scalar @$_ } values %cds;
    
16. 
    
17. print 'scaffold count = ', scalar keys %cds, "\n";
    
18. print 'cds count = ', $cdssum, "\n";

复制代码

回复 13# huang6894


   

huang6894

回复 16# pitonas


    => perl hh.pl
scaffold count = 41
cds count = 47536

然后咧？

pitonas

{:2_172:} 大神
sub bad_codon {
    my $bad = shift;
    print qq[Bad codon "$bad"!!\n];
    exit;   # will exit the script unless $code{$1}
}
108:

1.         $aa .= $code{$1} || bad_codon($1);

复制代码

exit the script unless $code{$1}
因为在中间退出,所以它不能打印完成
$aa .= $code{$1} || 'X';    # X => NO KEY IN %code;

回复 17# huang6894


   

pitonas

{:2_172:} 大神的代码
145, 146, 147

1.        return  "Bad codon \"$codon\"!!\n";  # return 返回
   
2.                 exit; # 因为返回return 所以它永远不会执行这条 exit
   
3.                 }
   
4. }

复制代码

huang6894

回复 18# pitonas


    ( ⊙o⊙ )哦，sorry，sorry。。。我没仔细研究好，问题解决了，谢谢大大谢谢大大