-
Notifications
You must be signed in to change notification settings - Fork 2
/
q20_filter.pl
63 lines (51 loc) · 1.05 KB
/
q20_filter.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#!/usr/bin/perl
use strict;
(my $usage = <<OUT) =~ s/\t+//g;
This script will extract high quality read q20:
OUT
die $usage unless @ARGV == 4;
my($f_db,$f_fa,$f_fq,$f_out)=@ARGV;
open(DB,"<$f_db");
open(IN1,"<$f_fa");
open(IN2,"<$f_fq");
open(OUT,">$f_out");
my %qscore=();
my %quality=();
while(<DB>)
{
my $line=$_;
chomp($line);
my @temp=split("\t",$line);
$qscore{$temp[0]}=$temp[1];
}
my $rd_id;
while(<IN2>)
{
my $line=$_;
chomp($line);
if($line=~/^\@/) { $rd_id=$line; $rd_id=~s/^\@//g; next; }
else { $quality{$rd_id}=$line; }
}
close IN2;
while(<IN1>)
{
my $line=$_;
chomp($line);
if($line=~/^\>/) { $rd_id=$line; $rd_id=~s/^>//g; }
else {
my $print_out=1;
my $qstr=$quality{$rd_id};
for(my $i=0;$i<length($qstr);$i++)
{
my $q=substr($qstr,$i,1);
if(defined $qscore{$q} && $qscore{$q}<20)
{
$print_out=0;
last;
}
}
if($print_out==1) { print OUT ">$rd_id\n"; print OUT $line,"\n"; }
}
}
close IN1;
close OUT;