Browse code

update bam_tally source for fix to quality tabulation

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@85550 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 15/01/2014 21:39:44
Showing 2 changed files

... ...
@@ -10,7 +10,7 @@ Description: GSNAP and GMAP are a pair of tools to align short-read
10 10
     to work with GMAP and GSNAP from within R. In addition, it provides 
11 11
     methods to tally alignment results on a per-nucleotide basis using 
12 12
     the bam_tally tool.
13
-Version: 1.5.4
13
+Version: 1.5.5
14 14
 Depends: R (>= 2.15.0), methods, GenomicRanges
15 15
 Imports: IRanges, Rsamtools (>= 1.7.4), rtracklayer (>= 1.17.15),
16 16
          GenomicFeatures, Biostrings, VariantAnnotation (>= 1.9.4), tools,
... ...
@@ -1,4 +1,4 @@
1
-static char rcsid[] = "$Id: bamtally.c 123356 2014-01-14 00:00:53Z twu $";
1
+static char rcsid[] = "$Id: bamtally.c 123611 2014-01-15 21:22:04Z twu $";
2 2
 #ifdef HAVE_CONFIG_H
3 3
 #include <config.h>
4 4
 #endif
... ...
@@ -701,8 +701,8 @@ Tally_clear (T this) {
701 701
   this->n_fromleft_minus = 0;
702 702
 
703 703
   if (this->use_array_p == true) {
704
-#if 0
705
-    /* print procedures clear each entry */
704
+#if 1
705
+    /* Note: these memset instructions are necessary to get correct results */
706 706
     memset((void *) this->matches_byshift_plus,0,this->n_matches_byshift_plus * sizeof(int));
707 707
     memset((void *) this->matches_byshift_minus,0,this->n_matches_byshift_minus * sizeof(int));
708 708
     memset((void *) this->matches_byquality,0,this->n_matches_byquality * sizeof(int));
... ...
@@ -2283,6 +2283,7 @@ iit_block (List_T *intervallist, List_T *labellist, List_T *datalist,
2283 2283
       debug2(printf("data for chrpos %u:\n",chrpos));
2284 2284
 
2285 2285
       /* Handle insertions at this position */
2286
+      debug2(printf("Pointers for insertions:\n"));
2286 2287
       pointers = push_int(&ignore,pointers,nbytes);
2287 2288
       if (this->insertions_byshift != NULL) {
2288 2289
 	unique_insertions_byshift = NULL;
... ...
@@ -2319,6 +2320,7 @@ iit_block (List_T *intervallist, List_T *labellist, List_T *datalist,
2319 2320
 	qsort(ins_array,ninsertions,sizeof(Insertion_T),Insertion_count_cmp);
2320 2321
 
2321 2322
 	/* Total number of different insertions at this position */
2323
+	debug2(printf("Number of insertions:\n"));
2322 2324
 	bytes = push_int(&nbytes,bytes,ninsertions);
2323 2325
 	for (i = 0; i < ninsertions; i++) {
2324 2326
 	  ins0 = ins_array[i];
... ...
@@ -2355,6 +2357,7 @@ iit_block (List_T *intervallist, List_T *labellist, List_T *datalist,
2355 2357
       }
2356 2358
 
2357 2359
       /* Handle deletions at this position */
2360
+      debug2(printf("Pointers for deletions:\n"));
2358 2361
       pointers = push_int(&ignore,pointers,nbytes);
2359 2362
       if (this->deletions_byshift != NULL) {
2360 2363
 	unique_deletions_byshift = NULL;
... ...
@@ -2391,6 +2394,7 @@ iit_block (List_T *intervallist, List_T *labellist, List_T *datalist,
2391 2394
 	qsort(del_array,ndeletions,sizeof(Deletion_T),Deletion_count_cmp);
2392 2395
 
2393 2396
 	/* Total number of different deletions at this position */
2397
+	debug2(printf("Number of deletions:\n"));
2394 2398
 	bytes = push_int(&nbytes,bytes,ndeletions);
2395 2399
 	for (i = 0; i < ndeletions; i++) {
2396 2400
 	  del0 = del_array[i];
... ...
@@ -2428,6 +2432,7 @@ iit_block (List_T *intervallist, List_T *labellist, List_T *datalist,
2428 2432
       
2429 2433
 
2430 2434
       /* Handle allele counts at this position */
2435
+      debug2(printf("Pointers for allele counts:\n"));
2431 2436
       pointers = push_int(&ignore,pointers,nbytes);
2432 2437
       if (pass_variant_filter_p(this->nmatches,this->mismatches_byshift,min_depth,variant_strands) == true) {
2433 2438
 	total_matches_plus = total_matches_minus = 0;
... ...
@@ -2461,6 +2466,7 @@ iit_block (List_T *intervallist, List_T *labellist, List_T *datalist,
2461 2466
 	}
2462 2467
 
2463 2468
 	/* Total signed counts at this position */
2469
+	debug2(printf("Total signed counts:\n"));
2464 2470
 	bytes = push_int(&nbytes,bytes,total_matches_plus + total_mismatches_plus);
2465 2471
 	bytes = push_int(&nbytes,bytes,total_matches_minus + total_mismatches_minus);
2466 2472
 	
... ...
@@ -2884,8 +2890,8 @@ revise_position (char querynt, char genomicnt, int mapq, int quality, int signed
2884 2890
       }
2885 2891
       
2886 2892
       if (max_shift_plus < this->n_matches_byshift_plus) {
2887
-	/* Clear existing array */
2888
-	memset(this->matches_byshift_plus,0,this->n_matches_byshift_plus * sizeof(int));
2893
+	/* Clear existing array.  Not necessary here if Tally_clear is doing the memset */
2894
+	/* memset(this->matches_byshift_plus,0,this->n_matches_byshift_plus * sizeof(int)); */
2889 2895
       } else {
2890 2896
 	/* Resize array */
2891 2897
 	oldsize = this->n_matches_byshift_plus;
... ...
@@ -2899,8 +2905,8 @@ revise_position (char querynt, char genomicnt, int mapq, int quality, int signed
2899 2905
       }
2900 2906
 
2901 2907
       if (max_shift_minus < this->n_matches_byshift_minus) {
2902
-	/* Clear existing array */
2903
-	memset(this->matches_byshift_minus,0,this->n_matches_byshift_minus * sizeof(int));
2908
+	/* Clear existing array.  Not necessary here if Tally_clear is doing the memset */
2909
+	/* memset(this->matches_byshift_minus,0,this->n_matches_byshift_minus * sizeof(int)); */
2904 2910
       } else {
2905 2911
 	/* Resize array */
2906 2912
 	oldsize = this->n_matches_byshift_minus;