From 6ac64a4bc8eb738cb7e1f7c67a9841afd7147fd2 Mon Sep 17 00:00:00 2001
From: Tim Reddy Tim <treddy@hudsonalpha.org>
Date: Tue, 16 Dec 2008 17:55:32 +0000
Subject: [PATCH] Forgot to add count_reads.pm script

---
 htswanalysis/scripts/count_reads.pm | 31 +++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100755 htswanalysis/scripts/count_reads.pm

diff --git a/htswanalysis/scripts/count_reads.pm b/htswanalysis/scripts/count_reads.pm
new file mode 100755
index 0000000..85b203d
--- /dev/null
+++ b/htswanalysis/scripts/count_reads.pm
@@ -0,0 +1,31 @@
+#!/usr/bin/perl -w
+
+use strict;
+use warnings;
+
+my $align_file = shift;
+my $all_gz_file = shift;
+
+if(!defined($all_gz_file)) { print STDERR "Usage: count_reads.pm [align file] [all file.gz]\n"; exit(1); }
+
+my $pf_count = 0;
+my $adapter_count = 0;
+my $align_count = 0;
+
+open(ALIGN,$align_file);
+while(<ALIGN>) {
+  chomp;
+  my(@a) = split(/ /,$_);
+  $pf_count++;
+
+  if(scalar(@a) > 3) {
+    if($_ =~ /contam/) { $adapter_count++; } else { $align_count++; }
+  }
+}
+close(ALIGN);
+
+my $all_count = `cat $all_gz_file | gzip -d | wc -l`;
+chomp $all_count; $all_count =~ s/\s//;
+
+print "#N_reads N_pf N_adapter N_align\n";
+print "$all_count\t$pf_count\t$adapter_count\t$align_count\n";
-- 
2.30.2