diff --git a/htsutils.c b/htsutils.c index fc055a0..a8ebdc6 100644 --- a/htsutils.c +++ b/htsutils.c @@ -14,6 +14,10 @@ int parseBAM(Settings *props) bam_hdr_t *header = sam_hdr_read(fh_in); bam1_t *bam = bam_init1(); + int absent; + khiter_t k; + khash_t(pass_hash) *hptr = kh_init(pass_hash); + uint32_t lines = 0; uint32_t poly = 0; @@ -47,23 +51,48 @@ int parseBAM(Settings *props) if ((polyTailSize >= props->sizepoly) && (polyReadSize <= (bam->core.l_qseq - props->sizepoly))) { poly++; + + // add polyBasePosition to hash + k = kh_put(pass_hash, hptr, polyBasePosition, &absent); + if (absent) + kh_value(hptr, k) = 1; + else + kh_value(hptr, k)++; + } + + //printf("Tail: %d|%d\tHead: %d|%d\tSeq:: ", polyTailA,polyTailT,polyHeadA,polyHeadT); //dumpBamSeq(bam); //printf("\n"); lines++; - //if (lines == 10) + //if (lines == 1) // break; } + fprintf(stderr, "BAM Lines: %d / %d\n", poly, lines); + + // loop through hash and print values + fprintf(stderr, "HASH Size: %d\n", kh_size(hptr)); + for (k = kh_begin(hptr); k != kh_end(hptr); ++k) + { + if (kh_exist(hptr, k)) + { + int key = kh_key(hptr, k); + int val = kh_val(hptr, k); + printf("%d\t%d\n", key, val); + } + } + + kh_destroy(pass_hash, hptr); bam_destroy1(bam); bam_hdr_destroy(header); sam_close(fh_in); - printf("BAM Lines: %d / %d\n", poly, lines); + return 0; } diff --git a/htsutils.h b/htsutils.h index 9dd989d..7a7a3f3 100644 --- a/htsutils.h +++ b/htsutils.h @@ -12,12 +12,17 @@ #include #include #include + #include +#include + #include "config.h" #define POLY_HEAD 0 #define POLY_TAIL 1 +KHASH_MAP_INIT_INT(pass_hash, int); + int parseBAM(Settings *); void dumpBamSeq(bam1_t *); diff --git a/main.c b/main.c index b2c0225..956f523 100644 --- a/main.c +++ b/main.c @@ -31,7 +31,7 @@ int main(int argc, const char * argv[]) double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; - printf("Time: %.4f sec\n", time_spent); + fprintf(stderr, "Time: %.4f sec\n", time_spent); return 0; }