哪位帮我看看，bioperl处理fastq文件

super_two

这是我写的去除质量值低于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.        

复制代码

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";

复制代码

super_two

我的文件是这样的，那个质量值是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
   

复制代码

回复 2# b114213903


   

b114213903

回复 3# super_two


    BioPerl模块自动就处理了！

super_two

原来是这样啊，谢谢


回复 4# b114213903


   

super_two

不能运行成功，报错是：Can't locate object method "length" via package "AAAAAATTTTTTAGATGAGCAGATTATGTGCTTCTTCCTCATTTACTTTTTTCGATGAA。。。。。。“，这个$seq->length()不能用吗？



回复 4# b114213903


   

b114213903

回复 6# super_two


    你是不是改脚本内容了？我这儿都运行正常！

super_two

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

$IN=Bio::SeqIO->new(-file=>"in.fastq",-format=>"fastq");中in.fastq  相应的变成了$infile.我的文件是.gz的压缩文件

回复 7# b114213903


   

b114213903

回复 8# super_two


    不能直接读压缩文档，只能读fastq格式的文档，可以通过管道，边解压边处理。

gaoyang636

lz本着学习Perl，BioPerl的心态来写的话，的确很好。但如果仅仅是想找个工具来处理数据，其实现在有大把现成的工具可以用哦，不但节省时间和精力，而且Performance还好的。比如你这个需求就可以使用Fastq toolkit完成。