+sub p2q_post_process {
+ my ($chr, $seq, $qual, $gaps, $l) = @_;
+ &p2q_filter_gaps($seq, $gaps, $l);
+ print "\@$chr\n"; &p2q_print_str($seq);
+ print "+\n"; &p2q_print_str($qual);
+}
+
+sub p2q_filter_gaps {
+ my ($seq, $gaps, $l) = @_;
+ for my $g (@$gaps) {
+ my $x = $g > $l? $g - $l : 0;
+ substr($$seq, $x, $l + $l) = lc(substr($$seq, $x, $l + $l));
+ }
+}
+
+sub p2q_print_str {
+ my ($s) = @_;
+ my $l = length($$s);
+ for (my $i = 0; $i < $l; $i += 60) {
+ print substr($$s, $i, 60), "\n";
+ }
+}
+
+#
+# varStats
+#
+
+sub varStats {
+ my %opts = (d=>'', c=>5);
+ getopts('d:c:', \%opts);
+ die("Usage: samtools.pl varStats [-d dbSNP.snp] [-c $opts{c}] <in.plp.snp>\n") if (@ARGV == 0 && -t STDIN);
+ my (@cnt, %hash);
+ my $col = $opts{c} - 1;
+ while (<>) {
+ my @t = split;
+ if ($t[2] eq '*') {
+ } else {
+ my $q = $t[$col];
+ $q = 99 if ($q > 99);
+ $q = int($q/10);
+ my $is_het = ($t[3] =~ /^[ACGT]$/)? 0 : 1;
+ ++$cnt[$q][$is_het];
+ $hash{$t[0],$t[1]} = $q;
+ }
+ }
+}
+
+#
+# Usage
+#
+