1 |
#!/usr/bin/perl |
2 |
use strict; |
3 |
my ($minq, $maxq)=(255,0); |
4 |
while (<>) { |
5 |
chomp; |
6 |
my @t = split(/\t/); |
7 |
next unless $t[10]==1; |
8 |
print '@'."$t[0]:$t[2]:$t[3]:$t[4]:$t[5]#$t[6]/$t[7]\n"; |
9 |
die ("Error: uneven length !\n") unless length($t[9])==length($t[8]); |
10 |
$t[8]=~tr/./N/; |
11 |
$t[8]=~s/^N+//; |
12 |
my $d=length($t[9])-length($t[8]); |
13 |
$t[9]=substr($t[9],$d) if $d; |
14 |
$t[8]=~s/N+$//; |
15 |
$d=length($t[9])-length($t[8]); |
16 |
$t[9]=substr($t[9],0,-$d) if $d; |
17 |
die ("Error: uneven length after trimming!\n") |
18 |
unless length($t[9])==length($t[8]); |
19 |
my $q=qv($t[9]); |
20 |
print "$t[8]\n"; |
21 |
print "+\n"; |
22 |
print "$q\n"; |
23 |
} |
24 |
|
25 |
print STDERR "minq=$minq('".chr($minq)."'), maxq=$maxq('".chr($maxq)."')\n"; |
26 |
$minq-=31; |
27 |
$maxq-=31; |
28 |
print STDERR " After conversion:\n"; |
29 |
print STDERR "minq=$minq('".chr($minq)."'), maxq=$maxq('".chr($maxq)."')\n"; |
30 |
|
31 |
sub qv { |
32 |
my @s=unpack('C*', $_[0]); |
33 |
# assuming phred64 data was given, convert it to phred33 |
34 |
my $qual=''; |
35 |
foreach my $c (@s) { |
36 |
#my $newc=$c-26; |
37 |
my $newc=$c-31; |
38 |
$qual .= chr($newc); |
39 |
$minq=$c if $c<$minq; |
40 |
$maxq=$c if $c>$maxq; |
41 |
} |
42 |
return $qual; |
43 |
} |