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

huang6894

是这样的：
两个文本，第一个文本格式：

1. >scaffold16 36.5
   
2. AATTTTCGCTGGGTTGGGTGTGTCAATTTTTGGCCCAAACACAAACCCACCGTGAGGGTACCCCCGAGTTGCCTTGAGAC
   
3. >scaffold35 33.2
   
4. GTGTCTTTATATACATATATATGCCGGCTATACTAGATCCTCTTTTATTTTTCGAAAGAAATAGTCTAGCCGTTCGGCTA
   
5. 
   

复制代码

第二个文本格式：

1. scaffold16     GLEAN   mRNA    1862    8809    0.997087        -       .
   
2. scaffold16     GLEAN   CDS     1862    2248    .       -       0       Parent=P
   
3. scaffold16     GLEAN   CDS     2709    2861    .       -       0       Parent=P

复制代码

我现在想根据第二个文本提取第一个文本的字符串，建立哈希使第一个文本“>"行（key）对应下面非">"的字符串（value），然后根据哈希对第二个文本进行处理，第二个文本第一列是key，第四列是起始位置，第五列是终止位置，想得到的是第三列是CDS，对应value的起始位置到终止位置的字符串，存成一个文档后，根据指定的规则转换成新的字符串，每60个字符换行。。。
所以有了下面这个程序：

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

复制代码

可是非常非常慢。。对于第一个文本65M，第二个文本3.8M（57021行）的情况下，花了非常久的时间和内存。。。
所以。。。如果。。。大神们可以的话，可以提供一下优化帮助吗？
谢谢谢谢！

q1208c

看你的说明, 我进了雾里, 看你的代码, 我掉沟里了.

我现在终于明白有些不喜欢perl而喜欢python的原因了.

给楼主个建议, 第一, 尽可能少用或不用 $_ 这个变量. 虽然没什么不对, 但程序会更可读.
第二, 变量名不要太长, 且尽可能有点意义.
第三, 良好的缩进有助于代码的理解.
第四, 过于专业的内容, 请尽量用功能类似的简单的数据来解释, 对于非专业人士来说, 没多少人能理解你那一串串 "碱基".

huang6894

回复 2# q1208c


    好吧。。又让大神们迷惑了。。。。 语言是算命老师教的。。实在抱歉。。。
    其实后面的就是一个hash，我在说明里用“规则转换”定义了。。。sorry。。。and，谢谢大神的建议。。
   现在一个问题是，我真心不知道怎么优化它让它快一点内存小一点

q1208c

回复 3# huang6894

我没细看你的代码.

你试着把 在函数中那两个 OPEN  file 挪到外面. 在那个 while 之外打开这两个文件, 函数中只管写入. 在 while 之后 close 这两个 文件. 这样, 可能会好一点.


   

huang6894

回复 4# q1208c


    大神v5。。果然快了一点。。。但是内存还是占用太大了。。。而且，感觉处理这么小的文件不应该这么慢才对。。。谢谢大神的建议，谢谢谢谢

pitonas

{:2_172:}大神! 这个, 可能会快一点. 但是内存还是占用大了。而且，更大了。{:2_172:}{:2_172:}
不花内存, 内存小一点
scaffold16, scaffold35.... 存成一个, 一个文档
这样, 慢了. 非常久的时间

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

复制代码

q1208c

如果内存还是占用过大, 可能要考虑调整代码的逻辑了.

因为我没办法看懂楼主的实际需求, 所以, 没办法提出调整的建议.

pitonas

{:2_172:}  %seqId => big hash (文本65M)

回复 7# q1208c


   

huang6894

回复 7# q1208c


    sorry。。。其实我真心不知道怎么说需求。。。谢谢你

huang6894

回复 6# pitonas


    非常非常感谢土司大神，学到了一些新的东西，谢谢谢谢