| 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 |
} |
