ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/qseq2fastq.pl
Revision: 24
Committed: Tue Jul 26 21:46:39 2011 UTC (8 years, 2 months ago) by gpertea
File size: 1131 byte(s)
Log Message:
Line User Rev File contents
1 gpertea 23 #!/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     }

Properties

Name Value
svn:executable *