哪位幫我看看，bioperl處理fastq文件

super_two

這是我寫的去除質(zhì)量值低于8的堿基占reads總長度40%的reads序列的腳本，誰能幫我看看有什么問題。謝謝~

1. #！/usr/bin/perl
   
2. use warnings;
   
3. use strict;
   
4. use Bio::SeqIO::fastq;
   
5. 
   
6. my $infile=$ARGV[0];
   
7. open (W,"<$infile") || die $!;
   
8. my %fileseq;
   
9. my ($id,$seq,$quality);
   
10. while(<W>){
    
11. if (/^(\S+)\s+(\S+)$/){
    
12. $fileseq{$id}{$seq}=$quality;
    
13.         }
    
14. }
    
15. close W;
    
16. 
    
17. open (M,">$infile.out")||die $!;
    
18. my $all=Bio::SeqIO::fastq->new(-file=>"$infile",-format=>'fastq');
    
19. while(my $seq=$all->next_seq()){
    
20. my $id=$seq->id;
    
21. my $seq=$seq->seq;
    
22. my $qual=join('',@{$seq->$qual});
    
23. my @quality=split(//,$qual);
    
24. 
    
25. my @quality1;
    
26. for my$a (@quality){
    
27. 
    
28.         my$b=chr($i-33);
    
29.         push (@quality1,$b);
    
30.         }
    
31. $qual=join ('',@quality1);
    
32.        
    
33.         my $base_8=0;
    
34.         foreach my$base (@quality1){
    
35.         if($base<=8){
    
36.         my$tatal_low=$base_8++;}
    
37.         my $total_quality=$fileseq{$id}{$seq}++;
    
38.         if($tatal_low/$total_quality>=0.4){
    
39.         s/$fileseq{$id}//;
    
40.         s/$id//g;
    
41.         s/$fileseq{$id}{$seq}//g;
    
42.                         }
    
43.                 }
    
44.        
    
45.         print M ("$id\n$seq\n\+\n$quality");
    
46.         }
    
47.         close M;
    
48.        
    
49.        
    
50.        

復(fù)制代碼

b114213903

1. #!/usr/bin/perl
   
2. use strict;
   
3. use Bio::SeqIO;
   
4. 
   
5. my $IN=Bio::SeqIO->new(-file=>"in.fastq",-format=>"fastq");
   
6. my $OUT=Bio::SeqIO->new(-file=>">out.fastq",-format=>'fastq');
   
7. while(my $seq=$IN->next_seq()){
   
8.         my $id=$seq->id;
   
9.         print "Now is:\t$id\n";
   
10.        
    
11.         my $len=$seq->length();
    
12.         my $quals=$seq->qual();
    
13.        
    
14.         my $n=0;
    
15.         foreach my $qual(@{$quals}){
    
16.                 if($qual<8){
    
17.                         $n++;
    
18.                 }
    
19.         }
    
20.         unless ($n/$len >= 0.4){
    
21.                 $OUT->write_seq($seq);
    
22.         }
    
23. }
    
24. $OUT->close();
    
25. $IN->close();
    
26. 
    
27. print "Finished!\n";

復(fù)制代碼

super_two

我的文件是這樣的，那個質(zhì)量值是ASCII碼，不用換算過來嗎？illuminate1.8+，文件如下

1. @SRR1.1 HWI-ST1:337:D26KUACXX:1:1101:1351:1997 length=51
   
2. CCTCTTTTGCCGTGGACCCCACTTGATCTGTCTCTTATACACATCTAGATG
   
3. +SRR1.1 HWI-ST1:337:D26KUACXX:1:1101:1351:1997 length=51
   
4. ??8B?;,=AFFF1EF1CEFE8FG@F3:?CGDG4?*B0D*B9F39??BD*/B
   

復(fù)制代碼

回復(fù) 2# b114213903


   

b114213903

回復(fù) 3# super_two


    BioPerl模塊自動就處理了！

super_two

原來是這樣啊，謝謝


回復(fù) 4# b114213903


   

super_two

不能運(yùn)行成功，報錯是：Can't locate object method "length" via package "AAAAAATTTTTTAGATGAGCAGATTATGTGCTTCTTCCTCATTTACTTTTTTCGATGAA。。。。。�！�，這個$seq->length()不能用嗎？



回復(fù) 4# b114213903


   

b114213903

回復(fù) 6# super_two


    你是不是改腳本內(nèi)容了？我這兒都運(yùn)行正常！

super_two

只是把輸入文件my $infile=$ARGV[0];后面“

$IN=Bio::SeqIO->new(-file=>"in.fastq",-format=>"fastq");中in.fastq  相應(yīng)的變成了$infile.我的文件是.gz的壓縮文件

回復(fù) 7# b114213903


   

b114213903

回復(fù) 8# super_two


    不能直接讀壓縮文檔，只能讀fastq格式的文檔，可以通過管道，邊解壓邊處理。

gaoyang636

lz本著學(xué)習(xí)Perl，BioPerl的心態(tài)來寫的話，的確很好。但如果僅僅是想找個工具來處理數(shù)據(jù)，其實(shí)現(xiàn)在有大把現(xiàn)成的工具可以用哦，不但節(jié)省時間和精力，而且Performance還好的。比如你這個需求就可以使用Fastq toolkit完成。